fft,ntt

在被两题卡了常数之后，花了很久优化了自己的模板

现在的一般来说任意模数求逆1s跑3e5,exp跑1e5是没啥问题的（自己电脑，可能比luogu慢一倍）

当模数是$998244353,1004535809,9985661441$的时候（这$3$个的原根都是$3$）

我们会使用$ntt$来求解

$ntt$的模板本身常数不大优化效果不明显

const int mo=;

const int G=;

IL int fsp(int x,int y)

{

    ll now=;

    while (y)

    {

        if (y&) now=now*x%mo;

        x=1ll*x*x%mo;

        y>>=;

    }

    return now;

}

IL void ntt_init()

{

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

    for (int i=;i<n;i++) r[i]=(r[i/]/)|((i&)<<(l-));

}

IL void clear()

{

    for (int i=;i<=n;i++) a[i]=b[i]=;

}

void ntt(int *a,int o)

{

    for (int i=;i<n;i++) if (i>r[i]) swap(a[i],a[r[i]]);

    for (int i=;i<n;i<<=)

    {

        int wn=fsp(G,(mo-)/(i*)); w[]=;

        rep(j,,i-) w[j]=(1ll*w[j-]*wn)%mo;

        for (int j=;j<n;j+=(i*))

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

          {

              int x=a[j+k],y=1ll*a[i+j+k]*w[k]%mo;

              a[j+k]=(x+y)%mo; a[i+j+k]=(x-y)%mo;

          }

    }

    if (o==-)

    {

      reverse(&a[],&a[n]);

      for (int i=,inv=fsp(n,mo-);i<n;i++)

        a[i]=1ll*a[i]*inv%mo;

    }

}

IL void getcj(int *A,int *B,int len)

{

    m=len*2; ntt_init();

    for (int i=0;i<len;i++) a[i]=A[i],b[i]=B[i];

    ntt(a,1); ntt(b,1);

    for(int i=0;i<n;i++) a[i]=1ll*a[i]*b[i]%mo;

    ntt(a,-1);

    for (int i=0;i<len;i++) B[i]=a[i];

    clear();

}

当模数不为这$3$个，我们就需要$mtt$来实现

而$mtt$的实现为用$mx$的方法将数的实部虚部分别放$x \& 65536,x(>>15)$

另外一个重要的地方是要预处理出$w$，我们采用指针来存，避免使用vector

代码$p$的初始值为$2*n$

所有数组大小为$4*n$

$getcj$的时候要先把数组中的负数变正

IL void clear()

{

    for (int i=;i<=n;i++) a[i].a=a[i].b=b[i].a=b[i].b=c[i].a=c[i].b=d[i].a=d[i].b=;

}

cp *w[N],tmp[N*];

int p;

IL void init()

{

    cp *now=tmp;

    for (int i=;i<=p;i<<=)

    {

        w[i]=now;

        for (int j=;j<i;j++) w[i][j]=(cp){cos(pi*j/i),sin(pi*j/i)};

        now+=i;

    }

}

IL void fft_init()

{

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

    for (int i=;i<n;i++) r[i]=(r[i/]/)|((i&)<<(l-));

}

void fft(cp *a,int o)

{

    for (int i=;i<n;i++) if (i>r[i]) swap(a[i],a[r[i]]);

    for (int i=;i<n;i<<=)

        for (int j=;j<n;j+=(i*))

        {

          cp *x1=a+j,*x2=a+i+j,*W=w[i];

          for (int k=;k<i;k++,x1++,x2++,W++)

          {

              cp x=*x1,y=(cp){(*W).a,(*W).b*o}*(*x2);

              *x1=x+y,*x2=x-y;

          }

        }

    if (o==-) for(int i=;i<n;i++) a[i].a/=n;

}

IL void getcj(int *A,int *B,int len)

{

    rep(i,,len)

    {

        A[i]=(A[i]+mo)%mo,B[i]=(B[i]+mo)%mo;

    }

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

    {

       a[i]=(cp){A[i]&,A[i]>>};

       b[i]=(cp){B[i]&,B[i]>>};

    }

    m=len*; fft_init();

    fft(a,); fft(b,);

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

    {

        int j=(n-)&(n-i);

        c[j]=(cp){0.5*(a[i].a+a[j].a),0.5*(a[i].b-a[j].b)}*b[i];

        d[j]=(cp){0.5*(a[i].b+a[j].b),0.5*(a[j].a-a[i].a)}*b[i];

    }

    fft(c,); fft(d,);

    double inv=ee/n;

    rep(i,,n) c[i].a*=inv,c[i].b*=inv;

    rep(i,,n) d[i].a*=inv,d[i].b*=inv;

    rep(i,,len)

    {

        ll a1=c[i].a+0.5,a2=c[i].b+0.5;

        ll a3=d[i].a+0.5,a4=d[i].b+0.5;

        B[i]=(a1+((a2+a3)%mo<<)+((a4%mo)<<))%mo;

    }

    clear();

}

对于其他的多项式函数

用$fft$还是$ntt$是差不多的（除了数组类型）

巴特西

fft,ntt

最新文章

热门文章