irrational base discrete weighted transform

前回の記事から少し調べたところ、Lucas-Lehmer testの乗算ではSchönhage-Strassen法は使われていないと知りました。代わりにirrational base discrete weighted transform (IBDWT)というアルゴリズムが使われているそうです。
(cf. GpuOwL - Prime-Wiki)

Schönhage-Strassen法では、N桁の整数x,yに対して適切なnとkを選び、それらを $2^{n}$ 進数 $2^{k}$ 桁(kは4とか16とか)で表して、FFT→要素積→IFFTします。x,yが巨大すぎてnも巨大になる場合は再帰的に処理します。これの時間計算量は桁数Nに対して $O(N log(N) log(log(N)))$ ですが、最後の $log(log(N))$ はこの再帰的に処理するところに由来します。また、Schönhage-Strassen法はxyそのものを計算するので、xとyは実際の桁数の2倍の桁あるかのように計算します。

一方でLucas-Lehmer testの乗算では、メルセンヌ数 $M=2^{q}-1$ 未満の数xを2乗した直後にmod Mしたいので、剰余環 $\mathbb{Z}/M\mathbb{Z}$ の上で計算できれば理想的です。（"2倍の桁あるかのように"が不要になるので）IBDWTではそれが可能で、かつ一度のFFT→要素積→IFFTで済ませられます。
(cf. Irrational base discrete weighted transform - Prime-Wiki)
このページのreferenceの1個目が元論文です。ちなみに3つ目の論文の"generalized IBDWT"は、メルセンヌ数以外を法とする剰余環の上でも同じように計算できるという意味でのgeneralizedです。あと3つ目の論文では数値誤差の評価も詳しくやっています。

実装

FFTW3を使ってIBDWTを実装してみました。前回のようにGMPで愚直にやるバージョンも用意して、計算結果の値の下32桁が同じであることを確認しました。
github.com
結果として、シングルスレッドで1イテレーションあたり0.7秒くらいで計算できました。前回の記事で考えていた値に近くて満足しました。手元の環境ではFFTW3の並列化効率がなぜか悪いとか、繰り上がり・繰り下がりの処理を並列化するのって面倒そうだなとか色々ありますが、果てしないのは仕方ないものです。

ハマったところ

GitHub上のコードに以下の記述があります。無理数の重みベクトルaを計算しています。

for (int64_t j = 0; j < N; ++j) {
	const int64_t index = ceiling(q, j, N) * N - (q * j);
	a[j] = std::pow(2.0, double(index) / double(N));
	assert(1.0 <= a[j] && a[j] <= 2.0);
}

上のコードを下のコードに置き換えると、計算結果の数値誤差がめちゃくちゃ増えてしまいます。std::pow関数自体は割と高精度なのですが、それでも引数に誤差が混ざっていては仕方ありません。

for (int64_t j = 0; j < N; ++j) {
	a[j] = std::pow(2.0, double(ceiling(q, j, N)) - double(q * j) / double(N));
	assert(1.0 <= a[j] && a[j] <= 2.0);
}

そもそも元論文での数値誤差に関する議論では、すべての無理数が最近傍の浮動小数点数に丸められていることを仮定していました。それを読んで、私も最初はSymPyとか使って厳密な最近傍値を計算してテーブル引きしたのですが、 $N=2^{25}$ とかだとSymPyの計算に1日以上かかったりして、それでいてstd::powも（上のコードでは）厳密な値とほとんど変わりませんでした。なので現在の実装に至りました。あと、FFTで使う $e^{\frac{2\pi i j}{N}}$ も最近傍値であるべきですが、FFTW3で高速に計算することにしたのでというのもありました。