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

最新推荐文章于 2022-03-14 19:03:25 发布

SSL_TJH

最新推荐文章于 2022-03-14 19:03:25 发布

阅读量136

点赞数

分类专栏：模板题 # 数学或数论文章标签：拆系数FFT MTT 多项式

要转发的话记得附上链接哦~（不过也不会有人转发的啦）

本文链接：https://blog.csdn.net/weixin_43346722/article/details/121155688

版权

数学或数论同时被 2 个专栏收录

285 篇文章 3 订阅

订阅专栏

模板题

154 篇文章 2 订阅

订阅专栏

该博客介绍了如何使用拆系数FFT（MTT）方法解决任意模数多项式乘法问题，通过将数拆分成两部分进行运算，减少计算次数。文章详细阐述了优化过程，包括复数运算、DFT与IDFT的应用，以及最终如何将结果整合。代码实现中展示了C++的处理过程，强调了预处理单位根以进一步优化效率。

摘要由CSDN通过智能技术生成

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

题目链接：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 + b i) * (c + d i) = a c - b d + (a d + b c) i$
$(a - b i) * (c + d i) = a c + b d + (a d - b c) i$

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

然后你考虑加在一起，就是 $2 (a c + a d * 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;
}