数値計算のリスク

精度保証付き数値計算についてのQiitaのAdvent Calendarの記事が話題になっていた。数値計算はなかなか奥が深いので困る。

[Rumpの例題]
$f(a,b) = 333.75 b^6 + a^2 (11 a^2 b^2 − b^6−121 b^4−2) + 5.5 b^8 +\dfrac{a}{2b}$
に $(a,b)=(77617.0, 33096.0)$ を代入した値は？

Mathematicaの場合

In[1]:= f[a_, b_] :=  333.75 b^6 + a^2 (11 a^2 b^2 - b^6 - 121 b^4 - 2) + 5.5 b^8 + a/(2 b)
In[2]:= f[77617.0, 33096.0]
Out[2]= -1.18059*10^21

In[3]:= g[a_, b_] :=  1335/4 b^6 + a^2 (11 a^2 b^2 - b^6 - 121 b^4 - 2) + 11/2 b^8 + a/(2 b)
In[4]:= N[g[77617, 33096], 20]
Out[4]= -0.82739605994682136814

Juliaの場合

[1] function g(a::BigInt,b::BigInt)
    x = 1335/4*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 11/2*b^8 + a/(2*b)
end
[2] g(BigInt(77617), BigInt(33096))
[2] -0.8273960599468213681411650954798162919990331157843848199178148416727096930142628

[3] function g128(a::Int128,b::Int128)
   x = 1335/4*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 11/2*b^8 + a/(2*b)
end
[4] g128(Int128(77617), Int128(33096))
[4] 1.1805916207174113e21

[5] function f(a::BigFloat,b::BigFloat)
    x = 1335/4*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 11/2*b^8 + a/(2*b)
end
[6] f(BigFloat(77617), BigFloat(33096))
[6] -0.8273960599468213681411650954798162919990331157843848199178148416727096930142628

[7] function f64(a::Float64,b::Float64)
    x = 1335/4*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 11/2*b^8 + a/(2*b)
end
[8] f(Float64(77617), Float64(33096))
[8] -1.1805916207174113e21