多倍長整数の掛け算への応用を前提とした高速フーリエ変換高速化に関する検証

大きな桁数の多倍長整数同士の乗算速度は、高速フーリエ変換を用いた畳み込みの速度に大きく依存します。今回は、高速フーリエ変換（FFT）を用いた乗算について、アルゴリズムや実装の工夫による実行速度の改善度合いを検証してまとめてみることにしました。本記事は個人的な記録であり、決して高速な実装ではない点にお気をつけください。本取り組みのきっかけは約1年ほど前に公開した以下の記事です。

任意精度で計算してみた。

実装の違いによるFFT実行速度の検証

系列長$N$がどのような数か（素数、合成数、2の冪、など）によってFFTにはさまざまなバリエーションがあります。FFTやそのバリエーションの仕組みに関する解説は素晴らしい資料が多数存在するため、本記事ではあまり踏み込まないこととします。

同じFFTアルゴリズムでも実装方法によって実際の計算時間には大きな違いが生じます。本記事の主な目的は、FFTの実装方法及び乗算の計算方法によって実行速度がどの程度異なるのかを実測することです。今回は以下の実装を比較実験に使用します。

再帰関数を用いたFFT

系列を奇数成分と偶数成分の2つに再帰的に分割しながら計算する方法です。再帰関数を用いた実装はわかりやすく、要素の並べ替えも必要ありません。一方で、重複した処理が含まれること、再帰関数がそこまで高速ではないことから計算速度に難点があります。

ループを用いたFFT

再帰関数をループとして展開することで計算の重複を省き、より効率的に処理することができます。本記事ではCooley-Turleyのアルゴリズムを用いて検証します。

キャッシュ効率を考慮したFFT

上記アルゴリズムのループ順序を入れ替え、配列へのアクセスがなるべく連続的になるように、つまり、隣接した要素にアクセスするように修正することでキャッシュメモリのヒット率が向上し、実行速度の改善につながります。

乗算高速化に向けた検証

準備

多倍長整数を表現するためのデータ構造

16ビット整数型の配列に4桁ずつ数値を詰めることで任意の桁数の数を表現しています。

多倍長整数の乗算

簡単のために自然数に限定して議論します。$m$進数で$n$桁の数$A$は、各桁の数字を$a_i$として$A=\sum_{i=0}^{n}a_{i}m^{i}$と表現できます。例えば10進数で$a_5a_4a_3a_2a_1a_0$と記述される数は、$\sum_{i=0}^{5}a_{i}10^{i}$と表現することができます。

従って、二つの数$A$と$B$の積$C$を求めるには、各桁の数を表す整数列$c_{k}=\sum_{i+j=k}a_ib_j$がわかれば良い、ということになります。筆算方式は二重ループを使う最も直感的な計算方法です。この方式の計算量はAとBの桁数を$N$として、$O(N^2)$になります。

大きな桁数の数値の積を高速に計算したい場合に登場するのが高速フーリエ変換 (FFT)です。FFTを用いることで、畳み込みと呼ばれる操作を系列長$N$に対して$O(NlogN)$の計算量で行うことができます。

実数列を対象としたFFTの効率化

多倍長整数の乗算を前提とする場合、実数列の畳み込み演算ができれば十分です。通常のFFTは複素数を対象としているため、実数列を複素数列に変換する必要があります。最も単純な変換は、実数列を実部にもつような複素数列への変換です。この方法では虚部が活用されていないため、効率化の余地があります。

2本の実数列をまとめてFFT

一方の実数列$\{a_k\}_{k=1}^{n}$を実部に、もう一方の実数列$\{b_k\}_{k=1}^{n}$を虚部に持つような複素数列$\{c_k\}_{k=1}^{n}=\{a_k+i b_k\}_{k=1}^{n}$を考えます。$\{c_k\}$を入力としてFFTを行ったのち、後処理を行うことで$\{a_k\}, \{b_k\}$をそれぞれ独立にFFTした場合と同じ計算結果を得ることができます。

後処理の方法は以下の通りです。

$$
\hat{a}_k=\dfrac{\hat{c}_k+\bar{\hat{c}}_{n-k}}{2}, \hat{b}_k=\dfrac{\hat{c}_k-\bar{\hat{c}}_{n-k}}{2i}
$$
ここで、$\hat{a}, \hat{b}, \hat{c}$はフーリエ変換後の数列を、$\bar{c}$は$c$の複素共役を表します。

実数列を分割してFFT長を半分にする

虚部を有効活用するというアイデアは2本の実数列をまとめたFFTと同様です。まずは1本の実数列$\{a_k\}_{k=1}^{n}$を添え字の偶奇に従って二分割します。分割後の実数列をもとに複素数列を構成します。実部を$\{a_{2l}\}$、虚部を$\{a_{2l+1}\}$とすることでFFTへの入力長を半分にすることができます。この複素数列を入力としてFFTを行い、最後に後処理を行うことで元の実数列に対する計算結果を復元します。後処理はこちらの記事に記載の方法で行います。

実験

実装の違いによる計算速度の違い、入力長の違いによる計算速度の違いを評価するため、ベンチマークを行いました。

Rustを用いて実装し、v1.86の環境でCriterionを用いてベンチマークを行いました。使用した計算機はMacbook Air 2022 (Apple M2 chip, 16GB)です。OSはSequoia 15.4.1です。実装の都合上、入力長は2の冪に限定します。

recursiveが再帰関数を用いたFFTの実装、radix2がループを用いたFFT実装、radix2ditがキャッシュ効率を考慮したFFT実装です。radix2_real, radix2dit_realは実数列を分割してFFT長を半分にする実装を表します。

どの実装も計算量のオーダーは同じですが、定数倍の違いで速度が大きく異なります。特に再帰関数を用いた実装とそれ以外では100倍以上の差があります。系列長が短い場合、ループを用いたFFTの実行速度に大きな違いはありません。一方で系列長が長くなるとキャッシュ効率を考慮した実装の優位性が目立つようになります。

FFT長を半分にする工夫の効果はキャッシュ効率の考慮ほどは大きくありませんが、虚部を有効活用しない実装よりも一貫して高速に動いています。