【题解】任意模数 NTT MTT luoguP4245 —Qiuly

$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;
}