uoj34: 多项式乘法

题目链接
模板题目，写的是递归+std::complex
以后会补上迭代版本。

#include <cmath>
#include <cstdio>
#include <complex>
#include <iostream>
#include <algorithm>
using namespace std;
const double pi = acos(-1.0);
const int N = (2<<17) + 10;
typedef complex<double> cp;

int n, m;
cp A[N], B[N];

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

int main() {
	cin >> 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 <cmath>
#include <cstdio>
#include <complex>
#include <iostream>
#include <algorithm>
using namespace std;
const int N = 1<<18;
typedef complex<double> cp;
const double pi = acos(-1.0);
char buf;
inline short read() {
	do buf = getchar(); while (! isdigit(buf)); return buf - '0';
}

bool flag;
int n, m, bit, tar[N];
short a[N], b[N];
cp A[N], B[N];

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 omega_n(cos(2 * pi / fac), sin(2 * pi / fac));
		if (flag) omega_n.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)], u = y[j + k];
				y[j + k] = u + t; y[j + k + (fac >> 1)] = u - t;
				omega *= omega_n;
			}
		}
	}
}

int main() {
	scanf("%d%d", &n, &m); n++; m++;
	for (int i = 0; i < n; i++) a[i] = read();
	for (int i = 0; i < m; i++) b[i] = read();
	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);
	return puts(""), 0;
}