「任意模数多项式乘法」

「任意模数多项式乘法」

前置知识

多项式乘法

基本问题

给定一个 \(n\) 次多项式 \(F(x)\) 和一个 \(m\) 次多项式,求出

\[F(x)\times G(x) \]

系数对 \(p\) 取模,且不保证 \(p\) 可以分解成 \(p=2^ka+1\) 之形式,\(0\leq a_i,b_i\leq 10^9\),\(2\leq p\leq 10^9+9\)

考虑直接用 \(FFT\),但是值域太大,\(long\; double\) 都炸了,精度也无法保证

直接用 \(NTT\),但是是任意模数,根本用不了

处理这种问题,我们常有两种做法:

三模NTT

都 \(1202\) 年了,不会还有人写三模 \(NTT\) 吧

好吧,其实是我不会

另一种比较常用的做法是

MTT

既然 \(FFT\) 处理不了值域很大的情况,我们就从问题入手,将值域缩小

不妨将两个多项式拆成:

\[F(x)=M\times A(x)+B(x) \]

\[G(x)=M\times C(x)+D(x) \]

当 \(M=2^{15}\) 时,可以完美避免炸 \(double\) 的问题

现在问题就转化为了

\[(M\times A(x)+B(x))\times (M\times C(x)+D(x)) \]

\[M^2A(x)C(x)+M(B(x)C(x)+A(x)D(x))+B(x)D(x) \]

所以直接 \(7\) 次 \(FFT\) 就解决啦

\(PS\):可以省到 \(4\) 次 \(FFT\),但是我还不会

代码
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>

typedef long long ll;
typedef unsigned long long ull;

using namespace std;

const int maxn = 3e5 + 50, INF = 0x3f3f3f3f;
const double pi = acos (-1);

inline int read () {
	register int x = 0, w = 1;
	register char ch = getchar ();
	for (; ch < '0' || ch > '9'; ch = getchar ()) if (ch == '-') w = -1;
	for (; ch >= '0' && ch <= '9'; ch = getchar ()) x = x * 10 + ch - '0';
	return x * w;
}

inline void write (register int x) {
	if (x / 10) write (x / 10);
	putchar (x % 10 + '0');
}

int n, m, M, mod, len = 1, bit;
int rev[maxn], f[maxn];

struct Complex {
	double x, y;
	Complex () {}
	Complex (register double a, register double b) { x = a, y = b; }
	inline Complex operator + (const Complex &a) const { return Complex (x + a.x, y + a.y); }
	inline Complex operator - (const Complex &a) const { return Complex (x - a.x, y - a.y); }
	inline Complex operator * (const Complex &a) const { return Complex (x * a.y + y * a.x, y * a.y - x * a.x); }
} g[maxn], a[maxn], b[maxn], c[maxn], d[maxn], omega[maxn];

inline void FFT (register int len, register Complex * a, register int opt) {
	for (register int i = 1; i < len; i ++) if (i < rev[i]) swap (a[i], a[rev[i]]);
	for (register int d = 1; d < len; d <<= 1) {
		for (register int i = 0; i < len; i += d << 1) {
			for (register int j = 0; j < d; j ++) {
				register Complex w = omega[len / (d << 1) * j]; w.x *= opt;
				register Complex x = a[i + j], y = w * a[i + j + d];
				a[i + j] = x + y, a[i + j + d] = x - y;
			}
		}
	}
}

int main () {
	n = read(), m = read(), mod = read(), M = 1 << 15;
	while (len <= n + m) len <<= 1, bit ++;
	for (register int i = 0; i < len; i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << bit - 1), omega[i] = Complex (sin (2 * pi * i / len), cos (2 * pi * i / len));
	for (register int i = 0, x; i <= n; i ++) x = read(), a[i].y = x / M, b[i].y = x % M;
	for (register int i = 0, x; i <= m; i ++) x = read(), c[i].y = x / M, d[i].y = x % M;
	FFT (len, a, 1), FFT (len, b, 1), FFT (len, c, 1), FFT (len, d, 1);
	for (register int i = 0; i < len; i ++) g[i] = a[i] * c[i]; FFT (len, g, -1);
	for (register int i = 0; i <= n + m; i ++) f[i] = (f[i] + (ll) (g[i].y / len + 0.5) % mod * M % mod * M % mod) % mod;
	for (register int i = 0; i < len; i ++) g[i] = a[i] * d[i] + b[i] * c[i]; FFT (len, g, -1);
	for (register int i = 0; i <= n + m; i ++) f[i] = (f[i] + (ll) (g[i].y / len + 0.5) % mod * M % mod) % mod;
	for (register int i = 0; i < len; i ++) g[i] = b[i] * d[i]; FFT (len, g, -1);
	for (register int i = 0; i <= n + m; i ++) f[i] = (f[i] + (ll) (g[i].y / len + 0.5) % mod) % mod;
	for (register int i = 0; i <= n + m; i ++) printf ("%d ", f[i]); putchar ('\n');
	return 0;
}
上一篇:FFT模板


下一篇:matlab之FFT(离散傅里叶)分析原理