log2

C

トラックバック元のコードを読んでいて、アレ？と思った。ショートコードに挑戦している人たちでもやってしまうものか？
さて、以下のコードの問題点は何でしょう。

/**
 * Returns binary logarithm of n.
 */
int intlog2(int n)
{
    return log(1.0 * n) / log(2.0);
}

遅いとかそういう問題ではないです。あしからず。

追記：簡単な解説

分からない人は以下のコードをコンパイルして実行してみましょう。

#include <stdio.h>
#include <math.h>
#include <limits.h>

int main()
{
    unsigned int i;

    for (i = 0; i < CHAR_BIT * sizeof(unsigned int); i++) {
        unsigned int n = 0x01 << i;
        double log2 = log((double) n) / log(2.0);
        double diff = log2 - i;
        printf("%u\t%u\t%.2f\t%.2e\n", i, (unsigned int) log2, log2, diff);
    }

    return 0;
}

手元の環境(gcc 4.1 on coLinux) ならこうなる。

0       0       0.00    0.00e+00
1       1       1.00    0.00e+00
2       2       2.00    0.00e+00
3       2       3.00    -1.60e-16
4       4       4.00    0.00e+00
5       5       5.00    1.60e-16
6       5       6.00    -3.20e-16
7       6       7.00    -1.60e-16
8       8       8.00    0.00e+00
9       9       9.00    1.60e-16
10      10      10.00   3.20e-16
11      11      11.00   4.81e-16
12      11      12.00   -6.41e-16
13      12      13.00   -4.81e-16
14      13      14.00   -3.20e-16
15      14      15.00   -1.60e-16
16      16      16.00   0.00e+00
17      17      17.00   1.60e-16
18      18      18.00   3.21e-16
19      19      19.00   4.81e-16
20      20      20.00   6.40e-16
21      21      21.00   8.01e-16
22      22      22.00   9.61e-16
23      23      23.00   1.12e-15
24      23      24.00   -1.28e-15
25      25      25.00   1.44e-15
26      25      26.00   -9.61e-16
27      27      27.00   1.76e-15
28      27      28.00   -6.40e-16
29      29      29.00   2.08e-15
30      29      30.00   -3.21e-16
31      31      31.00   2.40e-15

当然、コンピュータによる計算だから誤差が含まれる。浮動小数点のまま扱うならこの程度の誤差はこの程度の誤差でしかないわけだけど、6.0 - 3.20e-16 を int にキャストすると 5 になってしまって誤差がひどく拡大してしまってひどい目にあうという話。
で、この場合、誤差が負であることが問題で、入出力は所詮整数値なので、

log(n + 0.5) / log(2.0)

として適当に誤差が正になるように offset してやるという workaround もある。
ちなみに、MSVC8 だと誤差がすべて正の値になって特に問題なかったりするけど、基本的に int にキャストする場合は気をつけないといけない。