任意模数ntt_【模板篇】NTT和三模数NTT

之前写过FFT的笔记. 我们知道FFT是在复数域上进行的变换.

而且经过数学家的证明, DFT是复数域上唯一满足循环卷积性质的变换.

而我们在OI中, 经常遇到对xxxx取模的题目, 这就启发我们可不可以在模运算的意义下找一个这样的变换.

然后我们发现有个神奇的东西, 原根\(g\), 这东西在模意义下相当于单位复根\(-e^{\frac{2\pi i}{n}}\).

所以我们预处理一下\(g\)的幂和逆元, 然后改一下fft的代码就出现了快速数论变换ntt

懒得写了 直接上代码:

void getwn(){ //预处理原根的幂和逆元

int x=qpow(3,p-2);

for(int i=0;i<20;++i){

wn[i]=qpow(3,(p-1)/(1<

inv[i]=qpow(x,(p-1)/(1<

}

}

void ntt(int *y,bool f){ rev(y); //翻转代码和fft无异

for(int m=2,id=1;m<=n;m<<=1,++id){ //id用来记录转到第几下了

for(int k=0;k

int w=1,wm=f?wn[id]:inv[id]; //如果是dft就用幂, idft就用幂的逆元

for(int j=0;j

//这里跟fft一样, 不过要对p取模

int u=y[k+j]%p,t=1ll*w*y[k+j+m/2]%p;

y[k+j]=u+t; if(y[k+j]>p) y[k+j]-=p;

y[k+j+m/2]=u-t; if(y[k+j+m/2]<0) y[k+j+m/2]+=p;

w=1ll*w*wm%p;

}

}

}

if(!f){

int x=qpow(n,p-2);

for(int i=0;i

y[i]=1ll*y[i]*x%p;

}

}

好像差不多呢~ 不过这样就要求我们找一个原根好求的数. 比如著名的uoj数: 998244353 还有1004535809和469762049等, 这三个数原根都是3~

好像因为当时一看到模数不是1e9+7一般就会想到ntt, vfk为了防止这一点, 模数统一采用998244353, 现在看看收效不错.

不过 有些丧心病狂的人就是要用1e9+7作为ntt的模数, 甚至还出现了可以不模质数的情况!

那我们怎么解决任意模数ntt呢? 我们可以采用拆系数ntt或者三模数ntt. 这里介绍一下三模数ntt.

对于一般的数据范围, \(n\leq10^5, a_i\leq10^9\), 这样可能会到10^5*10^{9^2}=10^{23}级别.

所以我们可以找三个乘积\(>10^{23}\)的ntt-friendly的数, 然后分别ntt再想办法合并.

我们假如答案是ans, 那我们做三次ntt后就能得到如下三个柿子.

\[

\left\{\begin{matrix}

ans\equiv a_1(\mod m_1)\\

ans\equiv a_2(\mod m_2)\\

ans\equiv a_3(\mod m_3)

\end{matrix}\right.

\]

我们把前两个柿子通过中国剩余定理合并, 就可以得到

\[

\left\{\begin{matrix}

ans\equiv A(\mod M)\\

ans\equiv a_3(\mod m_3)

\end{matrix}\right.

\]

其中, \(M=m_1*m_2\)

这样我们设\(ans=kM+A\),

\[

kM+A\equiv a_3(\mod m_3) \k=(a_3-A)*M^{-1} (\mod m_3)

\]

这样我们求出\(k\)然后代回到\(ans=kM+A\)就可以求对任意模数取模的结果了.

中国剩余定理合并的时候直接乘是可以爆long long的, 所以我们要用到\(O(1)\)快速乘~

下面上一波代码: luogu4245 【模板】MTT

哎呀觉得自己码风有点丑啊qwq

#include

#include

#include

typedef long long LL;

const int N=600020,p0=469762049,p1=998244353,p2=1004535809;

const LL M=1ll*p0*p1;

int wn[20],nw[20],rev[N],n,lg,p;

int qpow(int a,int b,int p,int s=1){

for(;b;b>>=1,a=1ll*a*a%p)

if(b&1) s=1ll*s*a%p;

return s;

}

LL mul(LL a,LL b,LL p){ a%=p; b%=p;

return (a*b-(LL)((long double)a*b/p)*p+p)%p;

}

void calcw(int p){

int x=qpow(3,p-2,p);

for(int i=0;i<20;++i){

wn[i]=qpow(3,(p-1)/(1<

nw[i]=qpow(x,(p-1)/(1<

}

}

void init(){

for(int i=0;i

rev[i]=(rev[i>>1]>>1)|((i&1)<

}

void ntt(int *y,bool f,int p){ calcw(p);

for(int i=0;i

for(int m=2,id=1;m<=n;m<<=1,++id){

for(int k=0;k

int w=1,wm=f?wn[id]:nw[id];

for(int j=0;j>1;++j){

int &a=y[k+j]; int &b=y[k+j+m/2];

int u=a%p,t=1ll*w*b%p;

a=u+t; if(a>p) a-=p;

b=u-t; if(b<0) b+=p;

w=1ll*w*wm%p;

}

}

} int x=qpow(n,p-2,p);

if(!f) for(int i=0;i

}

char c1[N],c2[N]; int a[N],b[N],c[N],d[N],ans[3][N];

int main(){

int l1,l2; scanf("%d%d%d",&l1,&l2,&p);

for(int i=0;i<=l1;++i) scanf("%d",&a[i]),a[i]%=p;

for(int i=0;i<=l2;++i) scanf("%d",&b[i]),b[i]%=p;

for(n=1;n

std::copy(a,a+n,c); std::copy(b,b+n,d);

ntt(c,1,p0); ntt(d,1,p0);

for(int i=0;i

std::copy(a,a+n,c); std::copy(b,b+n,d);

ntt(c,1,p1); ntt(d,1,p1);

for(int i=0;i

std::copy(a,a+n,c); std::copy(b,b+n,d);

ntt(c,1,p2); ntt(d,1,p2);

for(int i=0;i

ntt(ans[0],0,p0); ntt(ans[1],0,p1); ntt(ans[2],0,p2);

for(int i=0;i

LL A=mul(1ll*ans[0][i]*p1%M,qpow(p1%p0,p0-2,p0),M)

+mul(1ll*ans[1][i]*p0%M,qpow(p0%p1,p1-2,p1),M);

if(A>M) A-=M;

LL k=((ans[2][i]-A)%p2+p2)%p2*qpow(M%p2,p2-2,p2)%p2;

a[i]=1ll*(k%p)*(M%p)%p+A%p;

if(a[i]>p) a[i]-=p;

}

for(int i=0;i<=l1+l2;++i) printf("%d",a[i]);

}

原文：https://www.cnblogs.com/enzymii/p/8925992.html