【BZOJ 4332】 4332: JSOI2012 分零食 （FFT+快速幂）
4332: JSOI2012 分零食
Time Limit: 10 Sec Memory Limit: 256 MB
Submit: 119 Solved: 66Description
这里是欢乐的进香河，这里是欢乐的幼儿园。今天是2月14日，星期二。在这个特殊的日子里，老师带着同学们欢乐地跳着，笑着。校长从幼儿园旁边的小吃店买了大量的零食决定分给同学们。听到这个消息，所有同学都安安静静地排好了队，大家都知道，校长不喜欢调皮的孩子。同学们依次排成了一列，其中有A位小朋友，有三个共同的欢乐系数O，S和U。如果有一位小朋友得到了x个糖果，那么她的欢乐程度就是f(x)=O*x2+S*x+U。现在校长开始分糖果了，一共有M个糖果。有些小朋友可能得不到糖果，对于那些得不到糖果的小朋友来说，欢乐程度就是1。如果一位小朋友得不到糖果，那么在她身后的小朋友们也都得不到糖果。（即这一列得不到糖果的小朋友一定是最后的连续若干位）所有分糖果的方案都是等概率的。现在问题是：期望情况下，所有小朋友的欢乐程度的乘积是多少？呆呆同学很快就有了一个思路，只要知道总的方案个数T和所有方案下欢乐程度乘积的总和S，就可以得到答案Ans=S/T。现在他已经求出来了T的答案，但是S怎么求呢？他就不知道了。你能告诉他么？因为答案很大，你只需要告诉他S对P取模后的结果。后记：虽然大家都知道，即便知道了T，知道了S对P取模后的结果，也没有办法知道期望情况下，所有小朋友欢乐程度的乘积。但是，当呆呆想到这一点的时候，已经彻底绝望了。Input
第一行有2个整数，分别是M和P。第二行有一个整数A，第三行有一个整数O。第四行有一个整数S，第五行有一个整数U。Output
一个整数S，因为答案可能很大，你只需要输出S 对P取模后的结果。Sample Input
4 100
4
1
0
0Sample Output
63
样例说明
函数f(x)=x^2。一共有4份零食，4位同学。如果只有第一个同学得到，欢乐程度为16，若前两位同学得到，欢乐程度的所有可能依次为9，9，16，若有三位同学得到，欢乐程度有4，4，4，最后一种情况，每一个同学都得到了零食，欢乐程度为1。相加后得到S=63。
应上传者要求，此题不公开，如有异议，请提出.HINT
对于100%的数据，M<=10000，P<=255，A<=108，O<=4，S<=300，U<=100。
Source
【分析】
O(n^2)做法：【实测可以过全部
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cstring> 4 #include<iostream> 5 #include<algorithm> 6 using namespace std; 7 int Mod; 8 int f[2][10010]; 9 10 int main() 11 { 12 int m,n,a,b,c; 13 scanf("%d%d%d%d%d%d",&m,&Mod,&n,&a,&b,&c); 14 if(n>m) n=m; 15 memset(f,0,sizeof(f)); 16 f[0][0]=1; 17 int ans=0,p=0,ss; 18 for(int j=1;j<=n;j++) 19 { 20 int nw=0,pp=p^1; 21 memset(f[pp],0,sizeof(f[pp])); 22 ss=0; 23 for(int i=1;i<=m;i++) 24 { 25 if(i>=3) ss=(ss+f[p][i-3])%Mod; 26 f[pp][i]=f[pp][i-1]; 27 if(i>=3) nw=(nw+ss*2*a)%Mod; 28 if(i>=2) nw=(nw+f[p][i-2]*(2*a+a+b))%Mod; 29 f[pp][i]=(f[pp][i]+nw+f[p][i-1]*(a+b+c))%Mod; 30 }p^=1; 31 ans=(ans+f[p][m])%Mod; 32 } 33 printf("%d\n",ans); 34 return 0; 35 }
具体解释见这里
O（n^2 logn）
【实测还不如暴力，不知道是不是我常数大。
【就是后面两个for改成卷积形式
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cstring> 4 #include<iostream> 5 #include<algorithm> 6 #include<cmath> 7 using namespace std; 8 const double pi=acos(-1); 9 #define Maxn 10010*4 10 int Mod; 11 int f[Maxn],g[Maxn]; 12 13 struct P 14 { 15 double x,y; 16 P() {x=y=0;} 17 P(double x,double y):x(x),y(y){} 18 friend P operator + (P x,P y) {return P(x.x+y.x,x.y+y.y);} 19 friend P operator - (P x,P y) {return P(x.x-y.x,x.y-y.y);} 20 friend P operator * (P x,P y) {return P(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);} 21 }a[Maxn],b[Maxn]; 22 23 int R[Maxn],nn,m; 24 void dft(P *a,int t) 25 { 26 for(int i=0;i<=nn;i++) if(i<R[i]) swap(a[i],a[R[i]]); 27 for(int i=1;i<nn;i<<=1) 28 { 29 P wn(cos(pi/i),t*sin(pi/i)); 30 for(int j=0;j<nn;j+=(i<<1)) 31 { 32 P w(1,0); 33 for(int k=0;k<i;k++,w=w*wn) 34 { 35 P x=a[j+k],y=w*a[j+k+i]; 36 a[j+k]=x+y;a[j+k+i]=x-y; 37 } 38 } 39 } 40 if(t==-1) 41 { 42 for(int i=0;i<=nn;i++) a[i].x/=nn; 43 } 44 } 45 46 void fft(int *X,int *Y) 47 { 48 for(int i=0;i<=nn;i++) 49 { 50 a[i].x=X[i];a[i].y=0; 51 b[i].x=Y[i];b[i].y=0; 52 } 53 dft(a,1);dft(b,1); 54 for(int i=0;i<=nn;i++) a[i]=a[i]*b[i]; 55 dft(a,-1); 56 for(int i=0;i<=m;i++) X[i]=(int)(a[i].x+0.5)%Mod; 57 } 58 59 int main() 60 { 61 int n,a,b,c; 62 scanf("%d%d%d%d%d%d",&m,&Mod,&n,&a,&b,&c); 63 if(n>m) n=m; 64 f[0]=g[0]=0; 65 for(int i=1;i<=m;i++) f[i]=g[i]=(a*i*i+b*i+c)%Mod; 66 int ans=0; ans+=f[m]; 67 68 nn=1;int ll=0; 69 while(nn<=2*m) nn<<=1,ll++; 70 for(int i=0;i<=nn;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(ll-1)); 71 for(int i=2;i<=n;i++) 72 { 73 fft(f,g); 74 ans=(ans+f[m])%Mod; 75 } 76 printf("%d\n",ans); 77 return 0; 78 }
O（nlogn^2）做法：
把上面的第一个循环改成快速幂
从7点搞到了现在。。
理解奥爷爷的代码好久啊。。。但是打的真心短。。。
下面那个qpow是一个矩阵乘法的快速幂！！【我傻啊看了很久才看出来。。
所以求$G^1+G^2+...G^n$
呵呵【我看这个看了好久
【呵呵
其实我觉得这个卷积有点迷
但是不管了，明天再说
关于$a^1+a^2+a^3...+a^n$，当然如果你只是要求数的话，可以直接套等比数列求和公式。但如果不是数呢？或者mod下没有逆元怎么求分母呢？
可以用 矩阵乘法快速幂 ，就类似上面的那个方法的。
打成数值的形式是这样：
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cstring> 4 #include<iostream> 5 #include<algorithm> 6 using namespace std; 7 8 int main() 9 { 10 int x,y; 11 scanf("%d%d",&x,&y); 12 int ans=0,xx=x; 13 while(y) 14 { 15 if(y&1) 16 { 17 ans=ans*x+xx*1; 18 } 19 xx=xx*x+xx*1; 20 x=x*x; 21 y>>=1; 22 } 23 printf("%d\n",ans); 24 }
上面的求卷积的幂的和，就是这样子做的。
然后是O（n logn^2）的代码
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cstring> 4 #include<iostream> 5 #include<algorithm> 6 #include<cmath> 7 using namespace std; 8 #define Maxn 10010*4 9 const double pi=acos(-1); 10 int Mod; 11 12 struct P 13 { 14 double x,y; 15 P() {x=y=0;} 16 P(double x,double y):x(x),y(y){} 17 friend P operator + (P x,P y) {return P(x.x+y.x,x.y+y.y);} 18 friend P operator - (P x,P y) {return P(x.x-y.x,x.y-y.y);} 19 friend P operator * (P x,P y) {return P(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);} 20 }a[Maxn],b[Maxn]; 21 22 int R[Maxn],nn,m; 23 void dft(P *a,int f) 24 { 25 for(int i=0;i<nn;i++) if(i<R[i]) swap(a[i],a[R[i]]); 26 for(int i=1;i<nn;i<<=1) 27 { 28 P wn(cos(pi/i),f*sin(pi/i)); 29 for(int j=0;j<nn;j+=i<<1) 30 { 31 P w(1,0); 32 for(int k=0;k<i;k++,w=w*wn) 33 { 34 P x=a[j+k],y=w*a[j+k+i]; 35 a[j+k]=x+y;a[j+k+i]=x-y; 36 } 37 } 38 } 39 if(f==-1) 40 { 41 for(int i=0;i<=nn;i++) a[i].x/=nn,a[i].y/=nn; 42 } 43 } 44 45 int A[Maxn],B[Maxn],C[Maxn],nw[Maxn]; 46 int aa,bb,cc; 47 void fft(int *A,int *B) 48 { 49 for(int i=0;i<nn;i++) 50 { 51 a[i].x=A[i];a[i].y=0; 52 b[i].x=B[i];b[i].y=0; 53 } 54 dft(a,1);dft(b,1); 55 for(int i=0;i<=nn;i++) a[i]=a[i]*b[i]; 56 dft(a,-1); 57 for(int i=1;i<=m;i++) A[i]=((int)(a[i].x+0.5)%Mod); 58 } 59 60 void add(int *A,int *B) 61 { 62 for(int i=1;i<=m;i++) A[i]=(A[i]+B[i])%Mod; 63 } 64 65 void qpow(int k) 66 { 67 for(int i=0;i<=m;i++) A[i]=0; 68 for(int i=1;i<=m;i++) C[i]=B[i]=(aa*i*i+bb*i+cc)%Mod; 69 while(k) 70 { 71 if(k&1) 72 { 73 fft(A,B); 74 add(A,C); 75 } 76 for(int i=1;i<=m;i++) nw[i]=C[i]; 77 fft(C,B); 78 add(C,nw); 79 fft(B,B); 80 k>>=1; 81 } 82 } 83 84 int main() 85 { 86 int n; 87 scanf("%d%d%d%d%d%d",&m,&Mod,&n,&aa,&bb,&cc); 88 if(n>m) n=m; 89 nn=1;int ll=0; 90 while(nn<=2*m) nn<<=1,ll++; 91 for(int i=0;i<=nn;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(ll-1)); 92 qpow(n); 93 printf("%d\n",A[m]); 94 return 0; 95 }
2017-04-14 21:46:24
