Rumpの例題とは、
\[ f = 333.75b^6 + a^2(11a^2b^2 - b^6 - 121b^4 - 2) + 5.5b^8 + \frac{a}{2b} \]を $a = 77617$、$b = 33096$ のときに計算する問題です（Rump's Example Revisited）。
とりあえずやってみます：
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)
1.1805916207174113e+21
f(77617.0, 33096.0)
-1.1805916207174113e+21
とんでもなく大きな値で、しかも引数を整数で与えたときと浮動小数点数で与えたときで正負が逆になってしまいました。どうやらとても難しい計算のようです。ちなみに正解は
\[ f = -0.827396059946821368141165095479816\ldots \]です。
多倍長演算ライブラリ mpmath を使って、何桁の精度なら正しい答えが得られるか、調べてみましょう：
from mpmath import mp
def f(a, b):
return (mp.mpf("333.75") * b**6 + a**2 * (11 * a**2 * b**2 - b**6 - 121 * b**4 - 2)
+ mp.mpf("5.5") * b**8 + a / (2 * b))
for mp.dps in range(30, 51):
print(mp.dps, f(mp.mpf("77617"), mp.mpf("33096")))
30 1.17260394005317863185883490452 31 1.17260394005317863185883490452 32 1.1726039400531786318588349045202 33 1.17260394005317863185883490452018 34 1.172603940053178631858834904520184 35 1.1726039400531786318588349045201837 36 -0.827396059946821368141165095479816292 37 -0.827396059946821368141165095479816292 38 -0.827396059946821368141165095479816292 39 -0.827396059946821368141165095479816291999 40 -0.827396059946821368141165095479816291999 41 -0.82739605994682136814116509547981629199903 42 -0.827396059946821368141165095479816291999033 43 -0.8273960599468213681411650954798162919990331 44 -0.82739605994682136814116509547981629199903312 45 -0.827396059946821368141165095479816291999033116 46 -0.8273960599468213681411650954798162919990331158 47 -0.82739605994682136814116509547981629199903311579 48 -0.827396059946821368141165095479816291999033115784 49 -0.8273960599468213681411650954798162919990331157844 50 -0.82739605994682136814116509547981629199903311578438
35桁まで試せば 1.1726… になりそうなのですが、その後36桁で急に値が変わるのがおもしろいところです。
333.75 と 5.5 はたまたま2進で正確な値になるので mp.mpf("...") の形にしないでも大丈夫のようです。
ところで、多倍長演算を使わなくても Python は整数の計算は正しいので、4倍して4で割る変形をすると、
def f(a, b):
return (1335 * b**6 + 4 * a**2 * (11 * a**2 * b**2 - b**6 - 121 * b**4 - 2)
+ 22 * b**8 + 2 * a / b) / 4
f(77617, 33096)
-0.8273960599468213
完全に正しい答えがえられました。
実は、割り算 / だけは浮動小数点演算になるのですが、最後の / 4 は全体を割るので精度に影響しませんし、これを除けば小数点以下が出るのは 2 * a / b だけですので、これが桁落ちを生じるような大きな値でないので、うまくいったのです。
Last modified: