FFT代碼詳解

關於FFT原理部分的介紹，在網上已經有很多了，所以在此只講代碼實現部分的內容。

原理可以參考https://www.cnblogs.com/RabbitHu/p/FFT.html

推薦看完它的原理解釋再來看這里的代碼解釋

廢話不多說，上代碼（多項式乘法）

#include <iostream>
#include <cstdio>
#include <cmath>
#define N 4000001
using namespace std;
struct cp//手寫復數類可以卡常
{
	double real,imag;
};
cp operator +(cp a,cp b)
{
	return (cp){a.real+b.real,a.imag+b.imag};
}
cp operator -(cp a,cp b)
{
	return (cp){a.real-b.real,a.imag-b.imag};
}

復數乘法：設$R_{a}$表示$a$的實部系數，$I_{a}$表示$a$的虛部系數

則$a*b$

$=(R_{a}+I_{a})*(R_{b}+I_{b})$

$=R_{a}*R_{b}+R_{a}*I_{b}+R_{b}*I_{a}+I_{a}*I_{b}$

因為$i^2=-1$

所以結果的實部為$R_{a}*R_{b}-I_{a}*I_{b}$

虛部為$R_{a}*I_{b}+R_{b}*I_{a}$

cp operator *(cp a,cp b)
{
　　return (cp){a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real};
}
double pi=acos(-1.0);
int lim,rev[N],len;
cp w[N],inv[N],a[N],b[N];
void get_w()
{
	for(int i=0;i<=lim;i++)
	{
		double angle=(double)i*2*pi/lim;
		w[i].imag=sin(angle);
		w[i].real=cos(angle);
		inv[i]=(cp){w[i].real,-w[i].imag};
	}
}

fft參數解釋

$arr:$系數數組，在$fft$后變為點值數組，$arr_{i}$表示將$w^i_n$帶入多項式后求得的值

$w:$預處理好的w單位根，在$fft$的時候正常帶入即可，在$idft$的時候帶入單位根的倒數（具體參見$idft$）void fft(cp *arr,cp *w)

{
	for(int i=0;i<lim;i++)
	{
		//處理每一個系數在分治過程中的實際位置；
		//if是因為只需交換一次，所以選擇由小的一方來執行
		if(i<rev[i]) swap(arr[i],arr[rev[i]]);
	}
	for(int i=2;i<=lim;i*=2)//枚舉區間長度
	{
		int l=i/2;
		for(int j=0;j<lim;j+=i)//枚舉區間位置，這些區間是互不相交的
		{
			//枚舉帶入的單位根w(k,l)，k>=l的單位根也可以在這里一並求出
			for(int k=0;k<l;k++)
			{

　　

意義變更

在這里$arr$的意義從系數變為$w^k_i$的點值，$a_{j,j+i-1}$分別表示將$w^{0,i-1}_i$的點值

下面的的t相當於文首博客中的$w^k_n * A_2(w^k_{n \over 2})$

				cp t=arr[j+k+l]*w[lim/i*k];//w(k,i)=w(k/i,1)=w(n*k/i,n)
				arr[j+k+l]=arr[j+k]-t;
				arr[j+k]=arr[j+k]+t;
			}
		}
	}
}
int main()
{
	int n,m;
	cin>>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);
	lim=1;
	while(lim<=n+m)	len++,lim<<=1;//這樣會比用cmath的log要快？
	for(int i=0;i<lim;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	get_w();
	fft(a,w);
	fft(b,w);
	for(int i=0;i<=lim;i++) a[i]=a[i]*b[i];
	fft(a,inv);
	for(int i=0;i<=n+m;i++) printf("%d ", (int)(a[i].real/lim+0.5));
　　//除以lim的原因具體參見idft，0.5是為了四舍五入
}