演算誤差と戦う
浮動小数点を使ったときの誤差についてはよく話題になる。 https://togetter.com/li/1289806
「0.01を10000回足したら100.003になる」といわれているが、実際にそうなるのか試してみよう。python3で0.01を10000回足す処理を書いてみる。
sum = 0 for i in range(10000): sum += 0.01 print(sum) print(sum == 100)
結果はこうなった。
100.00000000001425 False
100.003とは違うが、とにかく誤差は発生している。誤差が違うのは、pythonの浮動小数点がdouble型だからであろう。
計算機で整数以外を扱うときには数値表現の誤差・演算誤差に気をつけないといけない。 演算誤差を避けるための工夫について、pythonモジュールを使いながら見ていきたい。
decimal型
decimal型はいわゆる10進数の小数型である。 decimal型を用いれば、先ほどの計算で誤差が発生しない。 実際にやってみる。
pythonでは標準モジュールにdecimalが入っているのでそれを使う。 https://docs.python.org/3.6/library/decimal.html
from decimal import Decimal sum = Decimal(0) for i in range(10000): sum += Decimal("0.01") print(sum) print(sum == 100)
100.00 True
0.01を10000回足すと、100と一致することが確認できた。
decimalを使えばあらゆる誤差をさけられるかというと、そうでもない。 例えば以下のような例を考えてみる。
a = Decimal("1") / Decimal("3") b = Decimal("3") c = a * b print(a) print(c) print(c == 1)
0.3333333333333333333333333333 0.9999999999999999999999999999 False
(1/3)*3 を計算すると 0.9999999999999999999999999999 になってしまった。 有限桁の10進小数で 1/3 を正確に表現することは不可能なのである。 「3進法ライブラリ」のようなものを使えばたぶんいけるだろうが、もう少し汎用的にやりたい。
有理数
https://ja.wikipedia.org/wiki/%E6%9C%89%E7%90%86%E6%95%B0
有理数(ゆうりすう、英: rational number) とは、二つの整数 a, b (ただし b は 0 でない)をもちいて a/b という分数で表せる数のことをいう。
整数の分数は有理数である。 分母と分子を別々の整数として持っておけば、有理数の+-×÷の演算を整数の+-×で実行でき、誤差なく演算ができる(通分や約分の処理が厄介だが)。 (1/3)*3も正確に1としてくれるであろう。
pythonでは、有理数計算を行うfractionsが標準ライブラリに含まれる。 https://docs.python.org/3.6/library/fractions.html 実際に使ってみよう。
from fractions import Fraction a = Fraction(1, 3) # 2つの引数は (分子,分母)を表す b = Fraction(3) # 引数を1つだけ入力すると、自動で整数になる c = a * b print(a) print(c) print(c == 1)
1/3 1 True
(1/3)*3 が正確に 1 と一致することが確かめられた。また、print(a)が "1/3" という表示になっており、分子と分母を分けて保持していることが推察される。
ちなみに、Fractionは約分機能も持っている。
a = Fraction(4, 6) b = Fraction(2, 3) print(a) print(a == b)
2/3 True
4/6を2/3に約分してくれている。
さて、有理数が正確に計算できた。次はもちろん無理数へと進んでいく。
代数的数
無理数の中でも比較的扱いやすい、代数的数が誤差なく計算できるかどうか検証する。なお、この記事では実数のみを扱う。
https://ja.wikipedia.org/wiki/%E4%BB%A3%E6%95%B0%E7%9A%84%E6%95%B0
代数的数(だいすうてきすう、英: algebraic number)とは、 複素数であって、有理数係数(あるいは同じことだが、分母を払って、 整数係数)の 0 でない一変数多項式の根 (すなわち多項式の値が 0 になるような値)となるものをいう。 すべての整数や有理数は代数的数であり、またすべての整数の冪根も代数的数である。
まずは、Fractionでは が成り立たないことを確認する。
import math a = math.sqrt(Fraction(2)) b = a * a print(a) print(b) print(b == 2)
1.4142135623730951 2.0000000000000004 False
当然、有理数の範囲では√2は扱えない。そもそも math.sqrt() は平方根を近似値で求める関数であろう。
pythonで無理数を含めた代数計算を行うライブラリとして SymPy がある。 https://www.sympy.org/en/index.html
$ pip install sympy
さて、sympyでは が計算できるだろうか。
import sympy a = sympy.sqrt(2) b = a * a print(a) print(b) print(b == 2)
sqrt(2) 2 True
正しく計算できた。 print(a)が "sqrt(2)" という表示になっていることから「平方根として持っている」のだと思われる。
sympyの代数的数で遊んでみる。
- としたとき、が成り立つことを確認。
x = sympy.sqrt(2) + sympy.sqrt(3) y = x * x * x * x - 10 * x * x + 1 print(x) print(y) print(y.equals(0)) # 複雑な式のときは .equals()メソッドで等値判定をする
sqrt(2) + sqrt(3) -(sqrt(2) + sqrt(3))*(10*sqrt(2) + 10*sqrt(3)) + 1 + (sqrt(2) + sqrt(3))**4 True
- としたとき、が成り立つことを確認。
x = sympy.AlgebraicNumber("2**(2/3)") y = x * x * x print(float(x)) print(y) print(y.equals(4))
1.5874010519681994 2**(2/3)**3 True
楽しい。
超越数
https://ja.wikipedia.org/wiki/%E8%B6%85%E8%B6%8A%E6%95%B0
超越数(ちょうえつすう、英: transcendental number)とは、代数的数でない数、すなわちどんな有理係数の代数方程式 x^{n}+a_{n-1}x^{n-1} + ... + a_0 = 0 (n ≥ 1 かつ、各 ai は有理数)の解にもならないような複素数のことである。...よく知られた超越数にネイピア数(自然対数の底)や円周率がある。
さて最後は超越数である。超越数の計算誤差を実用上問題にする人はいるのか...? とは思うものの、興味本位で検証してみる。
高校数学で出てくる という超越数πを含む等式がある。 まずは浮動小数では上手くいかないことを確認する。
a = math.pi / 6 b = math.sin(a) print(a) print(b) print(b == 0.5) # 0.5は浮動小数の範囲で誤差なく表現できる
0.5235987755982988 0.49999999999999994 False
0.49999999999999994になってしまった。
sympyはこの等式を上手く扱えるだろうか。
a = sympy.pi / 6 b = sympy.sin(a) print(a) print(b) print(b.equals(0.5))
pi/6 1/2 True
すごい!計算機は無限精度の計算が可能だった!
(実際のところ、わけがわからない。特殊処理を入れているんだろうか?)
結論
sympyを使えば誤差で困ることはまずないと思われる(ただし、パフォーマンスは考えないものとする)