瞎讲：任意模数MTT

转载

mb5fe94c9f04536 2020-05-19 12:10:00

文章标签 #include #define 原根取模卷积 文章分类 代码人生

瞎讲一下。

三模数\(NTT\)

大概思路就是用三个满足\(a*2^b+1\)形式的质数来做\(NTT\)，

然后用数论方法搞出它的具体值（当长度为\(10^5\)级别时，卷积之后数字最多为\(10^{23}\)，所有同余的数中只有最小的那个在范围内）。

一般选\(469762049,998244353,1004535809\)，原根都是\(3\)

调用\(9\)次\(DFT\)，常数极大。

二模数\(NTT\)

和上面的那个思路差不多，只不过用一个大质数和一个小质数来搞。

long long相乘取模用强制转long double来解决。

质数取\(29*2^{57}+1\)和\(998244353\)，原根都是\(3\)。

拆分\(FFT\)

考虑将每个数拆成\(aW+b\)的形式，\(W\)取\(\sqrt P\)

结果长这样：\(acW^2+(ad+bc)W+bd\)

（这里\(a,b,c,d\)都应该理解成函数）

这样一个位置上的数字最大是\(10^{14}\)级别，精度好就可以接受。

调用\(7\)次\(DFT\)。

拆分\(FFT\)优化

分别计算：

\[(a+bi)(c+di)=(ac-bd)+(ad+bc)i \\(a-bi)(c+di)=(ac+bd)+(ad-bc)i\]

将实部和虚部拆开来，就可以分别求出所需的\(ac\)、\(bd\)、\((ad+bc)\)。

于是只需求\((a+bi)\)和\((a-bi)\)和\((c+di)\)的\(DFT\)，以及两个乘积要做\(IDFT\)。

这样就可以优化到\(5\)次\(DFT\)。

利用FFT三次变两次中提到的性质（也就是共轭的两个多项式，求出其中一个的\(DFT\)之后可以\(O(n)\)地求出另一个的\(DFT\)）。

于是\((a-bi)\)的\(DFT\)可以通过\((a+bi)\)的\(DFT\)求。

所以只需要用\(4\)次\(DFT\)。

贴个代码

using namespace std;
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cassert>
#define N 262144
#define db double
#define ll long long
const db PI=acos(-1);
struct com{db a,b;};
inline com operator+(const com &x,const com &y){return {x.a+y.a,x.b+y.b};}
inline com operator-(const com &x,const com &y){return {x.a-y.a,x.b-y.b};}
inline com operator*(const com &x,const com &y){return {x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a};} 
inline com operator*(const com &x,db y){return {x.a*y,x.b*y};}
int re[N];
int n,m,nN,p,W;
int a[N],b[N],c[N];
com wnk[18][N];
void dft(com A[],int flag=1){
  for (int i=0;i<nN;++i)
    if (i<re[i])
      swap(A[i],A[re[i]]);
  int cnt=0;
  for (int i=1;i<nN;i<<=1,++cnt){
//    com wn={cos(PI/i),flag*sin(PI/i)};
    for (int j=0;j<nN;j+=i<<1){
//      com wnk={1,0};
      for (int k=j;k<j+i;++k/*,wnk=wnk*wn*/){
        com tmp=wnk[cnt][k-j];
        tmp.b*=flag;
        com x=A[k],y=tmp*A[k+i];
        A[k]=x+y;
        A[k+i]=x-y;
      }
    }
  }
  if (flag==-1){
    db inv=(db)1/nN;
    for (int i=0;i<nN;++i)
      A[i]=A[i]*inv;
  }
}
com A[N],B[N],C[N];
void multi(int a[],int b[],int n,int m){
  int bit=0;
  for (nN=1;nN<=n+m;nN<<=1,++bit);
  re[0]=0;
  for (int i=1;i<nN;++i)
    re[i]=re[i>>1]>>1|(i&1)<<bit-1;
  for (int i=0;i<bit;++i)
    for (int j=0;j<1<<i;++j)
      wnk[i][j]={cos(PI*j/(1<<i)),sin(PI*j/(1<<i))};
  /*
  for (int i=0;i<nN;++i){
    A[i]={a[i],0};
    B[i]={b[i],0};
  }
  dft(A),dft(B);
  for (int i=0;i<nN;++i)
    A[i]=A[i]*B[i];
  dft(A,-1);
  for (int i=0;i<nN;++i)
    c[i]=((ll)(A[i].a+0.5)%p+p)%p;
  */
  W=sqrt(p);
  for (int i=0;i<nN;++i)
    A[i]={a[i]/W,a[i]%W};
  dft(A);
  /*
  for (int i=0;i<nN;++i)
    B[i]={a[i]/W,-a[i]%W};
  dft(B);
  for (int i=0;i<nN;++i){
    printf("(%lf,%lf)=(%lf,%lf)\n",B[i].a,B[i].b,A[(nN-i)%nN].a,-A[(nN-i)%nN].b);
    assert(abs(B[i].a-A[(nN-i)%nN].a)<1e-8);
    assert(abs(B[i].b+A[(nN-i)%nN].b)<1e-8);
  }
  */
  B[0]={A[0].a,-A[0].b};
  for (int i=1;i<nN;++i)
    B[i]={A[nN-i].a,-A[nN-i].b};
  
  for (int i=0;i<nN;++i)
    C[i]={b[i]/W,b[i]%W};    
  dft(C);
  for (int i=0;i<nN;++i){
    A[i]=A[i]*C[i];
    B[i]=B[i]*C[i];
  }
  dft(A,-1),dft(B,-1);
  for (int i=0;i<nN;++i){
    ll Aa=round(A[i].a),Ab=round(A[i].b),Ba=round(B[i].a);
    ll x=((ll)((Aa+Ba)/2)%p+p)%p;
    ll y=((ll)((Ba-Aa)/2)%p+p)%p;
    ll z=((ll)(Ab)%p+p)%p;
    c[i]=(W*W*x+W*z+y)%p;
//    printf("(%.16lf->%lld,%.16lf->%lld) (%.16lf->%lld,/) %x=%lld y=%lld z=%lld c[i]=%lld\n",A[i].a,Aa,A[i].b,Ab,B[i].a,Ba,x,y,z,c[i]);
  }
}
int main(){
//  freopen("in.txt","r",stdin);
//  freopen("out.txt","w",stdout);
  scanf("%d%d%d",&n,&m,&p);
  for (int i=0;i<=n;++i)
    scanf("%d",&a[i]),a[i]%=p;
  for (int i=0;i<=m;++i)
    scanf("%d",&b[i]),b[i]%=p;
  multi(a,b,n,m);
  for (int i=0;i<=n+m;++i)
    printf("%d ",c[i]);
  return 0;
}