快速傅立叶变换(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)

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注