快速傅里葉變換（FFT）求解多項式乘法

在我還會FFT的時候趕快寫下一篇博客留着以后看。。。。。。

FFT是用來求解多項式乘法，那么首先我們要知道多項式是啥。

\[A(x) = a_0+a_1x^1+a_2x^2+···+a_{n-1}x^{n-1} \]

這是個n-1次多項式（最高項是\(x^{n-1}\)），\(a_0,a_1,···a_{n-1}\)是它各項的系數，該多項式可以寫成：

\[A(x) = \sum_{j=0}^{n-1}a_jx^j \]

一個多項式可以通過一組系數所確定，而這組系數所組成的向量也叫做系數向量（如\(A(x)\)的系數向量為：[\(a_0 \ \ \ a_1\ ··· \ a_{n-1}\)]），通過系數向量表示一個多項式的方式叫做系數表示法

多項式乘法就是，我們已知兩個多項式：\(A(x),B(x)\)，分別是n-1次和m-1次多項式。

\[A(x) = \sum_{j = 0}^{n - 1}a_jx^j\\ B(x) = \sum_{j = 0}^{m-1}b_jx^j \]

現在講這兩個多項式相乘，得到一個最高n+m-1次多項式：

\[C(x) = \sum_{j = 0}^{n+m-1}c_jx^j \]

那么求解\(C(x)\)也就是求出它的系數向量：[\(c_0 \ \ c_1 \ \ ··· \ \ c_j \ \ ··· \ \ c_n+m-1\)]

如果直接遍歷\(A(x),B(x)\)的每一項的話是\(O(n*m)\)，我們需要更快的算法，就是我們所討論的：快速傅里葉變換(FFT)

前置芝士：

點值表示法

對於一個已知的多項式\(A(x)\)，將\(x_0\)代入可以得到一個確定的\(y_0\)，並且，可以將\((x_0,y_0)\)看做是坐標系上的一個點。進一步可以將任意多的互不相等的自變量代入\((x_1,x_2,···)\)，從而得到更多的點：\((x_1,y_1),(x_2,y_2),···\)

我們知道，通過兩個（不同的）點，可以解析一次多項式\((ax+b)\)，通過三個不同的點可以解析二次多項式\((ax^2+bx+c)\)，所以可以進一步的推出，通過n+1個不同的點，可以解出n次多項式

通過n+1個不同的點\(\left\{ (x_1,y_1),(x_2,y_2),··· \right\}\)，唯一確定一個n次多項式，該方式也叫做多項式的點值表示法

我們驚奇的發現兩個用點值表示的多項式相乘，復雜度是\(O(n)\)的！具體方法：\(C(x_i)=A(x_i)×B(x_i)\)，所以O\((n)\)枚舉\(x_i\)即可。

要是我們把兩個多項式轉換成點值表示，再相乘，再把新的點值表示轉換成多項式豈不就可以\(O(n)\)解決多項式乘法了！

……很遺憾，顯然，把多項式轉換成點值表示的朴素算法是\(O(n^2)\)的。另外，即使你可能不會——把點值表示轉換為多項式的朴素“插值算法”也是\(O(n^2)\)的。

但是可以發現，大整數乘法復雜度的瓶頸可能在“多項式轉換成點值表示”這一步（以及其反向操作），只要完成這一步就可以\(O(n)\)求答案了。如果能優化這一步就可以了。

於是乎介紹另外兩個前置芝士

離散傅里葉變換（DFT）以及離散傅里葉變換逆變換（IDFT），公式如下：

\[X_k = \sum_{j = 0}^{n-1}x_je^{\frac{2\pi i}{n}kj} (DFT)\\ X_j = \frac{1}{n}\sum_{k = 0}^{n-1}X_ke^{\frac{-2\pi i}{n}jk}(IDFT) \]

上面的兩個式子中\(i\)是虛根（\(i^2 = -1\)），關於這兩個式子的證明在文章最后會給出。

我們設置一些變量，將這兩個公式代入到我們之前的問題里。

設：\(w = e^{\frac{2\pi i}{N}}\)，其中\(N = n+m\)，關於變量\(w\)的性質先不討論，之后會有，再令

\[c_j = f(j),C_k = C(w^k) \]

將這些變量代入到DFT的公式里可得到：

\[C_k = C(w^k)=\sum_{j = 0}^{N-1}c_j(w^k)^j= \sum_{j=0}^{N-1}c_je_{\frac{2\pi i}{N}kj} \]

通過上式講DFT與我們目標多項式\(C(x)\)關聯了起來，自然也就可以通過IDFT來計算多項式的系數\(c_j\)：

\[c_j = \frac{1}{N}\sum_{k=0}^{N-1}C_k(w^{-j})^k=\frac{1}{N}\sum_{k=0}^{N-1}C_ke^{-\frac{2\pi i}{N}jk} \]

但是，由於不知道\(c_j\)的具體值（知道就完結撒花了），所以無法用DFT計算\(C(w^k)\)

不過，我們知道：\(C(x)=A(x)*B(x)\)，而且多項式\(A(x)B(x)\)的系數是已知的，所以可以“曲線救國”：

\[C(w^k)= A(w^k)*B(w^k) \]

通過該式計算\(C(w^k)\)，然后再使用IDFT計算\(c_j\)。

以上是FFT的概括，下面進行展開討論

一個芝士：N次單位根

為了更好的討論\(w = e^{\frac{2\pi i}{N}}\)，我們還是從歐拉公式說起

\[e^{i\theta} = \cos\theta+i\sin\theta \]

該公式的證明也在文章最后給出。

\(e^{i\theta}\)是一個復數，他的實部的值為\(\cos\theta\)，虛部的值為\(\sin\theta\)，所以，\(e^{i\theta}\)對應復數平面上的點與原點之間的距離恆為1，如圖

因為\(\sin\theta,\cos\theta\)是周期函數（周期為\(2\pi\)），所以可以得出：

\[e^{i(\theta+2\pi)=e^{i\theta}} \]

並且無論\(\theta\)取何值，\(e^{i\theta}\)都是分布在復平面的單位圓上，由於是單位元，所以\(\theta\)此時就等於是從1（起點）到\(e^{i\theta}\)在圓上的弧長。

上面簡要介紹了歐拉公式，下面回到我們定義的變量：\(w = e^{\frac{2 \pi i}{n}}\)，前面的概述中為了方便，\(w\)沒有加上下標，在此為更清楚的表示會給\(w\)加一個下標：\(w_n=e^{\frac{2\pi i}{n}}\)

從面歐拉公式的討論中可以指導，\(e^{2\pi i} = 1\)，相當於在單位圓上繞了一圈又回到了起點，而\(\frac{2\pi}{n}\)相當於是這個單位圓分成了\(n\)等分，而\(w_n^k\)就對應單位圓上第\(k\)個等分點（\(k=0,1,2,···，n-1\)）

這是個\(n=18\)的例子，並標出了兩個值\(\left\{ w_n^k|k=2,10 \right\}\)，通過前面的討論以及這個栗子，應該對\(w_n^k\)有了一個比較直觀的認識，說了這么多，應該能猜到N次單位根就是指的它了，關於N次單位根介紹幾個性質：

\[1. \ \ w_n^{k+\frac{n}{2}}=e^{\frac{2\pi i}{n}(k+\frac{n}{w})}=e^{\frac{2\pi i}{n}k}e^{i\pi} =w_n^k(\cos\pi+i\sin\pi)=-w_n^k \\2. \ \ w_{2n}^{2k} = w_n^k \\3. \ \ w_n^{-k} = 1·e^{\frac{i2\pi}{n}(-k)} = e^{2\pi i}e^{\frac{2 \pi i}{n}(-k)}=e^{\frac{2\pi i}{n}(n-k)}=w_n^{n-k} \]

性質1是消去引理，性質2是折半引理，性質3我忘了叫啥。。。

計算\(C(w_N^k)\)

首先，定義\(N = m+n\)（n,m分別是多項式\(A(x),B(x)\)的系數數量），前面已經說過無法直接使用DFT計算\(C(w_n^k)\)，需要先計算\(A(w_N^k),B(w_N^k)\)

單數如果直接把\(\left\{ w_N^k|k=0,1,2,···，N-1 \right\}\)簡單粗暴的代進去，時間復雜度依舊是\(O(N^2)\)，無論后面如何計算，也無法降低整個算法的復雜度，因為，需要使用的一些方法，來降低\(A(w_N^k),B(w_N^k)\) 的時間

for example，我們有個多項式：

\[H(x) = h_0+h_1x^1+h_2x^2+h_3x^3 \]

我們可以將它的奇數項與偶數項分開：

\[H(x)=(h_0+h_2x^2)+(h_1x^1+h_3x^3)=(h_0+h_2x^2)+x(h_1+h_3x^2)\\= H1(x^2)+xH2(x^2) \]

並且\(H1(x),H2(x)\)還可以遞歸的計算下去，該方法用到\(A(x),B(x)\)上就是：

\[A(x) = A1(x^2)+xA2(x^2)B(x) = B1(x^2)+xB2(x^2) \]

時間復雜度就變成了：

\[T(n) = 2T(\frac{n}{2})+O(n)=O(nlogn) \]

除此之外，還可以利用前面討論的單位根的性質來降低計算量，例如，\(0 \leq k < \frac{n}{2}\)：

\[A(w_n^k) = A1(w_n^{2k})+w_n^kA2(w_n^{2k})=A1(w_{\frac{n}{2}}^k)+w_n^kA2(w_{\frac{n}{2}}^k)\\A(w_n^{k+\frac{n}{2}}) = A1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}A2(w_n^{2k+n}) = A1(w_n^{2k})-w_n^kA2(w_n^{2k}) \\= A1(w_{\frac{n}{2}}^k)-w_n^kA2(w_{\frac{n}{2}}^k)) \]

這樣，計算\(k\)一般的取值范圍就可以得到整個取值范圍的結果，所以，將\(A(w_N^k),B(w_N^k)\)計算完畢之后，就可以通過簡單的乘法來計算\(C(w_N^k)\)：

\[C(w_N^k) = A(w_N^k)B(w_N^k),k=0,1,2,···，N-1 \]

求解\(C(x)\)系數

計算 這一步基本已經完成一大半了，只需要使用IDFT公式來求解\(C(x)\)的各項系數就可以了：

\[c_j = \frac{1}{N}\sum_{k=0}^{N-1}C_ke^{-\frac{2\pi i}{N}jk} = \frac{1}{N}\sum_{k=0}^{N-1}C_k(e^{-\frac{2\pi i}{N}j})^k=\frac{1}{N}\sum_{k=0}^{N-1}C_k(w_N^{-j})^k \]

上面式子中\(C_k = C(w_N^k)\)，另外，為了方面表示，通常會構造一個N-1次多項式\(D(x)\)，它的系數就是各個\(C_k\)的值：

\[D(x) = C_0 + C_1x^1+C_2x^2+···+C_{N-1}x^{N-1} = \sum_{k=0}^{N-1}C_kx^k \]

所以，最終\(c_j\)的計算公式可以寫成：

\[c_j = \frac{1}{N}\sum_{k=0}^{N-1}C_k(w_N^{-j})^k=\frac{1}{N}\sum_{k=0}^{N-1}C_k(w_N^{N-j})^k=\frac{1}{N}D(w_N^{N-j}) \]

上式第二部使用了單位根的性質3：

\[c_0 = \frac{1}{N}D(w_N^{N-0}) = \frac{1}{N}D(w_N^N)=\frac{1}{N}D(w_N^0) \]

總結

1.取N個單位根：[\(w_N^0 \ \ w_N^1 \ \ w_N^2 \ \ ··· \ \ w_N^{N-1}\)]

2.將它們分別代入到多項式\(A(x),B(x)\)中

3.計算\(C(w_N^k) = A(w_N^k)B(w_N^k)\)

4.構造一個新的多項式\(D(x) = C_0+C_1x^1+C_2x^2+···+C_{N-1}x^{N-1}\)

5.計算\(D(w_N^k)\)

6.計算多項式\(C(x)\)的系數\(c_j=\frac{1}{N}D(w_N^{N-j}),j=0,1,2,···，N-1\)

at last（未完待續）

• 歐拉公式證明
• DFT與IDTF證明