Luogu P4191 [CTSC2010]性能优化

Link
我们知道做长度为(n)的循环卷积其实就是做长度为(n)的DFT，这可以用单位根反演证明。
注意到这里(n=2^a3^b5^c7^d)，所以并不能够像平时做长度为(2^k)的DFT那样直接分半，而是分为(n)的最小质因子块。
剩下的就跟平时的DFT差不多了，处理rev数组可能会比较麻烦。

#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
const int N=500007;
int n,m,k,P,p[20],g[N],f[N],h[N];char ibuf[(1<<24)+1],*iS=ibuf,obuf[(1<<23)+1],*oS=obuf;
int read(){int x=0;while(isspace(*iS))++iS;while(isdigit(*iS))(x*=10)+=*iS++&15;return x;}
void output(int x){if(!x)return ;output(x/10),*oS++=x%10+'0';}
void write(int x){x? output(x),0:*oS++='0',*oS++='
';}
void mod(int&x){x+=x>>31&P;}
int mul(int a,int b){return 1ll*a*b%P;}
int pow(int a,int k){int r=1;for(;k;k>>=1,a=mul(a,a))if(k&1)r=mul(a,r);return r;}
int nroot(int x){for(int i=1;i<=m;++i)if(pow(x,n/p[i])==1)return 1;return 0;}
int getroot(int x)
{
    int g=2;
    for(int i=2;i<=7;++i) for(;!(x%i);p[++m]=i,x/=i);
    for(g=2;nroot(g);++g);
    return g;
}
void reverse(int*a)
{
    static int t[N];
    for(int i=m,len=n;i;memcpy(a,t,n<<2),len/=p[i--]) for(int id=0,j=0;j<n;j+=len) for(int k=0;k<p[i];++k) for(int l=0;l<len;l+=p[i]) t[id++]=a[k+j+l];
}
void NTT(int*a)
{
    static int t[N];
    reverse(a);
    for(int i=1,len=1;i<=m;memcpy(a,t,n<<2),++i)
    {
	memset(t,0,n<<2);
        for(int mid=len,wi=g[n/(len*=p[i])],j=0;j<n;j+=len) for(int k=0,wk=1;k<len;wk=mul(wk,wi),++k) for(int l=k%mid,w=1;l<len;w=mul(w,wk),l+=mid) mod(t[j+k]+=mul(w,a[j+l])-P);
    }
}
int main()
{
    fread(ibuf,1,1<<24,stdin);
    n=read(),P=n+1,k=read(),g[0]=1,g[1]=getroot(n);
    for(int i=2;i<n;++i) g[i]=mul(g[i-1],g[1]);
    for(int i=0;i<n;++i) f[i]=read();
    for(int i=0;i<n;++i) h[i]=read();
    NTT(f),NTT(h);
    for(int i=0;i<n;++i) f[i]=mul(f[i],pow(h[i],k));
    NTT(f),std::reverse(f+1,f+n);
    for(int i=0;i<n;++i) write(mul(f[i],n));
    fwrite(obuf,1,oS-obuf,stdout);
}