模板 P4245题解

Prean · 2022-06-29 19:16:49 · 题解

这是一个 5 次 FFT 的做法，但是可以把 4 次 FFT 吊起来锤。

我给同学讲了这种做法，他说这个做法吊打别的 5 次 FFT，建议我发篇题解，那就发吧。。。

这种 5 次 FFT 的做法不需要根据系数的关系解方程什么的，只需要设复多项式然后大力卷就行了。。。重要的是 4 次 FFT 不能用转置但是这种 5 次的可以。

假设我们要计算 F(x)\times G(x)。

设 F(x)=A(x)+c\times B(x),G(x)=C(x)+c\times D(x)，那么我们要求的就是：

对比一下：

设 T(x)=C(x)+i\times D(x)，那么我们可以在 A(x)T(x) 和 B(x)T(x) 中找到我们需要的所有项。

需要对 A(x),B(x),T(x) 做 DFT，A(x)T(x) 和 B(x)T(x) 做 IDFT，跑起来还是相当快的。

并且我使用这种 MTT 跑长度为 2^{21} 的卷积时仍然没掉精度。

#include<cstdio>
#include<cctype>
#include<cmath>
#define IMP(lim,act) for(int qwq=(lim),i=0;i^qwq;++i)act
typedef double db;
const int M=1<<18|5;
const db Pi=acos(-1);
int n,m,P,F[M],G[M],H[M];
struct Barrett{
    typedef unsigned long long ull;
    typedef __uint128_t LL;
    ull B,m;
    Barrett(const ull&m=2):m(m),B((LL(1)<<64)/m){}
    friend inline ull operator%(const ull&a,const Barrett&mod){
        ull r=a-mod.m*(LL(mod.B)*a>>64);return r>=mod.m?r-mod.m:r;
    }
}mod;
struct complex{
    db x,y;
    complex(const db&x=0,const db&y=0):x(x),y(y){}
    inline complex operator+(const complex&it)const{
        return complex(x+it.x,y+it.y);
    }
    inline complex operator-(const complex&it)const{
        return complex(x-it.x,y-it.y);
    }
    inline complex operator*(const complex&it)const{
        return complex(x*it.x-y*it.y,x*it.y+y*it.x);
    }
}buf[M<<1],*w[20];
inline int Getlen(const int&n){
    int len(0);while((1<<len)<n)++len;return len;
}
inline void swap(complex&a,complex&b){
    complex c=a;a=b;b=c;
}
inline int Get(const db&x){
    return((long long)(x+.5))%mod;
}
inline void init(const int&n){
    const int&m=Getlen(n);complex*now=buf;w[m]=now;now+=1<<m;
    IMP(1<<m,w[m][i]=complex(std::cos(i*Pi/(1<<m)),std::sin(i*Pi/(1<<m))));
    for(int k=m-1;k>=0&&(w[k]=now,now+=1<<k);--k)IMP(1<<k,w[k][i]=w[k+1][i<<1]);
}
inline void DFT(complex*f,const int&M){
    const int&n=1<<M;
    for(int len=n>>1,d=M-1;d>=0;--d,len>>=1)for(int k=0;k^n;k+=len<<1){
        complex*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=*R)),*L++=(x+y),*R++=*W++*(x-y);
    }
}
inline void IDFT(complex*f,const int&M){
    const int&n=1<<M;
    for(int len=1,d=0;d^M;++d,len<<=1)for(int k=0;k^n;k+=len<<1){
        complex*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=*W++**R)),*L++=(x+y),*R++=(x-y);
    }
    IMP(n,(f[i].x/=n,f[i].y/=n));for(int i=1;(i<<1)<n;++i)swap(f[i],f[n-i]);
}
inline void MTT(int*f,int*g,int*h,const int&n,const int&m,const int&LEN){
    static complex Q[M],P[M],T[M];const int&len=Getlen(n+m-1);
    IMP(n,(Q[i].x=f[i]&32767,P[i].x=f[i]>>15));IMP(m,T[i]=complex(g[i]&32767,g[i]>>15));
    DFT(Q,len);DFT(P,len);DFT(T,len);IMP(1<<len,(Q[i]=Q[i]*T[i],P[i]=P[i]*T[i]));IDFT(Q,len);IDFT(P,len);
    IMP(LEN,h[i]=(Get(Q[i].x)+(1ll*(Get(Q[i].y)+Get(P[i].x))<<15)+(1ll*Get(P[i].y)<<30))%mod);
    IMP(1<<len,Q[i]=P[i]=T[i]=complex());
}
inline int read(){
    int n(0);char s;while(!isdigit(s=getchar()));while(n=n*10+(s&15),isdigit(s=getchar()));return n;
}
inline void write(int n){
    static char s[15];int top(0);while(s[++top]=n%10^48,n/=10);while(putchar(s[top]),--top);putchar(' ');
}
signed main(){
    n=read()+1;m=read()+1;mod=Barrett(P=read());init(n+m-1);
    IMP(n,F[i]=read());IMP(m,G[i]=read());MTT(F,G,H,n,m,n+m-1);IMP(n+m-1,write(H[i]));
}