任意模数多项式乘法（洛谷4245）

代码：

#define _CRT_SECURE_NO_WARNINGS
#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<utility>
#include<stack>
#include<queue>
#include<vector>
#include<set>
#include<map>
#include<unordered_map>
using namespace std;
#define LL  long long
const int N = 5e5 + 10;
const double PI = acos(-1);
LL p,res[4][N], mod[5] = {0,998244353,1004535809,2281701377},gg=3,gi;
int n,m,tot,bit;
int rel[N];
LL f[N], g[N];
LL mul(LL a, LL b, LL mod)
{
    LL ans = (a * b - (LL)((long double)a / mod * b + 1e-8) * mod);
    return ans < 0 ? ans + mod : ans;
}
LL quick(LL a, LL b, LL mod)
{
    LL ans = 1;
    while (b)
    {
        if (b & 1) ans = mul(ans,a,mod);
        b = b >> 1;
        a = mul(a,a,mod);
    }
    return ans;
}
void ntt(LL a[], LL tot, int op,LL mod)
{
    for (int i = 0; i < tot; i++) if (i < rel[i]) swap(a[i], a[rel[i]]);
    for (int m = 2; m <= tot; m <<= 1)
    {
        LL g1 = quick(((op == 1) ? gg : gi), (mod - 1) / m, mod);
        for (int i = 0; i < tot; i += m)
        {
            LL gk = 1;
            for (int j = 0; j < m / 2; j++)
            {
                LL x = a[i + j], y =a[i + j + m / 2]*gk%mod;
                a[i + j] = (x+y)%mod;
                a[i + j + m / 2] = ((x - y) % mod + mod) % mod;
                gk = gk*g1%mod;
            }
         }
    }
    if (op != 1)
    {
        LL gk = quick(tot, mod - 2, mod);
        for (int i = 0; i < tot; i++)
            a[i] = a[i]*gk%mod;
    }
}
void ntt(LL a[], LL b[], LL ans[], LL tot, LL mod)
{
    for (int i = 0; i < tot; i++) f[i] = a[i], g[i] = b[i];
    gi = quick(gg, mod - 2, mod);
    ntt(f, tot, 1,mod);
    ntt(g, tot, 1,mod);
    for (int i = 0; i < tot; i++) ans[i] = f[i]*g[i]%mod;
    ntt(ans, tot, -1,mod);
}
LL exgcd(LL a, LL b, LL& x, LL &y,LL mod)
{
    if (!b)
    {
        x = 1, y = 0;
        return a;
    }
    LL d = exgcd(b, a % b, y, x,mod);
    y = y -(a/b)*x;
    return d;
}
LL china(int op)
{     
    LL pp = 1, ans=0;
    for (int i = 1; i <= 2; i++) pp *= mod[i];
    for (int i = 1; i <= 2; i++)
    {
        LL mm = pp / mod[i];
        LL x, y;
        exgcd(mm, mod[i], x, y,pp);
        x = (x % pp + pp) % pp;
        ans = (ans + mul(res[i][op], mul(mm, x, pp), pp)%pp)%pp;
    }
    return ans;
}
LL a[N], b[N];
int main() {
    cin >> n >> m >> p;
    for (int i = 0; i <= n; i++) scanf("%lld",&a[i]);
    for (int i = 0; i <= m; i++) scanf("%lld", &b[i]);
    while ((1 << bit) < n + m + 1) bit++;
    tot = 1 << bit;
    for (int i = 0; i < tot; i++) rel[i] = rel[i / 2] / 2 + ((i & 1) ? tot / 2 : 0);
    for (int i = 1; i <= 3; i++)
        ntt(a, b, res[i], tot, mod[i]);
    for (int i = 0; i <= n + m; i++)
    {
        LL a4 = china(i);
        LL a5 = mul(res[3][i] - a4, quick(mod[1], mod[3] - 2, mod[3]), mod[3]);
            a5=mul(a5, quick(mod[2], mod[3] - 2, mod[3]),mod[3]);
            a5 = mul(a5, mod[1], p);
            a5 = mul(a5, mod[2], p);
        printf("%lld ", ((a5+ a4) % p + p) % p);
    }
    return 0;
}