演算誤差と戦う

浮動小数点を使ったときの誤差についてはよく話題になる。 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では \sqrt{2}^2 = 2 が成り立たないことを確認する。

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では \sqrt{2}^2 = 2 が計算できるだろうか。

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=\sqrt{2} + \sqrt{3}としたとき、x^{4} - 10x^{2} + 1 = 0が成り立つことを確認。
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=2^{\frac{2}{3}}としたとき、x^{3} = 4が成り立つことを確認。
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 は有理数)の解にもならないような複素数のことである。...よく知られた超越数ネイピア数自然対数の底)や円周率がある。

さて最後は超越数である。超越数の計算誤差を実用上問題にする人はいるのか...? とは思うものの、興味本位で検証してみる。

高校数学で出てくる sin(\frac{\pi}{6}) = \frac{1}{2} という超越数πを含む等式がある。 まずは浮動小数では上手くいかないことを確認する。

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を使えば誤差で困ることはまずないと思われる(ただし、パフォーマンスは考えないものとする)