快速傅立叶变换(FFT) 多项式乘法模板

#include <algorithm>
#include <stdexcept>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>

const int N = 1e5 + 5;
const long double PI = acos(-1.0);
const long double EPS = 1E-8;
class complex {
  public:
    long double re, im;
    complex() {}
    complex(long double re, long double im) : re(re), im(im) {}
    complex operator + (const complex &x) { return complex(re + x.re, im + x.im); }
    complex operator - (const complex &x) { return complex(re - x.re, im - x.im); }
    complex operator * (const complex &x) { return complex(re * x.re - im * x.im, im * x.re + re * x.im); }
    complex operator / (const complex &x) {
        return complex((re * x.re + im * x.im) / (x.re * x.re + x.im * x.im),
                       (im * x.re - re * x.im) / (x.re * x.re + x.im * x.im));
    }
};

void print(complex x, char c) {
    printf("\t%.10lf", fabs(x.re) < EPS ? 0.0 : (double)x.re);
    printf("\t\t(%.10lf)", fabs(x.im) < EPS ? 0.0 : (double)x.im);
    printf("%c", c);
}

int n, rev[N];
complex F[2][N], w[N];

int bitReverse(int x, int y, int r = 0) {
    for (int i = 0; i < y; i++) r = (r << 1) ^ (x & 1), x >>= 1;
    return r;
}

void FFT(complex * F, int n, int oft) {
    for (int i = 0; i < n; i++) if (rev[i] > i) std::swap(F[i], F[rev[i]]);
    for (int i = 2; i <= n; i <<= 1) {
        complex wi(cos(oft * 2 * PI / i), sin(oft * 2 * PI / i));
        for (int j = 0; j < n; j += i) {
            complex w(1, 0);
            for (int k = j, h = i >> 1; k < j + h; k++) {
                complex t = w * F[k + h], u = F[k];
                F[k] = u + t;
                F[k + h] = u - t;
                w = w * wi;
            }
        }
    }
    if (oft == -1) for (int i = 0; i < n; i++) F[i].re = F[i].re / n;
}

complex ans[N];

int main() {

    int mx = 0;
    for (int j = 0; j < 2; j++) {
        scanf("%d", &n);
        mx = std::max(mx, n);
        for (int i = 0; i < n; i++) {
            double x;
            scanf("%lf", &x);
            F[j][i].re = x;
        }
    }
    n = 1 << (int)ceil(log2(mx));
    int BIT = log2(2 * n);
    for (int i = 0; i < 2 * n; i++) rev[bitReverse(i, BIT)] = i;

    FFT(F[0], 2 * n, 1);
    FFT(F[1], 2 * n, 1);
    for (int i = 0; i < 2 * n; i++) ans[i] = F[0][i] * F[1][i];
    FFT(ans, 2 * n, -1);
    for (int i = 0; i < 2 * n; i++) print(ans[i], '\n');

    return 0;
}


/*

(4 - x + 9x^2 - 2x^3) * (2x - 5x^2)
= 0 + 8x - 22x^2 + 23x^3 - 49x^4 + 10x^5

4
4 -1 9 -2
3
0 2 -5

*/

Code language: PHP (php)