離散フーリエ変換 (DFT) および数論変換 (NTT) の原理、そしてそれらのプログラミングにおける実装方法について記述します。
各用語の定義の違いを明確にするため、それぞれについての回を分割します。
前回は高速フーリエ変換 (FFT) の仕組みとそのプログラミングについて説明しました。
定理で記述した漸化式をそのまま実装すればよいため仕組みがわかりやすいという利点があるものの、この実装では各層で配列を作成して値をコピーしているため、空間計算量も
そこで、今回はビット反転置換 (bit-reversal permutation) を利用して、空間計算量を
前回の FFT に関する定理において、
と2つの多項式に分割していきました。
つまり、その階層において偶数番目と奇数番目の値の列に分割する操作であり、
例えば
空間計算量が
そこで、初めから最下層のインデックスである
以下では、この置換を求める方法について考察します。
上の図を観察すると、各インデックスを
このようなパターンを見ることで、次の定理が得られます。
(i)
(ii)
とし、
とする。
このとき、ある整数を
すなわち、
が成り立つ。
この
例えば、
です。
この置換
static int[] BitReversal_Original(int n)
{
var b = new int[n];
for (int u = 1, d = n >> 1; u < n; u <<= 1, d >>= 1)
for (int k = 0; k < n; ++k)
if ((k & u) != 0)
b[k] |= d;
return b;
}
それぞれのビットごとに
これでも十分なのですが、次のようにさらに
static int[] BitReversal(int n)
{
var b = new int[n];
for (int p = 1, d = n >> 1; p < n; p <<= 1, d >>= 1)
for (int k = 0; k < p; ++k)
b[k | p] = b[k] | d;
return b;
}
これは、
という手順になっています。
入力のベクトルにビット反転置換を施し、前回に実装した FFT の再帰関数を書き換えてみましょう。
なお、置換前と置換後には対称性があるため、「コピー先のインデックス」と考えても「コピー元のインデックス」と考えても同じです。
// c の長さは 2 の冪とします。
// h: 更新対象の長さの半分
static void TransformRecursive(Complex[] c, int l, int h)
{
if (h == 0) return;
TransformRecursive(c, l, h >> 1);
TransformRecursive(c, l + h, h >> 1);
for (int k = 0; k < h; ++k)
{
var v0 = c[l + k];
var v1 = c[l + k + h] * NthRoot(h << 1, k);
c[l + k] = v0 + v1;
c[l + k + h] = v0 - v1;
}
}
// 戻り値の長さは 2 の冪となります。
public static Complex[] Transform(Complex[] c, bool inverse)
{
if (c == null) throw new ArgumentNullException(nameof(c));
var n = ToPowerOf2(c.Length);
var br = BitReversal(n);
var t = new Complex[n];
for (int k = 0; k < c.Length; ++k)
t[br[k]] = c[k];
TransformRecursive(t, 0, n >> 1);
if (inverse && n > 1)
{
Array.Reverse(t, 1, n - 1);
for (int k = 0; k < n; ++k) t[k] /= n;
}
return t;
}
これで空間計算量を
上記と同等の処理を、再帰のない関数として次のように実装することもできます。
再帰関数をトップダウンで順次呼び出していく方法とは逆に、ループを利用してボトムアップで各階層を処理します。
// 戻り値の長さは 2 の冪となります。
public static Complex[] Transform(Complex[] c, bool inverse)
{
if (c == null) throw new ArgumentNullException(nameof(c));
var n = ToPowerOf2(c.Length);
var br = BitReversal(n);
var t = new Complex[n];
for (int k = 0; k < c.Length; ++k)
t[br[k]] = c[k];
for (int h = 1; h < n; h <<= 1)
{
for (int l = 0; l < n; l += h << 1)
{
for (int k = 0; k < h; ++k)
{
var v0 = t[l + k];
var v1 = t[l + k + h] * NthRoot(h << 1, k);
t[l + k] = v0 + v1;
t[l + k + h] = v0 - v1;
}
}
}
if (inverse && n > 1)
{
Array.Reverse(t, 1, n - 1);
for (int k = 0; k < n; ++k) t[k] /= n;
}
return t;
}
ただし、再帰関数の呼び出しは
以上で、より効率のよい FFT を実装できました。
さて、次回以降では数論変換 (NTT) を扱います。
前回:
(3) 高速フーリエ変換 (FFT)
次回:
(5) 数論変換 (NTT)