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

現在、かなりシビアな精度での数値計算を実行中。零割が頻発して、計算が止まる事象が続発。問題は差分計算。微分公式 lim_d→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という表現自体は使用できるようだ。しかし、４倍精度の三角関数が無いんでは、実際問題プログラミングするのはかなり難しいだろう。g95に取り付けるライブラリとか調べてみたが、どうもまだなさそう。バージョンアップを待つしかない。

しかたないので、Fedora12に移り、そこでifortで試してみることにした。すると、今度は一発でコンパイルが通った！そして、その答えは少なくとも３３桁まで正確に計算されていた。

π=3.14159265358979323846264338329750

あっぱれ。これからは、quadratic precisionの時代に突入することとなろう。