はじめに
しばらく前に自作のEasy-ISLispにFFTによる乗算を組み込んでみました。巨大数になると計算がうまく行われず複素数を使ってるから誤差かな？と誤差のないNTTに取り組みました。FFTでうまくいかなかった原因は誤差ではなく単なる私の勘違いによるバグでした。ビットリバースで12bitで計算してました。そりゃ、ダメでしょうよ。ともあれそんなきっかけで剰余環を利用したNTTによる高速乗算に取り組み、動作するところまでこぎつけました。後日のために記録しておこうと思います。
初等整数論のお勉強
数十年前に買った高木貞二先生の「初等整数論講義」を引っ張り出してきて、原始根のあたりを読み直しました。この本には剰余環のことは書いてないようですが、どうもNTTというのは剰余環という整数論の理論を応用しているようです。下記のページの記事がとても勉強になりました。実装例も示されています。
離散フーリエ変換と数論変換 (6) NTT の高速化
https://mathlog.info/articles/2661
うんうん、なんとなくわかったような気になれました。要するに複素数でやっていた回転因子の計算を整数の剰余でやればいいのでしょう。それなら普通のFFTでやったバタフライ計算のコードに手を加えればできるのでは？とやってみました。それほど甘くはありませんでした。どうもうまくいきません。バタフライ計算の順序が先日勉強したあのわかりやすい動画の順序と逆になっているのでした。
【フーリエ解析05】高速フーリエ変換(FFT)ってなに？内側のアルゴリズムを解説！
https://www.youtube.com/watch?v=zWkQX58nXiw&t=39s
手計算であれこれと試行錯誤してみたのですがどうもわかりません。複素数のときには単位円でπ回転すると-1で都合がよかったのですが、剰余を使った１の原始ｎ乗根ではそれは使えません。上記の説明のとおりにやってみることとしました。
理屈はわかったけども
理解したつもりで実装を作り始めたのですが、やってみてハマったことが多々ありました。以下、列挙しますと・・・
桁そろえ
入力データであるベクトルを表す配列にどのようにデータを配置するのか？ 半分にわけて前半に乗数のデータ、後半は０で埋めないと正しく計算されません。
空いた領域は？
乗数を前半部分に置いたとして、空いている領域はどうするのか？ 単純に０で埋めておけばOKでした。
遅いんだけど
実装が素朴すぎました。回転因子に相当する計算をその都度mod計算で冪乗を計算してました。計算前にあらかじめ計算しておいたら速くなりました。
桁数が多くなると変な値がでるよ
最後はこれに悩まされました。なんのことはないビットリバースの計算をたったの12bit分だけしか計算していませんでした。最初は小さいデータで計算実験をしていたので32bitフルに使うなんて考えていませんでした。
性能評価
どうやらまともに動き始めたようでしたのでベンチマークをとってみました。ntt*という関数名でFFT乗算を計算しています。従来の筆算アルゴリズムの＊関数と比較してみました。
Easy-ISLisp Ver2.47
> (load "tests/big.lsp")
T
> (time (bigtest1 100))
Elapsed Time(second)=1.274482
<undef>
> (time (bigtest2 100))
Elapsed Time(second)=31.525318
<undef>
>
(defglobal a (expt 8 70000))
(defglobal b (expt 9 70000))
(defun bigtest1 (n)
(if (= n 0)
t
(progn (ntt* a b) (bigtest1 (- n 1)))))
(defun bigtest2 (n)
(if (= n 0)
t
(progn (* a b) (bigtest2 (- n 1)))))
おお！だいぶ高速になっています。数万桁どうしの乗算だと25倍ほど高速になっていました。
正しく計算できてる？
計算結果がおかしいのでは高速になっても意味がありません。brent-salaminのアルゴリズムで円周率1万桁を計算してみました。BIGNUMの桁数の和が100以上のときにNTT乗算を起動するようにしました。
Easy-ISLisp Ver2.48
> (load "tests/pi.lsp")
T
> (pi-brent-salamin 10000)
314159265358979323846264338327950288419...
...
22955249887275846101264836999892256959688159205600101655256375574
>
This is SBCL 2.0.1.debian, an implementation of ANSI Common Lisp.
* (pi-brent-salamin 10000)
314159265358979323846264338327950288419...
...
22955249887275846101264836999892256959688159205600101655256375574
最初と最後は合ってます。たぶん大丈夫だよ。（大雑把な性格）
GMPよ、あんたはすごいぜ！
Easy-ISLispでは1万行の円周率計算で5秒ほどかかっていました。NTT乗算を組み込みましたら3秒程度まで短縮されていました。ところがSBCLは1秒もかかりません。さすがです。GMPは世界一です。
それでも作ってみたい
車輪の再発明といわれようがともかく自分でも作ってみたいのです。その昔、真空管でオーディオアンプを作ったことを思い出します。そりゃメーカー製のものを買ってきた方がコスパが圧倒的にいいのですけれども。あれから40年以上経って、相変わらず自作に拘っています。仕組みを理解したいのです。GMPを組み込んで・・・なんてのは自称Lisperの私のプライドに反するのです。
終わりに
最後まで駄文にお付き合いいただきまして誠にありがとうございました。実装はgithubに置いてあります。bignum.cの最後の方にNTT関係がまとめられています。コメントも英語でできるだけ入れておきました。後日のデバッグに備えるためです。よかったら実装の参考にしてみてください。