超越関数？

　超越関数にちゃんとした（広く通じる）定義はあるのか、という話。

　なんとなく「代数的数係数の多項式で表せない関数」ぐらいに考えていたけれど、そうすると、

は超越関数ということに？　それはちょっと違うような気が……。

　日本語版ウィキペディアから調べてみます。
超越関数 - Wikipedia

超越関数（ちょうえつかんすう、英: Transcendental Function）とは、係数が多項式であるような多項式で表せない関数である。より正確に言えば、1変数の関数が超越的であるのは、その変数について代数的独立性がある場合である。

　うーん。「係数が多項式であるような多項式」の意味が分かりません。どういう人向けの説明なんだろう？

　続いて英語版へ。
Transcendental function - Wikipedia, the free encyclopedia

A transcendental function is an analytic function that does not satisfy a polynomial equation, in contrast to an algebraic function. (The polynomials are sometimes required to have rational coefficients.) In other words, a transcendental function "transcends" algebra in that it cannot be expressed in terms of a finite sequence of the algebraic operations of addition, multiplication, and root extraction.

　こちらの方は明快！　超越関数とは多項式を満たさない解析関数で、時々その多項式は有理数係数のものに限定される。
　この定義だと、(1)は文脈によって超越関数に分類してもそうでないとしてもよさそうです。

　次はWolfram Mathworld。
Transcendental Function -- from Wolfram MathWorld

A function which is not an algebraic function. In other words, a function which "transcends," i.e., cannot be expressed in terms of, algebra. Examples of transcendental functions include the exponential function, the trigonometric functions, and the inverse functions of both.

　超越関数とは代数関数でない関数のことで、言い換えると、超越的な関数は代数的に表せない。……これは同語反復じゃないの？？
　ちなみに代数関数は、
Algebraic Function -- from Wolfram MathWorld

An algebraic function is a function f(x) which satisfies p(x,f(x))=0, where p(x,y) is a polynomial in x and y with integer coefficients. Functions that can be constructed using only a finite number of elementary operations together with the inverses of functions capable of being so constructed are examples of algebraic functions. Nonalgebraic functions are called transcendental functions.

　代数関数f(x)は、整数係数の多項式p(x, y)によって、p(x, f(x)) = 0を満たす。
　整数係数は意味的に有理数係数と等価で、もともとの考えに近いです。そして、(1)がはっきり超越関数の側に分類されます。

　どうなんだろうね、これは……。

GCD(0, 0) = ?

　前回わざと触れなかったところ。

　最大公約数は普通正の数に対して考えますが、実は片方が0でもOKです。
　「すべての整数は0の約数（0を割り切る）」と考えれば、0でないnと0の最大公約数はnになります。
　では、両方とも0の場合はどうすればいいでしょうか？

• 案A：∞とする

　すべての整数に最大値はないので、∞。または定義しない。
　論理的にはこれが正しいような……。

• 案B：0とする

　こちらは実用上の便利さを考慮すると捨てがたい案。
　3つ以上の数（うち1つ以上は0でない）の最大公約数を、例えばGCD(a, b, c) = GCD(GCD(a, b), c)として求めることができます。

　案Aは元の定義を重視して、案Bは結合的な二項演算の再定義という感じ？
　コンピュータ用なら、案Bを採用するのが無難かなーと思います。
　（参考：JavaのBigInteger、Commons Math、Guavaなどは案B。）

Binary GCD

　GCDは最大公約数。計算方法としてはユークリッドの互除法が有名です。
　別なやり方もあって、例えばBinary GCDアルゴリズムなんかおもしろいと思うのですが、なぜか日本語での解説が少ない……。
　なら自分で書いて残しておこうかと思ったり。

　Binary GCDアルゴリズムでは次の3つの定理が基本になります。

1. nとmが両方とも偶数なら、GCD(n, m) = 2*GCD(n/2, m/2)
2. nが偶数でmが奇数なら、GCD(n, m) = GCD(n/2, m)
   同じようにnが奇数でmが偶数なら、GCD(n, m) = GCD(n, m/2)
3. nとmが共に奇数でn > mなら、GCD(n, m) = GCD(n - m, m)
   同じようにnとmが共に奇数でn < mなら、GCD(n, m) = GCD(n, m - n)

　証明は略。（面倒……

　アルゴリズムの本体を文章に直して書くと次の通り。

1. nとmのどちらかが奇数になるまで2で割る（回数kを記録しておく）
2. 残りの側も奇数になるまで2で割る（定理2）
3. 以下、n = mになるまでくり返し
   • n > mのとき、
   1. nにn - mを代入（定理3）
   2. nが奇数になるまで2で割る（定理2）
   • n < mのとき、
   1. mにm - nを代入（定理3）
   2. mが奇数になるまで2で割る（定理2）
4. 最後にnに2^kを掛けて終わり！（定理1）

　互除法と比べると、特に任意精度整数（BigInteger）で除算を避けられる*1のがメリットらしいです。（それなら、CPUで除算できる普通の整数型では大差ないのかも？）

　Javaで書くとこうなります。ちょっとトリッキーな部分あり。

public static int gcd(int n, int m) {
    // 負の数はダメ
    if (n < 0 || m < 0) throw new ArithmeticException();
    // GCD(0, m) = m; GCD(0, 0) = 0
    // GCD(n, 0) = n
    if (n == 0 || m == 0) return n + m;
    // ここから
    int k = Integer.numberOfTrailingZeros(n | m);
    n >>= Integer.numberOfTrailingZeros(n);
    m >>= Integer.numberOfTrailingZeros(m);
    while (n != m) {
        if (n > m) {
            n -= m;
            n >>= Integer.numberOfTrailingZeros(n);
        }
        else {
            m -= n;
            m >>= Integer.numberOfTrailingZeros(m);
        }
    }
    return n << k;
}

（最初から最後までnumberOfTrailingZerosだらけ！）

*1:2で割るのはシフト演算で済むから。

罠が潜んでいる！

　ワンライナーの陥穽。

　剰余（n rem m）演算とモジュロ（n mod m）演算はnが負のとき結果が違ってきます。端的に言うとモジュロ演算で結果が負の数になることはありません。（詳細は英語版ウィキペディアのモジュロ演算。）
　例を挙げた方がいいか……。

n	-7	-6	-5	-4	-3	-2	-1	0	1	2	3	4
rem 3	-1	0	-2	-1	0	-2	-1	0	1	2	0	1
mod 3	2	0	1	2	0	1	2	0	1	2	0	1

　C言語やJavaの%演算子は剰余演算なので、モジュロ演算の方は（必要なら）自分で書くことになります。

public static int mod(int n, int m) {
    if (m <= 0) {
        throw new ArithmeticException();
    }
    int rem = n % m;
    if (n < 0) {
        return rem + m;
    }
    else {
        return rem;
    }
}

　一部ではif文抜きでモジュロ演算する裏技（？）として、

(n % m + m) % m

が紹介されています。が、これは途中で足し算がオーバーフローすると悲惨なことになるので避けましょう。

    int m = Integer.MAX_VALUE;
    System.out.println(mod(10, m));
    System.out.println((10 % m + m) % m);

10
-2147483639

　一見よさそうに見えて罠が潜んでいる、という話でした。了。

log2(x)

　JavaのMathクラスにはlog2(x)関数が見当たりません。
　他のプログラミング言語でもlog2(x)がない場合には役立つ話かも……？

　さて、下のやり方で何か問題があるでしょうか？

Math.log(x) / Math.log(2)

　対数の底の変換は、

なので数学的には自然というか素朴な方法ですね。

　しかし、浮動小数点演算では予想外の結果になってしまうことがあります。

    System.out.println(Math.log(1 << 29) / Math.log(2));

29.000000000000004

　xがぴったり2のn乗なら、正確にnを返してほしいと思いませんか？（思わなければ続きは読まなくて大丈夫です。）
　かといって2のn乗を特別扱いすると、単調性（ならば）が保たれるか分からないので、それはそれでうーんな感じです。
　ということで、なるべく誤差を発生させない計算方法を考えます。

　まず、

と表します。（kは整数で、sの範囲を制限すると決まる。）

　ここで、の絶対値が小さいとき、浮動小数点数の性質（0に近いほど細かく値の区別ができる）から誤差も小さく計算できます。

　ここではsの範囲を、

とします。（絶対にこうすべきというわけでもない。）
　あとは非正規化数に注意しつつコード化すると次のようになります。

public static double log2(double x) {
    // 特殊な値
    if (Double.isNaN(x) || x < 0.0) return Double.NaN;
    if (x == Double.POSITIVE_INFINITY) return Double.POSITIVE_INFINITY;
    if (x == 0.0) return Double.NEGATIVE_INFINITY;
    // ここから
    int k = Math.getExponent(x);
    if (k < Double.MIN_EXPONENT) {
        // 非正規化数は取扱い注意！
        k = Math.getExponent(x * 0x1.0p52) - 52;
    }
    if (k < 0) {
        k++;
    }
    double s = Math.scalb(x, -k);
    final double LOG2_E = 1.4426950408889634;
    return k + LOG2_E * Math.log(s);
}

　比較テスト。（）

    int j = 0, k = 0;
    for (int i = -1074; i <= 1023; i++) {
        double x = Math.scalb(1.0, i);
        if (log2(x) != i) j++;
        if (Math.log(x) / Math.log(2) != i) k++;
    }
    System.out.println("log2(x) function errors: " + j);
    System.out.println("Math.log(x) / Math.log(2) errors: " + k);

log2(x) function errors: 0
Math.log(x) / Math.log(2) errors: 441

　ライブラリでも作るのでなければこれで十分かと……。