離散フーリエ変換と数論変換 (2) 数列の畳み込み

畳み込み,DFT,離散フーリエ変換

838

この連載について

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

本文

今回は、離散フーリエ変換における数列の畳み込み (convolution) と、そのプログラミングについて記述します。計算の高速化については次回とし、原理の記述に焦点を当てます。

前回に引き続き、多項式をベクトルとして考えることにします。
すなわち、 $n - 1$ 次多項式は $C^{n}$ に属します。
また、ベクトル $f$ の第 $k$ 成分を $f_{k}$ で表します (つまり多項式における $x^{k}$ の係数)。

数列の畳み込み

畳み込みには多くの種類があります。これは、畳み込みという用語が特定の式を指すのではなく、同じ性質を持つものすべてに対して定義できるためです。
この連載では多項式を取り扱っている文脈から、数列の畳み込み (この用語は文献 [2] による。単に畳み込みでもよいはず) を定義します。

数列の畳み込み・多項式の積

$n, m$ を正の整数とする。
$n - 1$ 次多項式 $f \in C^{n}$ および $m - 1$ 次多項式 $g \in C^{m}$ に対し、 $n + m - 2$ 次多項式 $f * g \in C^{n + m - 1}$ を
$(f * g) (x) := f (x) g (x)$
で定義する。
この $f * g$ を、 $f$ と $g$ の畳み込み (convolution) と呼ぶ。

これでは畳み込みの定義になっていないように見えるかもしれませんが、
$f (x) = \sum_{0 \leq k < n + m - 1} a_{k} x^{k} (a_{k} \in C) g (x) = \sum_{0 \leq k < n + m - 1} b_{k} x^{k} (b_{k} \in C) a_{k} = 0 (n \leq k < n + m - 1) b_{k} = 0 (m \leq k < n + m - 1)$
であるとし、多項式の乗算をしてみると、分配法則により、
$f (x) g (x) = \sum_{0 \leq k < n + m - 1} (\sum_{0 \leq j \leq k} a_{j} b_{k - j}) x^{k}$
となります。

この式に現れた
$c_{k} := \sum_{0 \leq j \leq k} a_{j} b_{k - j} (0 \leq k < n + m - 1)$
により定まる数列あるいはベクトル $(c_{k}) \in C^{n + m - 1}$ を畳み込みの定義としてもよいです。
こちらのほうが畳み込みらしい形をしており原義通りですが、どちらも同じものです。
このように、数列の畳み込みは2つの多項式の積そのものとして定義できます。

畳み込み定理

$n, m_{1}, m_{2}$ を正の整数とし、 $m_{1} + m_{2} - 1 \leq n$ とする。
$m_{1} - 1$ 次多項式 $f \in C^{m_{1}}$ および $m_{2} - 1$ 次多項式 $g \in C^{m_{2}}$ に対し、
$F_{n} (f * g)_{k} = F_{n} (f)_{k} F_{n} (g)_{k} (0 \leq k < n)$
すなわち、
$F_{n} (f * g) = (F_{n} (f)_{k} F_{n} (g)_{k}) \in C^{n}$
が成り立つ。

DFT の定義および畳み込みの定義から、 $0 \leq k < n$ に対し、
$F_{n} (f * g)_{k} = (f * g) (ω_{n}^{k}) = f (ω_{n}^{k}) g (ω_{n}^{k}) = F_{n} (f)_{k} F_{n} (g)_{k}$

制約 $m_{1} + m_{2} - 1 \leq n$ は、多項式 $f$ と $g$ の積が $n - 1$ 次以下であることを保証するために設定しています。

この定理の中で、 $f * g, f, g$ の3つに対する DFT が現れます。これらの $F_{n}$ の $n$ を一致させておかないと、標本点に使われる $1$ の $n$ 乗根 $ω_{n}$ の値が変わってしまいます ( $m_{1}$ や $m_{2}$ としないように)。

この結果を少し言い換えて、次の系が得られます。

畳み込み定理

$n, m_{1}, m_{2}$ を正の整数とし、 $m_{1} + m_{2} - 1 \leq n$ とする。
$m_{1} - 1$ 次多項式 $f \in C^{m_{1}}$ および $m_{2} - 1$ 次多項式 $g \in C^{m_{2}}$ が既知であるとき、
多項式の積 $f * g \in C^{n}$ を DFT および IDFT により求めることができる。

具体的な手順としては、 $F_{n} (f * g)$ の第 $k$ 成分を $F_{n} (f)_{k} F_{n} (g)_{k}$ として計算し、
この $F_{n} (f * g)$ に $F_{n}^{- 1}$ を作用させればよい。

多項式 $f, g$ からそれらの積 $f * g$ を求めることは、単純な分配法則により時間計算量 $O (m_{1} m_{2})$ で可能です。
したがって、上記の変換手順で $O (n^{2})$ となるのでは意味がありませんが、後述する高速フーリエ変換 (FFT) により $O (n \log n)$ で実行できるようになります。

また、やはりこの時点においても、 $n$ が $2$ の冪である必要はありません。
FFT では $2$ の冪の性質を利用することになります。

プログラミング

今回は数列の畳み込みを前節で整理した手順に従って実装します。
DFT および IDFT は前回に実装したものを使います。

前回に引き続き、プログラミング言語は C# です。
配列 a, b の長さはそれぞれ $m_{1}, m_{2}$ に相当します。
DFT を実行する前に、配列の長さを $n$ に変更しています (拡張された部分の係数は $0$ )。

      using System;
using System.Numerics;

public static class DFT
{
    // 前回の Transform メソッドは省略。
    // public static Complex[] Transform(Complex[] c, bool inverse)

    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 = 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);
    }
}

前回実装した DFT および IDFT を実行するため、時間計算量は $O (n^{2})$ です。FFT ではこれを $O (n \log n)$ にすることができます。
ではこの関数を呼び出してみましょう。

      using System;
using System.Numerics;

static class Program
{
    static void Main()
    {
        var a = new Complex[] { 2, 1, 1 };
        var b = new Complex[] { -1, -1, 1 };

        // { -2, -3, 0, 0, 1 }
        var c = DFT.Convolution(a, b);
    }
}

この例は $(x^{2} + x + 2) (x^{2} - x - 1) = x^{4} - 3 x - 2$ を表しています。
$m_{1} = 3, m_{2} = 3$ としているため、 $n = 5$ となります。
今回も浮動小数点数を利用しているため、誤差が生じることがあります。計算結果が整数とわかっている場合は丸めればよいです。