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取模后的结果,也没有办法知道期望情况下,所有小朋友欢乐程度的乘积。但是,当呆呆想到这一点的时候,已经彻底绝望了。
对于100%的数据,M<=10000,P<=255,A<=108,O<=4,S<=300,U<=100。
Solution
设f[i,j]表示前i个人分j颗糖的答案,那么
f[i,j]=∑j−1k=1f[i−1,j−k]∗g[k]
f
[
i
,
j
]
=
∑
k
=
1
j
−
1
f
[
i
−
1
,
j
−
k
]
∗
g
[
k
]
这是个卷积,其中
g[x]=Ox2+Sx+U
g
[
x
]
=
O
x
2
+
S
x
+
U
看起来好像快速幂一下答案就出来了,但是题目要求恰好m个的答案,那么再设一个
p[n,m]=∑ni=1f[i][m]
p
[
n
,
m
]
=
∑
i
=
1
n
f
[
i
]
[
m
]
p[n,m]=∑ni=1f[i][m]=p[n2][m]+∑n2i=1f[i+n2][m]
p
[
n
,
m
]
=
∑
i
=
1
n
f
[
i
]
[
m
]
=
p
[
n
2
]
[
m
]
+
∑
i
=
1
n
2
f
[
i
+
n
2
]
[
m
]
=∑ni=1p[n2][m]+∑n2i=1∑m−1j=1f[i][j]∗f[n2][m−j]
=
∑
i
=
1
n
p
[
n
2
]
[
m
]
+
∑
i
=
1
n
2
∑
j
=
1
m
−
1
f
[
i
]
[
j
]
∗
f
[
n
2
]
[
m
−
j
]
=∑ni=1p[n2][m]+∑m−1j=1f[n2][m−j]∗∑n2i=1f[i][j]
=
∑
i
=
1
n
p
[
n
2
]
[
m
]
+
∑
j
=
1
m
−
1
f
[
n
2
]
[
m
−
j
]
∗
∑
i
=
1
n
2
f
[
i
]
[
j
]
=∑ni=1p[n2][m]+∑m−1j=1f[n2][m−j]∗p[n2][j]
=
∑
i
=
1
n
p
[
n
2
]
[
m
]
+
∑
j
=
1
m
−
1
f
[
n
2
]
[
m
−
j
]
∗
p
[
n
2
]
[
j
]
可以发现后面那部分又是一个卷积,那么递归求就ok
终于打了一波蝴蝶变换,感觉还行
Code
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <complex>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
typedef std:: complex <double> com;
const double pi=acos(-1);
const double eps=0.1;
const int N=50005;
com a[N],b[N],c[N];
int g[N],f[N],p[N];
int rev[N],len,lg,MOD;
int A,B,C,n,m;
void FFT(com *a,int f) {
for (int i=0;i<len;i++) if (i<rev[i]) std:: swap(a[i],a[rev[i]]);
for (int i=1;i<len;i<<=1) {
com wn(cos(pi/(double)i),(double)f*sin(pi/(double)i));
for (int j=0;j<len;j+=i*2) {
com w(1,0);
for (int k=0;k<i;k++) {
com u=a[j+k],v=w*a[j+k+i];
a[j+k]=u+v; a[j+k+i]=u-v;
w*=wn;
}
}
}
if (f==-1) for (int i=0;i<len;i++) a[i]/=len;
}
void solve1() {
rep(i,0,len-1) a[i]=f[i];
rep(i,0,len-1) b[i]=g[i];
FFT(a,1); FFT(b,1);
rep(i,0,len-1) c[i]=a[i]*b[i];
FFT(c,-1);
rep(i,1,m) f[i]=(int)(c[i].real()+eps)%MOD,p[i]=(p[i]+f[i])%MOD;
}
void solve2() {
rep(i,0,len-1) a[i]=f[i];
rep(i,0,len-1) b[i]=p[i];
FFT(a,1); FFT(b,1);
rep(i,0,len-1) c[i]=a[i]*b[i];
FFT(c,-1);
rep(i,1,m) p[i]=(p[i]+(int)(c[i].real()+eps)%MOD)%MOD;
rep(i,0,len-1) c[i]=a[i]*a[i];
FFT(c,-1);
rep(i,1,m) f[i]=(int)(c[i].real()+eps)%MOD;
}
int main(void) {
freopen("data.in","r",stdin);
freopen("myp.out","w",stdout);
scanf("%d%d%d",&m,&MOD,&n);
scanf("%d%d%d",&A,&B,&C);
for (len=1,lg=0;len<=m*2;len*=2,lg++);
for (int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
rep(i,1,m) p[i]=g[i]=f[i]=(A%MOD*i%MOD*i%MOD+B*i%MOD+C)%MOD;
int now=0,w=1; n=tmp;
while (w*2<=n) w*=2,now++;
while (now--) {
solve2();
if (n&(1<<now)) solve1();
}
printf("%d\n", p[m]);
return 0;
}