[JSOI2012]分零食

VIII.[JSOI2012]分零食

首先，暴力的DP是非常轻松的。设$f_{i,j}$为（前$i$个人，分了$j$颗糖）的种数，再设$g_i$为（一个人拿到$i$颗糖的快乐度），即$(Oi^2+Sx+U)$。

首先，我们可以将人数$A$与糖数$M$取$\min$，因为反正最多只有前$M$个人可以拿到糖，后面的人都是跑来撑场子的，可以直接忽略不计。

则有$f_{i,j}=\sum_{k=i-1}^{j-1}f_{i-1,k}*g_{j-k}$。直接暴力DP，复杂度$O(M^3)$，期望得分$40$。代码：

#include<bits/stdc++.h>
using namespace std;
int m,P,n,O,S,U,f[2][10100],res;
int calc(int x){
	x%=P;
	return (O*x*x+S*x+U)%P;
}
int main(){
	scanf("%d%d%d%d%d%d",&m,&P,&n,&O,&S,&U),n=min(n,m);
	for(int i=1;i<=m;i++)f[1][i]=calc(i);
	res=f[1][m];
	for(int i=2;i<=n;i++){
		for(int j=i;j<=m;j++){
			f[i&1][j]=0;
			for(int k=i-1;k<j;k++)f[i&1][j]=(f[i&1][j]+f[!(i&1)][k]*calc(j-k))%P;
		}
		res=(res+f[i&1][m])%P;
		for(int j=1;j<=m;j++)printf("%d ",f[i&1][j]);puts("");
	}
	printf("%d\n",res);
	return 0;
}

观察转移式子$f_{i,j}=\sum_{k=i-1}^{j-1}f_{i-1,k}*g_{j-k}$，可以发现后面的东西是一个卷积的形式。我们设一个函数$f(x)$，表示原来$g$的含义。则有$ans=\sum_{i=1}^A f^i(M)$。直接暴力FFT，复杂度$O(M^2\log M)$，期望得分$60$。代码：

#include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1);
int m,P,n,O,S,U,res,lim=1,lg,rev[1<<16],f[1<<16],g[1<<16],bs[1<<16];
int calc(int x){
	x%=P;
	return (O*x*x+S*x+U)%P;
}
struct cp{
	double x,y;
	cp(double u=0,double v=0){x=u,y=v;}
	friend cp operator +(const cp &u,const cp &v){return cp(u.x+v.x,u.y+v.y);}
	friend cp operator -(const cp &u,const cp &v){return cp(u.x-v.x,u.y-v.y);}
	friend cp operator *(const cp &u,const cp &v){return cp(u.x*v.x-u.y*v.y,u.x*v.y+u.y*v.x);}
}p[1<<16],q[1<<16];
void FFT(cp *a,int tp){
	for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int md=1;md<lim;md<<=1){
		cp rt=cp(cos(pi/md),tp*sin(pi/md));
		for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
			cp w=cp(1,0);
			for(int i=0;i<md;i++,w=w*rt){
				cp x=a[pos+i],y=w*a[pos+md+i];
				a[pos+i]=x+y;
				a[pos+md+i]=x-y;
			}
		}
	}
}
void mul(int *a,int *b,int *c){
	for(int i=0;i<lim;i++)p[i]=cp(a[i],0);
	for(int i=0;i<lim;i++)q[i]=cp(b[i],0);
	FFT(p,1),FFT(q,1);
	for(int i=0;i<lim;i++)p[i]=p[i]*q[i];
	FFT(p,-1);
	for(int i=0;i<lim;i++)c[i]=((int)(p[i].x/lim+0.5))%P;
}
int main(){
	scanf("%d%d%d%d%d%d",&m,&P,&n,&O,&S,&U),n=min(n,m);
	while(lim<=(2*m))lim<<=1,lg++;
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1)); 
	for(int i=1;i<=m;i++)bs[i]=calc(i);
	f[0]=1;
	for(int i=1;i<=n;i++){
		mul(f,bs,f),res=(res+f[m])%P;
		for(int j=0;j<i;j++)f[j]=0;
		for(int j=m+1;j<lim;j++)f[j]=0;
	}
	printf("%d\n",res);
	return 0;
}

如果答案不是求$f^1$到$f^A$所有东西的和的话，我们完全可以直接上多项式快速幂。考虑令一个函数$g_x=\sum_{i=1}^x f^i$，则答案为$g_A(M)$。

当$2|x$时，可以有$g_x=g_{x/2}+g_{x/2}*f^{x/2}$，因为

$\begin{aligned}&g_{x/2}+g_{x/2}*f^{x/2}\\=&(\sum_{i=1}^{x/2}f^i)+(\sum_{i=1}^{x/2}f^i)f^{x/2}\\=&\sum_{i=1}^{x/2}f^i+\sum_{i=1}^{x/2}f^{i+x/2}\\=&\sum_{i=1}^xf^i\\=&g_x\end{aligned}$

而有$f^x=f^{x/2}*f^{x/2}$。

当$2\nmid x$时，有$g_x=g_{x-1}+f^x,f^x=f^{x-1}*f^1$。

然后就直接像快速幂一样处理即可。复杂度$O(M\log^2 M)$，期望得分$100$。代码：

#include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1);
int m,P,n,O,S,U,res,lim=1,lg,rev[1<<16],f[1<<16],g[1<<16],bs[1<<16],ass[1<<16];
int calc(int x){
	x%=P;
	return (O*x*x+S*x+U)%P;
}
struct cp{
	double x,y;
	cp(double u=0,double v=0){x=u,y=v;}
	friend cp operator +(const cp &u,const cp &v){return cp(u.x+v.x,u.y+v.y);}
	friend cp operator -(const cp &u,const cp &v){return cp(u.x-v.x,u.y-v.y);}
	friend cp operator *(const cp &u,const cp &v){return cp(u.x*v.x-u.y*v.y,u.x*v.y+u.y*v.x);}
}p[1<<16],q[1<<16];
void FFT(cp *a,int tp){
	for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int md=1;md<lim;md<<=1){
		cp rt=cp(cos(pi/md),tp*sin(pi/md));
		for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
			cp w=cp(1,0);
			for(int i=0;i<md;i++,w=w*rt){
				cp x=a[pos+i],y=w*a[pos+md+i];
				a[pos+i]=x+y;
				a[pos+md+i]=x-y;
			}
		}
	}
}
void mul(int *a,int *b,int *c){
	for(int i=0;i<lim;i++)p[i]=cp(a[i],0);
	for(int i=0;i<lim;i++)q[i]=cp(b[i],0);
	FFT(p,1),FFT(q,1);
	for(int i=0;i<lim;i++)p[i]=p[i]*q[i];
	FFT(p,-1);
	for(int i=1;i<=m;i++)c[i]=((int)(p[i].x/lim+0.5))%P;
}
void ksm(int x){
	if(x==1){memcpy(f,bs,sizeof(bs)),memcpy(g,bs,sizeof(bs));return;}
	ksm(x>>1);
	mul(f,g,ass);
	for(int i=0;i<lim;i++)g[i]=(g[i]+ass[i])%P;
	mul(f,f,f);
	if(x&1){
		mul(f,bs,f);
		for(int i=0;i<lim;i++)g[i]=(g[i]+f[i])%P;
	}
}
int main(){
	scanf("%d%d%d%d%d%d",&m,&P,&n,&O,&S,&U),n=min(n,m);
	while(lim<=(2*m))lim<<=1,lg++;
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1)); 
	for(int i=1;i<=m;i++)bs[i]=calc(i);
	ksm(n);
	printf("%d\n",g[m]);
	return 0;
}