【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;
} 

__EOF__

本文作者：あおいSakura
本文链接：https://www.cnblogs.com/Sakura-TJH/p/luogu_P4245.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！