病的な計算式?(『ソフト・エッジ』から)

4621053833中島 震・みわ よしこ『ソフト・エッジ: ソフトウェア開発の科学を求めて』(丸善出版, 2013)で、数値計算の面白い例が紹介されていました。(参考文献リストあり。索引なし。)

精度を向上させれば近似はよくなるというのが、自然な考え方でしょう。でも、この直観に当てはまらない「病的な計算式」の例があります。(p.132)

ということで、f(a, b)=(333.75-a**2)*b**6+a**2*(11*a**2*b**2-121*b**4-2)+5.5*b**8+a/(2*b)としたときの、f(77617, 33096)が紹介されています。単精度・倍精度・4倍精度の浮動小数点数での計算結果は約1.17で、正しい値(約-0.827)とは符号も違ってしまうとのことです。

これだけ読むと、精度を向上させても近似はよくならないと思うかもしれませんが、そんなことはないはずです。実際、精度をある程度高くすれば、正しく計算できます。

Pythonで試します(ideone.comでも動きます)。

from decimal import getcontext, Decimal
for n in range(16, 40):
  getcontext().prec = n
  a = Decimal("77617")
  b = Decimal("33096")
  print(n, (Decimal("333.75")-a**2)*b**6+a**2*(11*a**2*b**2-121*b**4-2)+Decimal("5.5")*b**8+a/(2*b))
16 1.000000000000000E+21
17 -4.0000000000000000E+20
18 -2.00000000000000000E+19
19 2000000000000000001
20 -99999999999999998.827
21 1.17260394005317863186
22 -999999999999998.8273961
23 100000000000001.17260394
24 10000000000001.1726039401
25 -3999999999998.827396059947
26 -99999999998.827396059946821
27 -19999999998.8273960599468214
28 -999999998.8273960599468213681
29 100000001.17260394005317863186
30 20000001.1726039400531786318588
31 -999998.8273960599468213681411651
32 300001.17260394005317863185883490
33 -9998.82739605994682136814116509548
34 -1998.827396059946821368141165095480
35 -198.82739605994682136814116509547982
36 21.1726039400531786318588349045201837
37 -0.827396059946821368141165095479816292
38 -0.8273960599468213681411650954798162920
39 -0.82739605994682136814116509547981629200

このようにPythonでは、getcontext().precの値が37以上なら適切な結果が得られます。

Rubyでも、ちょっと工夫すれば適切な結果が得られます(ideone.comで確認)

Windows付属の電卓では、何の工夫も要りません(WIndows 7で確認)。xyのショートカットキーが「Y」なので、「(333.75-77617Y2)*33096Y6+77617Y2*(11*77617Y2*33096Y2-121*33096Y4-2)+5.5*33096Y8+77617/(2*33096)=」と打てば、適切な結果が得られます。

Mac OS Xでも計算できます(10.8.4で確認。できるようになったのは比較的最近のことだと思いますが)。「計算機」で計算する方法はよくわかりませんが、Spotlightに「(333.75-77617^2)*33096^6+77617^2*(11*77617^2*33096^2-121*33096^4-2)+5.5*33096^8+77617/(2*33096)」と入力すれば、適切な結果が得られます。

Windowsの電卓は32桁までの入力にしか対応していないので、230928922463004537535516678003369の平方根を求めたいというようなときに不便ですが、Spotlightなら、sqrt(230928922463004537535516678003369)とすればとりあえず計算できます。

参考:電卓に求められるコト

本稿執筆時点のWolframAlphaでは、まったく違う結果が2つ返ってきて、片方はあっていました

Googleはだめでした

Googleの例は別として、精度を上げてもうまくいかないようにするには、もっと「病的」でなければなりません。