洛谷 4245 【模板】任意模数NTT——三模数NTT / 拆系数FFT

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

三模数NTT：

　　大概是用3个模数分别做一遍，用中国剩余定理合并。

　　前两个合并起来变成一个 long long 的模数，再要和第三个合并的话就爆 long long ，所以可以用一种让两个模数的乘积不出现的方法：https://blog.csdn.net/qq_35950004/article/details/79477797

　　x*m1+a1 = -y*m2 + a2 <==> x*m1+y*m2 = a2-a1 <==> x*m1 = a2-a1 (mod m2) <==> x=(a2-a1)*m1^{-1} (mod m2)

　　然后根据该博客里的证明，在mod m2意义下算出来的 x 就是真的 x 。这样的话答案就是 x*m1+a1 ，可以在快速乘的过程中对题目中给的模数取模，就不会爆 long long 啦。

　　注意输入的 a[ ] 和 b[ ] 不能 ntt( ,0, ) 之后再 ntt( ,1, ) 回来，因为值已经模了刚才那个模数了；所以要多开一些数组。

　　注意输入进 mul 里的 a 和 b 应该是正的，不然没法 b>>=1 之类的。

#include<iostream>

#include<cstdio>

#include<cstring>

#include<algorithm>

#define ll long long

using namespace std;

const int N=1e5+;

int m[]={,,};

int n0,n1,mod,len,r[N<<],a[][N<<],b[][N<<],c[][N<<];

ll M=(ll)m[]*m[],d[N<<];

int rdn()

{

  int ret=;bool fx=;char ch=getchar();

  while(ch>''||ch<''){if(ch=='-')fx=;ch=getchar();}

  while(ch>=''&&ch<='') ret=ret*+ch-'',ch=getchar();

  return fx?ret:-ret;

}

void upd(ll &x,ll md){x>=md?x-=md:;}

void upd(int &x,ll md){x>=md?x-=md:;}

ll mul(ll a,ll b,ll md)

{

  a%=md; b%=md;//

  ll ret=;while(b){if(b&1ll)ret+=a,upd(ret,md);a+=a;upd(a,md);b>>=1ll;}return ret;

}

ll pw(ll x,ll k,ll md)

{ll ret=;while(k){if(k&1ll)ret=mul(ret,x,md);x=mul(x,x,md);k>>=1ll;}return ret;}

void ntt(int *a,bool fx,int md)

{

  for(int i=;i<len;i++)

    if(i<r[i])swap(a[i],a[r[i]]);

  for(int R=;R<=len;R<<=)

    {

      int m=R>>;

      int Wn=pw(,(md-)/R,md);

      fx?Wn=pw(Wn,md-,md):;

      for(int i=;i<len;i+=R)

    for(int j=,w=;j<m;j++,w=(ll)w*Wn%md)

      {

        int tmp=(ll)w*a[i+m+j]%md;

        a[i+m+j]=a[i+j]+md-tmp; upd(a[i+m+j],md);

        a[i+j]=a[i+j]+tmp; upd(a[i+j],md);

      }

    }

  if(!fx)return; int inv=pw(len,md-,md);

  for(int i=;i<len;i++) a[i]=(ll)a[i]*inv%md;

}

int main()

{

  n0=rdn()+; n1=rdn()+; mod=rdn();

  for(int i=;i<n0;i++)a[][i]=a[][i]=a[][i]=rdn();

  for(int i=;i<n1;i++)b[][i]=b[][i]=b[][i]=rdn();

  for(len=;len<=n0+n1;len<<=);

  for(int i=;i<len;i++)r[i]=(r[i>>]>>)+((i&)?len>>:);

  for(int i=;i<;i++)//don't ntt(a,1,m[i]) for it can't return(already mod)

    {

      ntt(a[i],,m[i]); ntt(b[i],,m[i]);

      for(int j=;j<len;j++)c[i][j]=(ll)a[i][j]*b[i][j]%m[i];

      ntt(c[i],,m[i]);

    }

  ll inv=pw(m[],m[]-,m[]),t;

  for(int i=,lm=n0+n1-;i<lm;i++)

    {

      t=mul((c[][i]-c[][i])%m[]+m[],inv,m[]);

      d[i]=(mul(t,m[],M)+c[][i])%M;

    }

  inv=pw(M,m[]-,m[]);

  for(int i=,lm=n0+n1-;i<lm;i++)

    {

      t=mul((c[][i]-d[i])%m[]+m[],inv,m[]);

      d[i]=(mul(t,M,mod)+d[i])%mod;

      printf("%lld ",d[i]);

    }

  puts(""); return ;

}

拆系数FFT：

　　参考材料：https://blog.csdn.net/lvzelong2014/article/details/80156989

　　因为数值最大可能 10^(9+9+5) ，double 的精度可能很不好。所以把 a[ i ] 拆成 k*m + b ，其中 m 大约是 sqrt(mod) ；

　　这样的话2个多项式变成了4个多项式，卷积的时候只卷积 k 和 b ，不用带上 m ，数值最大就是 sqrt(mod) 级别的，精度就能行了。

　　但 ( k1 + b1 )*( k2 + b2 ) = k1*k2 + b1*k2 + k1*b2 + b1*b2 ，需要4次卷积，做8次DFT（似乎可以弄成7次），太慢了。

　　所以就像那个博客里说的那样：

　　令 \( P(x) = A(x)+i*B(x) , Q(x) = A(x)-i*B(x) \) （写的时候就是把 a[ i ] 放在 p[ i ] 的实部、把 b[ i ] 放在 p[ i ] 的虚部）

　　则 \( Q(w_{n}^{k}) = conj(P(w_{n}^{-k})) \)（这里已经代入了点值。即把 P DFT之后，可以根据它得到 Q 的点值）（证明见那个博客）

　　然后由 P(x) 和 Q(x) 的定义，得

　　\( A(w_{n}^{k}) = \frac{P(w_{n}^{k})+Q(w_{n}^{k})}{2} \)（点值）

　　\( i*B(w_{n}^{k}) = \frac{P(w_{n}^{k})-Q(w_{n}^{k})}{2} \) 即 \( B(w_{n}^{k}) = i*\frac{Q(w_{n}^{k})-P(w_{n}^{k})}{2} \)

　　把向量写成 ( 实部+ i * 虚部 ) 的样子，就能把 i 代进去，从而得到 A(x) 和 B(x) 的点值向量。这种方法利用 P(x) 和 Q(x) 的关联，只 DFT 了一下P(x)，就求出了A(x)和B(x)两个多项式的点值！

　　所以 DFT 两次，求出 k1、b1、k2、b2 的点值。相乘一番，得到 k1*k2 、b1*k2 、k1*b2 、b1*b2 的点值。

　　考虑 iDFT 回去。还是设 \( P(w_{n}^{k}) = A(w_{n}^{k}) + i*B(w_{n}^{k})\)（点值），写的时候还是把向量写开，就能把 i 代进去，得到 \( P(w_{n}^{k}) \) 的向量。

　　这时的 \( P(w_{n}^{k}) \) 是一组点值， iDFT 回去后直接是系数的 \( A(x) + i*B(x) \) ！也就是 iDFT 后 p[ i ] 的实部是 A(x) 的系数 a[ i ] ，虚部是 b[ i ] 。

　　　　这个的证明可以看那个博客，大概就是如果系数是( A 的系数 + B 的系数 )的话，点值也应该是 ( A 的点值 + B 的点值 )；倒着推一下得到点值是 ( A 的点值 + B 的点值 ) 的话，iDFT 回去的系数也应该是 ( A 的系数 + B 的系数 ) 。

　　这里的A(x)和B(x)就是 k1*k2 、b1*k2 、k1*b2 、b1*b2 中的两个。两次 iDFT 后就得到了这4个东西的系数表达！

　　然后统计答案就行了。注意乘上 m² 或者 m 。

　　这道题得开 long double 。

#include<iostream>

#include<cstdio>

#include<cstring>

#include<algorithm>

#include<cmath>

#define db long double

#define ll long long

using namespace std;

const int N=1e5+,M=(<<)+; const db pi=acos(-);

int mod,a[N],b[N],len,r[M],ans[M];

struct cpl{db x,y;}pa[M],pb[M],Ta[M],Tb[M],Tc[M],Td[M],I;

cpl operator+ (cpl a,cpl b){return (cpl){a.x+b.x,a.y+b.y};}

cpl operator- (cpl a,cpl b){return (cpl){a.x-b.x,a.y-b.y};}

cpl operator* (cpl a,cpl b){return (cpl){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}

cpl cnj(cpl a){return (cpl){a.x,-a.y};}

void upd(int &x){x>=mod?x-=mod:;}

int rdn()

{

  int ret=;bool fx=;char ch=getchar();

  while(ch>''||ch<''){if(ch=='-')fx=;ch=getchar();}

  while(ch>=''&&ch<='') ret=ret*+ch-'',ch=getchar();

  return fx?ret:-ret;

}

void fft(cpl *a,bool fx)

{

  for(int i=;i<len;i++)

    if(i<r[i])swap(a[i],a[r[i]]);

  for(int R=;R<=len;R<<=)

    {

      int m=R>>;

      cpl Wn=(cpl){ cos(pi/m),fx?-sin(pi/m):sin(pi/m) };

      for(int i=;i<len;i+=R)

    {

      cpl w=I;

      for(int j=;j<m;j++,w=w*Wn)

        {

          cpl x=a[i+j],y=w*a[i+m+j];

          a[i+j]=x+y;  a[i+m+j]=x-y;

        }

    }

    }

  if(!fx)return;

  for(int i=;i<len;i++)a[i].x/=len,a[i].y/=len;//y!! for use y

}

void solve(int n,int *a,int m,int *b)

{

  int bin=(<<)-; cpl ta,tb,tc,td;

  for(int i=;i<=n;i++)pa[i]=(cpl){a[i]>>,a[i]&bin};//pa=A+iB(a=km+b,A=k,B=b)

  for(int i=;i<=m;i++)pb[i]=(cpl){b[i]>>,b[i]&bin};//<=m!!!

  for(len=;len<=n+m;len<<=);

  for(int i=;i<len;i++)r[i]=(r[i>>]>>)+((i&)?len>>:);

  fft(pa,);  fft(pb,);

  pa[len]=pa[]; pb[len]=pb[];

  for(int i=,j=len;i<len;i++,j--)//j=-i+len(if i=0 then j=0),qa[i]=cnj(pa[j])

    {

      //ta:point value of A

      ta=(pa[i]+cnj(pa[j]))*(cpl){0.5,};//ta={(pa[i].x+pa[j].x)/2,(pa[i].y+pa[j].y)/2}

      tb=(pa[i]-cnj(pa[j]))*(cpl){,-0.5};//tb=i(Q-P)/2=i*((Qx-Px)+i(Qy-Py))/2={(Py-Qy)/2,(Qx-Px)/2}

      tc=(pb[i]+cnj(pb[j]))*(cpl){0.5,};

      td=(pb[i]-cnj(pb[j]))*(cpl){,-0.5};

      Ta[i]=ta*tc;Tb[i]=ta*td;Tc[i]=tb*tc;Td[i]=tb*td;

    }

  pa[len]=pb[len]=(cpl){,};

  for(int i=;i<len;i++)pa[i]=Ta[i]+Tb[i]*(cpl){,};//pa=Ta+i*Tb={Tax-Tby,Tay+Tbx}

  for(int i=;i<len;i++)pb[i]=Tc[i]+Td[i]*(cpl){,};

  fft(pa,);  fft(pb,);//pa.x=Ta,pa.y=Tb

  for(int i=,j=n+m;i<=j;i++)

    {

      int Da=(ll)(pa[i].x+0.5)%mod;

      int Db=(ll)(pa[i].y+0.5)%mod;

      int Dc=(ll)(pb[i].x+0.5)%mod;

      int Dd=(ll)(pb[i].y+0.5)%mod;

      ans[i]=(((ll)Da<<) + ((ll)(Db+Dc)<<) + Dd)%mod+mod; upd(ans[i]);

    }

}

int main()

{

  int n,m; I.x=;

  n=rdn();m=rdn();mod=rdn();

  for(int i=;i<=n;i++)a[i]=rdn()%mod+mod,upd(a[i]);

  for(int i=;i<=m;i++)b[i]=rdn()%mod+mod,upd(b[i]);

  solve(n,a,m,b);

  for(int i=,j=n+m;i<=j;i++)printf("%d ",ans[i]);puts("");

  return ;

}

巴特西

洛谷 4245 【模板】任意模数NTT——三模数NTT / 拆系数FFT

最新文章

热门文章