洛谷P4245:【模板】MTT (CRT+三模数NTT)

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


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

具体做法是这样:先选取三个乘积在1023以上的便于使用NTT的模数。在这里我选的是m1=998244353=223119+1m2=1004535809=221479+1m3=469762049=2267+1。选这三个模数的好处在于它们的原根都是3

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

ansc1(modm1)

ansc2(modm2)

ansc3(modm3)

虽然这三条式子可以在1023以内唯一固定ans的值,但问题也随之而来:m1m2m3很大,无法直接用long long存下,而用long double之类的则会丢失精度,所以无法用普通的CRT。难道要写高精度?不,有一种很妙的方法可以解决这个问题。

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

ans(c1m2Inv(m2,m1)+c2m1Inv(m1,m2))(modm1m2)

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

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

为了方便,把上式化成这样的形式:

ansC(modM)

然后设:

ans=xM+C=ym3+c3

接下来的部分才是精髓。我们求出xmodm3意义下的值:

xMc3C(modm3)

modm3意义下,ym3被消掉了。然后有:

x(c3C)M1(modm3)

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

ans=(km3+q)M+C=km3M+qM+C

也就是说:

ans=km1m2m3+qM+C

而由于ans[0,m1m2m3),所以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;
}
发布了160 篇原创文章 · 获赞 76 · 访问量 10万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 技术黑板 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览