1. float, Decimal の有効桁数を確認しておく
これまでに 「float よりも精度の高い Decimal」 、「Decimal を文字列として出力するときの str() と repr() の違い」について見てきた。
- Python で正確な小数の計算 (1) - Decimal モジュールを使う
- Python で正確な小数の計算 (2) - Decimal の値を文字列にして、eval () で評価するときは repr 関数を使う
今回は、float と Decimal 有効桁数、str() と repr() の有効桁数について試しておく。こういうのは、一度自分で動作を確認しておかないと気持ちが悪い。 (+_+)
2. 事前準備
毎回 print 文で、式と結果を書くのは面倒。予め、与えられた式の文字列を受けとると、str() と repr() の結果を表示する関数を作成しておいた。
def p(E):
print E, " = ", str(eval(E)), " : ", repr(eval(E))
print 文は、str() を呼出していたと思うけれど、念のため…。
3. 有効桁数を確認する
str(), repr() の丸め
最初に、思い出しておくべき、基本事項は、
残念なことに、ほとんどの小数は 2 進法の分数として正確に表わすことができません。その結果、一般に、入力した 10 進の浮動小数点数は、 2 進法の浮動小数点数で近似された後、実際にマシンに記憶されます。
(B. 浮動小数点演算、その問題と制限 より)
float で小数は近似値であるということ。
str() と repr() についての違いを確認しておく。、
str() 関数は有効数字 12 桁しか生成しません
repr(float)は真の 10 進値を有効数字 17 桁で丸め（B. 浮動小数点演算、その問題と制限 より）
「真の」値を丸めているということを忘れずに。
ハードウェアが浮動小数点数を記憶するのに用いているビット数がマシンによって異なり、Python は単にマシンに 2 進で記憶されている、真の 10 進の値を近似した値を、されに 10 進で近似して出力する
(同上より)
念のため原文で確認。
because the number of bits used by the hardware to store floating-point values can vary across machines, and Python only prints a decimal approximation to the true decimal value of the binary approximation stored by the machine.
(B. Floating Point Arithmetic: Issues and Limitations より、太字は引用者による。)
以上を整理して、図にすると、
float
境界を試してみる。
p('0.111111111019 ') # 12 桁 ここまで
p('0.1111111110119') # 13 桁
p('0.11111111101111119 ') # 17 桁 ここまで
p('0.111111111011111119') # 18 桁
結果は、
0.111111111019 = 0.111111111019 : 0.11111111101899999 0.1111111110119 = 0.111111111012 : 0.11111111101189999
0.11111111101111119 = 0.111111111011 : 0.11111111101111119 0.111111111011111119 = 0.111111111011 : 0.11111111101111112
なぜ 17 桁なのかと言うと、
repr(float)が有効数字 17桁 の値を生成するのは、この値が (ほとんどのマシン上で) 、全ての有限の浮動小数点数 x についてeval(repr(x)) == xが成り立つのに十分で、かつ有効数字 16 桁に丸めると成り立たないからです。これは 2 進法の浮動小数点の性質です: Python のバグでも、ソースコードのバグでもなく、浮動小数点演算を扱えるハードウェア上の、すべての言語で同じ類の現象が発生します
(同上より、太字は引用者による)
Decimal
Decimal は文字列として与えると、演算が行われるまでそのまま保持しているようなので 1 を掛けて結果を見た。
p('Decimal("0.1111111110111111111011111119") * 1') # 28 桁 ここまで
p('Decimal("0.11111111101111111110111111119") * 1') # 29 桁
結果は、
Decimal("0.1111111110111111111011111119") * 1 = 0.1111111110111111111011111119 : Decimal("0.1111111110111111111011111119")
Decimal("0.11111111101111111110111111119") * 1 = 0.1111111110111111111011111112 : Decimal("0.1111111110111111111011111112")
4. 今回ハマったところ
1 / 3 * 3
を float と Decimal で試してみると、
p('1. / 3. * 3.')
p('1. / 3. * 3. == 1')
p('Decimal(1) / Decimal(3) * Decimal(3)')
p('Decimal(1) / Decimal(3) * Decimal(3) == 1')
結果は、
1. / 3. * 3. = 1.0 : 1.0 1. / 3. * 3. == 1 = True : True
Decimal(1) / Decimal(3) * Decimal(3) = 0.9999999999999999999999999999 : Decimal("0.9999999999999999999999999999")
Decimal(1) / Decimal(3) * Decimal(3) == 1 = False : False
表示にイコール = があるので見にくいけれど… ^^;
上記では、一見正確でないはずの float の方が正確な答えを出していて、精度が高いと思っていた Decimal が間違っているように見える。なぜだろう。 (@_@;) ？
数
まずは、「数」について、
整数は０と０から＋１または、－１ずつしていき得られる数のこと。負の整数、０、正の整数があります。→-3,0,7など
自然数は整数のうち、正の整数のものをさします。→1,7,13など
有理数は分数で現されるものをさします。いい替えると、整数および有限小数と無限循環小数です。→0.2,1/3,0.125など
無理数は分数で現されない小数です。つまり、循環しない無限小数のことです。→円周率（3.14…），√2(1.414…)
実数は虚数でない数。つまり無理数、有理数、整数、自然数をすべて含めた数です。(実数、無理数、有理数、整数、自然数の関係を分かりやすく教えてください！！ それ... - Yahoo!知恵袋 より、太字は引用者による)
念のためもう一つの説明を、
自然数【1、2、・・・】 ⇒ ものの個数など、数の基礎。
↓
整数【・・・、-2、-1、0、1、2、・・・】 ⇒ マイナスという概念とゼロという概念が加わったもの。
↓
有理数【p/q】 ⇒ 割り算（分数）という概念によるもの。2/3、-3/4など。
↓
無理数 ⇒ 有理数では表せないもの。3.14159…など。
↓
実数 ⇒ 上記の有理数（自然数や整数を含む）・無理数を合体したもの。(数学Aの事で質問です 素数、自然数、有理数、無理数、実数、整数の違いがごちゃ... - Yahoo!知恵袋 より、太字は引用者による)
つまり、 1 / 3 = 0.3333333333… とずっとつづく無限小数。
float による計算
では、float による計算の過程を追うと、
1. / 3. = 0.333333333333 : 0.33333333333333331
上記のように、str() による出力は13 桁目が丸められ、 repr() の出力は 18 桁目が丸められる。その後ろがどうなっているかは表示に表われない。
次に、正確ではないけれど (真の値ではないという意味) 、上記の答えに 2 を掛けたのと同じ状態にしてやると、
0.33333333333333331 + 0.33333333333333331 = 0.666666666667 : 0.66666666666666663
str() の出力は、13 桁目が丸められて 12 桁目が 7 となる。 repr() の出力は、最後が 31 + 31 で 62 となると思いきや 63 。これは、31 以降にあると考えらえる真の値の数が丸められて 3 になったということなのかな。
3 を掛けたときの真の値がどのようなものかわからないけれど、少しずつ 17 桁の値を変化させて見てみると、
0.99999999999999990 = 1.0 : 0.99999999999999989 0.99999999999999991 = 1.0 : 0.99999999999999989 0.99999999999999992 = 1.0 : 0.99999999999999989 0.99999999999999993 = 1.0 : 0.99999999999999989 0.99999999999999994 = 1.0 : 0.99999999999999989 0.99999999999999995 = 1.0 : 1.0 0.99999999999999996 = 1.0 : 1.0 0.99999999999999997 = 1.0 : 1.0 0.99999999999999998 = 1.0 : 1.0
更にもう少し桁数を多くして見ていくと、
0.99999999999999994454010 = 1.0 : 1.0 0.99999999999999994454011 = 1.0 : 0.99999999999999989 0.99999999999999994454012 = 1.0 : 0.99999999999999989 0.99999999999999994454013 = 1.0 : 0.99999999999999989 0.99999999999999994454014 = 1.0 : 0.99999999999999989 0.99999999999999994454015 = 1.0 : 0.99999999999999989 0.99999999999999994454016 = 1.0 : 1.0 0.99999999999999994454017 = 1.0 : 1.0 0.99999999999999994454018 = 1.0 : 1.0 0.99999999999999994454019 = 1.0 : 1.0
あれ？
0.99999999999999994454010 < 0.99999999999999994454011 = False : False 0.99999999999999994454015 < 0.99999999999999994454016 = True : True
こんなこともあるのかぁ … (@_@;)
それは横に置いておいて、上記を見ると 1.0 に丸めらる境界の辺りが何となくわかる。 1 / 3 の結果が丸められた値に 3 に極めて近い数値を掛けてみると、
0.33333333333333331 * 2.99999999999999977 = 1.0 : 0.99999999999999978 0.33333333333333331 * 2.99999999999999978 = 1.0 : 1.0
どうやらこの辺りが境界のようだ。
結局、 float の方は、値が丸められたことによって、 1 / 3 * 3 の結果が 1 になったということ。
Decimal による計算過程
まずは 1 / 3 の値は、
Decimal(1) / Decimal(3) = 0.3333333333333333333333333333 : Decimal("0.3333333333333333333333333333")
28 桁でピタッと丸められている。これに 3 を掛けるで、
Decimal(1) / Decimal(3) * Decimal(3) = 0.9999999999999999999999999999 : Decimal("0.9999999999999999999999999999")
ちょうどその 3 倍の値が正確に返ってくる。
5.5 decimal -- 10進浮動小数点数の算術演算 によると、
ハードウェアによる 2 進浮動小数点表現と違い、decimal モジュールでは計算精度をユーザが指定できます(デフォルトでは 28 桁です)。
デフォルトの設定のまま 29 桁の 9 を並べると、丸められて 1 となる。
Decimal("0.99999999999999999999999999999") * 1 = 1.000000000000000000000000000 : Decimal("1.000000000000000000000000000")
上記で試した float において 1 に丸められてしまった値も 1 になることはない。
Decimal("0.99999999999999994454016") * 1 = 0.99999999999999994454016 : Decimal("0.99999999999999994454016")
できるだけ誤差を少なく
結局、もし Decimal を使う場合で、計算結果が循環小数になる可能性があるなら、なるべく誤差を少なくするためには、乗算を除算よりも先にしておく。
Decimal(1) * Decimal(3) / Decimal(3) = 1 : Decimal("1")
5. 正確に計算するには
2008-09-09 - ひとり勉強会 には、RealLib を使う方法が紹介されている。
また、192:有理数を使う には、
Python に有理数型はありませんが、サンプルソースコード Rat.py にはクラスの拡張方法の例として有理数の実装が収められています。
そして、PEP 239 -- Adding a Rational Type to Python には、
Python has no numeric type with the semantics of an unboundedly
precise rational number. This proposal explains the semantics of
such a type, and suggests builtin functions and literals to
support such a type.
早く組込み関数で使えるようになるといいなぁ～。
0コメント:
コメントを投稿