FFT

紅樓首席船師發表於2024-10-04
FFT

前置芝士：

尤拉說：$e^{i\theta}=cos(\theta)+isin(\theta)$

定義單位根$\omega_{n}{k}=e{n}}$

記$\theta=2\pi/n$ ,則$\omega_{n}{k}=e=cos(k\theta)+isin(k\theta)$

$1.\omega_{n}{0}=e=1,\omega_{n}{n}=e=cos(2\pi)+isin(2\pi)=1$

$2.(\omega_{n}{k})2=(e^{i 2\pi \frac{k}{n}})2=e{n/2}}=\omega_{n/2}^{k}$

$3.\omega_{n}{k}=cos(k\theta)+isin(k\theta)=-cos(k\theta+\pi)-isin(k\theta+\pi)=-e=-e^{i(\frac{2\pi k}{n}+\pi)}= -e{i2\pi(\frac{k+n/2}{n})}=-\omega_{n}$

$4.\omega_{n}{a}*\omega_{n}=e^{i2\pi \frac{a+b}{n}}=\omega_{n}^{a+b}$

基本思路

對於一個最高次為$n$次的多項式，我們可以透過$n+1$組$(x_i,F(x_i))$來唯一確定這個多項式的係數，這種表示方式叫做點值表示法。

$FFT$：把一個多項式從係數表示轉化為點值表示法

$IFFT$：把一個多項式從點值表示法轉為係數表示法

為求解$G(x)=H(x)F(x)$ 的係數，我們的思路是先把$H(x)$和$F(x)$透過$FFT$轉為點值表示，分別為$(x_i,H(x_i))和(x_i,F(x_i))$，則$(x_i,H(x_i)F(x_i))$為$G(x)$的點值表示，再透過$IFFT$求出$G(x)$ 的係數

FFT

設$F(x)=a_0+a_1x+a_2x2...+a_{n-1}x$

我們可以利用單位根的一些性質進行快速計算，不妨將點值表示定為$(\omega_{n}{0}，F(\omega_{n})),(\omega_{n}{1}，F(\omega_{n})),...(\omega_{n}{n-1}，F(\omega_{n}))$

問題轉化為快速求解這$n$個點值

將$F(x)$的奇偶項分開（我們暫把$n$當作偶數）

$F(x)=(a_0+a_2x2+...+a_{n-2}x)+(a_1x+a_3x3+...+a_{n-1}x)$

$F(x)=(a_0+a_2x2+...+a_{n-2}x)+x(a_1+a_3x2+...+a_{n-1}x)$

記為$F(x)=A_0(x2)+xA_1(x2)$

代入 $x=\omega_{n}^{k}$

① $k<=n/2$

​ $F(\omega_{n}{k})=A_0(\omega_{n/2})+\omega_{n}{k}A_1(\omega_{n/2})$

②$k>n/2$

​ 記$k'=k-n/2$

​ $F(\omega_{n}{k})=A_0((-\omega_{n})2)+\omega_{n}A_1((-\omega_{n}{k'})2)$

​ $F(\omega_{n}{k})=A_0(\omega_{n/2})-\omega_{n}{k'}A_1(\omega_{n/2})$

那麼對於$k<=n/2$

​ $F(\omega_{n}{k})=A_0(\omega_{n/2})+\omega_{n}{k}A_1(\omega_{n/2})$

​ $F(\omega_{n}{k+n/2})=A_0(\omega_{n/2})-\omega_{n}{k}A_1(\omega_{n/2})$

發現求出$k<=n/2$部分的所有$A_0(x)$和$A_1(x)$，則可同時得到$k>n/2$部分的$F(x)$

那麼可以分治：

​ 我們最開始要求解$n$個長度為$n$的多項式$F(x)$

​ 而我們只要求$\frac{n}{2}$ 個長度為$\frac{n}{2}$ 的多項式$A_0(x)$和另$\frac{n}{2}$個長度為$\frac{n}{2}$的$A_1(x)$

​ ......以此類推

​ 最終只要求解共$n$ 個長度為 1 的多項式即可，長度為1意味著最高次是0，直接返回係數即可

則每一層的計算複雜度為O(n)，共log(n)層，複雜度為O(nlogn)

為了避免一些不必要的討論，且為了後續的一個最佳化，我們始終希望計算的每一個多項式的長度$n$是一個偶數，因此我們常把多項式的長度填充到$2^k$次為止，這也是前面直接把$n$當作偶數的原因。

IFFT

假設$G(x)=b_0+b_1x+b_2x2...+b_{n-1}x$

而我們現在已經求出了$G(x)$的點值表示$(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})$其中$x_i=\omega_{n}^{i}$

類似於中國剩餘定理，我們可以構造出一個多項式

令$\phi(x)=y_0+y_1x+y_2x2...+y_{n-1}x$

則$\phi(\omega_{n}{-i})=\displaystyle\sum_{j=0} y_j *(\omega_{n}{-i})j$

開啟$y_i$有：$\phi(\omega_{n}{-i})=\displaystyle\sum_{j=0}\displaystyle\sum_{k=0}{n-1}b_k*(\omega_{n})k*(\omega_{n})j=\displaystyle\sum_{j=0}\displaystyle\sum_{k=0}{n-1}b_k*\omega_{n}=\displaystyle\sum_{k=0}{n-1}\displaystyle\sum_{j=0}b_k\omega_{n}{j(k-i)}=\displaystyle\sum_{k=0}(b_k\displaystyle\sum_{j=0}{n-1}\omega_{n})$

這個形式跟中國剩餘定理的答案構造異曲同工，那我們進行分類

① $k=i$

​ $b_i\displaystyle\sum_{j=0}{n-1}\omega_{n}=b_in$

②$k\neq i$

​ $\displaystyle\sum_{j=0}{n-1}\omega_{n}=\displaystyle\sum_{j=0}{n-1}(\omega_{n})j=(\omega_{n})0+(\omega_{n})1...+(\omega_{n}){n-1}=\frac{(\omega_{n})n-1}{\omega_{n}-1}=\frac{(\omega_{n}{n})-1} {\omega_{n}^{k-i}-1} =\frac{1-1}{\omega_{n}^{k-i}-1}=0$

$\phi(\omega_{n}^{-i})=n*b_i$

$b_i=\frac{\phi(\omega_{n}^{-i})}{n}$

所以我們只需要對$\phi(x)$代入$\omega_{n}{-0},\omega_{n},...,\omega_{n}^{-(n-1)}$做一遍FFT，把他的每個點值除以$n$就得到了$G(x)$的係數

最佳化

直接進行分治需要預開大量空間，而且常數大，我們嘗試不要遞迴分治，而從最底下往上進行多項式的計算

所以我們需要把多項式長度擴充到$2^k$整，並且因此最底層的序數十分有規律,我們模擬一次最底層係數的分配過程

$(a_0,a_1,a_2,a_3, a_4,a_5,a_6,a_7)$

$(a_0,a_2,a_4,a_6),(a_1,a_3,a_5,a_7)$

$(a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)$

$(a_0),(a_4),(a_2),(a_6),(a_1),(a_5),(a_3),(a_7)$

如果把最高層當作第一層，則第$i$層相當於按照第$i$位的二進位制對序數進行分類，為$0$放左，為$1$放右

所以整個過程相當於從低位到高位進行了一次高位到低位的排序

其結果就是最底層的二進位制是反過來的

正常序$(000),(001),(010),(011),(100),(101),(110),(111)$

最底層$(000),(100),(010),(110),(001),(101),(011),(111)$

記$rev(i)$為$i$的逆序

顯然有遞推關係 $i=rev(rev(i/2))*2+i&1$

則有$rev(i)=rev(rev(rev(i/2))*2+i&1)$

注意到$rev(rev(i/2))*2$是一個偶數，則最低位是空出來的，與$i&1$不產生進位，所以我們可以直接拆開括號，有

$rev(i)=rev(rev(rev(i/2))*2)+rev(i&1)$

由於rev本身代表了逆序操作，開啟$rev$時要把位數倒過來，如下

$rev(i)=rev(rev(rev(i/2)))/2+(i&1)*2^{k-1}$

於是有$rev(i)=rev(i/2)/2+(i&1)*2^{k-1}$

#include<bits/stdc++.h>
using namespace std;

const int N = 10000000 + 5;
const double Pi = acos(-1.0); 

inline int read(){
	char ch=getchar();int x=0;
	while(ch<'0'||ch>'9') ch=getchar();
	for(;ch>='0'&&ch<='9';ch=getchar())
		x=(x<<3)+(x<<1)+(ch^48);
	return x;
}

complex <double> a[N],b[N];
int n,m,lim,l,rev[N];

void FFT(complex<double> *A,int opt){
	for(int i=0;i<lim;++i) 
		if(i<rev[i]) swap(A[i],A[rev[i]]);
	for(int m=2;m<=lim;m<<=1){int mid=m>>1;
		complex<double> Wn(cos(2.0*Pi/m),opt*sin(2.0*Pi/m));
		for(int i=0;i<lim;i+=m){
			complex <double> w(1,0);
			for(int k=0;k<mid;++k,w=w*Wn){
				complex <double> x=A[i+k],y=w*A[i+k+mid];
				A[i+k]=x+y;
				A[i+k+mid]=x-y;
			}
		}
	}
}

int main(){
	n=read(),m=read();
	for(int i=0;i<=n;++i) a[i].real(read());
	for(int i=0;i<=m;++i) b[i].real(read());
	lim=1,l=0;
	while(lim<n+m+1) lim<<=1,++l;
	for(int i=0;i<lim;++i)
		rev[i]=(rev[i>>1]>>1)+((i&1)<<(l-1));
	FFT(a,1);
	FFT(b,1);
	for(int i=0;i<lim;++i) a[i]=a[i]*b[i];
	FFT(a,-1);
	for(int i=0;i<=n+m;++i)
		printf("%d ",int(a[i].real()/lim+0.5));
	return 0;
}

相關文章