2011年6月13日月曜日

Intel Fortranを使うもう一つの理由:qatan

現在、かなりシビアな精度での数値計算を実行中。零割が頻発して、計算が止まる事象が続発。問題は差分計算。微分公式 limd→0 ( f(x+d) - f(x) ) / dは、数値計算で実行すると、ある程度dが小さいところまではよいのだが、あまりにもdが小さくなると、小さくなればなるほど答えが狂っていく....(微分の定義と正反対の結果)。もちろん、これは浮動小数点計算に特有のマシンイプシロン、つまり相対誤差εのせい。

そこで四倍精度にしてみることにした。

練習のため、まずはMacBook Airにインストールしたg95で円周率πを計算してみた。F77風にかくと、real*16 pi_quadと宣言し、
pi_quad = 4.0q0 * qatan(1.0q0)
とやってみた。するとqatanがありません、というエラーメッセージが出た....調べてみると、1.0q0という表現自体は使用できるようだ。しかし、4倍精度の三角関数が無いんでは、実際問題プログラミングするのはかなり難しいだろう。g95に取り付けるライブラリとか調べてみたが、どうもまだなさそう。バージョンアップを待つしかない。

しかたないので、Fedora12に移り、そこでifortで試してみることにした。すると、今度は一発でコンパイルが通った!そして、その答えは少なくとも33桁まで正確に計算されていた。
π=3.14159265358979323846264338329750
あっぱれ。これからは、quadratic precisionの時代に突入することとなろう。

0 件のコメント: