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

关于FFT

这个博客的讲解超级棒

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

算法导论上的讲解也不错

模板就是抄一抄别人的啦

首先是递归版本

#include <cstdio>
#include <complex>
#include <cmath>
using namespace std;

const double pi = acos(-1);
const int N = (2 << 17) + 10;
typedef complex<double> cp;
cp A[N], B[N];
int n, m;

void FFT(cp *y, int n, int o) {
    if (n == 1)    return ;
    cp l[n >> 1], r[n >> 1];
    for (int i = 0; i <= n; i++)
        if (i & 1)    r[i >> 1] = y[i]; 
        else    l[i >> 1] = y[i];
    FFT(l, n >> 1, o); FFT(r, n >> 1, o);
    cp omegan(cos(2 * pi / n), sin(2 * pi * o / n)), omega(1, 0);
    for (int i = 0; i < n >> 1; i++) {
        y[i] = l[i] + omega * r[i];
        y[i + (n >> 1)] = l[i] - omega * r[i];
        omega *= omegan;
    }
}

int main() {
    scanf("%d %d", &n, &m);
    for (int i = 0; i <= n; i++)
        scanf("%lf", &A[i].real());
    for (int i = 0; i <= m; i++)
        scanf("%lf", &B[i].real());
    m += n;
    for (n = 1; n <= m; n <<= 1);
    FFT(A, n, 1); FFT(B, n, 1);
    for (int i = 0; i <= n; i++)
        A[i] *= B[i];
    FFT(A, n, -1);
    for (int i = 0; i <= m; i++)
        printf("%d ", (int)(A[i].real() / n + 0.5));
    return 0;
}

迭代版本

#include <cstdio>
#include <cmath>
#include <complex>
#include <iostream>
using namespace std;

const int N = 1 << 18;
typedef complex<double> cp;
const double pi = acos(-1.0);
cp A[N], B[N];
bool flag;
int a[N], b[N], n, m, tar[N], bit;

inline void read(int &ans) {
    static char buf = getchar();
    ans = 0;
    for (; !isdigit(buf); buf = getchar());
    for (; isdigit(buf); buf = getchar())
        ans = ans * 10 + buf - '0';
}

inline int rev(int val) {
    int rst = 0;
    for (int i = 0; i < bit; i++) {
        rst <<= 1; rst |= val & 1; val >>= 1; 
    } return rst;        
}

inline void FFT(cp *y) {
    for (int i = 1; i <= bit; i++) {
        int fac = 1 << i;
        cp omegan(cos(2 * pi / fac), sin(2 * pi / fac));
        if (flag)    omegan.imag() *= -1;
        for (int j = 0; j < n; j += fac) {
            cp omega(1, 0);
            for (int k = 0; k < fac >> 1; k++) {
                cp t = omega * y[j + k + (fac >> 1)];
                cp u = y[j + k]; y[j + k] = u + t;
                y[j + k + (fac >> 1)] = u - t;
                omega *= omegan;
            }
        }
    }
}
int main() {
    read(n); read(m); n++; m++;
    for (int i = 0; i < n; i++)    read(a[i]);
    for (int i = 0; i < m; i++)    read(b[i]);
    m += n; for (n = 1; n < m; n <<= 1)    bit++;
    for (int i = 0; i < n; i++)    tar[i] = rev(i);
    for (int i = 0; i < n; i++)    A[i].real() = a[tar[i]];
    for (int i = 0; i < n; i++)    B[i].real() = b[tar[i]];
    FFT(A); FFT(B);
    for (int i = 0; i < n; i++)    A[i] *= B[i];
    for (int i = 0; i < n; i++)    if (i < tar[i])    swap(A[i], A[tar[i]]);
    flag = true; FFT(A);
    for (int i = 0; i < m - 1; i++) 
        printf("%.0lf ", 0.0001 + A[i].real() / n);
       puts("");
    return 0;
}