[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;
}