$$
\newcommand{vand}{\mathcal{V}}
\newcommand{set}[1]{\left\{ #1 \right\}}
\newcommand{finf}{\mathbb{F}}
$$
[TOC]
**キーワード:** 離散フーリエ変換、高速離散フーリエ変換、多倍長演算、1変数多項式
# はじめに: 内容と背景
離散フーリエ変換は、導出の部分やそこに絡む直観をすっ飛ばすと、やっていることはVandermonde行列という特殊行列を用いた行列積です。
* これは一般的な教科書にも書かれている知られた話です。
* 例えば、手元にある [この応用数学の本](https://www.morikita.co.jp/books/mid/005741) は（Vandermonde行列という固有名詞は出てこないのですが) このスタイルで説明されています。
一方で、離散フーリエ変換といえば、高速っしょ！！という名言（要出典）もある通り、高速フーリエ変換です。はい。[実験のセクション](#7-実装)でみる通り、相当速くなります。
ということで、高速化をVandermonde行列を起点に再導出する、という**備忘録**的な記事というのが、この文章の正体です。同じような説明は探せば出るので、新しい結果などは全く含んでいないです。
高速化の話をしたいものの、離散フーリエ変換が既に存在しているところから話を始めると面白くないのと、一方で連続フーリエ変換経由で離散フーリエ変換の導出をやる一般的な方法だと長くなりすぎます。
そこで本記事では一つの工夫として、多倍長演算の実装という、これまた結局Vandermonde行列積に行き着く別のストーリーを採用します。多倍長演算の実装に離散フーリエ変換を使うのは非常に枯れた方法で、この話の展開の仕方も新しいものではないです。
あくまでよく知られた話をまとめたメモ程度の記事です。
## 前提知識
* 多項式を知っていて、多項式の積が定義・実行できる。
* 行列積を定義・実行ができる。
* 逆行列の定義を知っている。
* Vandermonde行列というのをきいたことがある。
# 1. 主役: Vandermonde Matrix
この記事では、$n \times n$ サイズの [Vandermonde（ヴァンデルモンド）行列](https://en.wikipedia.org/wiki/Vandermonde_matrix)を、2つのパラメータ $\alpha$ と $n$ を持つ形で表すことにします:
$$
\vand_{\alpha, n} =
\begin{pmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \alpha & \alpha^2 & \cdots & \alpha^{n-1} \\
1 & \alpha^2 & (\alpha^2)^2 & \cdots & (\alpha^2)^{n-1} \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
1 & \alpha^{n-1} & (\alpha^{n-1})^2 & \cdots & (\alpha^{n-1})^{n-1}
\end{pmatrix}
$$
先に書いておいてなんですがVandermonde行列は、次の問題設定のセクションでは使いません。でも問題設定の後のセクションから使います。
**ノート:** 後で逆行列 $\vand^{-1}$ を使うことになるので、$\alpha$は何らかの体の値だと思ってください。具体的なイメージを持ちたいので使う体を先に教えてくれという方にお知らせすると、複素数体と有限体 $\finf_p$ を使います。
**疑問:** Vandermonde行列を既知の方は、この行列の中身には色んなバリエーションがある([Wikipediaのエントリ](https://en.wikipedia.org/wiki/Vandermonde_matrix))ことを知っているので、何故上の形を選ぶか？という *疑問* があるかと思います。これは解決されなければなりません。この要素の並びにしている理由は[セクション3](#セクション1で書いたVandermonde行列の形についての注意)ならびに[セクション6](#Vandermondeの行列式形に関する注意2)でお知らせします。
# 2. 問題設定: 多倍長整数とその積
多倍長整数型というのは、改めてなにかと言われるとやや説明がしづらいのですが、1バイト整数とか4バイト整数とか、ああいう最小値最大値がある数（数値型）ではなく、幾らでも大きい値が書けるような数値型です。
多倍長整数（型）では、例えば以下のような
$$
(12345678987654321)^2 =
152415789666209420210333789971041
$$
大きな数字や、その積を考えます。
この記事では、**多倍長整数型がない環境に、自力で新たに実装しようというストーリー**を考えます。PythonやRubyやHaskell([`Integer`型](https://downloads.haskell.org/~ghc/9.0.1/docs/html/libraries/ghc-bignum-1.0/GHC-Num-Integer.html)）やJava（[`BigInteger`型](https://docs.oracle.com/javase/jp/8/docs/api/java/math/BigInteger.html)）は言語組み込みの多倍長整数を持ちますが、一方でCやC++やRustだと多倍長整数型は組み込みの型ではないです。とはいえ提供ライブラリがあるのでそれを使えば良いので、別に困りはしないです。
現実的には自分で頑張るよりも、きっちり実装してくれているライブラリ、GNU GMP https://gmplib.org/ などを使うのが筋です。
ということは明記しつつ、それでも何とかやってみましょう。
## 多倍長整数の実現
多桁の数字をどう表現するか考えてすぐ思いつくのは、数を要素とするリストによって表現する方法かと思います。ここではそれと実質的に同じことをする「1変数多項式」による表現を採用します。つまり
$$
1234 \Rightarrow 1 x^3 + 2 x^2 + 3 x + 4
$$
のようにして、$x = 10$を代入すると元の数に戻るような1変数$x$上の多項式として記号的に表現する、ということです。
どうしてもデータ構造的に捉えたい方は、リストや辞書型で
```python
{ x^3: 1, x^2: 2, x^1: 3, x^0: 4 }
```
のように内部的に表現すると思ってください。
多項式での表現を前提にすると、乗算はどうなるか見てみます。
一例として、
$$56 \cdot 78 = 4368$$
という何の変哲もない計算を、多項式の乗算として処理します:
$$
\begin{array}{l}
(5 x + 6) \times (7 x + 8) = 35x^2 + 82x + 48.
\end{array}
$$
**記号の注意:** この記事では、多項式の積を $\times$ で表します。
結果の多項式を左から順番に読んで $35 82 48$ が答えかというと、そうではなくて、$x = 10$ によって評価してあげる必要があり、そうすると
$$
3500 + 820 + 48 = 4368
$$
として正しい結果を得ることができます。
いや、$x = 10$で評価するって言っても、評価した後の数字を表現する数値型がないという話じゃないの、というのはその通りです。
実際には、下の桁から上に向かって繰り上がり処理をして、多項式のままで「正しい形（係数が9以下の形)」に変形します。つまり:
1. $48 \to 4x + 8$にして、$35x^2 + (82 + 4)x + 8$にする。
2. $86x \to 8x^2 + 6x$として、$(35 + 8)x^2 + 6x + 8$にする。
3. $43x^2 \to 4x^3 + 3x^2$として、$4x^3 + 3x^2 + 6x + 8$にする。
こういう感じでやれば、積も実行できそうです。
**実装の注意:**
積を実行したい多項式の次数が巨大になるにつれ、結果の多項式の係数に「相当」巨大な数が出てきます。この係数が実装言語で表現可能な整数の最大値を超えていると、そもそも上のような正規化以前の問題として、正しい多項式が表現できなくなります。この辺をどう解決するかは実用上は重要なのですが、今回は実装はそこまで頑張りたくないのでこのトピックはskipします。
## 次のセクションへのつなぎ
多倍長整数型の積は定義できました。
今の実装は、多項式の積をナイーブに実行するものになっています。
この実装を高速化したいです。
というのは、２つの多項式の積は素朴にやると$O(N^2)$ (Nは次数)の計算になって、$N$が数十万とか数百万になると、しっかり遅くなるんですよね。
一方でこれから見る方法で $O(N \log N)$ まで持っていけて、[実験のセクション](#7-実装)でも見ますが、実行速度に顕著に差が出ます。
# 3. 高速化の準備: 畳み込み積とアダマール積
このセクションでは、高速化の準備として、２つの1変数$x$上の多項式 $p_1, p_2$ について、以下のことを示します：
$$
\vand^{-1} \bigl((\vand \cdot p_1) \ast (\vand \cdot p_2)\bigr) = p_1 \times p_2.
$$
失礼しました。記号が沢山出てきてしまいました。各記号は
* $\cdot$ はふつうの行列積 ← これは知っている
* $\ast$ は[Hadamard（アダマール）積](https://en.wikipedia.org/wiki/Hadamard_product_(matrices)) ← ???
* $\times$ は多項式の積、格好よく言うと[Convolution(畳み込み)積](https://en.wikipedia.org/wiki/Convolution) ← これも知っている
となっています。同じことですが、以下の図式を可換にしたいと言っています:
$$
\begin{array}{cccc}
\text{多項式世界} & p_1,\ p_2 & \overset{\times}{\longrightarrow} & p_1 \times p_2 \\
& \downarrow_{\vand} & & \uparrow_{\vand^{-1}} \\
\text{行列世界} & \vand \cdot p_1,\ \vand \cdot p_2 & \overset{\ast}{\longrightarrow} & (\vand \cdot p_1) \ast (\vand \cdot p_2)
\end{array}
$$
この定理には <font size="6">Convolution theorem</font> という格好良さげな名前がついているのですが（Link [1](https://en.wikipedia.org/wiki/Convolution_theorem) [2](https://ccrma.stanford.edu/~jos/sasp/Convolution_Theorem_DTFT.html))、Hadamard積が何者かわからないとなんとも言い難いので、とりあえず各積について見ていきます。
## ３つの積について
行列積、Hadamard積、Convolution積について一応述べます。
行列積は分かっているものとしますが、上の等式をみると、Vandermonde行列と多項式の積をとっています。いつのまに多項式が行列になったのかというと、これは記号をabuseしていまして、正確には多項式を次のように列ベクトルと同一視しています:
$$
k_3 x^3 + k_2 x^2 + k_1 x + k_0 \cong
\begin{pmatrix}
k_0 \\
k_1 \\
k_2 \\
k_3
\end{pmatrix}.
$$
さて問題の行列のHadamard積ですが、これはサイズの一致する２つの行列上の積で、point-wise（各点毎）に積を実行するタイプのものです:
$$
\begin{pmatrix}
a_1 & a_2 & a_3 \\
a_4 & a_5 & a_6
\end{pmatrix} \ast
\begin{pmatrix}
b_1 & b_2 & b_3 \\
b_4 & b_5 & b_6
\end{pmatrix} =
\begin{pmatrix}
a_1 b_1 & a_2 b_2 & a_3 b_3 \\
a_4 b_4 & a_5 b_5 & a_6 b_6
\end{pmatrix}
$$
これだけです。この積にHadamard積という格好良い名前がついているんだなあ以上の理解は、少なくとも今回は不要です。
**なんだ随分簡単だなあ** というとその感想そのものが一つのミソになっていて、Vandermonde行列をかけて行列の世界に移すと、（今まさに高速化しようとしている比較的重たい計算である）多項式の積が、要素別積に化けるというのがConvolution theoremのエッセンスです。
最後に、Convolution積ですが、これは多項式の積の格好良い呼び方です。以上なのですが、それだけだと寂しいので、一応よく見かけるようなformalな書き方をやっておくと:
$$
\sum^{i=K}_{i=0} n_i x^i \times \sum^{j=K}_{j=0} m_j x^j =
\sum_{0 \leq k \leq 2K}
\bigl(\sum_{i+j=k} n_i m_j\bigr) x^k
$$
となり、右辺は、２つの$K$次多項式（左辺）の畳み込み積を実行して得られる多項式での **$x^k$の係数** が
$$
n_k m_0 + n_{k-1} m_1 + \cdots + n_1 m_{k-1} + n_0 m_k
$$
になると言っています。
## Convolution theoremが 2次多項式の積 で成立することの確認
どのサイズのVandermonde行列でも成り立つのですが、具体的に話をみたいので、以下の２つの二次多項式について考えることにします:
$$
\begin{array}{l}
p_1 = a_2 x^2 + a_1 x + a_0 \cong \begin{pmatrix}
a_0 \\
a_1 \\
a_2 \\
0 \\
0
\end{pmatrix},\
p_2 = b_2 x^2 + b_1 x + b_0 \cong \begin{pmatrix}
b_0 \\
b_1 \\
b_2 \\
0 \\
0
\end{pmatrix}
\end{array}
$$
**待てい。どこから0が出てきた**ということについては急いで説明させてください。$p_1$と$p_2$の畳込み積の結果を格好つけて行列表示すると:
$$
p_1 \times p_2 =
\left(
\begin{array}{l}
a_0 b_0 \\
a_1 b_0 + a_0 b_1 \\
a_2 b_0 + a_1 b_1 + a_0 b_2 \\
a_2 b_1 + a_1 b_2 \\
a_2 b_2
\end{array}
\right)
$$
となります。ここのサイズと揃えたい、すなわち2次多項式の積の結果が4次多項式になるので、それに併せて$p_1$と$p_2$も行列で見るときには0でpaddingしています。そうしないと我々のconvolution theoremにおける行列の型が左辺と右辺で不一致するので……
あとはConvolution theoremが成立していることを確認したいのですが、実際には次の形の等式を示します:
$$
(\vand \cdot p_1) \ast (\vand \cdot p_2) = \vand \cdot (p_1 \times p_2) \tag{Conv}
$$
示すといっても、各行ごと $1 \leq \ell \leq 5$ に一致していることを確認するだけです。
### (Conv)の左辺の計算は多項式の畳み込み積
$$
\vand \times p_1 =
\begin{pmatrix}
1 & 1 & 1 & 1 & 1 \\
1 & \alpha & \alpha^2 & \alpha^3 & \alpha^4 \\
1 & \alpha^2 & (\alpha^2)^2 & (\alpha^2)^3 & (\alpha^2)^4 \\
1 & \alpha^3 & (\alpha^3)^2 & (\alpha^3)^3 & (\alpha^3)^4 \\
1 & \alpha^4 & (\alpha^4)^2 & (\alpha^4)^3 & (\alpha^4)^4 \\
\end{pmatrix}
\begin{pmatrix}
a_0 \\
a_1 \\
a_2 \\
0 \\
0
\end{pmatrix}
$$
です。$\lambda = \alpha^{\ell-1}$ と置くと、(Conv)の左辺の $1 \leq \ell \leq 5$ 行目は
$$
\begin{array}{l}
(
a_0 + a_1 \lambda^1 + a_2 \lambda^2 + 0 \lambda^3 + 0 \lambda^4
) \\
\qquad \times \\
(
b_0 + b_1 \lambda^1 + b_2 \lambda^2 + 0 \lambda^3 + 0 \lambda^4
) \\
\qquad = \\
a_0 b_0 + (a_1 b_0 + a_0 b_1) \lambda^1 + (a_2 b_0 + a_1 b_1 + a_0 b_2) \lambda^2 + (a_2 b_1 + b_1 a_2) \lambda^3 + (a_2 b_2) \lambda^4
\end{array}
$$
になります。
これは、一変数$\lambda$上の多項式の計算としてみることができます。
勿論、実際には$\lambda$は変数ではないです。気持ちの話です。
### (Conv)の右辺の計算は列ベクトル表現の多項式化
$$
a_0 b_0 + (a_1 b_0 + a_0 b_1) \lambda^1 + (a_2 b_0 + a_1 b_1 + a_0 b_2) \lambda^2 + (a_2 b_1 + b_1 a_2) \lambda^3 + (a_2 b_2) \lambda^4
$$
これは列ベクトルでの表現を、一つの多項式に読み替えただけになります。
### セクション1で書いたVandermonde行列の形についての注意
上の話は、次のVandermonde行列形でも成立します（確認してみてください):
$$
\begin{pmatrix}
1 & \alpha & \alpha^2 & \alpha^3 & \alpha^4 \\
1 & \alpha^2 & (\alpha^2)^2 & (\alpha^2)^3 & (\alpha^2)^4 \\
1 & \alpha^3 & (\alpha^3)^2 & (\alpha^3)^3 & (\alpha^3)^4 \\
1 & \alpha^4 & (\alpha^4)^2 & (\alpha^4)^3 & (\alpha^4)^4 \\
1 & \alpha^5 & (\alpha^5)^2 & (\alpha^5)^3 & (\alpha^5)^4 \\
\end{pmatrix} \tag{$\star$}
$$
一方で、
$$
\begin{pmatrix}
1 & 1 & 1 & 1 & 1 \\
\alpha & \alpha^2 & \alpha^3 & \alpha^4 & \alpha^5 \\
\alpha^2 & (\alpha^2)^2 & (\alpha^2)^3 & (\alpha^2)^4 & (\alpha^2)^5 \\
\alpha^3 & (\alpha^3)^2 & (\alpha^3)^3 & (\alpha^3)^4 & (\alpha^3)^5 \\
\alpha^4 & (\alpha^4)^2 & (\alpha^4)^3 & (\alpha^4)^4 & (\alpha^4)^5 \\
\end{pmatrix},
\begin{pmatrix}
\alpha & \alpha^2 & \alpha^3 & \alpha^4 & \alpha^5 \\
\alpha^2 & (\alpha^2)^2 & (\alpha^2)^3 & (\alpha^2)^4 & (\alpha^2)^5 \\
\alpha^3 & (\alpha^3)^2 & (\alpha^3)^3 & (\alpha^3)^4 & (\alpha^3)^5 \\
\alpha^4 & (\alpha^4)^2 & (\alpha^4)^3 & (\alpha^4)^4 & (\alpha^4)^5 \\
\alpha^5 & (\alpha^5)^2 & (\alpha^5)^3 & (\alpha^5)^4 & (\alpha^5)^5
\end{pmatrix}
$$
のように、左端が1で揃っていないVandermonde形では、左辺と右辺が一般には一致しません。確認してみてください。
ということは、セクション1のVandermonde行列の形ではなくて、$(\star)$でも良いのではないか?ということになります。Convolution theoremを成立させるところまでで言えばそうです。
しかし、$(\star)$は、以降で見ていく高速化のところで少し問題になります。具体的には[ここ](#Vandermondeの行列式形に関する注意2)で説明しますが、結局のところ、セクション1のVandermonde行列の形を選ぶと色々都合が良いということになります。
## 一般化
今回は二次多項式を具体的に使って計算に基づいて確認しましたが、任意次数でも明らかに成立する議論なので、任意サイズの行列でconvolution theoremが成立します。
フーリエ変換では複素数体上の話で同じことが展開されるので、どこまで一般化出来るか分かりづらいのですが、別に複素数体でなくても、どの体（または可換環）上のVandermonde行列を用いても、**(Conv)の形** は成立します。
注意点として、(Conv)の形はもともと示したかった convolution theorem:
$$
\vand^{-1} \bigl((\vand \cdot p_1) \ast (\vand \cdot p_2)\bigr) = p_1 \times p_2.
$$
とは違います。行列の世界に送る時に使う行列の逆行列がとれる必要があるのですが、その保証のためにVandermonde行列を使っていて、これが一つのポイントになります。
また、Vandermonde行列の旨味はそれだけに留まらず、以下で見るように分割統治法を用いた行列積の高速化の鍵ともなっています。
# 4. 高速Vandermonde行列積
ここまでで、convolution theoremを経由することで、多項式の積（畳み込み積）を実行する別の手段があることをみました。
すなわち、$\vand \cdot \hyphen$ によって行列の世界に送って、色々やったあとに$\vand^{-1} \cdot \hyphen$ で戻せばよいということです。
ここから高速フーリエ変換の肝である分割統治法を、Vandermonde行列積の高速実行という視点からみていきます。
まずは行列の世界に送る側の$\vand \cdot \hyphen$を考えますが、
具体的な計算も実行したいので、例として、$n = 8$のVandermonde行列 $\vand_{\alpha, n}$ を考えます。つまり
$$
\begin{pmatrix}
r_0 \\
r_1 \\
r_2 \\
\vdots \\
r_7
\end{pmatrix} =
\begin{pmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \alpha & \alpha^2 & \cdots & \alpha^7 \\
1 & \alpha^2 & (\alpha^2)^2 & \cdots & (\alpha^2)^7 \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
1 & \alpha^7 & (\alpha^7)^2 & \cdots & (\alpha^7)^7
\end{pmatrix} \cdot
\begin{pmatrix}
x_0 \\
x_1 \\
x_2 \\
\vdots \\
x_7
\end{pmatrix}.
$$
勿論これをそのまま計算しても良いのですが、$x$の偶数添字 ($x_0, x_2, x_4, x_6$) と奇数添字 ($x_1, x_3, x_5, x_7$) に分解して計算してみます（天下り的ではありますが）。
偶数添字部分の計算では、$x$について偶数部を取り出して、それに対応する列ベクトルたちを $\vand$ から取り出して、左からかけます:
$$
\begin{pmatrix}
e_0 \\
e_1 \\
e_2 \\
\vdots \\
e_7
\end{pmatrix} =
(\mathcal{E} =)
\begin{pmatrix}
1 & 1 & 1 & 1 \\
1 & \alpha^2 & \alpha^4 & \alpha^6 \\
1 & (\alpha^2)^2 & (\alpha^2)^4 & (\alpha^2)^6 \\
\vdots & \vdots & \vdots & \vdots \\
1 & (\alpha^7)^2 & (\alpha^7)^4 & (\alpha^7)^6 \\
\end{pmatrix}
\cdot
\begin{pmatrix}
x_0 \\
x_2 \\
x_4 \\
x_6
\end{pmatrix}
$$
同様にして奇数添字部分だけとってくると
$$
\begin{pmatrix}
o_0 \\
o_1 \\
o_2 \\
\vdots \\
o_7
\end{pmatrix} =
(\mathcal{O} =)
\begin{pmatrix}
1 & 1 & 1 & 1 \\
\alpha & \alpha^3 & \alpha^5 & \alpha^7 \\
\alpha^2 & (\alpha^2)^3 & (\alpha^2)^5 & (\alpha^2)^7 \\
\vdots & \vdots & \vdots & \vdots \\
\alpha^7 & (\alpha^7)^3 & (\alpha^7)^5 & (\alpha^7)^7 \\
\end{pmatrix}
\cdot
\begin{pmatrix}
x_1 \\
x_3 \\
x_5 \\
x_7
\end{pmatrix}
$$
とした上で、$r_i = e_i + o_i$が明らかに成立します。
ここから「２つ」やることがあります。
## 工夫1
$\mathcal{O}$はかなり$\mathcal{E}$に似た構造をしているので、この発想を元に、まず$\mathcal{E}$で計算して、$\alpha$が足りない部分をHadamard積で補うことで、$\mathcal{O}$行列を消去します:
$$
\begin{pmatrix}
o_0 \\
o_1 \\
\vdots \\
o_6 \\
o_7
\end{pmatrix} =
\begin{pmatrix}
1 \\
\alpha \\
\alpha^2 \\
\vdots \\
\alpha^7
\end{pmatrix}
\ast
\left(
\mathcal{E}
\cdot
\begin{pmatrix}
x_1 \\
x_3 \\
x_5 \\
x_7
\end{pmatrix}
\right)
$$
## 工夫2
工夫1を経て、計算のコアが
$$
\mathcal{E} \cdot \begin{pmatrix}
y_0 \\
y_1 \\
y_2 \\
y_3
\end{pmatrix}
$$
の形をしていることが分かったので、この計算ステップを減らす工夫をします。
$\beta = \alpha^2$としてこの式を書き直すと:
$$
\begin{pmatrix}
1 & 1 & 1 & 1 \\
1 & \beta & \beta^2 & \beta^3 \\
1 & \beta^2 & (\beta^2)^2 & (\beta^2)^3 \\
1 & \beta^3 & (\beta^3)^2 & (\beta^3)^3 \\\hline
1 & \beta^4 & (\beta^4)^2 & (\beta^4)^3 \\
1 & \beta^5 & (\beta^5)^2 & (\beta^5)^3 \\
1 & \beta^6 & (\beta^6)^2 & (\beta^6)^3 \\
1 & \beta^7 & (\beta^7)^2 & (\beta^7)^3
\end{pmatrix}
\cdot
\begin{pmatrix}
y_0 \\
y_1 \\
y_2 \\
y_3
\end{pmatrix}
$$
（横棒は単に視覚的なもので、値として深い意味があったりはしません。）
***ここで $\alpha^8 = \beta^4 = 1$ が成立する $\alpha$ を選んでおけば*** 嬉しいことが起こります。つまり
$$
\mathcal{E} = \begin{pmatrix}
\vand_{\beta, 4} \\\hline
\vand_{\beta, 4}
\end{pmatrix}
$$
が成立して、ブロック行列積により
$$
\mathcal{E} \cdot Y =
\begin{pmatrix}
\vand_{\beta, 4} \cdot Y \\
\vand_{\beta, 4} \cdot Y
\end{pmatrix}
\quad
\text{where}
\quad
Y = \begin{pmatrix}
y_0 \\
y_1 \\
y_2 \\
y_3
\end{pmatrix}
$$
となります。4x4のVandermonde行列積を一度やって、それをコピーすれば良いので、もとの8x8の行列積の計算を半分のサイズの行列積の計算に帰着できました。
まとめると、まず
$$
\begin{pmatrix}
v_0 \\
v_1 \\
v_2 \\
v_3
\end{pmatrix} =
\vand_{\alpha^2, 4} \begin{pmatrix}
x_0 \\
x_2 \\
x_4 \\
x_6
\end{pmatrix}
$$
の計算をして、次に
$$
\begin{pmatrix}
w_0 \\
w_1 \\
w_2 \\
w_3
\end{pmatrix} =
\vand_{\alpha^2, 4} \begin{pmatrix}
x_1 \\
x_3 \\
x_5 \\
x_7
\end{pmatrix}
$$
としたあとで、以下が成立します:
$$
\vand_{\alpha, 8} \cdot
\begin{pmatrix}
x_0 \\
x_1 \\
x_2 \\
x_3 \\
x_4 \\
x_5 \\
x_6 \\
x_7
\end{pmatrix} =
\begin{pmatrix}
v_0 \\
v_1 \\
v_2 \\
v_3 \\
v_0 \\
v_1 \\
v_2 \\
v_3
\end{pmatrix} +
\begin{pmatrix}
1 \\
\alpha \\
\alpha^2 \\
\alpha^3 \\
\alpha^4 \\
\alpha^5 \\
\alpha^6 \\
\alpha^7
\end{pmatrix}
\ast
\begin{pmatrix}
w_0 \\
w_1 \\
w_2 \\
w_3 \\
w_0 \\
w_1 \\
w_2 \\
w_3
\end{pmatrix}
$$
一般的な高速フーリエ変換を知っている方の頭の中には「あのどこからかやってくるマイナス1倍はどこ?」という疑問があるかもしれません。何を言っているか分からない方は、そこは自然にカバーできているので、安心して次のセクションまで読み飛ばしてください。
$\alpha^8 = 1$であり、かつ$\alpha, \alpha^2, \ldots, \alpha^7, \alpha^8$の全てが異なっている場合には（後でみますが、これはVandermonde行列が逆行列を持つために必要な条件です）、$\alpha^4 = -1$でなければならないので、以下が成立します:
$$
\begin{pmatrix}
1 \\
\alpha \\
\alpha^2 \\
\alpha^3 \\
\alpha^4 \\
\alpha^5 \\
\alpha^6 \\
\alpha^7
\end{pmatrix} =
\begin{pmatrix}
1 \\
\alpha \\
\alpha^2 \\
\alpha^3 \\
-1 \\
-\alpha \\
-\alpha^2 \\
-\alpha^3
\end{pmatrix}
$$
一般的な高速フーリエ変換の説明で$-1$が出てくるのは、まさにこの部分です。
## 再帰構造
既に明らかになっていますが、結局のところ$\vand_{\alpha, 8}$の計算が、うまいこと$\vand_{\alpha^2, 4}$という縦横半分サイズの行列積に置き換わりました。
ここで話を止める必要はなくて、この4x4の行列積すらも、再帰的に $\vand_{\alpha^4, 2}$ を用いた形に変形できます。
より大きい$n=32$から話をはじめると、$\vand_{\alpha, 32} \cdot \hyphen$ の形の計算は、$\vand_{\alpha^2, 16}$を用いた計算に書き換えることができて、この計算が更に$\vand_{\alpha^4, 8}$を用いた計算に書き換えることができて……という次第です。
# 5. Vandermonde行列上で要求される性質
高速Vandermonde変換では、行列のサイズを縦横半分にしていきたいので、スタート地点となる行列の形は、2のべき乗の
$$
\vand_{\alpha, 2^n} =
\begin{pmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \alpha & \alpha^2 & \cdots & \alpha^{2^n-1} \\
1 & \alpha^2 & (\alpha^2)^2 & \cdots & (\alpha^2)^{2^n-1} \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
1 & \alpha^{2^n-1} & (\alpha^{2^n-1})^2 & \cdots & (\alpha^{2^n-1})^{2^n-1}
\end{pmatrix}
$$
です。
また、$\alpha^{2^n} = 1$が成立する必要もあります。この条件を ***(V1)*** と呼びます。
Vandermonde行列が逆行列を持つためには、2列目の値が全て異なること、すなわち $1, \alpha, \alpha^2, \alpha^3, \ldots, \alpha^{2^n-1}$ の全てが異なる値を取らなければなりません。
[これはVandermonde行列の行列式が差積で与えられるという基本的な結果](https://en.wikipedia.org/wiki/Vandermonde_matrix)から明らかです。この条件を ***(V2)*** と呼びます。
それでは、両方の条件を満たす $\alpha$ を上手いこと作るための２つの例を見てみます。
## 複素数体による実現
１つ目の方法は、一般的なフーリエ変換で用いられるものですが、複素数体上で $\alpha = e^{-i \frac{2\pi}{n}}$ とする方法です。
[複素平面上でオイラーの公式を可視化するよくある図](https://ja.wikipedia.org/wiki/%E3%82%AA%E3%82%A4%E3%83%A9%E3%83%BC%E3%81%AE%E5%85%AC%E5%BC%8F)を思い出してもらいたいのですが、あの図のおかげで$\alpha^n=1$になりつつも、$\alpha, \alpha^2, \ldots, \alpha^n$は全て異なる値をとるというのが想起できます。よって簡単に (V1) と (V2) を満足できます。
この方法でVandermonde行列を具体化した場合に、上の方法で高速化したものが、離散高速フーリエ変換と呼ばれるものと一致しています。
じゃあ複素数体を使えば良いね、ということになるのですが、一般的なプログラミング言語の（組み込みの）複素数型は、`float型`とか`double型`とか`long double型`の、精度に限りある数値型が使われているため（参考 [1](https://en.cppreference.com/w/cpp/numeric/complex))、桁数に限りのない多倍長演算を実行するには誤差などの心配がつきまといます。
この点を考慮して、別の方法である有限体による実現もみておきます。
## 有限体による実現
別の方法による実現では、素数$p$で整数を割って得られるような[体](https://en.wikipedia.org/wiki/Finite_field) $\finf_p$があります。
$\finf_p$の定義は簡単で、
* この体には$\set{0, 1, 2, \ldots, p-1}$という$p$個の整数がおり
* 足し算や掛け算は整数として計算したあとに$p$で余りを求める
* 除算 $x / y$ については $y \neq 0$ であれば、$y \cdot z = 1 \pmod p$ となる $z$ が、$p$が素数であることより存在する(参考リンク [1](https://ja.wikipedia.org/wiki/%E3%83%A2%E3%82%B8%E3%83%A5%E3%83%A9%E9%80%86%E6%95%B0)) ので、それを使って$x \cdot z$をしたあとに$p$で余りを求める
とすれば良いです。
例えば $p = 17$ の場合は、$3 \cdot 9 = 10 \pmod{17}$ となり、$3 / 9 = 3 * 2 = 6$ になります。2が9の逆数なのは $2 \cdot 9 = 1 \pmod{17}$ から分かります。$(3 / 9) \cdot 9 = 6 * 9 = 3$ が成立するので、除算も出来ていそうです。
とくに $p=17=2^4+1$ のような $p = k2^n+1$ で表現されている素数については、フェルマーの小定理(参考 [1](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A7%E3%83%AB%E3%83%9E%E3%83%BC%E3%81%AE%E5%B0%8F%E5%AE%9A%E7%90%86) [2](https://manabitimes.jp/math/680)) から、任意の $1 < x < p$ について
$$
x^{k2^n} = (x^k)^{2^n} = 1 \pmod p
$$
が成立することが分かるので、$\alpha = x^k$ として採用すれば、(V1)を満たします。実際、$p = 17$のときには、$(3^1)^{16} = 1$が成立します。特にこの場合は
$$
\pmod{17}:
\begin{array}{llll}
3^1 = 3, & 3^2 = 9, & 3^3 =10, & 3^4 = 13, \\
3^5 = 5, & 3^5 = 15, & 3^6 = 11, & 3^7 = 16, \\
3^8 = 14, & 3^9 = 8, & 3^{10} = 7, & 3^{11} = 4, \\
3^{13} = 12, & 3^{14} = 2, & 3^{15} = 6, & 3^{16} = 1
\end{array}
$$
により(V2)も満たすので、高速Vandermonde行列積に使えます。
しかし$p=17$で$\alpha=3$として(V2)が成立するのは「偶然」です。
これを確認するために、別の例として、$p = 3 \cdot (2^5) + 1 = 97$を考えてみます。
条件(V1) $3^{3 \cdot 2^5} = {27}^{2^5} = 1 \pmod{97}$ がフェルマーの小定理から成立するのですが、一方で $27^{2^4 = 16} = 1 \pmod{97}$ が成立するので、(V2)が成立しません。
強調しますが、フェルマーの小定理は(V1)を保証してくれるだけです。
(V2)も満たす値の探索は、整数論の言葉で言えば、 $\pmod p$ における原始根を求める問題を解けば十分で、ということになると整数論に関する議論を使いたくなります。
このノートでは整数論の話は端折りますが、[代表的な効率の良い手法](https://en.wikipedia.org/wiki/Primitive_root_modulo_n#Finding_primitive_roots)は確立しているので、物凄く困ったことには一応なりません。
# 6. 高速逆Vandermonde行列積
ここまではVandermonde行列によって、多項式を行列の世界に送る側を高速化したのですが、最後に多項式の世界に戻ってくる際にはVandermonde行列の逆行列をかける必要があります。
逆行列についても高速化のための構造が備わっていなければ困るのですが、高速逆Vandermonde行列積が実行できることが、次のことから分かります:
$$
\vand_{\alpha, n}^{-1} =
\frac{1}{n}
\begin{pmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \overline{\alpha} & \overline{\alpha}^2 & \cdots & \overline{\alpha}^{n-1} \\
1 & \overline{\alpha}^2 & (\overline{\alpha}^2)^2 & \cdots & (\overline{\alpha}^2)^{n-1} \\
& & & \vdots \\
1 & \overline{\alpha}^{n-1} & (\overline{\alpha}^{n-1})^2 & \cdots & (\overline{\alpha}^{n-1})^{n-1}
\end{pmatrix}
$$
ここで $\overline{\alpha}$ は $\alpha$ の積逆元で、複素数体で $\alpha = e^{-i\frac{2 \pi}n}$ の場合は複素共役$\overline{\alpha} = e^{i\frac{2 \pi}n}$ を取れば良いですし、$\finf_p$ で $\alpha^n = 1$ ならば $\overline{\alpha} = \alpha^{n-1}$ です。
$\frac{1}{n}$がついているのがちょっと気になりますが、これがないと
$$
\vand \cdot \vand^{-1} = \begin{pmatrix}
n & 0 & 0 & \cdots \\
0 & n & 0 & \cdots \\
0 & 0 & n & \cdots \\
\vdots & \vdots & \vdots & \ddots \\
\end{pmatrix}
$$
となるので仕方がないところです。ですが分割統治法を行う上では問題になりません。
$\alpha$が(V1)と(V2)を満たす時、$\overline{\alpha}$もこれを満たすことに注意してください。すなわち、全く同じ方法で分割統治法により行列サイズを再帰的に半分にしていくことで、逆Vandermonde行列積も実行出来ます（最後に $\frac{1}{n}$ を全要素にかける必要はありますが）。
## Vandermondeの行列式形に関する注意2
もし我々がVandermonde行列の形として
$$
\mathcal{M} =
\begin{pmatrix}
1 & \alpha & \alpha^2 & \alpha^3 & \alpha^4 \\
1 & \alpha^2 & (\alpha^2)^2 & (\alpha^2)^3 & (\alpha^2)^4 \\
1 & \alpha^3 & (\alpha^3)^2 & (\alpha^3)^3 & (\alpha^3)^4 \\
1 & \alpha^4 & (\alpha^4)^2 & (\alpha^4)^3 & (\alpha^4)^4 \\
1 & \alpha^5 & (\alpha^5)^2 & (\alpha^5)^3 & (\alpha^5)^4 \\
\end{pmatrix}
$$
を採用した場合は、既にみたようにconvolution theoremは成立します。
ただし、この行列の逆行列は、各要素を逆元で置き換えた
$$
\mathcal{N} =
\begin{pmatrix}
1 & \overline{\alpha} & \overline{\alpha}^2 & \cdots & \overline{\alpha}^{n-1} \\
1 & \overline{\alpha}^2 & (\overline{\alpha}^2)^2 & \cdots & (\overline{\alpha}^2)^{n-1} \\
& & & \vdots \\
1 & \overline{\alpha}^n & ({\overline{\alpha}}^n)^2 & \cdots & ({\overline{\alpha}}^n)^{n-1}
\end{pmatrix}
$$
とはなりません。
例えば $\mathcal{M} \cdot \mathcal{N}$ の左上の要素を考えてみると、これは
$$
\begin{array}{l}
1 + {\alpha} + {\alpha}^2 + \cdots + {\alpha}^{n-1}
\end{array}
$$
となり、一方で
$$
\begin{array}{lll}
& (1 - \alpha)(1 + \alpha^2 + \alpha^3 + \cdots + \alpha^{n-1}) \\
= & 1 - \alpha^n & (\text{by calculus on polynomials}) \\
= & 0 & (\text{by $\alpha^n = 1$})
\end{array}
$$
で、$\alpha \neq 1$なので、左上の要素は0となります。
よって
$$
\mathcal{M} \cdot \mathcal{N} =
\begin{pmatrix}
0 & \cdots \\
\vdots & \ddots
\end{pmatrix}
$$
となり、明らかに単位行列ではなく、$\mathcal{N}$は$\mathcal{M}$の逆行列ではないです。
# 7. 実装
複素数体を使った実装はたくさんあると思うので、$\finf_p$ の方法でやってみます。こちらの実装もたくさんあるわけですが。
そこそこ大きな多項式の積を考えたいので、次の設定を使います:
$$
p = 7 \cdot 2^{20} + 1,\quad \alpha = 3^7 = 2187.
$$
この $p = 7340033$ が素数かどうかは、自分で頑張っても良いですし、[Wolframalpha](https://www.wolframalpha.com/input/?i=7340033&lang=ja)でも確認できます。
このときは、多項式積実行後の次数が$2^{20}$までなら表現できるので、今回は $\frac{2^{20}}{3} \sim{350k}$ 次数の多項式の積を実行してみます。
これはすなわち、約 35万桁 の整数同士の積を実行することになります。
C++で書いたプログラムはここにおいてあります: https://gist.github.com/yuezato/f965e51aac2ffd66924de9252b9bd897
コンパイルして実行してみると
```
$ g++ -std=c++11 -O3 fast_vt.cpp -Wall
$ ./a.out
Alpha(2187) = 3^7
P(7340033) = (7)*(2^20)+1
Multiply two polynomials whose degrees are 349525
Start calc
End calc. Elapsed: 377ms
Start naive
```
などが得られ、まずこのノートに書いた手法だと400ms秒足らずで自分の環境だと実行できていることが分かります。
その次の `Start naive` ですが、実装の正しさの確認のために、ナイーブな多項式積を実行して、これを高速化手法の結果と比較しようとするものの、数分待った程度では結果が全く返ってきません。
ということで十分高速化できているかなと思います。
念のためにもう少し小さいパラメータで実行してみると……
```
Alpha(243) = 3^5
P(163841) = (5)*(2^15)+1
Multiply two polynomials whose degrees are 10922
Start calc
End calc. Elapsed: 7ms
Start naive
End naiive. Elapsed: 1376ms
Ok
```
```
Alpha(125) = 5^3
P(786433) = (3)*(2^18)+1
Multiply two polynomials whose degrees are 87381
Start calc
End calc. Elapsed: 69ms
Start naive
End naiive. Elapsed: 92637ms
Ok
```
```
Alpha(177147) = 3^11
P(5767169) = (11)*(2^19)+1
Multiply two polynomials whose degrees are 174762
Start calc
End calc. Elapsed: 207ms
Start naive
End naiive. Elapsed: 530223ms
Ok
```
ということで相当高速化できていることが分かりました。
# おわりに
高速フーリエ変換は、フーリエ変換の式の導出が先立っていて、その出てきた式に対して「上手いこと」分割統治法が働くので、とにかく高速に実行できるという説明になりがちです。
それはそれで良いのですが
* 「なぜ」上手いこといくのか、その条件や構造が知りたい
* どのくらい脆い・あるいは奇跡的にうまくいく状況なのか?
* どういうケースであれば同じ様なことができるのか
といった類のことがどうしても気になるタイプの人達というのは一定居ると思っていて、自分がそれなのですが、その人達が納得するためには、もう少し一般的な立場で展開される構成が知りたいです。
そのような構成として、Vandermonde行列を中心にした説明はある程度うまくいくかなと思います。とはいえ、はじめにでも書いたのですが、このようなVandermonde行列で捉え直す話は新しい話では全くないです、ということは再度強調しておきます。
しかし一般化という観点からはあんまり踏み込めていないので、離散コサイン変換と関連があるChebyshev-Vandermonde変換(リンク [1 pdf](
http://www.icsi.berkeley.edu/ftp/global/pub/techreports/1993/tr-93-024.pdf))とその高速化まで綺麗に書けるのかなあ、なんかは調べてみたいのですが、結局使ったことがないので手がついていません。
他にやや気になっていることとしては、Convolution theoremまわりの証明を行列の言葉でもっと綺麗に書けないかなあ、というのがあります。
今回は多項式の畳み込み積を、行列の言葉や行列演算ではなく、普通の定義で済ませています。
一方で畳み込み積は、[Toeplitz行列](https://en.wikipedia.org/wiki/Toeplitz_matrix)を用いた行列積の形で定義できるので、こっちを経由するともう少し綺麗かつ一般的かつ行列の性質ベースの話で済んで、そうするとHadamard行列以外とのinteractionも自然に色々出せないかなあ〜、みたいなやつです。
何かご存知の方は教えてください。
Sign in with Wallet
Connect another wallet