15. 浮動小数点演算、その問題と制限

浮動小数点数はコンピューターのハードウェア上は2進数(binary)の分数で表されます。 たとえば、10進数 の分数では 0.625 は 6/10 + 2/100 + 5/1000 という値を持ち、2進数 の分数では 0.101 は 1/2 + 0/4 + 1/8 という値を持ちます。 この2つの分数はまったく同じ値を持ち、唯一異なる点は1つ目が10進数の分数で書かれており、2つ目は2進数の分数で書かれているということです。

残念なことに、ほとんどの小数は 2 進法の分数として正確に表わすことができません。その結果、一般に、入力した 10 進の浮動小数点数は、 2 進法の浮動小数点数で近似された後、実際にマシンに記憶されます。

最初は基数 10 を使うと問題を簡単に理解できます。分数 1/3 を考えてみましょう。分数 1/3 は、基数 10 の分数として、以下のように近似することができます:

0.3

さらに正確な近似は、

0.33

さらに正確な近似は、

0.333

となり、以後同様です。何個桁数を増やして書こうが、結果は決して厳密な 1/3 にはなりません。しかし、少しづつ正確な近似にはなっていくでしょう。

同様に、基数を 2 とした表現で何桁使おうとも、10 進数の 0.1 は基数を 2 とした小数で正確に表現することはできません。基数 2 では、1/10 は循環小数 (repeating fraction) となります

0.0001100110011001100110011001100110011001100110011...

どこか有限の桁で止めると、近似値を得ることになります。近年の殆どのコンピュータでは float 型は、最上位ビットから数えて最初の 53 ビットを分子、2 の冪乗を分母とした、二進小数で近似されます。1/10 の場合は、二進小数は 3602879701896397 / 2 ** 55 となります。これは、1/10 に近いですが、厳密に同じ値ではありません。

値が表示される方法のために、ほとんどのユーザは、近似に気づきません。Python はマシンに格納されている二進近似値の10進小数での近似値を表示するので、格納されている値が元の10進小数の近似値でしか無いことを忘れがちです。ほとんどのマシンで、もし Python が2進数で近似された 0.1 の近似値をそのまま10進数で表示していたら、その結果は次のようになったでしょう:

>>> 0.1
0.1000000000000000055511151231257827021181583404541015625

これは、ほとんどの人が必要と感じるよりも多すぎる桁数です。なので、Python は丸めた値を表示することで、桁数を扱いやすい範囲にとどめます:

>>> 1 / 10
0.1

表示された結果が正確に 1/10 であるように見えたとしても、実際に格納されている値は最も近く表現できる二進小数であるということだけは覚えておいてください。

幾つかの異なる10進数の値が、同じ2進有理数の近似値を共有しています。例えば、0.1 と 0.10000000000000001 と 0.1000000000000000055511151231257827021181583404541015625 はどれも 3602879701896397 / 2 ** 55 に近似されます。同じ近似値を共有しているので、どの10進数の値も eval(repr(x)) == x という条件を満たしたまま同じように表示されます。

昔の Python は、プロンプトと repr() ビルトイン関数は 17 桁の有効数字を持つ 0.10000000000000001 のような10進数の値を選んで表示していました。 Python 3.1 からは、ほとんどの場面で 0.1 のような最も短い桁数の10進数の値を選ぶようになりました。

この動作は2進数の浮動小数点にとってはごく自然なものです。これは Python のバグではありませんし、あなたのコードのバグでもありません。ハードウェアの浮動小数点演算をサポートしている全ての言語で同じ種類の問題を見つけることができます (いくつかの言語ではデフォルトの、あるいはどの出力モードを選んでも、この差を 表示 しないかもしれませんが)。

よりよい出力のために、文字列フォーマットを利用して有効桁数を制限した10進数表現を得ることができます:

>>> format(math.pi, '.12g')  # give 12 significant digits
'3.14159265359'

>>> format(math.pi, '.2f')   # give 2 digits after the point
'3.14'

>>> repr(math.pi)
'3.141592653589793'

これが、実際のコンピューター上の値の 表示 を丸めているだけの、いわば錯覚だということを認識しておいてください。

もう一つの錯覚を紹介します。例えば、0.1 が正確には 1/10 ではないために、それを3回足した値もまた正確には 0.3 ではありません:

>>> 0.1 + 0.1 + 0.1 == 0.3
False

0.1 はこれ以上 1/10 に近くなることができない値で、 0.3 もまた 3/10 に一番近い値なので、 round() 関数を使って計算前に丸めを行なっても意味がありません:

>>> round(0.1, 1) + round(0.1, 1) + round(0.1, 1) == round(0.3, 1)
False

数値を意図した正確な値に近づけることはできませんが、 math.isclose() 関数は不正確な値を比べるのに便利です:

>>> math.isclose(0.1 + 0.1 + 0.1, 0.3)
True

あるいは、 round() 関数を粗い近似値比較に使うこともできます:

>>> round(math.pi, ndigits=2) == round(22 / 7, ndigits=2)
True

このように2進数の浮動小数点の演算には多くの驚きがあります。「0.1」の問題について詳しい説明は、「表現エラー」セクションで行います。2進数の浮動小数点の仕組みと、実際によく遭遇する問題各種についての分かりやすい概要は、 Examples of Floating Point Problems を参照してください。その他よくある驚きの より詳細な説明は The Perils of Floating Point も参照してください。

究極的にいうと、"容易な答えはありません"。ですが、浮動小数点数のことを過度に警戒しないでください！ Python の float 型操作におけるエラーは浮動小数点処理ハードウェアから受けついたものであり、ほとんどのマシン上では一つの演算あたり高々 2**53 分の 1 です。この誤差はほとんどの作業で充分以上のものですが、浮動小数点演算は 10 進の演算ではなく、浮動小数点の演算を新たに行うと、新たな丸め誤差の影響を受けることを心にとどめておいてください。

異常なケースが存在する一方で、普段の浮動小数点演算の利用では、単に最終的な結果の値を必要な 10 進の桁数に丸めて表示するのなら、最終的には期待通りの結果を得ることになるでしょう。たいては str() で十分ですが、きめ細かな制御をしたければ、 書式指定文字列の文法 にある str.format() メソッドのフォーマット仕様を参照してください。

正確な10進数表現が必要となるような場合には、 decimal モジュールを利用してみてください。このモジュールは会計アプリケーションや高精度の計算が求められるアプリケーションに適した、10進数の計算を実装しています。

別の正確な計算方法として、 fractions モジュールが有理数に基づく計算を実装しています (1/3 のような数を正確に表すことができます)。

あなたが浮動小数点演算のヘビーユーザーなら、SciPy プロジェクトが提供している NumPy パッケージやその他の数学用パッケージを調べてみるべきです。 <https://scipy.org> を参照してください。

Python は 本当に float の正確な値が必要なレアケースに対応するためのツールを提供しています。 float.as_integer_ratio() メソッドは float の値を有理数として表現します:

>>> x = 3.14159
>>> x.as_integer_ratio()
(3537115888337719, 1125899906842624)

この分数は正確なので、元の値を完全に復元することができます:

>>> x == 3537115888337719 / 1125899906842624
True

float.hex() メソッドは float の値を16進数で表現します。この値もコンピューターが持っている正確な値を表現できます:

>>> x.hex()
'0x1.921f9f01b866ep+1'

この正確な16進数表現はもとの float 値を正確に復元するために使うことができます:

>>> x == float.fromhex('0x1.921f9f01b866ep+1')
True

この16進数表現は正確なので、値を (プラットフォームにも依存せず) バージョンの異なるPython 間でやり取りしたり、他のこのフォーマットをサポートした言語 (Java や C99 など) と正確にやり取りするのに利用することができます。

別の便利なツールとして、合計処理における精度のロスを緩和してくれる sum() 関数があります。これは累計加算中の丸めに拡張精度を使います。これにより、誤差が最終的な合計値に影響を与えるまで蓄積されなくなり、結果が改善されます:

>>> 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 == 1.0
False
>>> sum([0.1] * 10) == 1.0
True

The math.fsum() goes further and tracks all of the "lost digits" as values are added onto a running total so that the result has only a single rounding. This is slower than sum() but will be more accurate in uncommon cases where large magnitude inputs mostly cancel each other out leaving a final sum near zero:

>>> arr = [-0.10430216751806065, -266310978.67179024, 143401161448607.16,
...        -143401161400469.7, 266262841.31058735, -0.003244936839808227]
>>> float(sum(map(Fraction, arr)))   # Exact summation with single rounding
8.042173697819788e-13
>>> math.fsum(arr)                   # Single rounding
8.042173697819788e-13
>>> sum(arr)                         # Multiple roundings in extended precision
8.042178034628478e-13
>>> total = 0.0
>>> for x in arr:
...     total += x                   # Multiple roundings in standard precision
...
>>> total                            # Straight addition has no correct digits!
-0.0051575902860057365

15.1. 表現エラー

この章では、"0.1" の例について詳細に説明し、このようなケースに対してどのようにすれば正確な分析を自分で行えるかを示します。ここでは、 2 進法表現の浮動小数点数についての基礎的な知識があるものとして話を進めます。

表現エラー(Representation error)は、いくつかの (実際にはほとんどの) 10 進の小数が 2 進法 (基数 2)の分数として表現できないという事実に関係しています。これは Python (あるいは Perl, C, C++, Java, Fortran. およびその他多く) が期待通りの正確な 10 進数を表示できない主要な理由です。

なぜそうなるのでしょう？ 1/10 は２進法の小数で厳密に表現できません。少なくとも2000年以降、ほぼすべてのマシンは IEEE 754 2進数の浮動小数点演算を用いており、ほぼすべてのプラットフォームでは Python の浮動小数点を IEEE 754 binary64 "倍精度 (double precision)" 値に対応付けます。 IEEE 754 binary64 値は 53 ビットの精度を持つため、計算機に入力を行おうとすると、可能な限り 0.1 を最も近い値の分数に変換し、J/2**N の形式にしようと努力します。J はちょうど 53 ビットの精度の整数です。

1 / 10 ~= J / (2**N)

を書き直すと

J ~= 2**N / 10

となります。 J は厳密に 53 ビットの精度を持っている (>= 2**52 だが < 2**53 ) ことを思い出すと、 N として最適な値は 56 になります:

>>> 2**52 <=  2**56 // 10  < 2**53
True

すなわち、56 は J をちょうど 53 ビットの精度のままに保つ N の唯一の値です。J の取りえる値はその商を丸めたものです:

>>> q, r = divmod(2**56, 10)
>>> r
6

剰余が 10 の半分以上なので、最良の近似は切り上げて丸めたものになります。:

>>> q+1
7205759403792794

従って、IEEE 754の倍精度における 1/10 の取りえる最良の近似は:

7205759403792794 / 2 ** 56

分子と分母を2で割って分数を小さくします:

3602879701896397 / 2 ** 55

丸めたときに切り上げたので、この値は実際には 1/10 より少し大きいことに注目してください。 もし切り捨てをした場合は、商は 1/10 よりもわずかに小さくなります。どちらにしろ 厳密な 1/10 ではありません！

つまり、計算機は 1/10 を "理解する" ことは決してありません。計算機が理解できるのは、上記のような厳密な分数であり、IEEE 754 の倍精度浮動小数点数で得られるもっともよい近似は以下になります:

>>> 0.1 * 2 ** 55
3602879701896397.0

この分数に 10**55 を掛ければ、55 桁の十進数の値を見ることができます:

>>> 3602879701896397 * 10 ** 55 // 2 ** 55
1000000000000000055511151231257827021181583404541015625

これは、計算機が記憶している正確な数値が、10 進数値 0.1000000000000000055511151231257827021181583404541015625 にほぼ等しいということです。多くの言語 (古いバージョンの Python を含む) では、完全な 10 進値を表示するのではなく、結果を有効数字 17 桁に丸めます:

>>> format(0.1, '.17f')
'0.10000000000000001'

fractions モジュールと decimal モジュールを使うとこれらの計算を簡単に行えます:

>>> from decimal import Decimal
>>> from fractions import Fraction

>>> Fraction.from_float(0.1)
Fraction(3602879701896397, 36028797018963968)

>>> (0.1).as_integer_ratio()
(3602879701896397, 36028797018963968)

>>> Decimal.from_float(0.1)
Decimal('0.1000000000000000055511151231257827021181583404541015625')

>>> format(Decimal.from_float(0.1), '.17')
'0.10000000000000001'