FFT 学习笔记

例题：

https://www.luogu.org/problem/show?pid=1919

GDKOI2018 D1T3 基站

FFT：

OI竞赛中，快速傅里叶变换(fft) 可以在 O(nlogn)O(nlogn) 的时间内计算 多项式乘法 卷积

需要用到 复数（1的n次方虚数根）和分治

这个教程的公式推导比较详细

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-3

搬运结论：

A(ωkn)=∑i=0n2−1a2iωkin2+ωkn∑i=0n2−1a2i+1ωkin2A(ωnk)=∑i=0n2−1a2iωn2ki+ωnk∑i=0n2−1a2i+1ωn2ki

A(ωk+n2n)==∑i=0n2−1a2iωkin2+ωk+n2n∑i=0n2−1a2i+1ωkin2∑i=0n2−1a2iωkin2−ωkn∑i=0n2−1a2i+1ωkin2A(ωnk+n2)=∑i=0n2−1a2iωn2ki+ωnk+n2∑i=0n2−1a2i+1ωn2ki=∑i=0n2−1a2iωn2ki−ωnk∑i=0n2−1a2i+1ωn2ki

递归实现FFT（较慢）

void fft(comp *p,int n,int f) {
if (!(n>>=1)) return ;
comp p0
, p1
;
for (int i=0;i<n;i++) p0[i] = p[i*2], p1[i] = p[i*2+1];
fft(p0,n,f);  fft(p1,n,f);
comp w(1,0), x(cos(M_PI/n),f*sin(M_PI/n));
for (int i=0;i<n;i++,w *= x) p[i] = p0[i] + w * p1[i], p[i+n] = p0[i] - w * p1[i];
}

AC Code

#include <cstdio>
#include <cmath>
#include <complex>
#define N 200200
typedef std::complex<double> comp;

void fft(comp *p,int n,int f) {
if (!(n>>=1)) return ;
comp p0
, p1
;
for (int i=0;i<n;i++) p0[i] = p[i*2], p1[i] = p[i*2+1];
fft(p0,n,f);  fft(p1,n,f);
comp w(1,0), x(cos(M_PI/n),f*sin(M_PI/n));
for (int i=0;i<n;i++,w *= x) p[i] = p0[i] + w * p1[i], p[i+n] = p0[i] - w * p1[i];
}
inline int f(double x) {return (int)(x + 0.5);}

int n,L,i,x,ans
;
comp a
,b
;
char A
,B
;

int main() {
scanf("%d%s%s",&n,A,B);  --n;
for (i=0;i<=n;i++) a[i].real() = A[n-i] - '0';
for (i=0;i<=n;i++) b[i].real() = B[n-i] - '0';
for (L=1;L<=n*2;L<<=1);  fft(a,L,1);  fft(b,L,1);
for (i=0;i<=L;i++) a[i] *= b[i];     fft(a,L,-1);

for (i=0,x=0;i<=L;i++) if ((ans[i] = f(a[i].real()/L) + x) > 9) x = ans[i] / 10, ans[i] %= 10;  else x = 0;
for (i=L;~i;i--) if (ans[i]) break;
for (;~i;i--) printf("%d",ans[i]);  putchar('\n');
}
//a*b problem 要考虑前导0和进位，要注意细节，例如倒序存 a[i]

补充：

FFT非递归

inline void fft(comp a[],int f) {
for (i=0;i<N;i++) if (i < rev[i]) swap(a[i],a[rev[i]]);
for (i=1;i<N;i<<=1) {
comp wn(cos(M_PI/i),f*sin(M_PI/i)),w,X,Y;
for (j=0;j<N;j+=(i<<1)) {
w = comp(1,0);
for (k=0;k<i;k++,w = w * wn) {
X = a[j+k]; Y = w*a[j+i+k];
a[j+k] = X+Y;  a[j+i+k] = X-Y;
}
}
}
if (f == -1) for (i=0;i<N;i++) a[i].r

/= N;
}