十進多倍長計算

Excelに =0.3-0.2-0.1 と打ち込んでみてください。当然ながら(?)ぴったりゼロになります。

では今度は =(0.3-0.2-0.1) とカッコ付きで打ち込んでみてください。あれあれ! -2.77556E-17 になってしまいました。この E-17 は ×10-17 つまり1のあとに0が17個並ぶ数 100000000000000000 分の 1 を掛け算した値という意味です。ごくわずかですが誤差があります。

Excelの計算結果としては、この -2.77556E-17 のほうが正しいのです。Excelの中では10進表現ではなく2進表現が使われているので、0.1 のような10進で表現すると有限小数になる値でも、2進で表現すると無限に循環する小数になり、正確に表せません。このため、わずかな誤差が出るのが正しい計算です。

では最初の =0.3-0.2-0.1 はなぜぴったりゼロになったのかというと、Excelはカッコを付けないと10進でキリの良い値に丸めるようにプログラムされているのです。

これは数値計算としてはまずい仕様です。本当はゼロでない答えが出ても、ゼロにされてしまうのですから。たとえばExcelに =74724506/23785549-PI() と打ち込んでみてください。ぴったりゼロになります。円周率は有理数 74724506/23785549 と(Excelの精度で)ぴったり一致するのでしょうか。でも =(74724506/23785549-PI()) とカッコ付きで打ち込めば、答えは -2.22045E-15 になります。

ちなみに、74724506/23785549-π の15桁まで正しい値は -2.22712210211597E-15 です。

以上はExcelの話ですが、AppleのNumbersという表計算ソフトの最新版では =74724506/23785549-PI() と打ち込むとぴったり -2.22712210211597E-15 になります(標準では -0.00000... となってしまって最後まで見えないので、列の幅を広げるか、セルのデータフォーマットを「指数表示」に直してください)。

どうしてNumbersはこんなに正確に答えを出せるのかというと、Pages、Numbers、Keynote の計算の精度向上について(現在はリンク切れ)に書いてあるように、128ビット10進浮動小数点(DFP)演算を使うようになったためです。一方、Excelは64ビットIEEE 754にアドホックな修正を加えたものです(Wikipediaの Numeric precision in Microsoft Excel 参照)。なお、LibreOfficeについてはLibreOffice CalcとExcelの計算結果誤差をみなさんで検証していただいた結果まとめをご覧ください。

Pythonも標準では同じことが起きます:

0.3 - 0.2 - 0.1
-2.7755575615628914e-17

Pythonでは decimal というパッケージを使えば任意精度の10進演算ができます。

from decimal import Decimal, getcontext

getcontext().prec = 100  # 100桁精度にする(デフォルトは28桁)

Decimal('0.3') - Decimal('0.2') - Decimal('0.1')
Decimal('0.0')

2進でも正確に表せる数(整数など)なら文字列として与える必要はありません:

Decimal(2).sqrt()
Decimal('1.41421356237309504880168872420969807856967187537694...') # 略

使える関数は sqrt()exp()ln()log10() くらいのものです。pi()sin()cos() の定義例は decimal のページに載っています。これを使って 74724506 / 23785549 - π の値を求めてみましょう:

from decimal import Decimal, getcontext

def pi():
    """Compute Pi to the current precision.

    >>> print(pi())
    3.141592653589793238462643383

    """
    getcontext().prec += 2  # extra digits for intermediate steps
    three = Decimal(3)      # substitute "three=3.0" for regular floats
    lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
    while s != lasts:
        lasts = s
        n, na = n+na, na+8
        d, da = d+da, da+32
        t = (t * n) / d
        s += t
    getcontext().prec -= 2
    return +s               # unary plus applies the new precision

getcontext().prec = 100  # 100桁精度にする(デフォルトは28桁)
print(Decimal(74724506) / Decimal(23785549) - pi())

結果が -2.22712210211597472...E-15 のように出力されます。

ちなみに、74724506 / 23785549 がほぼ π になることは次のようにして探索しました:

import math

a = 1
for i in range(1, 100000000):
    j = round(i * math.pi)
    x = j / i - math.pi
    if abs(x) < a:
        a = abs(x)
        print(j, i, a)
        if a == 0:
            break

ループを使わないで1から100までの整数を列挙できるかという問題もやってみましょう(How can I print 1 to 100 in C++ without a loop, goto or recursion? にある Conner Davis の方法):

getcontext().prec = 298

1000 / Decimal(999)**2
Decimal('0.001002003004005006007008009010011012...(中略)...098099100')

Wikipediaの精度保証付き数値計算に載っているRumpの例題をやってみましょう。まず普通のコードです(行末の \ は行が続く意味です):

def f(a, b):
    return 333.75 * b**6 + a**2 * (11 * a**2 * b**2 - b**6 - 121 * b**4 - 2) \
        + 5.5 * b**8 + a / (2 * b);

これで f(77617, 33096)f(77617.0, 33096.0) を表示してみると、まったくわけのわからない値になります。そこで

getcontext().prec = 100

def decimalf(a, b):
    return Decimal('333.75') * b**6 + a**2 * (11 * a**2 * b**2 - b**6 - 121 * b**4 - 2) \
        + Decimal('5.5') * b**8 + a / (2 * b);

として decimalf(Decimal('77617'), Decimal('33096')) を求めると、正しい -0.827396... になります。

(おまけ)1 / 0.9899 を計算して、小数点以下を2桁ずつ区切って読んでみてください。フィボナッチ数列になっています(10000/9899 を計算するとなぜ商 (1.0102030508) にフィボナッチ数列が現れるのですか?