Rumpの例題のやらしさ

qiita.com
この記事を読んで,ちょっと計算してみようかとRumpの例題を計算し始めてしまったのが運の尽き。全然真値が計算できないと泣きそうになって夜を明かす羽目に。

Rumpの例題はこのページにもあるけど
a=77617.0, b=33096.0で,
f(a,b) = 333.75 b^{6}+a^{2}(11a^{2}b ^{2}-b^{6}-121b^{4}-2)+5.5b^{8}+ \displaystyle \frac{a}{2b}
を計算すると言うだけの話。

これをただ単にrubyに打ち込むと一見して巨大数と巨大数を足したあとにa/2bという小さい数を足してるから失敗するだろうなというのはわかる。ただ,その直感と違って,ruby(2.6.3p62)にこの計算をそのまま打ち込んでさせると以下の結果になる

$ ruby -e 'a=77617.0;b=33096; puts 333.75*b**6+a**2*(11*a**2*b**2-b**6-121*b**4-2)+5.5*b**8+a/2/b'
-1.1805916207174113e+21

元々の論文(https://tore.tuhh.de/bitstream/11420/318/1/Ru88a.pdf)や,直感的な狙い(手前の3項の足し算が打ち消し合って0になるので,a/2bだけが残る)と異なり,10^21の巨大な数字が残る。ここまで清々しく違うと逆に異常だと分かるので親切かもしれない。*1

で,これを高等で洗練された精度保証付き計算の技法ではなく, 1970年代から連綿と伝わる本当のプログラマの手法で解決しようと考える。まずは巨大数と巨大数の引き算を解析的な変形により除去することを考える。
特に,5.5b^2-a^2 = -1.0であることに気づいた上で,斉次の項をまとめるのを優先して変形すると以下のように変形できる事に気づく。
\left[ 333.75 + 5.5b^2-a^2 +  11  ( \frac{a}{b} )^2 \left(  ( \frac{a}{b} )^2-11 \right) \right] b^6-2a^2 +  \frac{a}{2b}
この内側の
 333.75 + 5.5b^2-a^2 +  11  ( \frac{a}{b} )^2 \left(  ( \frac{a}{b} )^2-11 \right)がぴったり0になるように見えるのがトラップで,別に巨大な数が出てきてないから問題ないじゃないかと,楽観的に思うと,
f(a,b) = -12048797376.827396
という謎の数字にたどり着いてしまう。
実際にrubyで単純に計算するとそうなる

$ ruby -e 'a=77617.0;b=33096.0; puts 333.75+5.5*b**2-a**2+11*(a/b)**2*((a/b)**2-11);'
0.0

ただ,これを細かく見てみるとa/bが絡む項は整数でないわけで,きっちり0になるのは不自然だと気づく。b^6はかなりの巨大数なのでものすごい小さい数(それこそ,10^-13とかでも)が残っても最終結果を崩壊させるくらいの巨大な数として残ってくるので,もっと精緻な取り扱いをしないと死ぬ。

ということで,色々と諦めてこれをすべて有理数として多倍長整数で計算することを考える。5.5b^2-a^2 = -1は整数の範囲で成立する式だから,これを前提として式を変形すると,以下のようになる。
 f(a,b) = \displaystyle \frac{(1335-4)b^6 + 44a^4b^2 - 4\cdot 121 a^2b^4 - 8a^2}{4} + \frac{a}{2b}

この1項目の分数の分子を構成する項をそれぞれ多倍長整数で計算すると以下のようになる。

$ ruby -e 'a=77617;b=33096; puts (1335-4)*b**6; puts 44*a**4*b**2 - 4*121*a**2*b**4; puts -8*a**2'
 1749166305248087785856979173376
-1749166305248087785808783983872
-48195189512

上で足したらピッタリ0になると喜んだのがこの最初の二つの項の和を4b^6で割った数字で,この二つは21桁目まで全部同じという悪魔的な数字なので,普通にfloatで足したらそりゃピッタリ0になって地獄を見る数字であることがわかる。
ということで,分子を計算すると,

$ ruby -e 'a=77617;b=33096; puts ((1335-4)*b**6+44*a**4*b**2 - 4*121*a**2*b**4-8*a**2)'
-8

なので,結果,

$ ruby -e 'a=77617;b=33096; puts -8/4.0 + a/2.0/b'
-0.8273960599468213

f(a,b) = \displaystyle \frac{-8}{4}+ \frac{a}{2b} = -0.827396
となる。

Rumpさん悪魔的にすごいというのが計算してみた感想。どうやったら,一見して同じに見えないに22桁目まで同じ値になる式を作れるのか全くわからない。

*1:ちなみに,python3でも似たような結果が出てきた