BZOJ 4332 4332: JSOI2012 分零食 （FFT+快速幂）

Posted konjak魔芋

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了BZOJ 4332 4332: JSOI2012 分零食 （FFT+快速幂）相关的知识，希望对你有一定的参考价值。

4332: JSOI2012 分零食

Time Limit: 10 Sec  Memory Limit: 256 MB
Submit: 119  Solved: 66

Description

这里是欢乐的进香河，这里是欢乐的幼儿园。 
今天是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
0

Sample 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 }
View Code

具体解释见这里

 

 

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 }
View Code

 

 

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 }
View Code

上面的求卷积的幂的和，就是这样子做的。

 

 然后是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 }
View Code

 

2017-04-14 21:46:24

 

以上是关于BZOJ 4332 4332: JSOI2012 分零食 （FFT+快速幂）的主要内容，如果未能解决你的问题，请参考以下文章

bzoj千题计划309：bzoj4332: JSOI2012 分零食（分治FFT）

BZOJ 4332 FFT+快速幂

bzoj 4332 FFT型的快速幂（需要强有力的推导公式能力）

洛谷 P4332 [SHOI2014]三叉神经树 题解

BZOJ2246 [SDOI2011]迷宫探险 记忆化搜索dp + 概率

想要唯一的记录——啥是 SQL？