BZOJ4332: JSOI2012 分零食

首先，因为糖果数M<=10000，最多分给前M个人，所以人数A>M的情况可以视为A=M来做
然后朴素的想法是一个DP，即f[i][j]表示前i个人得到j个糖果的欢乐程度乘积和，有$f[i][j]=\sum_{k=1}^{j-1}f[i-1][k]*g(j-k)$，其中$g(i)$表示一个人得到i个糖果的欢乐程度，这样直接DP是$O(n^3)$的
经观察可以发现，这是一个卷积的形式，可以用FFT优化，因为需要统计$\sum_{i=1}^nf[i][M]$，所以设转移的多项式为$g$，要求的即是$f[0]*(g^1+g^2+g^3+.....+g^n)$在T位置的值，这个东西可以分治处理，总复杂度$O(nlog^2n)$

code：

#include<set>
#include<map>
#include<deque>
#include<queue>
#include<stack>
#include<cmath>
#include<ctime>
#include<bitset>
#include<string>
#include<vector>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<complex>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;

const int maxn = 81000;
const int maxd = 18;
const double pi = acos(-1);

int T,P,n,O,S,U;

struct E
{
    double x,y;
    E(){x=y=0.0;}
    E(const double _x,const double _y){x=_x;y=_y;}
}zero,w[maxn],g[maxn],f1[maxn],f2[maxn],f[maxn],tt[maxn],t2[maxn];
int id[maxn],N,ln;
E operator +(E x,E y){return E(x.x+y.x,x.y+y.y);}
E operator -(E x,E y){return E(x.x-y.x,x.y-y.y);}
E operator *(E x,E y){return E(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}

void DFT(E *s,int sig)
{
    for(int i=0;i<N;i++) if(i<id[i]) swap(s[i],s[id[i]]);

    for(int m=2;m<=N;m<<=1)
    {
        int t=m>>1,tt=N/m;
        for(int i=0;i<t;i++)
        {
            E wn=sig==1?w[i*tt]:w[N-i*tt];
            /*E wn;
            if(sig==1) wn=E(cos(2*pi*i/m),sin(2*pi*i/m));
            else wn=E(cos(2*pi*i/m),sin(-2*pi*i/m));*/
            for(int j=i;j<N;j+=m)
            {
                E tx=s[j],ty=s[j+t]*wn;
                s[j]=tx+ty;
                s[j+t]=tx-ty;
            }
        }
    }
    if(sig==-1) for(int i=0;i<N;i++) s[i].x/=N;
}

void get_(const int u,const int k)
{
    if(u==1)
    {
        for(int i=0;i<N;i++) f[i]=tt[i]=f1[i];
        return ;
    }
    const int nex=k+1;
    get_(u>>1,nex);

    for(int i=0;i<N;i++) t2[i]=f[i];
    DFT(t2,1); DFT(tt,1);
    for(int i=0;i<N;i++) t2[i]=t2[i]*tt[i];
    DFT(t2,-1);
    for(int i=0;i<=T;i++) t2[i].x=((int)(t2[i].x+0.5))%P,t2[i].y=0;
    for(int i=T+1;i<N;i++) t2[i]=zero;

    for(int i=0;i<N;i++) f[i]=f[i]+t2[i];

    //DFT(tt,1);
    for(int i=0;i<N;i++) tt[i]=tt[i]*tt[i];
    DFT(tt,-1);
    for(int i=0;i<=T;i++) tt[i].x=((int)(tt[i].x+0.5))%P,tt[i].y=0;
    for(int i=T+1;i<N;i++) tt[i]=zero;
    if(u&1)
    {
        DFT(tt,1);
        for(int i=0;i<N;i++) tt[i]=tt[i]*f2[i];
        DFT(tt,-1);
        for(int i=0;i<=T;i++) tt[i].x=((int)(tt[i].x+0.5))%P,tt[i].y=0;
        for(int i=T+1;i<N;i++) tt[i]=zero;
    }

    if(u&1) for(int i=0;i<N;i++) f[i]=f[i]+tt[i];
}

int main()
{
    scanf("%d%d%d%d%d%d",&T,&P,&n,&O,&S,&U);
    O%=P; S%=P; U%=P;
    if(n>T)n=T;

    N=1,ln=0; while(N<=(T<<1)) N<<=1,ln++;
    for(int i=1;i<N;i++) id[i]=(id[i>>1]>>1)|(i&1?(1<<ln-1):0);

    for(int m=2;m<=N;m<<=1)
    {
        int t=m>>1,tt=N/m;
        for(int i=0;i<t;i++)
        {
            w[i*tt]=E(cos(2*pi*i/m),sin(2*pi*i/m));
            w[N-i*tt]=E(cos(2*pi*i/m),sin(-2*pi*i/m));
        }
    }

    g[0]=E(1,0);
    for(int i=1;i<=T;i++)
    {
        int fx=(i*i*O%P+i*S%P+U)%P;
        f1[i]=f2[i]=E(fx,0);
    }
    DFT(f2,1);

    get_(n,1);
    DFT(f,1); DFT(g,1);
    for(int i=0;i<N;i++) g[i]=g[i]*f[i];
    DFT(g,-1);

    int ans=(int)(g[T].x+0.5)%P;
    printf("%d\n",ans);

    return 0;
}