文献あり

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

1446

この連載について

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

本文

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

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

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

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

フーリエ和

関数 $f : [0, 1) \to C$ が次のような三角多項式の形式で表せるとする。
$f (x) = \sum_{k = - n}^{n} c_{k} \exp (2 π i k x) (c_{k} \in C)$

このとき、
$\hat{f} (k) = \int_{0}^{1} f (x) \exp (- 2 π i k x) d x$
とすると、係数 $c_{k}$ は
$c_{k} = \hat{f} (k)$
により定まる。

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

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

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

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

多項式

$n$ を正の整数とする。
$n - 1$ 次多項式 $f : C \to C$ は、次の形式で表せる写像とする。
$f (x) = \sum_{0 \leq k < n} c_{k} x^{k} (c_{k} \in C)$
すなわち、 $n - 1$ 次多項式は $n$ 個の係数 $c_{k}$ のみによって特徴づけられる。
また、必要であれば、 $k \geq n$ における係数を $c_{k} = 0$ とする。

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

離散フーリエ変換

$n$ を正の整数とする。
複素数体 $C$ における $1$ の原始 $n$ 乗根を、 $ω_{n} = e^{\frac{2 π i}{n}} = \exp \frac{2 π i}{n}$ で表す。

$n - 1$ 次多項式 $f : C \to C$ に対し、 $n - 1$ 次多項式 $\hat{f} : C \to C$ を次のように定義する。
$\hat{f} (x) := \sum_{0 \leq k < n} f (ω_{n}^{k}) x^{k}$

$f$ に $\hat{f}$ を対応させる写像を、離散フーリエ変換 (discrete Fourier transform, DFT) と呼び、 $F_{n}$ で表す。

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

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

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

$n$ を正の整数とする。
$1$ の原始 $n$ 乗根 $ω_{n} \in C$ について、次が成り立つ。
$\sum_{0 \leq k < n} ω_{n}^{j k} = {\begin{cases} n (j = 0) \\ 0 (0 < j < n) \end{cases}$

左辺を $S$ とすると、 $S = \sum_{0 \leq k < n} (ω_{n}^{j})^{k}$

(i) $j = 0$ のとき、 $ω_{n}^{j} = 1$ だから、
$S = \sum_{0 \leq k < n} 1 = n$
(ii) $0 < j < n$ のとき、 $ω_{n}^{j} \neq 1$ だから、
$S = \frac{1 - (ω_{n}^{j})^{n}}{1 - ω_{n}^{j}} = \frac{1 - (ω_{n}^{n})^{j}}{1 - ω_{n}^{j}} = 0$

逆離散フーリエ変換

$n$ を正の整数とする。
$n - 1$ 次多項式 $f : C \to C$ に対し、
$f (x) = \frac{1}{n} \sum_{0 \leq k < n} \hat{f} (ω_{n}^{- k}) x^{k}$
が成り立つ。

$\hat{f}$ に $f$ を対応させる写像を、逆離散フーリエ変換 (inverse discrete Fourier transform, IDFT) と呼ぶ。
これは $F_{n}$ の逆元であるため、 $F_{n}^{- 1}$ で表す。

$f (x) = \sum_{0 \leq j < n} c_{j} x^{j} (c_{j} \in C)$
とすると、 $0 \leq m < n$ に対して、
$\hat{f} (ω_{n}^{- m}) = \sum_{0 \leq k < n} (\sum_{0 \leq j < n} c_{j} ω_{n}^{k j}) ω_{n}^{- m k} = \sum_{0 \leq j < n} c_{j} (\sum_{0 \leq k < n} ω_{n}^{k (j - m)})$
補題より、
$\hat{f} (ω_{n}^{- m}) = n c_{m}$

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

逆離散フーリエ変換

$n$ を正の整数とする。
$n - 1$ 次多項式 $f : C \to C$ は、次の形式であるとする。
$f (x) = \sum_{0 \leq k < n} c_{k} x^{k} (c_{k} \in C)$
$n$ 個の係数 $c_{k}$ は未知だが、 $n$ 個の値 $f (ω_{n}^{k})$ は既知であるとする。
このとき、
$\hat{f} (x) := \sum_{0 \leq k < n} f (ω_{n}^{k}) x^{k}$
とすると、係数 $c_{k}$ は
$c_{k} = \frac{1}{n} \hat{f} (ω_{n}^{- k})$
により定まる。

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

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

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

$\hat{f} (x)$ の係数は、 $f$ を入力として $f (ω_{n}^{k})$ で定まる
$f (x)$ の係数は、 $\hat{f}$ を入力として $\frac{1}{n} \hat{f} (ω_{n}^{- k})$ で定まる

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

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

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

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

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

離散フーリエ変換

$n$ を正の整数とする。
$n - 1$ 次多項式 $f \in C^{n}$ に対し、 $n - 1$ 次多項式 $\hat{f} \in C^{n}$ を
$\hat{f} := (f (ω_{n}^{k}))$
で定義する。
写像 $F_{n} : C^{n} ∋ f \mapsto \hat{f} \in C^{n}$ を、離散フーリエ変換 (DFT) と呼ぶ。

逆離散フーリエ変換

$n$ を正の整数とする。
$n - 1$ 次多項式 $f \in C^{n}$ に対し、
$f = (\frac{1}{n} \hat{f} (ω_{n}^{- k}))$
が成り立つ。
写像 $F_{n}^{- 1} : C^{n} ∋ \hat{f} \mapsto f \in C^{n}$ を、逆離散フーリエ変換 (IDFT) と呼ぶ。

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

離散フーリエ変換 (DFT) $F_{n} : f \mapsto (f (ω_{n}^{k}))$
逆離散フーリエ変換 (IDFT) $F_{n}^{- 1} : f \mapsto (\frac{1}{n} f (ω_{n}^{- k}))$

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

プログラミング

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

さてここでは、前節で整理した DFT および IDFT を実装していきます。
ともに $C^{n}$ から $C^{n}$ への写像であり、 $n - 1$ 次多項式の係数 $(c_{k}) \in 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 + 4 x + 5 x^{2} + 6 x^{3} + 123 x^{4}$ を表しています。DFT, IDFT の順に変換すると、出力c2が入力c1とほぼ同じ値になることを確認できます。IDFT, DFT の順に実行しても同様です。
ただし、浮動小数点数を利用しているため、誤差が生じることがあります。計算結果が整数とわかっている場合は丸めればよいです。