この話のゴール
- Pythonの小数の表現について理解する
- 小数の表現による誤差、演算の誤差を知る
10回足してみる話
1を10回足すと10以上に（というか10に等しく）なる。
x = 0
for i in range(10):
x += 1
print(x < 10)
print(x == 10)
print(type(x))
当たり前ですよね。
% python int_sum.py
False
True
<type 'int'>
では、0.1を10回足すと？
x = 0.0
for i in range(10):
x += 0.1
print(x < 1)
print(type(x))
1以上に…なってない
% python float_sum.py
True
<type 'float'>
計算機上では数値を2進数で表現しますが、小数を有限桁の2進数で表現しているために、こういったことが起こります。
詳しい解説が
http://docs.python.jp/3/tutorial/floatingpoint.html
にありますが、ここではその解説の前提になっている「浮動小数点数」の表現方法について見ていきましょう。
2進数での数の表し方
本題に入る前に2進数による数の表し方をおさらいしておきましょう。
10進数による数の表し方は0〜9の数字と10のべき乗を使って、例えば
2015
= 2*1000 + 0*100 + 1*10 + 5*1
= 2*10^3 + 0*10^2 + 1*10^1 + 5*10^0
のように書いていました。
同じように2進数では0〜1の数字と2のべき乗で、
11001011
= 1*2^7 + 1*2^6 + 0*2^5 + 0*2^4 + 1*2^3 + 0*2^2 + 1*2^1 + 1*2^0
のように数を表します。なお、2進数表記は桁数が多くて書くのが大変なので、説明する時は16進数を使ってコンパクトな形で書くことが多いです。例えば2進数の11001011は16進数では0xcbと書きます。
浮動小数点数の表現
C言語の簡単なプログラムを使って、小数の表現が実際どのようになっているかを解析してみましょう。
#include <stdio.h>
int main(int argc, char *argv[])
{
double z;
if (argc != 2) {
fprintf(stderr, "usage: %s str\n", argv[0]);
return -1;
}
if (sscanf(argv[1], "%lf", &z) != 1) {
fprintf(stderr, "parse failed\n");
return -1;
}
printf("double: %f\n", z);
printf("hex: %016llx\n", *(unsigned long long *)&z);
return 0;
}
このプログラムでは
- 引数として与えられた数値を小数として読み取り
- 読み取った値を10進の小数としてプリント
- そのビット表現を16進数でプリント
ということをしています。
% gcc -o float float.c
でコンパイルし、引数を与えて実行してみると以下のようになります。
% ./float 1
double: 1.000000
hex: 3ff0000000000000
% ./float -1
double: -1.000000
hex: bff0000000000000
% ./float 0.5
double: 0.500000
hex: 3fe0000000000000
% ./float -0.5
double: -0.500000
hex: bfe0000000000000
% ./float 2
double: 2.000000
hex: 4000000000000000
1と-1、0.5と-0.5の16進数表記を見比べると、符号の違いが先頭が 3 か b かの違いと対応していそうです。それぞれ2進数で表記すると
- 0011
- 1011
なので、最上位のビットが0なら正の数、1なら負の数だと予想できます。
それが正しければ、
- 4 (16進) = 0100 (2進)
- c (16進) = 1100 (2進)
なので、-2 は c000000000000000 になるはずです。実際試してみましょう。
% ./float -2
double: -2.000000
hex: c000000000000000
予想が正しかったようです。次に符号ビットを除いた部分に着目すると
- 2 (2^1): 4000000000000000
- 1 (2^0): 3ff0000000000000
- 0.5 (2^(-1)): 3fe0000000000000
と、2倍すると先頭3文字部分(符号ビットを除いた11ビットの値)が1増えているように見えます。したがって 4 (2^2) は 4010000000000000 になると予想できますが、実際確認してみると
% ./float 4
double: 4.000000
hex: 4010000000000000
なりましたね。つまり符号ビットに続く11ビットで指数 (2の何乗をかけた数か)を表現していることがわかりました。
では、残りの13文字(52ビット)は何を表わしているのでしょうか？
2進数の 1.1、つまり 1 + 1/2 = 1.5 を与えてみると
% ./float 1.5
double: 1.500000
hex: 3ff8000000000000
上位に 0x8 (1000) が現われました。2進数の1.11 つまり、1 + 1/2 + 1/4 = 1.75 だと
% ./float 1.75
double: 1.750000
hex: 3ffc000000000000
この場合は 0xc (1100) から始まっています。つづけて 2進数の 1.111 つまり 1 + 1/2 + 1/4 + 1/8 = 1.875 だと
% ./float 1.875
double: 1.875000
hex: 3ffe000000000000
となり、ここまでをまとめると
- 1.0000: 0000
- 1.1000: 1000
- 1.1100: 1100
- 1.1110: 1110
左右を見比べると、左の値の小数点より上にある1が省略されたものが右の値になっているように思えます。
では、1.00001 だとどうなるでしょうか？ 1 + 1/32 = 1.03125 を引数に与えると
% ./float 1.03125
double: 1.031250
hex: 3ff0800000000000
0x08 (0000 1000) となりました。
- 1.00001: 00001
この場合も小数点より上にある1が省略されたということで説明できます。
最初の方で見ていた2のべき乗の場合は下位52ビット分が0になっていましたが、2のべき乗は2進数表記すれば最上位が1で残りが0となりますので、それも説明が付きます。
ここまで観察したことから、小数の表現を式で表すと
(-1)^(x >> 63) * 2^((0x7ff & (x >> 52)) - 0x3ff) * ((x & 0x000fffffffffffff) * 2^(-52) + 1)
となります。なお、0 はこの形だと表せませんが、掛け算をする時などに重要な役割を果たすので、
% ./float 0
double: 0.000000
hex: 0000000000000000
と、全てのビットを 0 にしたもので表すことになっています。
まとめると
- 最上位の1ビットで符号を表す
- 続く11ビットで指数部を表す
- 残り52ビットに、暗黙の最上位ビットを加えた53ビットで仮数（かすう）部を表す
- 0 は全てのビットが 0 で表す
という形で、小数が表現されていました。
足し算
最初に見ていた、0.1を10回足すという計算で、途中の値をプリントしてみましょう。
#include <stdio.h>
int main()
{
double z = 0.0;
for (int i = 0; i < 10; i++)
{
z += 0.1;
printf("hex: %016llx\n", *(unsigned long long *)&z);
}
return 0;
}
% gcc -o float_sum float_sum.c
% ./float_sum
hex: 3fb999999999999a
hex: 3fc999999999999a
hex: 3fd3333333333334
hex: 3fd999999999999a
hex: 3fe0000000000000
hex: 3fe3333333333333
hex: 3fe6666666666666
hex: 3fe9999999999999
hex: 3feccccccccccccc
hex: 3fefffffffffffff
まず、最初の行は10進数の 0.1 を表しているはずですが、
2^((0x3fb - 0x3ff) * (0x1999999999999a * 2^(-52))
= 0x1999999999999a * 2^(-56)
= 7205759403792794.0 / 72057594037927936.0
となっており、正確に0.1を表現できていません(0.1よりちょっと大きな数になっています)。10進数の0.1を2進数で表現すると循環小数になり、16進数で9が無限に続く形になりますが、それを有限桁数で表現するために誤差が入っていることを示しています。
次に、足し算の計算の仕方を見ていきましょう。
浮動小数点数は指数と仮数で表現されていましたが、その足し算は 10進数でいうところの
10 * 3 + 100 * 5 = 100 * (0.3 + 5) = 100 * 5.3
と同じようなことをします。実際、小数を (指数, 仮数) という形で書くことにすると、
3fb999999999999a + 3fb999999999999a
= (3fb, 0x1999999999999a) + (3fb, 0x1999999999999a)
= (3fb, 0x1999999999999a + 0x1999999999999a)
= (3fb, 0x33333333333334)
= (3fc, 0x1999999999999a)
= 3fc999999999999a
3fc999999999999a + 3fb999999999999a
= (3fc, 0x1999999999999a) + (3fb, 0x1999999999999a)
= (3fc, 0x1999999999999a) + (3fc, 0xccccccccccccd)
= (3fc, 0x1999999999999a + 0xccccccccccccd)
= (3fc, 0x26666666666667)
= (3fd, 0x13333333333334)
= 3fd3333333333334
ここで仮数のシフト時の丸めのルール
- シフトアウトしたビットの最上位が0なら切り捨て
- シフトアウトしたビットの最上位が1で、シフトアウトした最上位以外のビットに1があれば切り上げ
- シフトアウトしたビットの最上位が1で、シフトアウトした最上位以外のビットに1がなければ
- 残った数の最下位ビットが1なら切り上げ
- 残った数の最下位ビットが0なら切り捨て
を使っています
0x33333333333334 = ...00110100
1ビットシフト→ ...1010 = 0x1999999999999a
0x999999999999a = ...10011010
1ビットシフト→ ...1101 = 0xccccccccccccd
0x26666666666667 = ...01100111
1ビットシフト→ ...0100 = 0x13333333333334
ここでも53ビットで指数を表現するために丸めを行う必要があり、さらに誤差が混入していることがわかります。これらの誤差によって、0.1を10回足しても1に満たないという結果が出てしまいます。
Pythonの小数の表現
C言語のプログラムを使って浮動小数点数の表現について理解したところで、
最後に Python で実際に同じ表現が用いられていることを確認しておきましょう。
x = 0.0
for i in range(10):
x += 0.1
print(x.hex())
% python float_sum_hex.py
0x1.999999999999ap-4
0x1.999999999999ap-3
0x1.3333333333334p-2
0x1.999999999999ap-2
0x1.0000000000000p-1
0x1.3333333333333p-1
0x1.6666666666666p-1
0x1.9999999999999p-1
0x1.cccccccccccccp-1
0x1.fffffffffffffp-1
書き方が若干違いますが、それぞれC言語で見ていたものと同じ値が現れていますね。
実際 CPython の実装を確認すると、Pythonのfloatの値はC言語のdoubleで表現されています。
typedef struct {
PyObject_HEAD
double ob_fval;
} PyFloatObject;
まとめ
- Python(と C言語)の小数の表現 (IEEE 754) をいくつかの値を元に読み解いた
- Pythonのfloat (C言語の double) は64ビットを使って符号(1)、指数(11)、仮数(52+1)で表現される
- 固定のビット数で一定の精度の数を表現できる。
- 誤差により、0.1を10回足しても1にならないといったことが起こる
- 0.1 を有限桁の2進数で表現するために生じる誤差
- 足し算のために桁あわせをする際に生じる誤差
- 10進数の小数を正確に表すには→ decimal (http://docs.python.jp/3/library/decimal.html)