離散フーリエ変換と数論変換 (3) 高速フーリエ変換 (FFT)

高速フーリエ変換,DFT,離散フーリエ変換

864

この連載について

離散フーリエ変換 (DFT) および数論変換 (NTT) の原理、そしてそれらのプログラミングにおける実装方法について記述します。
各用語の定義の違いを明確にするため、それぞれについての回を分割します。

本文

今回は、離散フーリエ変換 (DFT) において演算を効率化する高速フーリエ変換 (fast Fourier transform, FFT) と、そのプログラミングについて記述します。

前回に引き続き、多項式をベクトルとして考えることにします。
すなわち、 $n - 1$ 次多項式は $C^{n}$ に属します。

DFT の定義を再確認しておきます。

離散フーリエ変換再確認

$n$ を正の整数とする。
$1$ の原始 $n$ 乗根を $ω_{n} = \exp \frac{2 π i}{n}$ で表す。

$n - 1$ 次多項式 $f = (c_{k}) \in C^{n}$ は、次の形式で表せる写像とする。
$f (x) = \sum_{0 \leq k < n} c_{k} x^{k}$
このとき、
$F_{n} (f) = (f (ω_{n}^{k})) \in C^{n}$
で表される変換 $F_{n}$ を離散フーリエ変換 (DFT) と呼ぶ。

実験

ここでは逆変換である IDFT については考えず、DFT の高速化を考えます。
すなわち、与えられた $(c_{k})$ に対して $(f (ω_{n}^{k}))$ を効率的に求める方法を考察していきます。
まず、 $n = 4$ として実験してみましょう。

$f (x) = c_{0} + c_{1} x + c_{2} x^{2} + c_{3} x^{3}$
として、各 $f (ω_{4}^{k})$ を変形してみると、
$F_{4} (f) = (\begin{matrix} f (1) \\ f (i) \\ f (- 1) \\ f (- i) \end{matrix}) = (\begin{matrix} (c_{0} + c_{2}) + 1 \cdot (c_{1} + c_{3}) \\ (c_{0} - c_{2}) + i \cdot (c_{1} - c_{3}) \\ (c_{0} + c_{2}) - 1 \cdot (c_{1} + c_{3}) \\ (c_{0} - c_{2}) - i \cdot (c_{1} - c_{3}) \end{matrix})$
この式を観察すると、例えば $c_{0} + c_{2}$ や $i \cdot (c_{1} - c_{3})$ など共通する部分を再利用できそうです。またこの部分が、長さが半分の DFT である $F_{2}$ の形となっていることから、効果的な短縮ができそうだと予想できます。

実際、
$g_{0} (x) := c_{0} + c_{2} x g_{1} (x) := c_{1} + c_{3} x$
とすると、
$f (x) = c_{0} + c_{2} x^{2} + x (c_{1} + c_{3} x^{2}) = g_{0} (x^{2}) + x g_{1} (x^{2})$
の形式となります。
また、 $(ω_{4}^{k})^{2} = ω_{2}^{k}$ であることから、 $F_{2} (g_{0})$ および $F_{2} (g_{1})$ の各成分が現れます。

高速フーリエ変換

前節における考察を一般化します。
$n$ を $2$ の冪とします。

$n > 1$ のとき、
$g_{0} (x) := c_{0} + c_{2} x + c_{4} x^{2} + \dots = \sum_{0 \leq k < n / 2} c_{2 k} x^{k} g_{1} (x) := c_{1} + c_{3} x + c_{5} x^{2} + \dots = \sum_{0 \leq k < n / 2} c_{2 k + 1} x^{k}$
とすると、 $g_{0}, g_{1} \in C^{n / 2}$ であり、
$f (x) = g_{0} (x^{2}) + x g_{1} (x^{2})$
と表せます。

次に、 $F_{n} (f) = (f (ω_{n}^{k}))$ のそれぞれの成分を求めていくと、 $0 \leq k < n / 2$ に対して、
$f (ω_{n}^{k}) = g_{0} (ω_{n / 2}^{k}) + ω_{n}^{k} g_{1} (ω_{n / 2}^{k}) f (ω_{n}^{k + n / 2}) = f (- ω_{n}^{k}) = g_{0} (ω_{n / 2}^{k}) - ω_{n}^{k} g_{1} (ω_{n / 2}^{k})$
と計算できます。
つまり、 $F_{n} (f)$ を前後半に分割して考えると、 $F_{n / 2} (g_{0}) = (g_{0} (ω_{n / 2}^{k}))$ および $F_{n / 2} (g_{1}) = (g_{1} (ω_{n / 2}^{k}))$ が既知であれば、 $F_{n} (f)$ を時間計算量 $O (n)$ で求められます。

上の式を書き換えて、
$F_{n} (f)_{k} = F_{n / 2} (g_{0})_{k} + ω_{n}^{k} F_{n / 2} (g_{1})_{k} F_{n} (f)_{k + n / 2} = F_{n / 2} (g_{0})_{k} - ω_{n}^{k} F_{n / 2} (g_{1})_{k}$
のようにも表現できます。

また、 $n = 1$ のときは $F_{1} (f) = c_{0} \in C$ とすればよいです。
時間計算量については、上記の操作を $n$ が $2$ の冪のときだけ実行するため、 $O (n \log n)$ となります。

以上の考察を定理の形でまとめておきます。

高速フーリエ変換

$n$ を $2$ の冪とする。
$n - 1$ 次多項式 $f = (c_{k}) \in C^{n}$ に対し、 $F_{n} (f) \in C^{n}$ を求めるための漸化式は、次で与えられる。

(i) $n = 1$ のとき、 $F_{1} (f) = c_{0} \in C$
(ii) $n > 1$ のとき、
${\begin{cases} g_{0} := (c_{2 k}) \in C^{n / 2} \\ g_{1} := (c_{2 k + 1}) \in C^{n / 2} \end{cases}$
とすると、
${\begin{cases} F_{n} (f)_{k} = F_{n / 2} (g_{0})_{k} + ω_{n}^{k} F_{n / 2} (g_{1})_{k} \\ F_{n} (f)_{k + n / 2} = F_{n / 2} (g_{0})_{k} - ω_{n}^{k} F_{n / 2} (g_{1})_{k} \end{cases} (0 \leq k < n / 2)$

また、時間計算量は $O (n \log n)$ である。
この漸化式による計算方法を、高速フーリエ変換 (fast Fourier transform, FFT) と呼ぶ。

上記の考察から、数学的帰納法により成り立つ。

このように、サイズの小さい同じ問題に分割して再帰的に解く方法は分割統治法 (divide-and-conquer algorithm) と呼ばれます。例えば、マージソートも同様の手法によります。

FFT の逆変換である逆高速フーリエ変換 (inverse fast Fourier transform, IFFT) についても、同様の考察から次の漸化式が得られます。

逆高速フーリエ変換

$n$ を $2$ の冪とする。
$n - 1$ 次多項式 $f = (c_{k}) \in C^{n}$ に対し、 $F_{n}^{- 1} (f) \in C^{n}$ を求めるための漸化式は、次で与えられる。

(i) $n = 1$ のとき、 $F_{1}^{- 1} (f) = c_{0} \in C$
(ii) $n > 1$ のとき、
${\begin{cases} g_{0} := (c_{2 k}) \in C^{n / 2} \\ g_{1} := (c_{2 k + 1}) \in C^{n / 2} \end{cases}$
とすると、
${\begin{cases} F_{n}^{- 1} (f)_{k} = \frac{1}{2} (F_{n / 2}^{- 1} (g_{0})_{k} + ω_{n}^{- k} F_{n / 2}^{- 1} (g_{1})_{k}) \\ F_{n}^{- 1} (f)_{k + n / 2} = \frac{1}{2} (F_{n / 2}^{- 1} (g_{0})_{k} - ω_{n}^{- k} F_{n / 2}^{- 1} (g_{1})_{k}) \end{cases} (0 \leq k < n / 2)$

FFT と異なる点は、 $ω_{n}^{k}$ の代わりに $ω_{n}^{- k}$ を使うことと、各階層ごとに計算結果を $2$ で割ることです。

また、次のように DFT の結果を利用して IDFT の結果を求める方法もあります。
これについては FFT である必要はなく、前回までの議論でも成り立ちます。

逆離散フーリエ変換別解

$n$ を正の整数とする。
$n - 1$ 次多項式 $f = (c_{k}) \in C^{n}$ に対し、 $F_{n}^{- 1} (f) \in C^{n}$ は、次の式で求められる。
$F_{n}^{- 1} (f)_{k} = {\begin{cases} \frac{1}{n} F_{n} (f)_{0} (k = 0) \\ \frac{1}{n} F_{n} (f)_{n - k} (0 < k < n) \end{cases}$

$F_{n}^{- 1} (f) = (\frac{1}{n} f (ω_{n}^{- k})) = (\frac{1}{n} f (ω_{n}^{n - k}))$
により成り立つ。

つまり、DFT の結果に対して、 $0 < k < n$ の範囲を逆順に並べ替え、 $n$ で割れば IDFT の結果になるということです。

以上で、DFT および IDFT が高速化され、時間計算量 $O (n \log n)$ で実行できることがわかりました。
これにより、畳み込み (多項式の積) も $O (n \log n)$ に高速化されます。

プログラミング

定理 1 および定理 3 で得られた式をそのまま実装します。
定理 1 の主張が再帰で表現されているため、この実装においても再帰関数を使います。

入力の配列の長さが $2$ の冪でない場合、それを超える最小の $2$ の冪を $n$ とします。戻り値の配列の長さもその $n$ となります。
Convolution 関数については前回の実装とほぼ同じで、 $n$ を $2$ の冪に変更する部分のみ異なります。

      using System;
using System.Numerics;

public static class FFT
{
    public static int ToPowerOf2(int n)
    {
        var p = 1;
        while (p < n) p <<= 1;
        return p;
    }

    // k 番目の 1 の n 乗根 (n-th root)
    static Complex NthRoot(int n, int k)
    {
        var t = 2 * Math.PI * k / n;
        return Complex.FromPolarCoordinates(1, t);
    }

    // c の長さは 2 の冪とします。
    static void TransformRecursive(Complex[] c)
    {
        var n = c.Length;
        if (n == 1) return;

        var n2 = n >> 1;
        var c0 = new Complex[n2];
        var c1 = new Complex[n2];
        for (int k = 0; k < n2; ++k)
        {
            c0[k] = c[2 * k];
            c1[k] = c[2 * k + 1];
        }

        TransformRecursive(c0);
        TransformRecursive(c1);

        for (int k = 0; k < n2; ++k)
        {
            var v0 = c0[k];
            var v1 = c1[k] * NthRoot(n, k);
            c[k] = v0 + v1;
            c[k + n2] = 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 t = new Complex[n];
        c.CopyTo(t, 0);

        TransformRecursive(t);

        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[] Convolution(Complex[] a, Complex[] b)
    {
        if (a == null) throw new ArgumentNullException(nameof(a));
        if (b == null) throw new ArgumentNullException(nameof(b));

        var n = ToPowerOf2(a.Length + b.Length - 1);
        Array.Resize(ref a, n);
        Array.Resize(ref b, n);

        var fa = Transform(a, false);
        var fb = Transform(b, false);

        for (int k = 0; k < n; ++k)
        {
            fa[k] *= fb[k];
        }
        return Transform(fa, true);
    }
}