BZOJ1919: [Ctsc2010]性能优化（FFT循环卷积）

传送门

题意：
求给出两个长度为n的整数序列$a[0..n−1]$,$b[0..n−1]$和非负整数$C$。
对于两个长度为$n$的整数序列，定义∗运算，结果为一个长度为n的整数序列，例如$f∗g=h$，则有$h[k]=∑_{i+j≡k(mod \,n)}f[i]⋅g[j]$。
求$a∗b∗b∗⋯∗b$每一位模$(n+1)$的值，其中有$C$个$∗$运算，$(n+1)$是质数，$n$的质因数大小均不超过$10$。
$n≤5⋅10^5,a[i],b[i],C≤10^9$

题解：
循环卷积裸题。。
首先卷积满足结合律，后面的部分可以求出点值之后快速幂。

考虑FFT如何求出循环的卷积：
由$Cooley-Tukey$算法中的IDFT可知，如果我们知道一个函数的n个点值则可以通过逆变换求出原函数，现在假设已经求出两个函数的n个点值。按照同样的方法，有：

$$A(\omega _n^k)=\sum_{i=0}^{n-1}a_iw_n^{ik},B(\omega _n^k)=\sum_{i=0}^{n-1}b_iw_n^{ik}$$

相乘可得：

$$C(\omega^k)=\sum_{i=0}^{n-1}a_iw_n^{ik}\sum_{j=0}^{n-1}b_jw_n^{jk}$$

因为有$\omega_n^{ik}·\omega_n^{jk}=\omega_n^{k·(i+j)}=\omega_n^{k·((i+j)mod\,n)}$

那么

$$C(w_n^k)=\sum_{(i+j\,mod \,n)=l}A_lw_n^{kl}$$

其中

$$A_l={a_i·b_j(i+j \mod n =l)}$$

发现$C$就是$DFT$后的标准形式，也就是说直接对这个点值进行$IDFT$就可以得到原函数，直接求点值就好了。

直接求点值？？$n$可能不是$2^k$的形式。

其实只需要对$n$进行质因数分解，每一层每一个数暴力选取上一层的结果。
因为$n$的最大质因数不超过7，那么每一层的每一个数选择的数不超过7个，复杂度是$O(7·n\log n)$。
推导过程可以参考:http://blog.csdn.net/skywalkert/article/details/51737272

#include<bits/stdc++.h>
using namespace std;
const int Maxn=5e5+50;
inline int read(){
    char ch=getchar();int i=0,f=1;
    while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
    while(isdigit(ch)){i=(i<<1)+(i<<3)+ch-'0';ch=getchar();}
    return i*f;
}
inline void W(int x){
    static int buf[50];
    if(!x){putchar('0');return;}
    while(x)buf[++buf[0]]=x%10,x/=10;
    while(buf[0])putchar(buf[buf[0]--]+'0');
}
int n,C,G,pw[Maxn],a[Maxn],b[Maxn],tp[Maxn],mod,pr[Maxn],tot,pos[Maxn];
inline int power(int a,int b){
    int res=1;
    for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)res=1ll*res*a%mod;
    return res; 
}
inline bool check(int x,int t){
    for(int i=1;i<=tot;i++)if(power(x,t/pr[i])==1)return true;
    return false;
}
inline void findori(int phi){
    for(int i=2;i<=10;i++)
        for(;!(phi%i);pr[++tot]=i,phi/=i);
    G=2;
    while(check(G,n))++G;
}
inline int getpos(int x,int s,int now,int al){
    if(now==tot+1)return s;
    int bl=al/pr[now],rs=x%pr[now];
    return getpos((x-rs)/pr[now],s+bl*rs,now+1,bl);
}
inline void dft(int *A){
    memcpy(tp,A,sizeof(int)*n);
    for(int i=0;i<n;i++)tp[pos[i]]=A[i];
    for(int i=0;i<n;i++)A[i]=tp[i];
    for(int bl=1,pos=tot;pos>=1;bl*=pr[pos],--pos){
        int bl_len=bl*pr[pos],ct=n/bl_len;
        for(int bg=0;bg<n;bg+=bl_len)
            for(int p=0;p<bl_len;p+=bl){
                for(int q=0;q<bl;++q){
                    int s=0,o=(p+q)*ct;
                    for(int r=0;r<pr[pos];++r)
                        (s+=1ll*pw[1ll*o*r%n]*A[bg+r*bl+q]%mod)%=mod; 
                    tp[bg+p+q]=s;
                }
            }
        for(int i=0;i<n;i++)A[i]=tp[i];
    }
} 
int main(){
    n=read(),C=read();mod=n+1;findori(n);C=(C-1)%n+1; 
    for(int i=0;i<n;i++)a[i]=read();
    for(int i=0;i<n;i++)b[i]=read();
    pw[0]=1;for(int i=1;i<n;i++)pw[i]=1ll*pw[i-1]*G%mod;
    for(int i=1;i<n;i++)pos[i]=getpos(i,0,1,n);
    dft(a),dft(b);
    for(int i=0;i<n;i++)a[i]=1ll*a[i]*power(b[i],C)%mod;
    dft(a);reverse(a+1,a+n);
    for(int i=0;i<n;i++)a[i]=1ll*a[i]*power(n,n-1)%mod;
    for(int i=0;i<n;i++)W(a[i]),putchar('\n');
}