3

fftとバイナリ分割法を使用して、非常に大きな n (最大 100 万)の階乗を計算するプログラムをCで実装しようとしています。

任意精度の integerを表す単純なライブラリを実装しました。fftとifftを計算するには、「C の数値レシピ」のtwofft.cとfour1.cルーチンを使用します。

特定の n まではすべてうまくいきますが、数値 (浮動配列) が大きすぎると、正規化と丸めの後のifft (four1 で計算) の値が正しくなくなります。

たとえば、40 個のゼロで終わる 2000 桁の 2 つの数値があり、それらを (fft を使用して) 互いに乗算する必要がある場合、ifft を計算すると、いくつかの末尾のゼロが「1」になります。これは、この「ゼロ」の 1 つ (たとえば 0,50009) を丸めたときに、「1」になったために発生します。さて、私の実装が間違っているのか、それともこの数値を別の方法で丸める必要があるのか​​ わかりません。バイナリ分割法と素因数分解の両方を使用しようとしましたが、n >= 9000の場合、結果は間違っています。

これを解決する方法はありますか？ご清聴ありがとうございました。下手な英語で申し訳ありません。

4

2 に答える 2

1

記録のために：

また、数値レシピの four1.c からの ifft によって間違った値が返されることにも気付きました。N=256 の複雑な値を入力としてテストしただけで、実際の時間領域信号のみが得られるように組み立てられました。

結果として得られる時間領域ベクトルは、他の実装の IFFT に対応するようにミラーリング (端から端へ、またはその逆...) し、1 つシフトする必要があります。(私は numpy.fft.ifft、オクターブの ifft、および逆離散フーリエ変換を最適化せずにテストしました。IDFT 式に基づいて、間違いなく正しいはずです)。

数値レシピによって提供されるバージョンには、基本的なアルゴリズムの誤りがなければなりません。彼らの本では、この問題に関連するものは何も説明されていません。

于 2015-12-10T14:35:42.200 に答える