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

浮動小数点数は、コンピューターのハードウェア上では 2 進小数で表されま
す。 例えば、 **10 進小数** の "0.625" は 6/10 + 2/100 + 5/1000 という
値を持ち、同様に **2 進小数** の "0.101" は 1/2 + 0/4 + 1/8 という値を
持ちます。 この二つの小数は同じ値を表しています。違いは、前者が 10 進
法、後者が 2 進法の小数表記で書かれていることだけです。

残念ながら、ほとんどの 10 進小数は 2 進小数として正確に表せません。そ
のため、入力した 10 進の浮動小数点数に対して、実際にコンピューターに格
納される 2 進の浮動小数点数は、一般に近似値にすぎません。

まず 10 進法で考えると、この問題を理解しやすくなります。分数 1/3 を考
えてみましょう。これは 10 進小数で次のように近似できます:

   0.3

さらに正確な近似は、

   0.33

さらに正確な近似は、

   0.333

というように続きます。何桁まで書いても、結果が正確に 1/3 になることは
ありませんが、桁数を増やすほど 1/3 に近い近似値になります。

同様に、2 進法で何桁使っても、10 進数の 0.1 を 2 進小数として正確に表
すことはできません。2 進法では、1/10 は次のように無限に繰り返す循環小
数になります:

   0.0001100110011001100110011001100110011001100110011...

有限のビット数で打ち切ると、得られるのは近似値です。現在のほとんどのコ
ンピューターでは、浮動小数点数は、最上位ビットから数えて最初の 53 ビッ
トを分子に使い、2 の冪乗を分母とする分数で近似されます。1/10 の場合、
この分数は "3602879701896397 / 2 ** 55" です。これは 1/10 に近い値です
が、正確に等しいわけではありません。

値の表示方式によって、ほとんどのユーザーはこの近似に気づきません。
Python が表示するのは、コンピューターに格納された 2 進の近似値を正確に
10 進数で表した値ではなく、それをさらに 10 進数で近似した値です。ほと
んどのコンピューターで、0.1 に対して格納された 2 進の近似値を Python
が正確に 10 進数で表示するとすれば、次のように表示する必要があります:

   >>> 0.1
   0.1000000000000000055511151231257827021181583404541015625

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

   >>> 1 / 10
   0.1

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

興味深いことに、異なる多くの 10 進数が、最も近い近似値として同じ 2 進
小数に対応します。例えば、"0.1" と "0.10000000000000001" と
"0.1000000000000000055511151231257827021181583404541015625" はどれも
"3602879701896397 / 2 ** 55" で近似されます。これらの 10 進数は同じ近
似値に対応するため、どれを選んで表示しても、 "eval(repr(x)) == x" とい
う条件を保てます。

以前は、Python のプロンプトと組み込み関数 "repr()" は、有効数字 17 桁
の "0.10000000000000001" を選んで表示していました。Python 3.1 以降は、
ほとんどのシステムで、これらのうち最も短い表現を選び、単に "0.1" と表
示できるようになりました。

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

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

   >>> format(math.pi, '.12g')  # 有効数字 12 桁で表示
   '3.14159265359'

   >>> format(math.pi, '.2f')   # 小数点以下 2 桁で表示
   '3.14'

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

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

一つの錯覚が別の錯覚を生むこともあります。例えば、0.1 は正確には 1/10
ではないため、0.1 を三つ足した結果も正確に 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 の浮動小数点演算
の誤差は浮動小数点演算を行うハードウェアに由来し、ほとんどのコンピュー
ターでは演算ごとに高々 2**53 分の 1 程度です。これはほとんどの用途で十
分な精度です。ただし、10 進の演算ではないこと、そして浮動小数点演算の
たびに新たな丸め誤差が生じる可能性があることを覚えておく必要があります
。

極端なケースはあるものの、普段の浮動小数点演算では最終結果の表示を必要
な 10 進の桁数に丸めるだけで、たいていは期待する結果が得られます。通常
は "str()" で十分ですが、より細かく制御するには、 書式文字列の構文 に
ある "str.format()" メソッドの書式指定子を参照してください。

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

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

浮動小数点演算を頻繁に使う場合は、NumPy パッケージや、SciPy プロジェク
トが提供するその他の数学・統計演算用の多くのパッケージを調べてみてくだ
さい。 <https://scipy.org> を参照してください。

浮動小数点数の正確な値を *本当に* 知りたいというまれな場合に役立つツー
ルも、Python は提供しています。 "float.as_integer_ratio()" メソッドは
、浮動小数点数の値を分数として表します:

   >>> 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

"math.fsum()" はさらに進んで、累計に値を加える際に「失われた桁」を全て
追跡し、結果を得るまでの丸めを一度だけに抑えます。これは "sum()" より
低速ですが、絶対値の大きい入力値が互いにほぼ打ち消し合い、最終的な合計
がゼロに近くなるようなまれな場合には、より正確な結果が得られます:

   >>> arr = [-0.10430216751806065, -266310978.67179024, 143401161448607.16,
   ...        -143401161400469.7, 266262841.31058735, -0.003244936839808227]
   >>> float(sum(map(Fraction, arr)))   # 正確に合計し、一度だけ丸める
   8.042173697819788e-13
   >>> math.fsum(arr)                   # 一度だけ丸める
   8.042173697819788e-13
   >>> sum(arr)                         # 拡張精度で複数回丸める
   8.042178034628478e-13
   >>> total = 0.0
   >>> for x in arr:
   ...     total += x                   # 標準精度で複数回丸める
   ...
   >>> total                            # 単純な加算では正しい桁が一つもない！
   -0.0051575902860057365


15.1. 表現誤差
==============

この節では「0.1」の例を詳しく説明し、このような場合に正確な分析を自分
で行う方法を示します。 2 進浮動小数点表現についての基本的な知識がある
ことを前提としています。

表現誤差 (*Representation error*) とは、一部の (実際にはほとんどの) 10
進小数を 2 進小数として正確に表せないことを指します。これが、Python (
あるいは Perl、C、C++、Java、Fortran など多くの言語) で、期待した正確
な 10 進数が表示されないことがよくある主な理由です。

なぜそうなるのでしょうか。1/10 は 2 進小数として正確に表せません。少な
くとも 2000 年以降、ほぼ全てのコンピューターは IEEE 754 の 2 進浮動小
数点演算を使い、ほぼ全てのプラットフォームで Python の浮動小数点数は
IEEE 754 binary64 の「倍精度」値に対応しています。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

つまり、 *J* がちょうど 53 ビットになる *N* の値は 56 だけです。このと
き、 *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 の倍精度で表せる、1/10 の最もよ
い近似値です:

   >>> 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'
