Python の小数

# この話のゴール

* Pythonの小数の表現について理解する
* 小数の表現による誤差、演算の誤差を知る

# 10回足してみる話

1を10回足すと10以上に（というか10に等しく）なる。

```py:int_sum.py
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回足すと？

```py:float_sum.py
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言語の簡単なプログラムを使って、小数の表現が実際どのようになっているかを解析してみましょう。

```c:float.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回足すという計算で、途中の値をプリントしてみましょう。

```c:float_sum.c
#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 で実際に同じ表現が用いられていることを確認しておきましょう。

```py:float_sum_hex.py
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で表現されています。

```c:cpython/Include/floatobject.h
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)