「任意模数多项式乘法」
前置知识
基本问题
给定一个 \(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;
}