循环卷积快速幂
两个注意点:
n+1不是2^k*P+1形式,任意模数又太慢?n=2^k1*3^k2*5^k3*7^k4
多路分治!深刻理解FFT运算本质:分治,推式子得到从下往上的迭代公式
最后求的是w_n^i的点值
快速幂:
循环卷积快速幂比较特殊,就是G*F,>=n的项的系数加到-n位置上
所以,由于w(n,p+n)=w(n,p),点值相乘直接得到G*F的点值表达
F,B点值,快速幂相乘即可
不用把n扩充系数等。
// luogu-judger-enable-o2 // luogu-judger-enable-o2 #include<bits/stdc++.h> #define reg register int #define il inline #define fi first #define se second #define mk(a,b) make_pair(a,b) #define numb (ch^‘0‘) using namespace std; typedef long long ll; template<class T>il void rd(T &x){ char ch;x=0;bool fl=false; while(!isdigit(ch=getchar()))(ch==‘-‘)&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb); (fl==true)&&(x=-x); } template<class T>il void output(T x){if(x/10)output(x/10);putchar(x%10+‘0‘);} template<class T>il void ot(T x){if(x<0) putchar(‘-‘),x=-x;output(x);putchar(‘ ‘);} template<class T>il void prt(T a[],int st,int nd){for(reg i=st;i<=nd;++i) ot(a[i]);putchar(‘\n‘);} namespace Miracle{ const int N=5e5+5; int pri[4]={2,3,5,7}; int n,mod,C; int G,GI; int ci[10]; int ad(int x,int y){ return x+y>=mod?x+y-mod:x+y; } int qm(int x,int y){ int ret=1; while(y){ if(y&1) ret=(ll)ret*x%mod; x=(ll)x*x%mod; y>>=1; } return ret; } int dep[33],cnt; void divi(){ int tmp=n; while(tmp%2==0) tmp/=2,dep[++cnt]=2,++ci[2]; while(tmp%3==0) tmp/=3,dep[++cnt]=3,++ci[3]; while(tmp%5==0) tmp/=5,dep[++cnt]=5,++ci[5]; while(tmp%7==0) tmp/=7,dep[++cnt]=7,++ci[7]; } bool che(int x){ for(reg i=0;i<4;++i){ if(ci[pri[i]]){ if(qm(x,n/pri[i])==1) return false; } } return true; } void fin(){ G=2; while(!che(G)) ++G; } int f[N]; int b[N]; int pos[N]; int getpos(int x){ // cout<<" getpos "<<x<<endl; int len=n,l=0; for(reg i=1;i<=cnt;++i){ // cout<<" i "<<i<<" x "<<x<<" ll "<<l<<endl; int be=(x-l)%dep[i],th=(x-l)/dep[i]; len=len/dep[i]; x=l+len*be+th; l=x-th; } // cout<<" bac "<<x<<endl; return x; } int g[N]; int pw[7*N][2]; int mem[10]; void FFT(int *f,int n,int c){ for(reg i=0;i<n;++i) g[i]=f[i]; for(reg i=0;i<n;++i) f[pos[i]]=g[i]; int las=1; for(reg i=cnt;i>=1;--i){ int p=dep[i]; int st=(mod-1)/(las*p); for(reg l=0;l<n;l+=las*p){ for(reg b=0;b<las;++b){ for(reg j=0;j<p;++j){ mem[j]=0; for(reg i=0;i<p;++i){ mem[j]=ad(mem[j],(ll)pw[st*(b+j*las)*i][c]*f[l+b+i*las]%mod); } } for(reg j=0;j<p;++j){ f[l+b+j*las]=mem[j]; } } } las*=p; } } int main(){ rd(n);rd(C);mod=n+1; divi();fin();GI=qm(G,mod-2); // cout<<" G "<<G<<" GI "<<GI<<endl;// pw[0][0]=pw[0][1]=1; for(reg i=1;i<=7*n;++i) pw[i][0]=(ll)pw[i-1][0]*GI%mod,pw[i][1]=(ll)pw[i-1][1]*G%mod; for(reg i=0;i<n;++i) rd(f[i]),f[i]%=mod; for(reg i=0;i<n;++i) rd(b[i]),b[i]%=mod; for(reg i=0;i<n;++i) pos[i]=getpos(i); // cout<<" pos "<<endl; // prt(pos,0,n-1); FFT(f,n,1); // cout<<" ff "<<endl; // prt(f,0,n-1); FFT(b,n,1); // cout<<" bb "<<endl; // prt(b,0,n-1); for(reg i=0;i<n;++i) f[i]=(ll)f[i]*qm(b[i],C)%mod; FFT(f,n,0); for(reg i=0;i<n;++i){ f[i]=(ll)f[i]*qm(n,mod-2)%mod; printf("%d\n",f[i]); } return 0; } } signed main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2019/3/19 18:42:47 */
算是对FFT本质的理解吧。
原文:https://www.cnblogs.com/Miracevin/p/10560913.html