FFT，快速傅立葉變換，蒟蒻看別人的題解都太深奧，看不懂，好不容易學會，以蒟蒻的理解寫給那些想學FFT卻又找不到合適的資料的OIer，蒟蒻理解有限，難免有許多錯誤，請大家多多包涵。

快速傅立葉變換

百度的各種講解都TM扯什麼頻率什麼的，蒟蒻完全看不懂，後來認真看了看算導，獲益匪淺，算導上講的真心不賴，有很多內容都來自算導。

1.多項式

        多項式的兩種表達方式：係數表達和點值表達

係數表達就是大家常用的表達方式，點值表達就像在這個多項式函式上取n個不同的點，這樣就可以確定原多項式。

比如說二次函式需要3個點就可以確定,一次函式需要2個點，一個n次多項式需要n個點(n次多項式意思是有0..n-1次冪的多項式)

A(x)=x^2+2*x-1可以被表達為{  ( 0 , -1 ) , ( 1 , 2 ) , ( 2 , 7 )  }

加法和乘法：

B(x)=x^2-x+2  { ( 0 , 2 ) , ( 1 , 2 ) , ( 2 , 4 ) }

C(x)=A(x)+B(x)=2x^2+x+1   { ( 0, 1) , ( 1 , 4 ) , ( 2, 11 ) }

     注意乘法需要2n個點 lz比較懶就不寫了……

於是我們得到一個計算多項式的方法：

 2.n次單位複數根

有關複數根的性質可以百度到，不再贅述

3.DFT&&FFT

使用單位根計算點值表示式叫DFT（離散傅立葉變換）複雜度n^2，FFT是其優化版複雜度nlogn

計算FFT的虛擬碼(好吧用的是python的高亮)

下劃線代表的是下標,括號代表上標，for 迴圈的range是左閉右開的

FFT(a):
	n=a.length()
	if n==1:
		return a
	w_n=e^(pi*i/n)=complex(cos(2*pi/n),sin(2*pi/n))
	w=1
	a(0)=[a0,a2,...a_n-2]
	a(1)=[a1,a3,...a_n-1]
	y(0)=FFT(a(0))
	y(1)=FFT(a(1))
	for k in range(0,n/2):
		y_k=y_k(0)+w*y_k(1)
		y_k+n/2=y_k(0)-w*y_k(1)
		w=w*w_n
	return y
 

FFT的for迴圈中有兩次w_n^k*y_k(1)的計算，於是可以改寫成這樣

for k in range(0,n/2):
	t=w*y_k(1)
	y_k=y_k(0)+t
	y_k+n/2=y_k(0)-t
	w=w*w_n
#這一過程被稱蝴蝶操作

觀察每次按照奇偶位置分割所形成的樹：

每個數和他二進位制相反的位置互換！！

虛擬碼（算導給的真是……）

BIT-REVERSE-COPY(a,A):
	n=a.length()
	for k in range(0,n):
		A[rev(k)]=a_k
#算導說rev函式很好寫，就沒寫……

FFT(a):
	BIT-REVERSE-COPY(a,A)
	n=a.length()
	for s in range(1,log2(n)+1):
		m=2^s
		w_m=e^(2*pi*i/m)=complex(cos(2*pi*m),sin(2*pi*m))
		for k in range(0,n,m):
			w=1
			for j in range(0,m/2):
				t=w*A[k+j+m/2]
				u=A[k+j]
				A[k+j]=u+t
				A[k+j+m/2]=u-t
				w=w*w_m
	return A

#include<bitset>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#define N 400005
#define pi acos(-1.0) // PI值
using namespace std;
struct complex
{
	double r,i;
	complex(double real=0.0,double image=0.0){
		r=real;	i=image;
	}
	// 以下為三種虛數運算的定義
	complex operator + (const complex o){
		return complex(r+o.r,i+o.i);
	}
	complex operator - (const complex o){
		return complex(r-o.r,i-o.i);
	}
	complex operator * (const complex o){
		return complex(r*o.r-i*o.i,r*o.i+i*o.r);
	}
}x1[N],x2[N];
char a[N/2],b[N/2];
int sum[N]; // 結果存在sum裡
int vis[N];
void brc(complex *a,int l){//原來神犇的二進位制平攤反轉置換太神看不懂，蒟蒻寫了一個O(n)的…… 
	memset(vis,0,sizeof(vis));//O(logn)的在後面 
	for(int i=1;i<l-1;i++){
		int x=i,y=0;
		int m=(int)log2(l)+0.1;
		if(vis[x])continue;
		while(m--){
			y<<=1;
			y|=(x&1);
			x>>=1;
		}
		vis[i]=vis[y]=1;
		swap(a[i],a[y]);
	}	
}
void fft(complex *y,int l,double on) // FFT O(nlogn)
				     		// 其中on==1時為DFT，on==-1為IDFT
{
	register int h,i,j,k;
	complex u,t; 
	brc(y,l); // 呼叫反轉置換
	for(h=2;h<=l;h<<=1) // 控制層數
	{
		// 初始化單位復根
		complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
		for(j=0;j<l;j+=h) // 控制起始下標
		{
			complex w(1,0); // 初始化螺旋因子
			for(k=j;k<j+h/2;k++) // 配對
			{
				u=y[k];
				t=w*y[k+h/2];
				y[k]=u+t;
				y[k+h/2]=u-t;
				w=w*wn; // 更新螺旋因子
			} // 據說上面的操作叫蝴蝶操作…
		}
	}
	if(on==-1)	for(i=0;i<l;i++)	y[i].r/=l; // IDFT
}
/* 
void fft2(complex *a,int s,int t){//蒟蒻自己寫的遞迴版FFT，不保證正確 ,程式碼內部有未定義變數 
	if((n>>t)==1)return;//s記錄起始,t記錄深度,呼叫時應從0開始 
	fft(a,s,t+1);
	fft(a,s+(1<<t),t+1);
	for(int i=0;i<(n>>(t+1));i++){
		p=(i<<(t+1))+s;
		wt=w[i<<t]*a[p+(1<<t)];
		tt[i]=a[p]+wt;
		tt[i+(n>>(t+1))]=a[p]-wt;
	}
	for(i=0;i<(n>>t);i++)a[(i<<t)+s]=tt[i];
}*/
int main(void)
{
	int l1,l2,l;
	register int i;
	while(scanf("%s%s",a,b)!=EOF)
	{
		l1=strlen(a);
		l2=strlen(b);
		l=1;
		while(l<l1*2 || l<l2*2)	l<<=1; // 將次數界變成2^n
					      				// 配合二分與反轉置換
		for(i=0;i<l1;i++) // 倒置存入
		{
			x1[i].r=a[l1-i-1]-'0';
			x1[i].i=0.0;
		}
		for(;i<l;i++)	x1[i].r=x1[i].i=0.0;
		// 將多餘次數界初始化為0
		for(i=0;i<l2;i++)
		{
			x2[i].r=b[l2-i-1]-'0';
			x2[i].i=0.0;
		}
		for(;i<l;i++)	x2[i].r=x2[i].i=0.0;
		fft(x1,l,1); // DFT(a)
		fft(x2,l,1); // DFT(b)
		for(i=0;i<l;i++)	x1[i]=x1[i]*x2[i]; // 點乘結果存入a
		fft(x1,l,-1); // IDFT(a*b)
		for(i=0;i<l;i++)	sum[i]=x1[i].r+0.5; // 四捨五入
		for(i=0;i<l;i++) // 進位
		{
			sum[i+1]+=sum[i]/10;
			sum[i]%=10;
		}
		l=l1+l2-1;
		while(sum[l]<=0 && l>0)	l--; // 檢索最高位
		for(i=l;i>=0;i--)	putchar(sum[i]+'0'); // 倒序輸出
		putchar('\n');
	}
	return 0;
}
/*void brc(complex *y,int l) // 二進位制平攤反轉置換 O(logn)
{
	register int i,j,k;
	for(i=1,j=l/2;i<l-1;i++)
	{
		if(i<j)	swap(y[i],y[j]); // 交換互為下標反轉的元素
								// i<j保證只交換一次
		k=l/2;
		while(j>=k) // 由最高位檢索，遇1變0，遇0變1，跳出
		{
			j-=k;
			k>>=1;
		}
		if(j<k)	j+=k;
	}
}*/

pyc神犇的寫法，bzoj3527，zjoi2014 力，無限YM

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
const int maxn=1000010;
int n,N,L;
int rev[maxn];
int dig[maxn];
double p[maxn];
struct cp{
	double r,i;
	cp(double _r=0,double _i=0):
		r(_r),i(_i){}
	cp operator+(cp x){return cp(r+x.r,i+x.i);}
	cp operator-(cp x){return cp(r-x.r,i-x.i);}
	cp operator*(cp x){return cp(r*x.r-i*x.i,r*x.i+i*x.r);}
};
cp a[maxn],b[maxn],c[maxn],A[maxn],x,y;
void FFT(cp a[],int flag){
	for(int i=0;i<N;i++)A[i]=a[rev[i]];
	for(int i=0;i<N;i++)a[i]=A[i];
	for(int i=2;i<=N;i<<=1){
		cp wn(cos(2*M_PI/i),flag*sin(2*M_PI/i));
		for(int k=0;k<N;k+=i){
			cp w(1,0);
			for(int j=0;j<i/2;j++){
				x=a[k+j];
				y=w*a[k+j+i/2];
				a[k+j]=x+y;
				a[k+j+i/2]=x-y;
				w=w*wn;  
			}
		}
	}
	if(flag==-1)for(int i=0;i<N;i++)a[i].r/=N;
}
double anss[maxn];
int main(){
	scanf("%d",&n);
	for(int i=0;i<n;i++)scanf("%lf",&p[i]);
	for(L=0,N=1;N<n;N<<=1,L++);L++;N<<=1;
	for(int i=0;i<N;i++){
		int len=0;
		for(int t=i;t;t>>=1)dig[len++]=t&1;
		for(int j=0;j<L;j++)rev[i]=rev[i]*2+dig[j];
	}
	for(int i=0;i<n;i++)a[i]=cp(p[i],0);
	for(int i=1;i<n;i++)b[i]=cp(1.0/i/i,0);
	FFT(a,1);FFT(b,1);
	for(int i=0;i<N;i++)c[i]=a[i]*b[i];
	FFT(c,-1);
	for(int i=0;i<n;i++)anss[i]=c[i].r;
	memset(a,0,sizeof(a));
	memset(b,0,sizeof(b));
	for(int i=0;i<n;i++)a[i]=cp(p[n-i-1],0);
	for(int i=1;i<n;i++)b[i]=cp(1.0/i/i,0);
	FFT(a,1);FFT(b,1);
	for(int i=0;i<N;i++)c[i]=a[i]*b[i];
	FFT(c,-1);
	for(int i=0;i<n;i++)anss[i]-=c[n-i-1].r;
	for(int i=0;i<n;i++)
		printf("%.9f\n",anss[i]);
	return 0;
}

重新過了一遍高精乘

#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
const int maxn=1e6+10;
struct cp{
	double r,i;
	cp(double _r=0,double _i=0):
		r(_r),i(_i){}
	cp operator+(cp x){return cp(r+x.r,i+x.i);}
	cp operator-(cp x){return cp(r-x.r,i-x.i);}
	cp operator*(cp x){return cp(r*x.r-i*x.i,r*x.i+i*x.r);}
};
cp a[maxn],b[maxn],A[maxn],x,y,c[maxn];
char s1[maxn],s2[maxn];
int sum[maxn],a1[maxn],a2[maxn],dig[maxn];
int len1,len2,rev[maxn],N,L;
void FFT(cp a[],int flag){
	for(int i=0;i<N;i++)A[i]=a[rev[i]];
	for(int i=0;i<N;i++)a[i]=A[i];
	for(int i=2;i<=N;i<<=1){
		cp wn(cos(2*M_PI/i),flag*sin(2*M_PI/i));
		for(int k=0;k<N;k+=i){
			cp w(1,0);
			for(int j=k;j<k+i/2;j++){
				x=a[j];
				y=a[j+i/2]*w;
				a[j]=x+y;
				a[j+i/2]=x-y;
				w=w*wn;
			}
		}
	}
	if(flag==-1)for(int i=0;i<N;i++)a[i].r/=N;
}
int main(){
	scanf("%s%s",s1,s2);
	len1=strlen(s1);
	len2=strlen(s2);
	for(N=1,L=0;N<max(len1,len2);N<<=1,L++);N<<=1;L++;
	for(int i=0;i<N;i++){
		int len=0;
		for(int t=i;t;t>>=1)dig[len++]=t&1;
		for(int j=0;j<L;j++)rev[i]=(rev[i]<<1)|dig[j];
	}
	for(int i=0;i<len1;i++)a1[len1-i-1]=s1[i]-'0';
	for(int i=0;i<len2;i++)a2[len2-i-1]=s2[i]-'0';
	for(int i=0;i<N;i++)a[i]=cp(a1[i]);
	for(int i=0;i<N;i++)b[i]=cp(a2[i]);
	FFT(a,1);FFT(b,1);
	for(int i=0;i<N;i++)c[i]=a[i]*b[i];
	FFT(c,-1);
	for(int i=0;i<N;i++)sum[i]=c[i].r+0.5;
	for(int i=0;i<N;i++){
		sum[i+1]+=sum[i]/10;
		sum[i]%=10;
	}
	int l=len1+len2-1;
	while(sum[l]==0&&l>0)l--;
	for(int i=l;i>=0;i--)
	putchar(sum[i]+'0');
	putchar('\n');
	return 0;
}