blog2019-12-26Rumpの例題のやらしさqiita.comこの記事を読んで,ちょっと計算してみようかとRumpの例題を計算し始めてしまったのが運の尽き。真値が計算できないと泣きそうになって夜を明かす羽目に。Rumpの例題はこのページにもあるけどで,を計算すると言うだけの話。これをただ単にrubyに打ち込むと一見して巨大数と巨大数を足したあとにa/2bという小さい数を足してるから失敗するだろうなというのはわかる。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の巨大な数字が残る。ここまで清々しく違うと逆に異常だと分かるので親切かもしれない。1970年代から連綿と伝わる本当のプログラマの手法で解決しようと考える。巨大数と巨大数の引き算を解析的な変形により除去することを考える。あることに気づいた上で,斉次の項をまとめるのを優先して変形すると以下のように変形できる事に気づく。この内側のがぴったり0になるように見えるのがトラップで,別に巨大な数が出てきてないから問題ないじゃないかと,楽観的に思うと,という謎の数字にたどり着いてしまう。333.75+5.5*b**2-a**2+11*(a/b)**2*((a/b)**2-11);'0.0ただ,これを細かく見てみるとa/bが絡む項は整数でないわけで,きっちり0になるのは不自然だと気づく。かなりの巨大数なのでものすごい小さい数(それこそ,10^-13とかでも)が残っても最終結果を崩壊させるくらいの巨大な数として残ってくるので,もっと精緻な取り扱いをしないと死ぬ。いうことで,色々と諦めてこれをすべて有理数として多倍長整数で計算することを考える。整数の範囲で成立する式だから,これを前提として式を変形すると,以下のようになる。この1項目の分数の分子を構成する項をそれぞれ多倍長整数で計算すると以下のようになる。1749166305248087785856979173376-1749166305248087785808783983872-48195189512上で足したらピッタリ0になると喜んだのがこの最初の二つの項の和を4b^6で割った数字で,この二つは21桁目まで全部同じという悪魔的な数字なので,普通にfloatで足したらそりゃピッタリ0になって地獄を見る数字であることがわかる。a/2.0/b'-0.8273960599468213となる。Rumpさん悪魔的にすごいというのが計算してみた感想。やったら,一見して同じに見えないに22桁目まで同じ値になる式を作れるのか全くわからない。12-14)月別アーカイブはてなブログをはじめよう!khwarizmiさんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。