.. _tut-fp-issues:

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

.. sectionauthor:: Tim Peters <tim_one@users.sourceforge.net>


浮動小数点数は、計算機ハードウェアの中では、基数を 2 とする (2進法の)
分数として表現されています。例えば、小数

::

   0.125

は、 1/10 + 2/100 + 5/1000 という値を持ちますが、これと同様に、 2 進法の分数

::

   0.001

は 0/2 + 0/4 + 1/8 という値になります。これら二つの分数は同じ値を持っていますが、
ただ一つ、最初の分数は基数 10 で記述されており、二番目の分数は基数 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...

どこか有限の桁で止めると、近似値を得ることになります。
これこそが、次のような事態に出くわす理由です。

::

   >>> 0.1
   0.10000000000000001

今日では、ほとんどのマシンでは、0.1 を Python のプロンプトから入力すると上の
ような結果を目にします。そうならないかもしれませんが、これはハードウェアが
浮動小数点数を記憶するのに用いているビット数がマシンによって異なり、Python は
単にマシンに 2 進で記憶されている、真の 10 進の値を近似した値を、されに 10 進で
近似して出力するだけだからです。ほとんどのマシンでは、Python が 0.1 を
記憶するために 2 進近似した真の値を10 進で表すと、以下のような出力になるでしょう！

::

   >>> 0.1
   0.1000000000000000055511151231257827021181583404541015625

Python プロンプトは、文字列表現を得るために何に対しても :func:`repr` を
使います。浮動小数点数の場合、 ``repr(float)`` は真の 10 進値を有効数字
17 桁で丸め、以下のような表示を行います。

::

   0.10000000000000001

``repr(float)`` が有効数字 17 桁の値を生成するのは、この値が (ほとんどの
マシン上で) 、全ての有限の浮動小数点数 *x* について ``eval(repr(x)) == x``
が成り立つのに十分で、かつ有効数字 16 桁に丸めると成り立たないからです。

これは 2 進法の浮動小数点の性質です: Python のバグでも、ソースコードの
バグでもなく、浮動小数点演算を扱えるハードウェア上の、すべての言語で
同じ類の現象が発生します (ただし、言語によっては、デフォルトのモードや
全ての出力モードでその差を *表示しない* かもしれません)。

Python の組み込みの :func:`str` 関数は有効数字 12 桁しか生成しません。
このため、この関数を代わりに使用したいと思うかもしれません。
この関数は ``eval(str(x))`` としたときに *x* を再現しないことが多いですが、
出力を目で見るには好ましいかもしれません。

::

   >>> print str(0.1)
   0.1

現実には、上の表示は錯覚であると気づくのは重要なことです。
マシン内の値は厳密に 1/10 ではなく、単に真のマシン内の  *表示される値*
を丸めているだけなのです。

まだ驚くべきことがあります。例えば、以下を見て、

::

   >>> 0.1
   0.10000000000000001

:func:`round` 関数を使って桁を切り捨て、期待する 1 桁にしたい誘惑に
かられたとします。しかし、結果は依然同じ値です。

::

   >>> round(0.1, 1)
   0.10000000000000001

問題は、"0.1" を表すために記憶されている 2 進表現の浮動小数点数の値は、
すでに 1/10 に対する最良の近似になっており、値を再度丸めようとしても
これ以上ましにはならないということです。
すでに値は、 :func:`round` で得られる値になっているというわけです。

もう一つの重要なことは、0.1 が正確に 1/10 ではないため、0.1 を 10 個
加算すると厳密に 1.0 にはならないこともある、ということです。

::

   >>> sum = 0.0
   >>> for i in range(10):
   ...     sum += 0.1
   ...
   >>> sum
   0.99999999999999989

2 進の浮動小数点数に対する算術演算は、このような意外性をたくさん持っています。
"0.1" に関する問題は、以下の "表現エラー" の章で詳細に説明します。
2 進法の浮動小数点演算にともなうその他のよく知られた意外な事象に関しては
`The Perils of Floating Point <http://www.lahey.com/float.htm>`_ を
参照してください。

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

異常なケースが存在する一方で、普段の浮動小数点演算の利用では、単に最終的な
結果の値を必要な 10 進の桁数に丸めて表示するのなら、最終的には期待通りの
結果を得ることになるでしょう。
こうした操作は普通 :func:`str` で事足りますし、よりきめ細かな制御をしたければ、
:ref:`formatstrings` にある :meth:`str.format` メソッドのフォーマット仕様を
参照してください。


.. _tut-fp-error:

表現エラー
==========

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

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

::

   >>> 0.1
   0.10000000000000001

なぜこうなるのでしょうか？ 1/10 は 2 進法の分数で厳密に表現することができません。
今日 (2000年11月) のマシンは、ほとんどすべて IEEE-754 浮動小数点演算を使用して
おり、ほとんどすべてのプラットフォームでは Python の浮動小数点を IEEE-754
における "倍精度(double precision)" に対応付けます。754 の double には 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
   4503599627370496L
   >>> 2**53
   9007199254740992L
   >>> 2**56/10
   7205759403792793L

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

::

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

残りは 10 の半分以上なので、最良の近似は丸め値を一つ増やした (round up)
ものになります。

::

   >>> q+1
   7205759403792794L

従って、754 倍精度における 1/10 の取りえる最良の近似は 2\*\*56 以上の値、
もしくは

::

   7205759403792794 / 72057594037927936

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

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

::

   >>> .1 * 2**56
   7205759403792794.0

となります。

この分数に 10\*\*30 を掛ければ、有効数字 30 桁の十進数の  (切り詰められた)
値を見ることができます。

::

   >>> 7205759403792794 * 10**30 / 2**56
   100000000000000005551115123125L

これは、計算機が記憶している正確な数値が、10 進数値
0.100000000000000005551115123125 にほぼ等しいということです。この値を
有効数字 17 桁で丸めると、Python が表示する値は 0.10000000000000001 に
なります (もちろんこのような値になるのは、 IEEE 754 に適合していて、C
ライブラリで可能な限り正確に値の入出力を行った場合だけです ---
読者の計算機ではそうではないかもしれません！)

