【luogu P4245】【模板】任意模数多项式乘法（拆系数FFT）（MTT）

【模板】任意模数多项式乘法

题目链接：luogu P4245

题目大意

给你两个多项式，要你求它们的卷积。
然后模数任意。

思路

这里用的是拆系数 FFT，即 MTT。

我们把一个数拆成 (a*2^{15}+b) 的形式（不一定是 (2^{15})，大概在 (sqrt{_{值域}}) 即可）

那两个数相乘：(c_1*c_2=(a_1*2^{15}+b_1)*(a_2*2^{15}+b_2)=a_1a_2*2^{30}+(a_1b_2+a_2b_1)*2^{15}+b_1b_2)。

那我们就只用做 (4) 次多项式乘法就可以啦~
（然而 (12) 次 FFT，直接炸飞）

而且就算你一开始的 FFT 只用 (4) 次，也要 (8) 次，难以接受。

考虑搞点优化，用一些式子来优化。
((a+bi)*(c+di)=ac-bd+(ad+bc)i)
((a-bi)*(c+di)=ac+bd+(ad-bc)i)

会发现我们让 (P=a_1+a_2i,P'=a_1-a_2i,Q=b_1+b_2i)。
那上面的两个式子就是 (PQ,P'Q)。

然后你考虑加在一起，就是 (2(ac+ad*i))，带入到这里就是 (2(a_1b_1+a_1b_2*i))。
那我们不难看出除二就可以求出这两个。

然后分别用上面的两个式子减去 (a_1b_1,a_1b_2) 就可以得到 (a_2b_2,a_2b_1) 了。
（用 (i) 项那边减）

然后带入回 (a_1a_2*2^{30}+(a_1b_2+a_2b_1)*2^{15}+b_1b_2) 就有答案了。

然后中间搞的过程会到 (10^{19})（IDFT 之前），所以要 long double。

不难看出现在就变成了 (5) 次 FFT（三次 DFT，两次 IDFT）。

其实还可以预处理单位根来再次优化。

代码

#include<cmath>
#include<cstdio>
#include<algorithm>
#define lim 32000
#define ll long long

using namespace std;

struct complex {
	long double x, y;
	complex (long double xx = 0, long double yy = 0) {
		x = xx; y = yy;
	}
}p1[100001 << 2], p2[100001 << 2], q[100001 << 2];
int n, m, limit, l_size, an[100001 << 2];
ll mo, x, ans[100001 << 2];
long double Pi = acos(-1.0);

complex operator +(complex x, complex y) {
	return (complex){x.x + y.x, x.y + y.y};
}

complex operator -(complex x, complex y) {
	return (complex){x.x - y.x, x.y - y.y};
}

complex operator *(complex x, complex y) {
	return (complex){x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x};
}

void FFT(complex *now, int op) {
	for (int i = 0; i < limit; i++)
		if (an[i] > i) swap(now[i], now[an[i]]);
	for (int mid = 1; mid < limit; mid <<= 1) {
		complex Wn(cos(Pi / mid), op * sin(Pi / mid));
		for (int R = (mid << 1), j = 0; j < limit; j += R) {
			complex w(1, 0);
			for (int k = 0; k < mid; k++, w = w * Wn) {
				complex x = now[j + k], y = w * now[j + mid + k];
				now[j + k] = x + y;
				now[j + mid + k] = x - y;
			}
		}
	}
}

int main() {
	scanf("%d %d %lld", &n, &m, &mo);
	for (int i = 0; i <= n; i++) {
		scanf("%lld", &x);
		p1[i] = (complex){x / lim, x % lim};//拆开
		p2[i] = (complex){x / lim, -x % lim};
	}
	for (int i = 0; i <= m; i++) {
		scanf("%lld", &x);
		q[i] = (complex){x / lim, x % lim};
	}
	
	limit = 1;
	while (limit <= n + m) {
		limit <<= 1;
		l_size++;
	}
	for (int i = 0; i < limit; i++)
		an[i] = (an[i >> 1] >> 1) | ((i & 1) << (l_size - 1));
	
	FFT(p1, 1); FFT(p2, 1); FFT(q, 1);
	for (int i = 0; i < limit; i++) q[i].x /= limit, q[i].y /= limit;//先除了后面就不用除
	for (int i = 0; i < limit; i++) p1[i] = p1[i] * q[i], p2[i] = p2[i] * q[i];
	FFT(p1, -1); FFT(p2, -1);
	
	for (int i = 0; i <= n + m; i++) {
		ll a1b1 = (ll)floor((p1[i].x + p2[i].x) / 2 + 0.5) % mo;//两个的和除二
		ll a1b2 = (ll)floor((p1[i].y + p2[i].y) / 2 + 0.5) % mo;
		ll a2b1 = ((ll)floor(p1[i].y + 0.5) - a1b2) % mo;
		ll a2b2 = ((ll)floor(p2[i].x + 0.5) - a1b1) % mo;
		ans[i] = (a1b1 * lim * lim + (a1b2 + a2b1) * lim + a2b2) % mo;
		ans[i] = (ans[i] + mo) % mo;//带进去式子里面算
		printf("%lld ", ans[i]);
	}
	
	return 0;
}