離散フーリエ変換と数論変換 (1) 離散フーリエ変換 (DFT)

この連載について

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

1. 離散フーリエ変換 (DFT) (本稿)
2. 数列の畳み込み
3. 高速フーリエ変換 (FFT)
4. ビット反転置換
5. 数論変換 (NTT)
6. NTT の高速化
7. 中国剰余定理による値の復元

本文

今回は、離散フーリエ変換 (discrete Fourier transform, DFT) の基本的な性質と、プログラミングにおける実装方法について記述します。畳み込みのような応用や、計算の高速化には触れず、まずは原理の記述に焦点を当てます。

この連載では虚数単位を$i$で表し、インデックス (添字) を$k, j$などで表します。
また、$e^{i \theta} = \cos \theta + i \sin \theta$を $\exp i\theta$ で表すことがあります。

(離散でない) フーリエ変換

(離散でない) フーリエ変換 (Fourier transform, FT) というものがあり、その一部として次の結果があります。

ここで、$f:[0, 1) \to \mathbf{C}$は半開区間$[0, 1)$から複素数全体への写像を表しています。

フーリエ変換では、このように積分を利用します。フーリエ和をさらにフーリエ級数に発展させると、それが収束すれば周期関数を三角関数の和 (正弦波の重ね合わせ) で表すことが可能になります。
それに対して、積分の代わりに和で扱うものが離散フーリエ変換と言えます。

離散フーリエ変換と逆離散フーリエ変換

離散フーリエ変換の定義および定理を記述します。
本稿では、三角多項式を用いる代わりに多項式を用いる定義を扱います。

この連載では$c_{n-1} \neq 0$という制約を設けないため、$c_{n-1} = 0$である多項式も便宜上$n-1$次と見なすことができます。文献 [1] では「次数上界が$n$の多項式 (polynomial of degree-bound $n$)」と表現されています。

DFT が$n$に依存していることを表現するため、記号$\mathcal{F}_n$を導入します。
ここでは、$n$は$2$の冪である必要はありません。
後述する高速フーリエ変換 (FFT) では、$n$は$2$の冪とし、その性質を利用することになります。

文献によっては、
$$ ~~ \widehat{f}(t) := \sum_{0 \leq k < n} f(k) \exp \left( -i \dfrac{2\pi k}{n} t \right) = \sum_{0 \leq k < n} f(k) \, \omega_n^{-kt} $$
のように三角多項式として定義されていることもありますが、この場合は$x = \omega_n^{-t}$とすれば本質的には同じものです。

次の定理の証明の核となる補題を用意します。

この結果を、冒頭のフーリエ変換の説明のように言い換えてみましょう。

$n$個の標本点における値をもとに$n$個の係数$c_k$を決定しているということは、つまり$n$元連立$1$次方程式を解いているのです。ここで、標本点として$1$の$n$乗根の集合 $\{ \omega_n^{~k} \}$ を採用すると、$c_k$が美しい形で求められるということです。

$f$の実際の次数が$n-1$より小さい場合でも、範囲外の係数を$0$として$n-1$次多項式と見なすことによりこれまでの議論が成り立ちます。逆に、$f$の次数が$n-1$より大きい場合はこの連立方程式を解けません。

もう一つの要点は、順変換と逆変換に次のような対称性が現れることです。

• $\widehat{f}(x)$の係数は、$f$を入力として$f(\omega_n^{~k})$で定まる
• $f(x)$の係数は、$\widehat{f}$を入力として$\dfrac{1}{n} \widehat{f}(\omega_n^{-k})$で定まる

なお、$\widehat{f}$の定義で$\dfrac{1}{\sqrt{n}}$倍しておくと、逆変換も$\dfrac{1}{\sqrt{n}}$倍となって統一できるのですが、離散フーリエ変換の分野では、逆変換を$\dfrac{1}{n}$倍とする定義が多いようです (おそらく実用上の都合)。

多項式をベクトルとして考える

今後、第$k$成分が$v_k$であるベクトルを$(v_k)$で表すことにします。

多項式の定義において、「$n-1$次多項式は$n$個の係数$c_k$のみによって特徴づけられる」と書きました。
したがって、$n-1$次多項式$f:\mathbf{C} \to \mathbf{C}$を、$n$次元ベクトル $(c_k) \in \mathbf{C}^{~n}$ と同一視できます。
これは、線形代数学でよくある考え方です。

ここでは、前節で見た定義や定理をベクトルで表現してみましょう。

IDFT の入力と出力だけに注目すると、$\widehat{f}$は任意の入力として考えてよい (記号を変えてよい) ため、まとめると次のようになります。

• 離散フーリエ変換 (DFT) $\mathcal{F}_n: f \mapsto (f(\omega_n^{~k}))$
• 逆離散フーリエ変換 (IDFT) $\mathcal{F}_n^{-1}: f \mapsto \left( \dfrac{1}{n} f(\omega_n^{-k}) \right)$

多項式$f$に対して、DFT により「標本点における$f$の値の集合」に変換し、類似の変換である IDFT により$f$を復元できることが明確になりました。
多項式をベクトル (あるいは数列) として考えることは必須ではありませんが、表記や入出力の観点でだいぶ見通しがよくなります。

プログラミング

プログラミングで実用されるのは後述する高速フーリエ変換 (FFT) ですが、原理の理解のために、まずは高速化せずにナイーブな実装を書くのもよいでしょう。ナイーブな実装は、FFT のテストを保証するためにも使うことができます。

さてここでは、前節で整理した DFT および IDFT を実装していきます。
ともに$\mathbf{C}^{~n}$から$\mathbf{C}^{~n}$への写像であり、$n-1$次多項式の係数$(c_k) \in \mathbf{C}^{~n}$を入力として受け取ります。
計算内容が非常に似ているため、実装する関数を共通化します。

C# では次のように実装できます。$n$次元ベクトルは、長さ$n$の配列で表せばよいです。

      using System;
using System.Numerics;

public static class DFT
{
    // 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);
    }

    // f(ω_n^k) の値
    static Complex f(int n, Complex[] c, int k)
    {
        Complex r = 0;
        for (int j = 0; j < c.Length; ++j)
            r += c[j] * NthRoot(n, k * j);
        return r;
    }

    // f の係数が整数のとき、f^ の係数も整数になるとは限りません。
    public static Complex[] Transform(Complex[] c, bool inverse)
    {
        if (c == null) throw new ArgumentNullException(nameof(c));

        var n = c.Length;
        var r = new Complex[n];
        for (int k = 0; k < n; ++k)
            r[k] = inverse ? f(n, c, -k) / n : f(n, c, k);
        return r;
    }
}
    

時間計算量は$O(n^2)$です。FFT ではこれを$O(n \log n)$にすることができます。
ではこの関数を呼び出してみましょう。

      using System;
using System.Numerics;

static class Program
{
    static void Main()
    {
        var c1 = new Complex[] { 3, 4, 5, 6, 123 };
        // DFT
        var t = DFT.Transform(c1, false);
        // IDFT
        var c2 = DFT.Transform(t, true);

        // c1 ≒ c2 (誤差あり)
    }
}
    

c1は、多項式$3 + 4x + 5x^2 + 6x^3 + 123x^4$を表しています。DFT, IDFT の順に変換すると、出力c2が入力c1とほぼ同じ値になることを確認できます。IDFT, DFT の順に実行しても同様です。
ただし、浮動小数点数を利用しているため、誤差が生じることがあります。計算結果が整数とわかっている場合は丸めればよいです。

作成したサンプル コード

• DftNttTest (GitHub)

次回: (2) 数列の畳み込み

改訂履歴

• 2021.09.03 記号$\mathcal{F_n}$を導入。
• 2021.09.13 $1$の原始$n$乗根に関する補題を追加。