luogu P4245 [模板]任意模数NTT

传送门

orz 远航之曲

因为答案最大可以取到\(10^{23}\),所以要选择三个乘积\(\ge 10^{23}\)的质数(见代码),分别算一遍,然后crt起来

注意这里要这样crt:

首先三个同余式为\(ans\equiv a1(mod\ m1)\)\(ans\equiv a2(mod\ m2)\)\(ans\equiv a3(mod\ m3)\)

记\(M=m1*m2,M_i=\frac{M}{m_i}\),\(t_i\)意义是\(t_i*M_i\equiv 1(mod\ m_i)\)

先将前两个合并,得到\(aa=a1t1m2+a2t2m1(mod M)\)

然后\(ans=bbM+aa=m3\ k+a3\)

然后\(bbM\equiv a3-aa(mod\ m3)\)

然后\(bb\equiv (a3-aa)M^{-1}(mod\ m3)\)

最后\(ans=bbM+aa(mod\ md)\)

#include<bits/stdc++.h>
#define LL long long
#define ldb long double
#define il inline
#define re register
  
using namespace std;
const int M=270000+10;
const LL m1=469762049,m2=998244353,m3=1004535809;
il int rd()
{
  int x=0,w=1;char ch=0;
  while(ch<'0'||ch>'9') {if(ch=='-') w=-1;ch=getchar();}
  while(ch>='0'&&ch<='9') {x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
  return x*w;
}
int n,m,md,l,rdr[M],a[M],b[M],c[M],d[M],an[3][M],inv;
il LL fmul(LL a,LL b,LL mod)
{
  LL an=a*b-(LL)((ldb)a*b/mod+0.5)*mod;
  return an<0?an+mod:an;
}
il LL fpow(LL a,LL b,LL mod)
{
  LL an=1;
  while(b){if(b&1) an=an*a%mod;a=a*a%mod,b>>=1;}
  return an;
}
void ntt(int *a,int op,int mod)
{
  int W,w,x,y;
  for(int i=0;i<n;++i) if(i<rdr[i]) swap(a[i],a[rdr[i]]);
  for(int i=1;i<n;i<<=1)
    {
      W=fpow(3,(mod-1)/(i<<1),mod);
      if(op<0) W=fpow(W,mod-2,mod);
      for(int j=0;j<n;j+=i<<1)
        {
          w=1;
          for(int k=0;k<i;++k,w=1ll*w*W%mod)
            {
              x=a[j+k],y=1ll*w*a[j+k+i]%mod;
              a[j+k]=(x+y)%mod,a[j+k+i]=(x-y+mod)%mod;
            }
        }
    }
    
}
  
int main()
{
  n=rd(),m=rd(),md=rd();
  for(int i=0;i<=n;++i) a[i]=rd();
  for(int i=0;i<=m;++i) b[i]=rd();
  m+=n;
  for(n=1;n<=m;n<<=1) ++l;
  for(int i=0;i<n;++i) rdr[i]=(rdr[i>>1]>>1)|((i&1)<<(l-1));//
  memcpy(c,a,4*(n+1)),memcpy(d,b,4*(n+1));
  ntt(c,1,m1),ntt(d,1,m1);
  for(int i=0;i<n;++i) an[0][i]=1ll*c[i]*d[i]%m1;
  ntt(an[0],-1,m1);
  inv=fpow(n,m1-2,m1);
  for(int i=0;i<=m;++i) an[0][i]=1ll*an[0][i]*inv%m1;//
  memcpy(c,a,4*(n+1)),memcpy(d,b,4*(n+1));
  ntt(c,1,m2),ntt(d,1,m2);
  for(int i=0;i<n;++i) an[1][i]=1ll*c[i]*d[i]%m2;
  ntt(an[1],-1,m2);
  inv=fpow(n,m2-2,m2);
  for(int i=0;i<=m;++i) an[1][i]=1ll*an[1][i]*inv%m2;//
  memcpy(c,a,4*(n+1)),memcpy(d,b,4*(n+1));
  ntt(c,1,m3),ntt(d,1,m3);
  for(int i=0;i<n;++i) an[2][i]=1ll*c[i]*d[i]%m3;
  ntt(an[2],-1,m3);
  inv=fpow(n,m3-2,m3);
  for(int i=0;i<=m;++i) an[2][i]=1ll*an[2][i]*inv%m3;//
  for(int i=0;i<=m;++i)
    {
      LL aa=(fmul(1ll*an[0][i]*m2%(m1*m2),fpow(m2%m1,m1-2,m1),m1*m2)+fmul(1ll*an[1][i]*m1%(m1*m2),fpow(m1%m2,m2-2,m2),m1*m2))%(m1*m2);
      LL bb=(((an[2][i]-aa)%m3+m3)%m3*fpow((m1*m2)%m3,m3-2,m3))%m3;
      printf("%lld ",((bb%md)*(m1*m2%md)+(aa%md))%md);
    }
  return 0;
}