关于FFT
这个博客的讲解超级棒
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform
算法导论上的讲解也不错
模板就是抄一抄别人的啦
首先是递归版本
1 #include <cstdio> 2 #include <complex> 3 #include <cmath> 4 using namespace std; 5 6 const double pi = acos(-1); 7 const int N = (2 << 17) + 10; 8 typedef complex<double> cp; 9 cp A[N], B[N]; 10 int n, m; 11 12 void FFT(cp *y, int n, int o) { 13 if (n == 1) return ; 14 cp l[n >> 1], r[n >> 1]; 15 for (int i = 0; i <= n; i++) 16 if (i & 1) r[i >> 1] = y[i]; 17 else l[i >> 1] = y[i]; 18 FFT(l, n >> 1, o); FFT(r, n >> 1, o); 19 cp omegan(cos(2 * pi / n), sin(2 * pi * o / n)), omega(1, 0); 20 for (int i = 0; i < n >> 1; i++) { 21 y[i] = l[i] + omega * r[i]; 22 y[i + (n >> 1)] = l[i] - omega * r[i]; 23 omega *= omegan; 24 } 25 } 26 27 int main() { 28 scanf("%d %d", &n, &m); 29 for (int i = 0; i <= n; i++) 30 scanf("%lf", &A[i].real()); 31 for (int i = 0; i <= m; i++) 32 scanf("%lf", &B[i].real()); 33 m += n; 34 for (n = 1; n <= m; n <<= 1); 35 FFT(A, n, 1); FFT(B, n, 1); 36 for (int i = 0; i <= n; i++) 37 A[i] *= B[i]; 38 FFT(A, n, -1); 39 for (int i = 0; i <= m; i++) 40 printf("%d ", (int)(A[i].real() / n + 0.5)); 41 return 0; 42 }
迭代版本
1 #include <cstdio> 2 #include <cmath> 3 #include <complex> 4 #include <iostream> 5 using namespace std; 6 7 const int N = 1 << 18; 8 typedef complex<double> cp; 9 const double pi = acos(-1.0); 10 cp A[N], B[N]; 11 bool flag; 12 int a[N], b[N], n, m, tar[N], bit; 13 14 inline void read(int &ans) { 15 static char buf = getchar(); 16 ans = 0; 17 for (; !isdigit(buf); buf = getchar()); 18 for (; isdigit(buf); buf = getchar()) 19 ans = ans * 10 + buf - '0'; 20 } 21 22 inline int rev(int val) { 23 int rst = 0; 24 for (int i = 0; i < bit; i++) { 25 rst <<= 1; rst |= val & 1; val >>= 1; 26 } return rst; 27 } 28 29 inline void FFT(cp *y) { 30 for (int i = 1; i <= bit; i++) { 31 int fac = 1 << i; 32 cp omegan(cos(2 * pi / fac), sin(2 * pi / fac)); 33 if (flag) omegan.imag() *= -1; 34 for (int j = 0; j < n; j += fac) { 35 cp omega(1, 0); 36 for (int k = 0; k < fac >> 1; k++) { 37 cp t = omega * y[j + k + (fac >> 1)]; 38 cp u = y[j + k]; y[j + k] = u + t; 39 y[j + k + (fac >> 1)] = u - t; 40 omega *= omegan; 41 } 42 } 43 } 44 } 45 int main() { 46 read(n); read(m); n++; m++; 47 for (int i = 0; i < n; i++) read(a[i]); 48 for (int i = 0; i < m; i++) read(b[i]); 49 m += n; for (n = 1; n < m; n <<= 1) bit++; 50 for (int i = 0; i < n; i++) tar[i] = rev(i); 51 for (int i = 0; i < n; i++) A[i].real() = a[tar[i]]; 52 for (int i = 0; i < n; i++) B[i].real() = b[tar[i]]; 53 FFT(A); FFT(B); 54 for (int i = 0; i < n; i++) A[i] *= B[i]; 55 for (int i = 0; i < n; i++) if (i < tar[i]) swap(A[i], A[tar[i]]); 56 flag = true; FFT(A); 57 for (int i = 0; i < m - 1; i++) 58 printf("%.0lf ", 0.0001 + A[i].real() / n); 59 puts(""); 60 return 0; 61 }