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

https://www.lydsy.com/JudgeOnline/problem.php?id=4332

因为如果一位小朋友得不到糖果，那么在她身后的小朋友们也都得不到糖果。

所以设g[i][j] 表示前i位小朋友，分到j个糖果，且前i位小朋友都分到糖果的方案数

令F(x) 表示分到x个糖果的欢乐程度

∴g[i][j] = ∑ g[i-1][j-k]*F(k)

记g[i]=g[i-1]*F，则 g[i]=F ^ i

但是要求的是 Σ g[i][m]

记f[n]=Σ g[i]  i∈[1,n] ，那么ans=f[n][m]

f[n]=Σ g[i]  i∈[1,n]

     =Σ f(n/2)+Σ g[i]  i∈[n/2+1,n]

     =Σ f(n/2)+Σ F^i  i∈[n/2+1,n]

     =Σ f(n/2)+Σ F^(n/2+i)  i∈[1,n/2]

     =Σ f(n/2)+F^(n/2) * Σ F^i  i∈[1,n/2]

     =Σ f(n/2)+g(n/2)*f(n/2)

然后可以分治解决

如果n是奇数，f(n)=f(n-1)+g[n]=f(n-1)+g(n-1)*f

边界条件：g[][0]=1

#include<cmath>
#include<cstdio>
#include<algorithm>

using namespace std;

const int M=1<<17;

#define N 10001

int m,mod;

int r[M+1];
int len;

const double pi=acos(-1);

struct Complex
{
    double x,y;
    Complex() { }
    Complex(double x_,double y_):x(x_),y(y_) { }
    Complex operator + (Complex p)
    {
        Complex C;
        C.x=x+p.x;
        C.y=y+p.y;
        return C;
    }
    Complex operator - (Complex p)
    {
        Complex C;
        C.x=x-p.x;
        C.y=y-p.y;
        return C;
    }
    Complex operator * (Complex p)
    {
        Complex C;
        C.x=x*p.x-y*p.y;    
        C.y=x*p.y+y*p.x;
        return C;
    }
    void clear()
    {
        x=y=0;
    }
};

typedef Complex E;

E F[M+1],f[M+1],g[M+1],tmp[M+1];

void FFT(E *a,int ty)
{
    for(int i=0;i<len;++i)
        if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=1;i<len;i<<=1)
    {
        E wn(cos(pi/i),ty*sin(pi/i));
        for(int p=i<<1,j=0;j<len;j+=p)
        {
            E w(1,0);
            for(int k=0;k<i;++k,w=w*wn)
            {
                E x=a[j+k],y=w*a[i+j+k];
                a[j+k]=x+y; a[i+j+k]=x-y;
            }
        }
    }
    if(ty==-1)
    {
        for(int i=0;i<len;++i) a[i].x=a[i].x/len,a[i].x=int(a[i].x+0.5)%mod,a[i].y=0;
    }
}

void solve(E *f,E *g,int n)
{
    if(!n)
    {
        g[0].x=1;
        return;
    }
    if(n&1)
    {
        solve(f,g,n-1);
        FFT(g,1);
        for(int i=0;i<len;++i) g[i]=g[i]*F[i];
        FFT(g,-1);
        for(int i=0;i<=m;++i) f[i]=f[i]+g[i];
        for(int i=0;i<=m;++i) f[i].x=int(f[i].x)%mod,f[i].y=0;
        for(int i=m+1;i<len;++i) f[i].clear(),g[i].clear();
    }
    else
    {
        solve(f,g,n/2);
        for(int i=0;i<len;++i) tmp[i]=f[i];
        FFT(tmp,1);
        FFT(g,1);
        for(int i=0;i<len;++i) tmp[i]=tmp[i]*g[i];
        FFT(tmp,-1);
        for(int i=0;i<len;++i) g[i]=g[i]*g[i];
        FFT(g,-1);
        for(int i=0;i<=m;++i) f[i]=f[i]+tmp[i];
        for(int i=0;i<=m;++i) f[i].x=int(f[i].x)%mod,f[i].y=0;
        for(int i=m+1;i<len;++i) f[i].clear(),g[i].clear();
    }
}

int main()
{
    int n,o,s,u;
    scanf("%d%d%d%d%d%d",&m,&mod,&n,&o,&s,&u);
    //F[0].x=1;
    for(int i=1;i<=m;++i) F[i].x=(o*i*i+s*i+u)%mod;
    int l=0;
    for(len=1;len<=m+m;len<<=1,l++);
    for(int i=0;i<len;++i) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
    FFT(F,1);
    solve(f,g,n);
    printf("%d",int(f[m].x));
    return 0;
}