快速傅里葉變換（FFT）略解

前言

如果我們能用一種時間上比 \(O(n^2)\) 更優秀的方法來計算大整數（函數）的乘法，那就好了。快速傅里葉變換（FFT） 可以幫我們在 \(O(n\log n)\) 的時間內解決問題。

函數乘積

計算兩個大整數之積時，我們發現

\[(2x+3)(4x+5)=8x^2+22x+15\quad...(*)\\ 23\times45=1035\]

而如果我們把 \((*)\) 式右邊的每一位的系數看做一個數每位上的數碼，正好得到了 \(1035\)。事實上，對於所有的多項式乘法，以上規律同樣成立。
證明： （提示）考慮豎式乘法的過程，和多項式乘法的過程，它們的本質都是一樣的。

這樣，我們就把問題轉換為：計算兩個已知函數之積的函數的解析式。

復平面、單位圓

考慮 \(\sqrt{-9}\) 的值。

\[\begin{aligned}\sqrt{-9}&=\sqrt{-1}\times\sqrt9=3\sqrt{-1}.\end{aligned} \]

類似地，\(\forall N\in \Z_-\) 我們都可以用類似的方法得到 $$\sqrt{N}=\sqrt{-N}\times\sqrt{-1}$$
引入虛數單位 \(\text{i}\)，使 \(\text{i}^2=-1.\) 這樣我們就重新認識了數的范圍，從實數擴充到復數。

一復數 \(a+b\text{i}\) 中的 \(a,b\in\R\)，\(a\) 是它的實數部分，\(b\text{i}\) 是虛數部分。若 \(b=0\)，則它是實數。復數服從實數的大部分運算法則。
若兩個復數，它們的實數部分相等，虛數部分之和為 \(0\)，我們稱它們互為 共軛復數。

我們知道，數軸上的每個點與每個實數一一對應。類似地，我們可以使用 復平面 上的點表示復數。復平面與平面直角坐標系類似，它的 \(x\) 軸單位長度為 \(1\)，\(y\) 軸單位長度為 \(\text{i}\)。復平面上的點之橫縱坐標表示的數之和即為該點表示的數。比如，\((1,2)\) 表示 \(1+2\text{i}\)。

以圓點為圓心，\(1\) 為半徑，在復平面上作圓，如圖所示。這個圓叫 單位圓。

在圓上任取一點 \(A\)，過此點作 \(x\) 軸的垂線段，垂足為 \(B\)。設 \(∠AOB=\theta\)。易知$$OB=OA\times\cos\theta=\cos\theta\
AB=OB\times\sin\theta=\sin\theta$$\(\therefore A(\cos\theta,\sin\theta).\) 這個 \(\theta\) 叫做 \(A\) 的 輻角。

單位根及其性質

1. \(\forall k,n\in\N^*\) 有 \(\omega_n^k\times\omega_n^1=\omega_n^{k+1}\)；
   證明：
   
   如圖所示，過點 \(A(1,0)\) 作圓內接 \(n\) 邊形。設 \(∠AOA'=\theta.\) 則

\[A'(\cos\theta,\sin\theta)\\A''(\cos2\theta,\sin2\theta)\\\ \\\begin{aligned}&\quad(\cos\theta+\sin\theta\times\text{i})\times(\cos2\theta+\sin2\theta\times\text{i})\\ &=\cos\theta\cos2\theta+(\cos\theta\sin2\theta+\sin\theta\cos2\theta)\times\text{i}-\sin\theta\sin2\theta\\ &=(\sin\theta\sin2\theta+\cos\theta\cos2\theta)+(\sin\theta\cos2\theta+\cos\theta\cos2\theta)\times\text{i}\end{aligned}\]

由三角函數恆等變換式

\[\sin(\alpha+\beta)=\sin\alpha\sin\beta-\cos\alpha\cos\beta,\\ \cos(\alpha+\beta)=\sin\alpha\cos\beta+\sin\beta\cos\alpha\]

知

\[\begin{aligned}原式&=\sin(\theta+2\theta)+\cos(\theta+2\theta)\times\text{i}\\ &=\sin3\theta+\cos3\theta\times\text{i}\end{aligned}\]

換句話說，\(A'\) 和 \(A''\) 表示的數相乘后得到了對應的 \(A'''\)。
設這個多邊形的頂點為 \(\{\omega_i^0\},\ i\in[0,n-1]\)。那么我們一定有$$\omega_n^k\times\omega_n^1=\omega_n^{k+1},\ k\in[0,n-1]$$
特別地，我們有 \(\omega_n^n=\omega_n^0\)。它們叫做 \(n\) 次 單位根。

3. \(\omega_n^0=\omega_n^n=1\)；
4. \(\omega_n^0,\omega_n^1,...,\omega_n^{n-1}\) 互不相同；
5. \(\forall k,n\in\N^*\) 有 \(\omega_n^k=\omega_{2n}^{2k}\)；

證明一： 感性理解。一個 \(n\) 邊形把單位圓分成了 \(n\) 個部分，取這 \(n\) 個圓弧的中點，順次連接這 \(2n\) 個點，得到一個新多邊形。則 \(\omega_n^k\) 即表示 \(n\) 變形的第 \(k\) 個單位根，也表示 \(2n\) 變形的第 \(2k\) 個點。

證明二：
\(\begin{aligned}\omega_{2n}^{2k}&=\cos2k\times\frac{2\pi}{2n}+\text{i}\times\sin2k\times\frac{2\pi}{2n}\\ &=\cos k\times\frac{2\pi}n+\text{i}\times\sin k\times\frac{2\pi}n\\ &=\omega_n^k.\end{aligned}\)

5. \(\forall k,n\in\N^*\) 有 \(\omega_n^{k+\frac2n}=-\omega_n^k\)；

證明一： 感性理解。乘以 \(\omega_n^2\) 的意思其實就是在單位圓逆時針轉半圈。單位圓上的一個點，逆時針轉了半圈后到達的點，與原來的點關於原點對稱。

證明二： $$\begin{aligned}\omega_n^{k+\frac2n}&=\omega_n^k\times\omega_n^{\frac2n}\
&=\omega_n^k\times(\cos\pi+\text{i}\times\sin\pi)\
&=\omega_n^k\times(-1+0)\
&=-\omega_n^k.\end{aligned}$$

6. \(\forall k,n\in\N^*\) 有 $$\sum_{i=0}^{n-1}{(\omega_n^k)^i}=\begin{cases}0,k\neq0\n,k=0\end{cases}$$
   證明： \((\text{I})\) 若 \(k\neq0\)，設 $$S=\sum_{i=0}^{n-1}{(\omega_n^k)^i}\quad...(1)$$ 則 $$\omega_n^k\times S=\sum_{i=1}^{n}{(\omega_n^k)^i}\quad...(2)$$
   \((2)-(1)\) 得$$\begin{aligned}(\omega_n^k-1)\times S&=(\omega_n^k)^n-1\
   S&=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}\
   &=\frac{1-1}{\omega_n^k-1}\
   &=0.\end{aligned}$$

\((\text{II})\) 當 \(k=0\) 時，

\[\begin{aligned}原式&=\sum_{i=0}^{n-1}{1}\\&=n.\end{aligned} \]

## 快速傅里葉變換 顯然地，**一個 $n$ 次多項式可以被 $n+1$ 個點唯一確定**。

那么，我們可以在單位圓上取 \(n+1\) 個單位根，代入已知的兩個函數，得到 \(n+1\) 對點，再把每對點相乘，得到結果函數上的 \(n+1\) 個點，再求出結果函數。

設已知函數（合並后）為 $$\begin{aligned}f(x)&=\sum_{i=0}^{n-1}{a_ix^i}\
&=(a_0+a_2x+...+a_{n-2}x^{n-2})\
&+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})\end{aligned}$$
設 $$\begin{aligned}g(x)&=a_0+a_2x+...+a_{n-2}x^{\frac{n}2-1}\
h(x)&=a_1+a_3x+...+a_{n-1}x^{\frac{n}2-1}\end{aligned}$$
則 $$f(x^2)=g(x)+x\times h(x^2)$$
令 \(x=\omega_n^k\)，由 單位根的性質1 得 $$\begin{aligned}f(x^2)&=g(\omega_n^{2k})+\omega_n^k\times h(\omega_n^{2k})\
&=g(\omega_{\frac{n}2}^k)+\omega_n^k\times h(\omega_{\frac{n}2}^k)\quad...(1)\end{aligned}$$
令 \(x=-\omega_n^k=\omega_n^{k+\frac{n}2}\) 得 $$\begin{aligned}f(x^2)&=g(\omega_n^{2k+n})+\omega_n^{k+\frac{n}2}\times h(\omega_n^{2k+n})\
&=g(\omega_n^{2k})+\omega_n^{k+\frac{n}2}\times h(\omega_n^{2k})\
&=g(\omega_{\frac{n}2}^k)-\omega_n^k\times h(\omega_{\frac{n}2}^k)\quad...(2)\end{aligned}$$
不難發現，\((1)\) 式與 \((2)\) 式 互軛！換句話說，當你求出 \((1)\) 式的值時，只需將虛數部分取相反數（時間復雜度為 \(O(1)\)）即得到了 \((2)\) 式。

這樣，我們就將問題的規模減為一半。同理，剩下的一半也可以使用類似方法，分為更小的兩半……沒錯，這就是 分治，這樣就將時間復雜度降為了 \(O(n\log n)\)。

快速傅里葉逆變換

設 \(\{y_i\}\) 是 \(\{a_i\}\) 的傅里葉變換，即在 \(\{\omega_n^i\}\) 處的值；
設 \(\{c_i\}\) 是 \(\{y_i\}\) 在 \(\{\omega_n^{-i}\}\) 的值。
則有 $$\begin{aligned}c_k&=\sum_{i=0}^{n-1}{y_i(\omega_n^{-k})^i}\
&=\sum_{i=0}^{n-1}{(\sum_{j=0}^{n-1}{a_j(\omega_n^k)^j})·(\omega_n^{-k})^i}\
&=\sum_{i=0}^{n-1}{a_j\times\sum_{j=0}^{n-1}{(\omega_n^{j-k})^i}}\end{aligned}$$
根據 單位根的性質6 ，有且僅有一個 \(j\in[0,n-1]\) 使得 \(j=k\)。
\(\therefore \forall j\) 有且僅有一個$${(\omega_n^{j-k})^i}=n$$有且僅有 \(n-1\) 個 $${(\omega_n^{j-k})^i}=0$$

\[\begin{aligned}\therefore\sum_{j=0}^{n-1}{(\omega_n^{j-k})^i}&=n,\\c_k&=\sum_{i=0}^{n-1}{a_j\times\sum_{j=0}^{n-1}{(\omega_n^{j-k})^i}}\\ &=a_k\times n,\\ a_k&=\frac{c_k}n.\end{aligned}\]

換句話說，經過快速傅里葉逆變換后的數組 \(\{c_i\}\)，除以 \(n\) 后就得到了結果 \(\{a_i\}\)。（這里的 \(n\) 是指 \(c\) 數組的長度）

小結

使用快速傅里葉變換（FFT）計算兩個函數之積的步驟如下：

1. 分別對兩個函數進行快速傅里葉變換；
2. 將兩組結果相乘；
3. 對結果進行快速傅里葉逆變換，並將結果除以 \(n\)。
   
   

迭代快速傅里葉變換

經過上文的學習，容易寫出以下代碼（感謝 linjiayang2016 大佬的代碼）

void fast_fast_tle(int limit,complex *a,int type){
    if(limit==1)
        return ;
    complex a1[limit>>1],a2[limit>>1];	//*
    for(int i=0;i<=limit;i+=2)
        a1[i>>1]=a[i],a2[i>>1]=a[i+1];
    fast_fast_tle(limit>>1,a1,type);
    fast_fast_tle(limit>>1,a2,type);
    complex Wn=complex(std::cos(2.0*Pi/limit),type*std::sin(2.0*Pi/limit)),w=complex(1,0);
    for(int i=0;i<(limit>>1);i++,w=w*Wn)
        a[i]=a1[i]+w*a2[i],a[i+(limit>>1)]=a1[i]-w*a2[i];
}

發現 \(*\) 行定義了兩個不小的數組，而 FFT 函數是被遞歸調用的，所以會造成爆棧。

原序列	1	2	3	4	5	6	7	8
快速傅里葉變換后序列	1	5	3	7	2	6	4	8
原下標的二進制	000	100	010	110	001	101	011	111
新下標的二進制	000	001	010	011	100	101	110	111

不難發現，快速傅里葉變換后的序列中，每個數的新下標的二進制，在數值上等於原下標二進制的反轉。

設原下標為 \(x\) 的數的新下標為 \(r[x]\)。

以十進制數 \(184\) 的二進制表示為例。設這個二進制數長度為 \(l\)。

\[\begin{array}{ccc}1&\ &0&\ &1&\ &1&\ &1&\ &0&\ &0&|&0\end{array} \]

我們把這個數分為兩部分，第一部分是前 \(l-1\) 位，第二部分是最后一位。不難得到，它的翻轉就是 在第二部分后面接上第一部分的翻轉形成的數。
即

	r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));

至此，快速傅里葉變換的相關內容已全部結束。

\(\text{luogu P3803}\)

給你兩個正整數 \(n,m\leq10^6\)，和 \(n+1\) 次多項式 \(f(x)\)、\(m+1\) 次多項式 \(g(x)\)。求 \(f(x)\) 和 \(g(x)\) 的卷積。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>

using namespace std;
#define reg register
const int MAXN=400010;
const double pi=acos(-1.0);

struct compex{
	double x,y;
	friend compex operator +(const compex a,const compex b){
		compex c;
		c.x=a.x+b.x;c.y=a.y+b.y;
		return c;
	}
	friend compex operator -(const compex a,const compex b){
		compex c;
		c.x=a.x-b.x;c.y=a.y-b.y;
		return c;
	}
	friend compex operator *(const compex a,const compex b){
		compex c;
		c.x=a.x*b.x-a.y*b.y;
		c.y=a.x*b.y+a.y*b.x;
		return c;
	}
}a[MAXN],b[MAXN];
int r[MAXN];
int n,m,l=0;
int limit=1;

void FFT(compex*t,int type){
	for(reg int i=0;i<limit;++i)
		if(i<r[i]) swap(t[i],t[r[i]]);
	for(reg int mid=1;mid<limit;mid<<=1){
		compex wn=(compex){cos(pi/mid),type*sin(pi/mid)};
		for(reg int j=0,R=(mid<<1);j<limit;j+=R){
			compex w={1.0,0};
			for(reg int k=0;k<mid;++k,w=w*wn){
				compex x=t[j+k],y=w*t[j+k+mid];
				t[j+k]=x+y;
				t[j+k+mid]=x-y;
			}
		}
	}
}
int main(){
	scanf("%d%d",&n,&m);
	for(reg int i=0;i<=n;++i)
		scanf("%lf",&a[i].x);
	for(reg int i=0;i<=m;++i)
		scanf("%lf",&b[i].x);
	while(limit<=n+m) limit<<=1,++l;
	for(reg int i=1;i<limit;++i)
		r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
	FFT(a,1);FFT(b,1);
	for(reg int i=0;i<limit;++i)
		a[i]=a[i]*b[i];
	FFT(a,-1);
	for(reg int i=0;i<=n+m;++i)
		printf("%d ",(int)(a[i].x/limit+0.5));
}