UOJ 34: 多项式乘法（FFT模板题）

关于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 }

 

转载于:https://www.cnblogs.com/cminus/p/7278928.html