[bzoj4332][JSOI2012]分零食

题目大意

把正整数M分解成至多N份且每份不为0（注意1+2+3与2+3+1是不一样的即存在顺序性），一份x的价值是f(x)=a2*x*x+a1*x+a0，总价值为每一份价值的乘积。求所有情况下总价值的和，答案模mo。

M<=10000，mo<=255，N<=10^8，a2<=4，a1<=300，a0<=100

DP

设g[i,j]表示把j分解成i份的总价值和。

显然g[i,j]=∑j−1k=1g[i−1,j−k]∗f[k]

发现这是标准卷积形式，也就是说

g[i+1]=g[i]∗f

假如我们只需要分成N份，那么只需要把g快速幂即可，但是有一个“至多”。

我们设p[i,j]=∑ik=1g[k,j]

因为g要倍增，那我们可不可以让p也倍增？

事实证明可以。

p[i,j]=p[i/2,j]+∑i/2k=1g[i/2+k,j]

p[i,j]=p[i/2,j]+∑i/2k=1∑j−1l=1g[k,l]∗g[i/2,j−l]

p[i,j]=p[i/2,j]+∑j−1l=1∑i/2k=1g[k,l]∗g[i/2,j−l]

p[i,j]=p[i/2,j]+∑j−1l=1p[i/2,l]∗g[i/2,j−l]

后面那部分是标准卷积形式！

参考程序

#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef double db;
struct node{
db x,y;
node friend operator +(node a,node b){
node c;
c.x=a.x+b.x;c.y=a.y+b.y;
return c;
}
node friend operator -(node a,node b){
node c;
c.x=a.x-b.x;c.y=a.y-b.y;
return c;
}
node friend operator *(node a,node b){
node c;
c.x=a.x*b.x-a.y*b.y;c.y=a.x*b.y+a.y*b.x;
return c;
}
};
const db pi=acos(-1);
const int maxm=10000*4;
int b[maxm],g[maxm],h[maxm],p[maxm];
node e[maxm],f[maxm],tt[maxm];
db ce;
int i,j,k,l,t,n,m,mo,len,a2,a1,a0;
void DFT(node *a,int sig){
fo(i,0,len-1){
int p=0;
for(int j=0,tp=i;j<ce;j++,tp/=2) p=(p<<1)+(tp%2);
tt[p]=a[i];
}
for(int m=2;m<=len;m*=2){
int half=m/2;
fo(i,0,half-1){
node w;
w.x=cos(i*sig*pi/half),w.y=sin(i*sig*pi/half);
for(int j=i;j<len;j+=m){
node u=tt[j],v=tt[j+half]*w;
tt[j]=u+v;
tt[j+half]=u-v;
}
}
}
if (sig==-1)
fo(i,0,len-1) tt[i].x/=len;
fo(i,0,len-1) a[i]=tt[i];
}
void FFT(int *a,int *b,int *c){
int i;
fo(i,0,len-1) e[i].x=a[i],f[i].x=b[i],e[i].y=f[i].y=0;
DFT(e,1);DFT(f,1);
fo(i,0,len-1) e[i]=e[i]*f[i];
DFT(e,-1);
fo(i,0,m) c[i]=round(e[i].x);
fo(i,0,m) c[i]%=mo;
}
void solve(int n){
int i;
if (n==1){
fo(i,0,m) p[i]=g[i]=b[i];
return;
}
solve(n/2);
FFT(p,g,h);
fo(i,0,m) p[i]=(p[i]+h[i])%mo;
FFT(g,g,g);
if (n%2){
FFT(g,b,g);
fo(i,0,m) p[i]=(p[i]+g[i])%mo;
}
}
int main(){
scanf("%d%d%d%d%d%d",&m,&mo,&n,&a2,&a1,&a0);
len=1;
while (len<m*2) len*=2;
ce=db(log(len)/log(2));
b[0]=0;
fo(i,1,m) b[i]=(a2*i%mo*i%mo+a1*i%mo+a0%mo)%mo;
solve(n);
printf("%d\n",p[m]);
}