int(0.5005*10000)が5004になるのはなぜですか。

awkについて語るスレ $2 - int(0.5005*10000)が5004になるのはなぜですか。からです。 この問題は以前、PHPの演算誤差で問題になったものですが、結論から言えば、「ある有限の有効桁数の範囲の 2 進数で 0.5005 を正確に表現できないため発生する誤差」です。 詳細はときどきの雑記帖 リターンズ - 7021 != 7021 に書かれてありますので熟読をしてみてください。

以下の内容は上記ページの 0.5005 で計算させたものとなっています。

上記ときどきの雑記帖 リターンズ - 7021 != 7021 にかかれてある awk のスクリプトを 0.5005 にして実行してみましょう。 (少し変更してあります)

BEGIN {
    val = 5005;
    limit = 1000;
    for (i = 0; val != 0 && i < limit; i++) {
        printf("%7.5f x 2 -> %7.5f ... ", val / 10000, val * 2 / 10000);
        val *= 2
        if (val >= 10000) {
            print "1";
            val -= 10000;
        }
        else {
            print "0";
        }
    }

    if (val != 0)
        print "...まだ続きます"
}

実行してみます。

$ nawk -f test_0.5005.awk
0.50050 x 2 -> 1.00100 ... 1
0.00100 x 2 -> 0.00200 ... 0
0.00200 x 2 -> 0.00400 ... 0
0.00400 x 2 -> 0.00800 ... 0
0.00800 x 2 -> 0.01600 ... 0
0.01600 x 2 -> 0.03200 ... 0
0.03200 x 2 -> 0.06400 ... 0
0.06400 x 2 -> 0.12800 ... 0
0.12800 x 2 -> 0.25600 ... 0
0.25600 x 2 -> 0.51200 ... 0
0.51200 x 2 -> 1.02400 ... 1
<snip>
0.67200 x 2 -> 1.34400 ... 1
0.34400 x 2 -> 0.68800 ... 0
...まだ続きます

というわけで、IEEE754 で示される倍精度 (nawk, gawk の場合 53 bit) では表現できません。

同様に上記ページに載っているビットパターンダンププログラムを使ってみることにします。

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

char *tbl[] = {
    "0000", "0001", "0010", "0011", "0100", "0101", "0110", "0111",
    "1000", "1001", "1010", "1011", "1100", "1101", "1110", "1111",
};

int
main()
{
    double val = 0.5005;
    double tval;
    char *p;
    int exponent;
    char mantissa[64] = {0};
    int i,j;

    p = (char *)&val;
    exponent = (p[7] & 0x7f) * 16 + ((p[6] >> 4) & 0xf) - 1023;

    strcpy(&mantissa[ 0], tbl[p[6]&0xf]);
    for (i=5, j=0; i>=0; i--) {
        strcpy(&mantissa[j+4], tbl[(p[i] >> 4)&0xf]);
        strcpy(&mantissa[j+8], tbl[p[i]&0xf]);
        j+=8;
    }

    printf("%5.4f = %c2^%d * 1.%s\n",
           val,
           (p[7] & 0x80) ? '-' : ' ',
           exponent, mantissa);

    for (tval=0, p=mantissa+strlen(mantissa); p>=mantissa; p--)
        tval = ((tval / 2) + (*p=='1' ? 1 : 0));

    tval = tval/2 + 1;
    tval *= pow(2, exponent);
    printf("%18.16f\n", tval);

    return 0;
}

コンパイル、実行してみましょう。

$ gcc -o bit_pattern_dump -O -lm bit_pattern_dump.c

$ ./bit_pattern_dump
0.5005 =  2^-1 * 1.0000000001000001100010010011011101001011110001101010
0.5004999999999999

つまり、0.5004999999999999 に 10000 をかけたものを整数化しているわけですから、結果的に 5004 になります。

なお、任意桁数での 2 進数での値は bc を用いると便利です。

$ echo 'scale=100; ibase=10; obase=2; 5005 / 10000' | bc
.1000000000100000110001001001101110100101111000110101001111110111110\
01110110110010001011010000111001010110000001000001100010010011011101\
00101111000110101001111110111110011101101100100010110100001110010101\
10000001000001100010010011011101001011110001101010011111101111100111\
01101100100010110100001110010101100000010000011000100100110111

さて、gawk の拡張版として xgawk がありますが、こちらは「MPFR の関数に限って基数が 10 で計算されます」。デフォルトでは 53 bit で計算されますが、こちらで計算してみましょう。

#! /usr/local/bin/xgawk -f
# mpfr_0.5005.awk

@load mpfr

BEGIN {
    print mpfr_rint(mpfr_mul(0.5005, 10000));
    printf("%d\n", mpfr_rint(mpfr_mul(0.5005, 10000)));
}

面倒なのが数値を扱うものは全て MPFR の関数で計算しないと通常の gawk の挙動、すなわち基数が 2 として扱われてしまうようです。 実行してみます。

$ xgawk -f mpfr_0.5005.awk
5.0050000000000000E3
5005

というわけで基数が 10 の世界では 0.5005 が表現できるので、5005 という答えを得ることができました。

tag_nawk.pngtag_nawk.pngtag_nawk.pngtag_nawk.png