TSperl 発言

名無しです。

先日、検証用に C でフーリエ変換のプログラムを作ったのですが、
16 点で動作確認すると問題なさそうなのに、13 万点だと結果が
予想とかけ離れていてなぜだろう？と調べたことがありました。
逆変換してみると確かに元の値に戻らない。

試しに perl の Math::FFT というモジュールを使って計算して見
ると、こっちは精確です。

原因は簡単に見つかって、N^2 ループの中の足し算の桁落ちでした。
伊理正夫・藤野和建の「数値計算の常識」という本の p18 に載っ
ている、桁落ちを集めて再度足しこむ方法を試すとバッチリでした。

Math::FFT で精確な値が出るのは、加算の回数が N*N の DFT と比
べて N * log(N) と少なくて済むからかなと思ったり（FFT を理解
してないので良くわかりませんが、ググるとそれらしい記事もあっ
た）。

これはネタだと思い、DFT のプログラムを perl で書き直して、
対策なし版，桁落ち対策版，bignum 版，Math::FFT 版の結果を並べ
てこんなに差が出たよー、と書こうと昨日の深夜からがんばったの
ですが、データ点数が少ないと結果に差が出ないし、データ点数を
多くすると一晩経っても計算が終わらないしで、結局書けません。

しょうがないので、もっと簡略化して差を出した例をば。

素直に計算した (1) と比べて桁落ち対策版の (2) は 100 倍、精
度が上がってます（三角関数を一周期ぶん足しているので結果は 0
になるのが正しい）。

[スクリプト]
use strict;
use Math::Trig;
my $n = shift or die;
my $sum;

$sum = 0;
foreach my $i (1 .. $n){
    my $x = sin(2 * pi * $i / $n);
    
    # ↓これだと桁落ちする
    $sum += $x;
}
printf "(1) %+.3lg\n", $sum;

$sum = 0;
my($t,$r);
foreach my $i (1 .. $n){
    my $x = sin(2 * pi * $i / $n);
    
    # これも $sum += $x を計算しているが、発生した桁落ちを $r に
    # かき集めて大きくしてから $sum に足し込んでいる。
    $r += $x; $t = $sum; $sum += $r; $t = $sum - $t; $r -= $t;
}
printf "(2) %+.3lg\n", $sum;

[実行結果]
D:$ perl gomi2.pl 100000
(1) +5.51e-013
(2) +3.99e-015

D:$
--