【算法笔记】多项式快速算法——卡常FFT

  • 本文纯属理性的愉悦，除非出题人松松松
  • 本文为博主原创，未经许可不得转载

共轭对称

在离散傅里叶变换（DFT）中，我们有如下公式：

$X(n)=\sum_{k=0}^{N-1}x_k w_{N}^{nk}$

其中$n\in [0,n-1]$，$X(n)$表示第$n$个单位根的点值表达

我们用$\overline{X(n)}$表示$X(n)$的共轭复数，则：

$\overline{X(n)}=\sum_{k=0}^{N-1}\overline{x_k}w_{N}^{-nk}=\sum_{k=0}^{N-1}x_k w_{N}^{(N-n)k}=X(N-n)$

所以我们说离散傅里叶变换$DFT[x(n)]$是共轭对称的

　　

构造离散复数

回想我们原先方法：先把$x(n),y(n)$做$DFT$，然后点值相乘后插值回来，一共做了三次$FFT$

接下来我们考虑构造复数$z(n)=x(n)+iy(n)$，如果能把$z(n)$做$DFT$，然后通过计算得到乘积，然后插值回去

那么我们就能只做$2$次$FFT$，成功卡掉$\frac{1}{3}$的常数

接下来就是推式子的愉(du)悦(liu)时间了

$Z(n)=DFT[x(n)+iy(n)]=X(n)+iY(n)$

注意这里的$X(n),Y(n)$都不是实数，所以它们并不是$Z(n)$的实部和虚部

设$Z(n)$的实部和虚部分别为$A(n),B(n)$，那么有

$Z(n)=\overline{Z(N-n)}=\overline{X(N-n)+iY(N-n)}=\overline{A(N-n)+iB(N-n)}$

$X(n)-iY(n)=A(N-n)-iB(N-n)$

$X(n)+iY(n)=A(n)+iB(n)$

可以得到：

$X(n)=\frac{A(n)+A(N-n)+i[B(n)-B(N-n)]}{2}$

$Y(n)=\frac{B(n)+B(N-n)-i[A(n)-A(N-n)]}{2}$

这样我们就完成了离散值$X(n),Y(n)$的计算，后面的步骤和普通FFT是一样的

　　

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 300010
const double pi=acos(-1);
int n,m,R[MAXN];
namespace INIT{
	char buf[1<<15],*fs,*ft;
	inline char getc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
	inline int read(){
		int x=0,f=1;  char ch=getc();
		while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getc();}
		while(isdigit(ch))  {x=x*10+ch-'0';  ch=getc();}
		return x*f;
	}
}using namespace INIT;
const int OUT=500000;
char Out[OUT],*ou=Out;int Outn[30],Outcnt;
inline void write(int x){
	if(!x)*ou++=48;
	else{
		for(Outcnt=0;x;x/=10)
			Outn[++Outcnt]=x%10+48;
		while(Outcnt)*ou++=Outn[Outcnt--];
	}
}
struct complex{
	double r,v;
	complex(double a=0,double b=0):r(a),v(b){}
	inline complex operator +(const complex &b){return complex(r+b.r,v+b.v);}
	inline complex operator -(const complex &b){return complex(r-b.r,v-b.v);}
	inline complex operator *(const complex &b){return complex(r*b.r-v*b.v,r*b.v+v*b.r);}
}A[MAXN],B[MAXN],w[MAXN];
inline void swap(complex &a,complex &b){complex t(a);a=b;b=t;}
void FFT(complex *a,int L,int f){
	for(int i=0;i<L;++i)  if(i<R[i])  swap(a[i],a[R[i]]);
	for(int len=2;len<=L;len<<=1){
		int l=(len>>1);  //complex wn(cos(pi/l),f*sin(pi/l));
		for(int i=1;i<l;++i)  w[i]=w[i-1]*wn;
		for(int st=0;st<L;st+=len)for(int k=0;k<l;++k){
			complex x=a[st+k],y=w[k]*a[st+k+l];
			a[st+k]=x+y;  a[st+k+l]=x-y;
		}
	}
}
void solve(){
	int L=1,H=0;  while(L<n+m)  L<<=1,++H; w[0]=1;
	for(int i=0;i<L;++i)  R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
	FFT(A,L,1);
	for(int i=1;i<L;++i){
		double x1=A[i].r,x2=A[L-i].r;
		double y1=A[i].v,y2=A[L-i].v;
		//complex a((x1+x2)*0.5,(y1-y2)*0.5);
		//complex b((y1+y2)*0.5,(x2-x1)*0.5);
		//为了代码高亮不得不注释掉
		B[i]=a*b;
	}B[0]=A[0].v*A[0].r;
	FFT(B,L,-1);
	for(int i=0;i<n+m-1;++i)  write(int(B[i].r/L+0.5)),*ou++=32;
	fwrite(Out,1,ou-Out,stdout);
}
void solve2(){
	for(int i=0;i<n+m-1;++i){
		int ret=0;
		for(int j=0;j<=i;++j) ret+=A[j].r*A[i-j].v;
		printf("%d ",ret);
	}
}
int main(){
	//freopen(FILE".in","r",stdin);
	//freopen(FILE".out","w",stdout);
	n=read()+1;  m=read()+1;
	for(int i=0;i<n;++i)  A[i].r=read();
	for(int i=0;i<m;++i)  A[i].v=read();
	if(n<=8&&m<=8)  {solve2(); return 0;}
	solve();
	return 0;
}

优化的效果也是显而易见的：

　　

参考文献

  • 《实序列离散傅里叶变换的快速算法》——Miskcoo传送门

本文标题:【算法笔记】多项式快速算法——卡常FFT

文章作者:chty

发布时间:2018年05月13日 - 17时05分

最后更新:2018年07月25日 - 21时33分

许可协议:本文为博主原创，未经许可不得转载。

文章目录
  1. 1. 共轭对称
  2. 2. 构造离散复数
  3. 3. 参考代码
  4. 4. 参考文献
,