[Luogu 4245] 任意模数NTT

Description

给定 \(2\) 个多项式 \(F(x), G(x)\)，请求出 \(F(x) * G(x)\)。

系数对 \(p\) 取模，且不保证 \(p\) 可以分解成 \(p = a \cdot 2^k + 1\) 之形式。

Solution

设 \(m_1=469762049,m_2=998244353,m_3=1004535809​\)，有

\[ \begin{cases} x\equiv c_1\pmod{m_1}\\ x\equiv c_2\pmod{m_2}\\ x\equiv c_3\pmod{m_3} \end{cases} \]

用中国剩余定理合并前两个同余式，得到

\[ \begin{cases} x\equiv c_3\pmod{m_3}\\ x\equiv c_4\pmod{m_1m_2} \end{cases} \]

设 \(x=km_1m_2+c_4\)，有

\[ km_1m_2+c_4\equiv c_3\pmod{m_3}\\ k\equiv (c_3-c_4)m_1^{-1}m_2^{-1}\pmod{m_3} \]

设 \(k=am_3+(c_3-c_4)m_1^{-1}m_2^{-1}​\)，有

\[ x=(am_3+(c_3-c_4)m_1^{-1}m_2^{-1})m_1m_2+c_4\\ x\equiv (c_3-c_4)m_1^{-1}m_2^{-1}m_1m_2+c_4\pmod{m_1m_2m_3} \]

其中 \(m_1^{-1}m_2^{-1}\) 是在模 \(m_3\) 意义下的。

注意 \(exgcd\) 的返回值可能是负数，要处理一下；在 \(ksm(a, b, p)\) 之前先将 \(a\) 对 \(p\) 取模。

Code

#include <cstdio>
#include <algorithm>

typedef long long LL; const int N = 262150;
int a[N], b[N], n, m, nn, mm, pp, R[N], L, p[N], g[N]; LL f[2][N];

int read() {
    int x = 0; char c = getchar();
    while (c < '0' || c > '9') c = getchar();
    while (c >= '0' && c <= '9') x = (x << 3) + (x << 1) + (c ^ 48), c = getchar();
    return x;
}
void exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) { x = 1, y = 0; return; }
    exgcd(b, a % b, y, x), y -= a / b * x;
}
LL ksm(LL a, LL b, LL p) {
    LL res = 1;
    for (; b; b >>= 1, a = 1LL * a * a % p)
        if (b & 1) res = 1LL * res * a % p;
    return res;
}
LL mul(LL a, LL b, LL p) {
    LL res = 0; int f = 1;
    if (a < 0) a = -a, f = -f; if (a >= p) a %= p;
    if (b < 0) b = -b, f = -f; if (b >= p) b %= p;
    for (; b; b >>= 1, a += a + a < p ? a : a - p)
        if ((b & 1) && (res += a) >= p) res -= p;
    return res * f; 
}
void NTT(LL *A, int f, int t) {
    for (int i = 0; i < n; ++i) if (i < R[i]) std::swap(A[i], A[R[i]]);
    for (int i = 1; i < n; i <<= 1) {
        int wn = ksm(f ? 3 : g[t], (p[t] - 1) / (i << 1), p[t]);
        for (int j = 0, r = i << 1; j < n; j += r) {
            int w = 1;
            for (int k = 0; k < i; ++k, w = 1LL * w * wn % p[t]) {
                int x = A[j + k], y = 1LL * w * A[i + j + k] % p[t];
                A[j + k] = (x + y) % p[t], A[i + j + k] = (x - y + p[t]) % p[t];
            }
        }
    }
}
void solve(int t, int k) {
    LL c[N] = {}, d[N] = {};
    for (int i = 0; i <= nn; ++i) c[i] = a[i];
    for (int i = 0; i <= mm; ++i) d[i] = b[i];
    NTT(c, 1, t), NTT(d, 1, t);
    for (int i = 0; i < n; ++i) f[k][i] = 1LL * c[i] * d[i] % p[t];
    NTT(f[k], 0, t);
    int inv = ksm(n, p[t] - 2, p[t]);
    for (int i = 0; i <= m; ++i) f[k][i] = 1LL * f[k][i] * inv % p[t];
}
int main() {
    nn = read(), mm = read(), pp = read();
    for (int i = 0; i <= nn; ++i) a[i] = read();
    for (int i = 0; i <= mm; ++i) b[i] = read();
    m = nn + mm; for (n = 1; n <= m; n <<= 1) ++L;
    for (int i = 0; i < n; ++i) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
    p[1] = 469762049, p[2] = 998244353, p[3] = 1004535809;
    g[1] = 156587350, g[2] = 332748118, g[3] = 334845270;
    solve(1, 0), solve(2, 1);
    LL mod = 1LL * p[1] * p[2], x, y;
    exgcd(p[1], p[2], x, y), x = (x % mod + mod) % mod;
    for (int i = 0; i <= m; ++i) f[0][i] = (f[0][i] + mul(mul(x, f[1][i] - f[0][i] + mod, mod), p[1], mod)) % mod;
    solve(3, 1);
    LL inv = ksm(mod % p[3], p[3] - 2, p[3]);
    for (int i = 0; i <= m; ++i) printf("%lld ", (mul(((f[1][i] - f[0][i]) % p[3] + p[3]) * inv % p[3], mod, pp) + f[0][i]) % pp);
    return 0;
}

转载于:https://www.cnblogs.com/fly-in-milkyway/p/10554396.html