logの計算 - Void 〜tomの雑記〜

たまたまlogの計算をするアルゴリズムを考えなくてはいけなかったので今回はその話。というか残しとかないと忘れそうなんで。
え？math libraryを使えば早いじゃないかって？いやいや、整数のみで演算がしたいんですよ。さいわい、底は2なのでそんなに難しくはない、と思っていたら、欲しいのは小数点１桁の精度。つまり10log2(x)が欲しい訳だ。これは困ったね。
とりあえずできたんで、ソース貼っておこう。

int floorLog2(unsigned int n) {
  if (n == 0)
    return -1;
 
  int pos = 0;
  if (n >= (1 <<16)) { n >>= 16; pos += 16; }
  if (n >= (1 << 8)) { n >>=  8; pos +=  8; }
  if (n >= (1 << 4)) { n >>=  4; pos +=  4; }
  if (n >= (1 << 2)) { n >>=  2; pos +=  2; }
  if (n >= (1 << 1)) {           pos +=  1; }
  return pos;
}

int calc10log2(unsigned int x)
{
    unsigned int y, z, i;
    int n = 0, sm = 0;
    int log = 0;

    n = floorLog2(x);
    if (n == -1)
        return -1;
    log = n << 10;
    if ((2 << (n-1)) == x)
        return ((log >> 6) + (log >> 8) + 1) >> 1;

    y = (x << 10) >> n;
    z = 1024;
    for(i=0;i<10;i++) {
        m = 0;
        while (y < 2048) {  /* 2048 = 1024*2 */
            y *= y;
            y >>= 10;
            sm++;
            z >>= 1;
        }
        y >>= 1;
        log += z;
        if (z == 0 || sm > 7) /* 2^7 = 128 */
            break;
    }
    return ((log >> 6) + (log >> 8) + 1) >> 1;
}

一見すると何やってるか不明なアルゴリズムだな、こりゃ。
元ネタはBinary logarithm -Wikipedia。
前半の関数はlog2(x)の整数部分を求めている。後半が、残りの部分。
さて、前半部分だが、log2(x)の整数部分が欲しい訳だ。つまりとなるような n を求めれば良い。２進数で考えると

だから、一番左に１がたっているビットの位置を求めれば良いことになる。という訳で、前半の部分はビットに1が立っているかを見ている。log2(x)の整数部分だけを求めたいのならここで終了で万事解決なのだが、小数部分を求めたいとなるととたんに難しくなる。

問題の後半が、小数部分を計算している。もちろん、 $x=2^n$ なら問題はない。それ以外のときどうしたら良いか。それがBinary logarithm -Wikipediaのreal numberの項目の部分になる。もちろんこのままやると小数演算になるので、全体を1000倍するなどして整数化し、最後に100で割ることで10logを計算する。あ、上のソースでは1000ではなく1024倍して整数化している。
この部分、ちょっとくだいておこう。

となるので、 $y=2^{-n}x$ とおけばlog2(x)=n+log2(y)となる。ここでであるので、全体を $2^{-n}$ 倍すれば、となる。つまりである。
さて、となる最小のmをMとするとの両辺にlogをとって
となるので

である。従って

となる。ここで y=z/2とおくと2<=z<4であるから 1<=y<2である。ここまでくると気づいた人もいるでしょう。そう、条件が同じなので上の手順がもう一回使える訳だ。ということは得られる結果もやはり

となるので、これを繰り返していけば

が得られる訳だ。さて、これを使って整数演算をすれば良いわけだが、上の式は無限級数なので、どこかで打ち切らないといけない。求めたい精度によるが、小数点第１位だったら1/100以下になる項で打ち切れば良い。
あとは、全体を1000倍して得られた結果を100で割れば整数演算のみで小数点第１位まで求められる。ソースでは、1024=2の10乗を使って、シフト演算で処理している。
ちなみにソースの方では四捨五入してるので、演算がややこしくなってます。

という感じで、アルゴリズムというよりは数学的要素が強いよね。
しっかし、よう考えるよな。上のような変形。結構トリッキーだよねぇ。