题解 P4191 【[CTSC2010]性能优化】
憨豆Beng
2018-08-29 20:25:42
题意：
求给出两个长度为n的整数序列a[0..n−1],b[0..n−1]和非负整数C。对于两个长度为n的整数序列，定义∗运算，结果为一个长度为n的整数序列，例如f∗g=h，则有h[k]=∑i+j≡k(modn)f[i]⋅g[j]。
求a∗b∗b∗⋯∗b每一位模(n+1)的值，其中有C个∗运算，(n+1)是质数，n的质因数大小均不超过10。
n≤5⋅105,a[i],b[i],C≤109
题解：
循环卷积裸题。首先卷积满足结合律，后面的部分可以求出点值之后快速幂。
考虑FFT如何求出循环的卷积：
由Cooley−Tukey算法中的IDFT可知，如果我们知道一个函数的n个点值则可以通过逆变换求出原函数，现在假设已经求出两个函数的n个点值。按照同样的方法，有：
A(ωkn)=∑i=0n−1aiwikn,B(ωkn)=∑i=0n−1biwikn
相乘可得：
C(ωk)=∑i=0n−1aiwikn∑j=0n−1bjwjkn
因为有ωikn⋅ωjkn=ωk⋅(i+j)n=ωk⋅((i+j)modn)n
那么
C(wkn)=∑(i+jmodn)=lAlwkln
其中
Al=ai⋅bj(i+jmodn=l)
发现C就是DFT后的标准形式，也就是说直接对这个点值进行IDFT就可以得到原函数，直接求点值就好了。
直接求点值？？n可能不是2k的形式。
其实只需要对n进行质因数分解，每一层每一个数暴力选取上一层的结果。
因为n的最大质因数不超过7，那么每一层的每一个数选择的数不超过7个，复杂度是O(7⋅nlogn)。
推导过程可以参考:[传送门](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');
}```