FFT

1 簡述

FFT是專門用來求解多項式乘法的一個高效演算法。 總所周知，樸素的多項式乘法的時間複雜度是\(O(n^2)\)，而FFT利用複數的知識做到了\(O(nlogn)\)。

2 點值表示式

設\(A(x)\)是一個n-1次方的多項式，那麼把n個不同的x代入，一定可以得到n個y，這n對(x,y)唯一確定了該多項式的係數。由多項式可以求出點值表示式，由點值表示式可以求出這個多項式。

這裡，我們把n作為2的冪次方存在。

我們發現，用點值表示式做多項式乘法是\(O(n)\)的，即\(A(x_i)=B(x_i)*C(x_i)\)。 所以把多項式先轉換為點值表示式，然後就可以\(O(n)\)解決多項式乘法問題。

樸素的，把多項式轉換為點值表示式也是\(O(n^2)\)的。

所以多項式乘法的瓶頸在於如果在更快的時間內把多項式轉換為點值表示式。

這個就需要傅立葉變換了。

3 複數

複數是數學上一個很常見的概念，複數的特點是對負數進行開方。負數相加的規則是實部（x軸）和實部相加，虛部（y軸）和虛部相加。 複數乘法的規則是模長相乘，幅角相加。 下面略作證明：

\[z_1=r_1*(cos\theta+i*sin\theta),z_2=r_2*(cos\xi+i*sin\xi) \]

\[z_1*z_2=r_1*r_2*(cos\theta*cos\xi+i*cos\theta*sin\xi+i*sin\theta*sin\xi-sin\theta*sin\xi) \]

\[z_1*z_2=r_1*r_2*(cos(\theta+\xi)+i*sin(\theta+\xi)) \]

c++裡面提供了複數的模板，可以直接進行加減乘除（相當於一個pair），當然我們也可以自己寫。

#include<complex>
complex<double> x;

我們要用到的複數都是模長為1的複數，這樣相乘的模長還是1，只有幅角進行了改變。 我們可以畫一個單位圓，把這個圓平均進行n等分，每個點都表示一個複數。

從（1，0）開始，我們逆時針旋轉n個點從0開始編號，第k個點的複數記作\(w^{k}\)，明顯的\(w^k=w*w*w...*w\)是k個w相乘，\(w^k\)

其他的一些性質：

1. \(w^{2k}_{2n}=w^k_n\)，這個很明顯，因為幅角是一樣的。
2. \(w^{k+\frac{n}{2}}=-w^k\)，這個也很明顯，因為差180度的複數剛好是相反數。

單位圓上的點有什麼特殊性質麼？

設\(y_0,y_1,...,y_{n-1}\)是多項式\(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\)的離散傅立葉變換，我們再設一個多項式\(B(x)=y_0+y_1x+y_2x^2+...+y_{n-1}x^{n-1}\)，現在我們把上面的n個單位根的倒數，即\(w^0_n,w^{-1}_n,w^{-2}_n,...\)作為x代入到B裡面去，得到一個新的離散傅立葉變換\((z_0,z_1,z_2,....)\)。

\[z_k=\sum_{i=0}^{n-1}y_i(w^{-k}_n)^i \]

\[z_k=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_jx^j_i)(w^{-k}_n)^i \]

\[z_k=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(w^i_n)^j)(w^{-k}_n)^i \]

\[z_k=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(w^{j-k})^i) \]

\[\sum_{i=0}^{n-1}(w^{j-k})^i,這個是可以求出來的。 \]

當\(j=k\)時，答案是n，如果\(j \ne k\)，那麼根據等比數列求和公式：\(\frac{1-(w^{j-k})^n}{1-w^{j-k}}=0\)

所以\(z_k=na_k\)，所以\(a_k=\frac{z_k}{n}\)

總結

• 把多項式A(x)的離散傅立葉變換的結果作為另一個多項式B(x)的係數，取單位根的倒數作為x代入B(x)，得到的點值再除以n，就是A(x)的各項係數。
• 從而實現了傅立葉變換的逆變換，把點值轉換為係數。
• 這個就是傅立葉變換神奇的性質。

4 離散傅立葉變換

樸素的傅立葉變換還是太慢了，所以我們要進行快速的傅立葉變換，藉助於分治的思想。

我們設：

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

把每個\(w^k_n\)代入，然後把A(x)按照下標分成奇偶兩半部分。

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

設有多項式：

\[A_1(x)=a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1} \]

\[A_2(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1} \]

於是：\(A(x)=A_1(x^2)+xA_2(x^2)\)

如果\(k<\frac{n}{2}\)，把\(w^k_n\)代入：

\[A(w^k_n)=A_1(w^{2k}_n)+w^k_nA_2(w^{2k}_n) \]

\[=A_1(w^k_{\frac{n}{2}})+w^k_nA_2(w^{k}_{\frac{n}{2}}) \]

那麼對於\(A(w^{k+\frac{n}{2}})\):

\[A(w^{k+\frac{n}{2}})=A_1(w^{2k+n}_n)+w^{k+\frac{n}{2}}A_2(w^{2k+n}_n) \]

\[=A_1(w^k_{\frac{n}{2}}\times w^n_n)-w^k_nA_2(w^k_{\frac{n}{2}}\times w^n_n) \]

\[=A_1(w^k_{\frac{n}{2}})-w^k_nA_2(w^k_{\frac{n}{2}}) \]

於是問題就變成我們只要先計算出\(\frac{n}{2}\)時候的答案就可以了。

#include<bits/stdc++.h>
using namespace std;  
int const N=2e5+10;  
double const PI=asin(1.0)*2; 
typedef complex<double> cp; 
cp tmp[N<<1],a[N<<1],b[N<<1];  
int n,m,ans[N<<1];  
cp c(int n,int k){
	return cp(cos(2*k*PI/n),sin(2*k*PI/n));  
}  
void fft(cp *a,int n,int inv){
	if(n==1) return;  
	int m=n/2;  
	for(int i=0;i<m;i++){
		tmp[i]=a[2*i];  
		tmp[i+m]=a[2*i+1];  
	}  
	for(int i=0;i<n;i++) a[i]=tmp[i];  
	fft(a,m,inv);  
	fft(a+m,m,inv);  
	for(int i=0;i<m;i++){
		cp x=c(n,i);  
		if(inv) x=conj(x);  
		tmp[i]=a[i]+x*a[i+m];  
		tmp[i+m]=a[i]-x*a[i+m];  
	}
	for(int i=0;i<n;i++)  
		a[i]=tmp[i];  
}
int main(){
	scanf("%d%d",&n,&m); 
	for(int i=0;i<=n;i++) 
		scanf("%lf",&a[i].real());  
	for(int i=0;i<=m;i++) 
		scanf("%lf",&b[i].real());  
	int k=1;  
	while (k<n+m+1) k*=2;  
	fft(a,k,0); 
	fft(b,k,0);  
	for(int i=0;i<k;i++)  
		a[i]*=b[i];  
	fft(a,k,1);  
	for(int i=0;i<=n+m;i++)  
		printf("%d ",int(a[i].real()/k+0.5)); 
	return 0; 
}