GGAutomaton's Blog
0%

FFT模板

发表于 分类于 ，

从多项式乘法到快速傅里叶变换，Miskcoo
FFT 学习笔记，Menci

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
#pragma GCC optimize(2)
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const int MAXN=2100000;//两倍空间
const double pi=acos(-1);
struct complex_t {
	double a,b;
	complex_t operator *(const complex_t &res) {
		return (complex_t){a*res.a-b*res.b,a*res.b+b*res.a};
	}
	complex_t operator +(const complex_t &res) {
		return (complex_t){a+res.a,b+res.b};
	}
	complex_t operator -(const complex_t &res) {
		return (complex_t){a-res.a,b-res.b};
	}
	void operator *=(const complex_t &res) {
		*this=(*this)*res;
	}
	void operator +=(const complex_t &res) {
		*this=(*this)+res;
	}
};
int n,m,step,r[MAXN];
complex_t a[MAXN],b[MAXN];
inline int read() {
	char ch=getchar();int sum=0;
	while (ch>'9' || ch<'0') ch=getchar();
	while (ch>='0' && ch<='9') sum=sum*10+ch-'0',ch=getchar();
	return sum;
}
void DFT(complex_t c[],int T) {//1:DFT -1:IDFT
	for (int i=0;i<n;i++) {
		if (i<r[i]) swap(c[i],c[r[i]]);
	}
	for (int i=1;i<n;i<<=1) {
		complex_t x=(complex_t){cos(pi/i),T*sin(pi/i)};
		for (int j=0;j<n;j+=(i<<1)) {
			complex_t y=(complex_t){1,0};
			for (int k=0;k<i;k++,y=y*x) {
				complex_t z=y*c[i+j+k];
				c[i+j+k]=c[j+k]-z;
				c[j+k]+=z;//顺序不可颠倒 
			}
		}
	}
}
void FFT() {
	m+=n;n=1;//n+m次多项式，n+m+1项 
	while (n<=m) n<<=1,step++;//必须>n+m，n=(1<<step) 
	for (int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(step-1));//计算元素最终的位置 
	DFT(a,1);DFT(b,1);//DFT
	for (int i=0;i<=n;i++) a[i]*=b[i];
	DFT(a,-1);//IDFT
	for (int i=0;i<=m;i++) a[i].a=a[i].a/n+0.5;//IDFT 1/n
}
int main() {
	n=read(),m=read();//n次多项式，n+1个项 
	for (int i=0;i<=n;i++) a[i].a=read();
	for (int i=0;i<=m;i++) b[i].a=read();
	FFT();
	for (int i=0;i<=m;i++) printf("%d ",(int)a[i].a);
	return 0;
}