洛谷P4245：【模板】MTT （CRT+三模数NTT）

题目传送门：https://www.luogu.org/problemnew/show/P4245

---

题目分析：一道任意模数多项式乘法的模板题。可以写拆项+FFT，或者三模数NTT。我暂时只写了后者。

具体做法是这样：先选取三个乘积在$10^{23}$以上的便于使用NTT的模数。在这里我选的是$m_1=998244353=2^{23}*119+1$，$m_2=1004535809=2^{21}*479+1$，$m_3=469762049=2^{26}*7+1$。选这三个模数的好处在于它们的原根都是$3$。

然后用这三个模数做NTT，可以得到以下三条式子：

$$ans\equiv c_1(\mod m_1)$$

$$ans\equiv c_2(\mod m_2)$$

$$ans\equiv c_3(\mod m_3)$$

虽然这三条式子可以在$10^{23}$以内唯一固定$ans$的值，但问题也随之而来：$m_1*m_2*m_3$很大，无法直接用long long存下，而用long double之类的则会丢失精度，所以无法用普通的CRT。难道要写高精度？不，有一种很妙的方法可以解决这个问题。

首先注意到这里只有三个模数，而且两个模数乘起来是不会爆long long的，所以可以先合并前两条式子。根据CRT，有：

$$ans\equiv (c_1m_2Inv(m_2,m_1)+c_2m_1Inv(m_1,m_2))(\mod m_1m_2)$$

其中$Inv(x,y)$表示x关于y的逆元。

这条式子涉及到两个很大的数相乘然后再取模，而直接相乘会爆long long。可以用$O(\log(m_1m_2))$的快速乘，或者$O(1)$转double后相乘。

为了方便，把上式化成这样的形式：

$$ans\equiv C(\mod M)$$

然后设：

$$ans=xM+C=ym_3+c_3$$

接下来的部分才是精髓。我们求出$x$在$\mod m_3$意义下的值：

$$xM\equiv c_3-C(\mod m_3)$$

在$\mod m_3$意义下，$ym_3$被消掉了。然后有：

$$x\equiv (c_3-C)M^{-1}(\mod m_3)$$

算出右半部分的值为$q$，则可令$x=km_3+q$。将其代入$ans=xM+C$：

$$ans=(km_3+q)M+C=km_3M+qM+C$$

也就是说：

$$ans=km_1m_2m_3+qM+C$$

而由于$ans\in[0,m_1m_2m_3)$，所以$k$必为0！也就是说$ans$就是$qM+C$！直接把这条式子对题面要求的模数取模即可！

参考文章：sro AntiLeaf

（其实dalao的blog已经写得很详细了，只不过为了加深记忆我还是再写了一遍QAQ）

---

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=1000000;
const long long M[3]={1004535809,998244353,469762049};
const long long g=3;
typedef long long LL;

LL val[3][maxn];

LL A[maxn];
LL B[maxn];

int Rev[maxn];
int N,Lg;

int F[maxn];
int G[maxn];

int n,m;
LL p;

LL Pow(LL x,LL y,LL Mod)
{
    if (!y) return 1LL;
    LL temp=Pow(x,y>>1,Mod);
    temp=temp*temp%Mod;
    if (y&1) temp=temp*x%Mod;
    return temp;
}

LL Inv0=Pow(M[1]%M[0],M[0]-2LL,M[0]);
LL Inv1=Pow(M[0]%M[1],M[1]-2LL,M[1]);
LL MM=M[0]*M[1];
LL Inv2=Pow(MM%M[2],M[2]-2LL,M[2]);

void DFT(LL *a,int f,LL Mod)
{
    for (int i=0; i<N; i++)
        if (i<Rev[i]) swap(a[i],a[ Rev[i] ]);

    for (int len=2; len<=N; len<<=1)
    {
        int mid=(len>>1);
        LL e=Pow(g,(Mod-1LL)/len,Mod);
        if (f==-1) e=Pow(e,Mod-2LL,Mod);

        for (LL *p=a; p!=a+N; p+=len)
        {
            LL wn=1LL;
            for (int i=0; i<mid; i++)
            {
                LL temp=wn*p[mid+i]%Mod;
                p[mid+i]=(p[i]-temp+Mod)%Mod;
                p[i]=(p[i]+temp)%Mod;
                wn=wn*e%Mod;
            }
        }
    }
}

void NTT(int x)
{
    for (int i=0; i<N; i++) A[i]=F[i]%M[x],B[i]=G[i]%M[x];
    DFT(A,1,M[x]);
    DFT(B,1,M[x]);
    for (int i=0; i<N; i++) A[i]=A[i]*B[i]%M[x];
    DFT(A,-1,M[x]);
    LL inv=Pow(N,M[x]-2LL,M[x]);
    for (int i=0; i<N; i++) val[x][i]=A[i]*inv%M[x];
}

LL Mul(LL x,LL y,LL Mod)
{
    if (!y) return 0LL;
    LL temp=Mul(x,y>>1,Mod);
    temp=(temp+temp)%Mod;
    if (y&1) temp=(temp+x)%Mod;
    return temp;
}

LL Get(LL v0,LL v1,LL v2)
{
    LL temp1=Mul(v0*M[1]%MM,Inv0,MM);
    LL temp2=Mul(v1*M[0]%MM,Inv1,MM);
    LL temp=(temp1+temp2)%MM,C=temp;
    temp=(v2-temp%M[2]+M[2])%M[2]*Inv2%M[2];
    temp=(MM%p*temp%p+C%p)%p;
    return temp;
}

int main()
{
    freopen("4245.in","r",stdin);
    freopen("4245.out","w",stdout);

    scanf("%d%d%I64d",&n,&m,&p);
    for (int i=0; i<=n; i++) scanf("%d",&F[i]),F[i]%=p;
    for (int i=0; i<=m; i++) scanf("%d",&G[i]),G[i]%=p;

    N=1,Lg=0;
    while (N<n+m+4) N<<=1,Lg++;
    for (int i=0; i<N; i++)
        for (int j=0; j<Lg; j++)
            if (i&(1<<j)) Rev[i]|=(1<<(Lg-j-1));

    for (int i=0; i<3; i++) NTT(i);
    for (int i=0; i<=n+m; i++) printf("%I64d ", Get(val[0][i],val[1][i],val[2][i]) );
    printf("\n");

    return 0;
}