自作LispにFFT乗算を組み込んだ話

はじめに

最近になって自作のEasy-ISLispのBIGNUMを大幅に改良しました。円周率1万桁程度の計算は可能になりました。しかし、遅いのです。SBCLなどが瞬時にやってのける計算に数秒かかっています。SBCLなど多くのLisp処理系ではGMPが採用されています。これは世界最強です。Lisperとしてはこれに挑戦、少しでも近づいてやろうという気持ちになります。離散高速フーリエ変換（FFT）を応用した乗算に取り組んでみました。カラツバは敢えて採用しませんでした。ロシアのウクライナ軍事進攻に対する私の抗議の気持ちが込められています。

この車輪の再発明には意味がある

フーリエ変換という数学的アイディアがベースになっています。フーリエという数学者は飲んだくれだったそうで私としては親近感を感じます。このフーリエ変換をコンピューターに応用できるようにしたのが離散的フーリエ変換であり、その高速計算のアルゴリズムがFFTというものです。その仕組みは再帰的な考え方になっておりLisperならずとも「エレガント！」と叫びたくなることでしょう。現代の画像圧縮処理など多くの場面でFFTは応用されています。FFTをしっかり理解し、マスターすることには大きな意味があると思います。

理論と数学のお勉強

QiitaにはFFTに関する数学的な説明、初学者向けの説明が多くありますので、説明はそれらをご覧ください。次の動画が私にとってとても理解の助けになりました。数式だけではピンとこなかった部分も平易な説明を聞いていると具体的なイメージが湧いてきます。

山口大学大学院　けんゆーさんによる解説動画
https://www.youtube.com/watch?v=zWkQX58nXiw&t=1242s

FFTの乗算への応用については下記の動画による講義を拝見して勉強しました。
https://www.youtube.com/watch?v=mIEhAPL0Ua8&t=557s

心臓部分

FFT計算、その逆変換をする関数を　FFT(X) IFFT(x) と表記しますと乗算は次のようになります。

乗数、被乗数のベクトル　X,Y

IFFT(FFT(x) * FFT(y))  -→正規化

手順

多倍長整数を複素数を要素とする配列にします。それをFFTによりフーリエ変換をします。その計算結果の各々の配列要素である複素数の積をとったものを計算します。さらにそれをIFFTで逆フーリエ変換をします。不思議なことにこれは乗算結果となっています。さらに桁上げの処理をして多倍長整数に変換します。

多倍長整数を複素数配列に変換する

Easy-ISLispでの多倍長整数は10^10進法になっています。１つの桁は32bit整数で実装されており０～999,999,999の整数が入っています。これを3桁つづ3つに分解して複素数に変換してfttxという配列にいれておきます。なぜ3桁づつに分解するのかといいますとオーバーフローを防ぐためです。999,999,999という整数を複素数に変換してさらに積をとるとオーバーフローを起こします。基数を10^3として分解することでこれを防ぎます。

pointer = get_pointer (x) - lenx + 1;	//LSB

  for (i = 0; i < lenx; i++)
    {
      //one bigcell separate to three FFT data.
      fftx[3 * i] = (double complex) (bigcell[pointer + i] % FFTBASE);
      fftx[3 * i + 1] =
	(double complex) ((bigcell[pointer + i] / FFTBASE) % FFTBASE);
      fftx[3 * i + 2] =
	(double complex) (bigcell[pointer + i] / (FFTBASE * FFTBASE));
    }

FFT IFFT 本体

再帰計算によっています。fttxという配列に変換された多倍長整数にFFT計算をしてfttx配列に戻します。複素数による回転因子を計算する必要があります。

//FFTの本体
void
fft1 (int n, int pos)
{
  double complex temp;

  if (n == 2)
    {
      temp = fftx[pos] + fftx[pos + 1];
      fftx[pos + 1] = fftx[pos] - fftx[pos + 1];
      fftx[pos] = temp;
    }
  else
    {
      int i, half;
      half = n / 2;
      for (i = 0; i < half; i++)
	{
	  temp = fftx[pos + i] + fftx[pos + half + i];
	  fftx[pos + half + i] =
	    w_factor (n, i) * (fftx[pos + i] - fftx[pos + half + i]);
	  fftx[pos + i] = temp;
	}
      //recursion
      fft1 (half, pos);
      fft1 (half, pos + half);
    }
}

//FFT
void
fft (int n)
{

  fft1 (n, 0);

}

//回転因子
double complex
w_factor (int n, int i)
{
  double complex z;

  z = cos (2.0 * M_PI / (double) n) - sin (2.0 * M_PI / (double) n) * I;
  return (cexpt (z, i));
}

逆離散フーリエ変換

これについて解説している文書が少ないのです。でもこの仕組みはとても単純です。回転因子の計算方法を変えるのと最後にNで割ることだけです。あとはFFTをまったく同じです。FFTでは回転因子は時計方向に回ります。IFFTでは時計と逆方向に回ります。複素関数の勉強が役に立ちました。オイラーの複素数の公式は美しいですね。

//FFT用
double complex
w_factor (int n, int i)
{
  double complex z;

  z = cos (2.0 * M_PI / (double) n) - sin (2.0 * M_PI / (double) n) * I;
  return (cexpt (z, i));
}
//逆FFT用　符号が違うだけ
// for inverse FFT
double complex
iw_factor (int n, int i)
{
  double complex z;

  z = cos (2.0 * M_PI / (double) n) + sin (2.0 * M_PI / (double) n) * I;
  return (cexpt (z, i));
}

ビットリバース

この方法に依った場合、計算された配列のインデックスにつき補正が必要になります。それを行うのがビットリバースです。

int
bit_reverse (int n, int bit)
{
  int binary[12], i;


  for (i = 0; i < bit; i++)
    {
      binary[i] = n % 2;
      n = n / 2;
    }

  n = 0;
  for (i = 0; i < bit; i++)
    {
      n = binary[i] + n * 2;
    }

  return (n);
}

正規化

乗数、被乗数のそれぞれについてFFT計算をし、その配列の各要素について積をとります。その結果につき逆離散フーリエ変換を施します。その結果について正規化をする必要があります。基数の10^3単位で桁上がりをしつつもとの10^10進法の整数に変換します。

// normalize
  int pool, carry;
  carry = 0;
  for (i = 0; i < ans_len; i++)
    {
      pool = (((int) round (creal (ffty[3 * i])) + carry) % FFTBASE);
      carry = ((int) round (creal (ffty[3 * i])) + carry) / FFTBASE;
      pool =
	pool +
	((((int) round (creal (ffty[3 * i + 1])) +
	   carry) % FFTBASE) * FFTBASE);
      carry = ((int) round (creal (ffty[3 * i + 1])) + carry) / FFTBASE;
      pool =
	pool +
	((((int) round (creal (ffty[3 * i + 2])) +
	   carry) % FFTBASE) * (FFTBASE * FFTBASE));
      carry = ((int) round (creal (ffty[3 * i + 2])) + carry) / FFTBASE;
      bigcell[big_pt0++] = pool;
    }

実装面での落とし穴

まずまず正しい答えがでたところで正規化でハマりました。複素数の実数部を整数に変換する部分です。当初、(int)でキャストをあてていました。これだと1つ小さな値になってしまいます。ceil関数で切り捨ててから(int)でキャストするとこれは防げました。しかし値が正しい計算結果を微妙にことなります。試行錯誤しつつroundで四捨五入をして(int)でキャストすると正しい結果になることがわかりました。

さて、それで計算は正しくなったものの、計算速度が遅いのです。筆算アルゴリズムよりも遅い結果になっています。原因を調べていくと原因は回転因子を計算するときに使ったcpow関数でした。複素数のべき乗を計算するものです。これは実数、複素数のべき乗にも対応するためテイラー展開をしているのでした。これでは遅いわけです。整数べきをとるなら単純に掛け算をすればいいのでした。SICPで紹介されているあのアルゴリズムで書き直しました。これにより筆算アルゴリズムよりも高速になりました。

(defun bigtest1 (n)
    (if (= n 0)
        t
        (progn (fft* a b) (bigtest1 (- n 1)))))

(defun bigtest2 (n)
    (if (= n 0)
        t
        (progn (* a b) (bigtest2 (- n 1)))))

Easy-ISLisp Ver2.44
> (load "tests/big.lsp")
T
> (time (bigtest1 100))
Elapsed Time(second)=0.136918
<undef>
> (time (bigtest2 100))
Elapsed Time(second)=0.236661
<undef>
> 

有効な計算の範囲

上記の実装により6000桁程度同士の乗算までなら正しい計算がされていることがわかりました。複素数の積をとるときにオーバーフローを起こしているものと思われます。FFTに与える要素は基数を10^3としました。最大で999が与えられます。FFTによりx0の要素には最大で999×配列要素数の値を実部にもつの複素数が入ります。これの積をとったときにオーバーフローしているものと思われます。複素数はc言語のdouble complex を利用しています。これですと実部、虚部はdouboleによっており15桁の有効桁を持っています。Easy－ISLISPの場合BIGNUMは10^10進法ですので1万桁の数は１００００÷９＝１１１１セルとなります。これをFFTに変換すると3倍の３３３３セルとなります。９９９×３３３３＝3329667です。これの二乗をとっても12桁です。15桁に収まるはずなのですが・・・・。その手前で計算がおかしくなります。まだまだ私の理解が不足しているようです。

追記
６０００桁以上で正しい答えが出ない理由は複素数計算における誤差だと思われます。初等整数論の原始根を応用したFFTに取り組む予定です。

追記（2022/6/11）
原因はbitreverseのバグでした。64ビットで計算するところをわずか12ビットで計算していたためでした。

終わりに

FFT乗算に初挑戦したのはおよそ10年ほど前でした。Kent教授のScheme本にあったFFTのコードで試していました。当時は資料も少なく理解不足でした。あれから中途半端にしていたのですが、この度一念発起して一応計算できるところまでたどり着きました。FFTの理屈はわかって実装したみたものの、いろいろとハマりどころがありました。実装してみようという人の参考になりましたら幸いです。
実装はGithubにてOSSとして公開しています。bignum.cというところにあるbigx_fft_multという関数がそれです。
https://github.com/sasagawa888/eisl