$NTT$ 兹磁值域更大的多项式相乘,但是不兹磁任意模数。

$FFT$ 兹磁任意模数,但是不自此值域很大的多项式相乘。

现在要求做一遍值域很大的多项式相乘,并且要求任意模数,怎么办?

有两种解决方法:

  • 将模数用扩展中国剩余定理拆开做 $NTT$ ,然后用扩展中国剩余定理合并答案 (不费)

  • 将多项式的值拆开,依次做 $FFT$ 最后将答案加起来。($MTT$)

那么我们来讨论一下 $MTT$ ,现在要求 $F(x)\cdot G(x)$,由于值域很大不方便跑 $FFT$ ,于是我们将 Ta 们拆开。

设一个常数 $M$ ,对于 $F(x)$ ,我们设 $A1(x)$ 为 $F(x)/M$ ,$B1(x)$ 为 $F(x)\% M$ ,那么显然有:

$$F(x)=A1(x)M+B1(x)$$

对于 $G(x)$ 同理,我们将其拆成 $A2$ 和 $B2$ ,那么最终可以得到:

$$F(x)\cdot G(x) =(A1(x)M+B1(x))\cdot (A2(x)M+B2(x))$$

$$=A1(x)A2(x)M^2+A1(x)B2(x)M+A2(x)B1(x)M+B1(x)B2(x)$$

$A1(x)A2(x)$ 的值域就不大了,可以直接跑 $FFT$ ,其余的一样,最后一共是跑 $4$ 次 $FFT$ ,当然也得跑 $4$ 次 $DFT$ 。

另外 $FFT$ 的精度也要提高,最好开 $long \ double$ ,然后以前我的 $\pi$ 是这样写的:

#define PI 3.1415926535898

这个精度太低了!跑到后面会 WA ,得换个精度更高的:

const long double PI=acos(-1);

这下就好了,然后还有就是 $FFT$ 中的 $sin,cos$ ,貌似 $std$ 中的 $sin,cos$ 精度更高一些。

另外 $M$ 的值最好取 $2^{15}$ ,因为这种方法一共要跑 $8$ 次 $FFT$ ,速度比较慢,如果被卡了的话可以尝试预处理 $FFT$ 中的单位根。

Code:

#include<cmath>
#include<cstdio>
#include<iostream>
#include<algorithm>

const int N=4e5+2;
const int BLOCK=32768;
const long double PI=acos(-1);

typedef long long ll;

int filp[N],Ans[N],limit=1,cnt=0,n,m,p;
struct complex{complex(long double a=0,long double b=0){x=a,y=b;}long double x,y;};
complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
complex A1[N],B1[N],A2[N],B2[N],X[N];

inline void FFT(complex *f,short inv){
    for(int i=0;i<limit;++i)if(i<filp[i])std::swap(f[i],f[filp[i]]);
    for(int p=2;p<=limit;p<<=1){
        int len=p>>1;
        complex tmp=complex(std::cos(PI/len),inv*std::sin(PI/len));
        for(int k=0;k<limit;k+=p){
            complex buf=complex(1,0);
            for(int l=k;l<k+len;++l,buf=buf*tmp){
                complex t=f[l+len]*buf;
                f[l+len]=f[l]-t,f[l]=f[l]+t;
            }           
        }
    }return;
}

inline void solve(complex *a,complex *b,int res){
    for(int i=0;i<limit;++i)X[i]=a[i]*b[i];
    FFT(X,-1);
    for(int i=0;i<=n+m;++i)(Ans[i]+=(ll)(X[i].x/limit+0.5)%p*res%p)%=p;
}
inline void MTT(complex *a,complex *b,complex *c,complex *d){
    FFT(a,1),FFT(c,1),FFT(b,1),FFT(d,1);
    solve(a,c,BLOCK*BLOCK%p);solve(a,d,BLOCK%p);
    solve(c,b,BLOCK%p);solve(b,d,1);
}

int main(){
    scanf("%d%d%d",&n,&m,&p);
    for(int i=0,x;i<=n;++i)
        scanf("%d",&x),A1[i]=x/BLOCK,B1[i]=x%BLOCK;
    for(int i=0,x;i<=m;++i)
        scanf("%d",&x),A2[i]=x/BLOCK,B2[i]=x%BLOCK;
    for(limit=1;limit<=n+m;limit<<=1)++cnt;--cnt;
    for(int i=0;i<limit;++i)filp[i]=(filp[i>>1]>>1)|((i&1)<<cnt);
    MTT(A1,B1,A2,B2);
    for(int i=0;i<=n+m;++i)
        printf("%d ",Ans[i]);
    printf("\n");
    return 0;
}
分类: 文章

Qiuly

井戸の底の空にはまだかすかな希望の光がある……

2 条评论

Qiuly · 2019年3月20日 3:48 下午

嗯,谢谢啦

victor · 2019年3月19日 8:29 下午

推荐看看毛啸 16 年的论文,可以把 8 次 DFT 优化成 4 次,常数小一半

发表评论

电子邮件地址不会被公开。 必填项已用*标注