bzoj4332 JSOI2012 分零食

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]=\sum_{k=1}^{j-1}f[i-1,j-k]*g[k]$
这是个卷积，其中$g[x]=Ox^2+Sx+U$
看起来好像快速幂一下答案就出来了，但是题目要求恰好m个的答案，那么再设一个$p[n,m]=\sum_{i=1}^{n}f[i][m]$
$p[n,m]=\sum_{i=1}^{n}f[i][m]=p[\frac{n}{2}][m]+\sum_{i=1}^{\frac{n}{2}}f[i+\frac{n}{2}][m]$
$\ \ \ \ \ \ \ \ \ \ \ \ =\sum_{i=1}^{n}p[\frac{n}{2}][m]+\sum_{i=1}^{\frac{n}{2}}\sum_{j=1}^{m-1}f[i][j]*f[\frac{n}{2}][m-j]$
$\ \ \ \ \ \ \ \ \ \ \ \ =\sum_{i=1}^{n}p[\frac{n}{2}][m]+\sum_{j=1}^{m-1}f[\frac{n}{2}][m-j]*\sum_{i=1}^{\frac{n}{2}}f[i][j]$
$\ \ \ \ \ \ \ \ \ \ \ \ =\sum_{i=1}^{n}p[\frac{n}{2}][m]+\sum_{j=1}^{m-1}f[\frac{n}{2}][m-j]*p[\frac{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;
}