P5075 [JSOI2012]分零食

$Pro$

你可以把一个数 $m$ 分成不超过 $n$ 个数 $a_i$ ，每一种方案对答案的贡献是 $\prod f(a_i)$ ，求最后的贡献值。

$Sol$

很容易列出 $O(nm^2)$ 的 $dp$ 方程

$dp[i][y]=\sum dp[i-1][y-x]+f[x]$

其中 $f[x]=f(x)=Ox^2+Sx+U$

（ $f[x]$ 代表多项式 $f$ 的第 $x$ 项代表的值， $f(x)$ 代表题目所说的函数）

注意到 $dp[i-1][y-x]+f[x]$ 是一个卷积的形式，可以化成 $dp[i]=dp[i-1] * f$

于是发现 $dp_i=f^i$ ，

但是发现题目要求的是 $\sum _{i=1} ^n dp_i$ ，于是不能暴力 $O(nm\log m)$ 枚举，怎么办呢。

先看一道简化版的题目：

给你 $a,n,p$ ，求 $\sum _{i=1}^n a^i \mod p$ 的值。（ $n \le 10^{18}$ ）

不要用数学方法乱搞。

考虑类似于快速幂，使用倍增，不同的是，我们在倍增过程中记录两个变量 $F,G$

其中 $F=\sum _{i=1}^n a^i \mod p,G=a^n \mod p$

发现 $F'=F+F \times G,G'=G \times G$

当 $n$ 为奇数的时候还有处理一下边角， $G=G \times a,F=F+G$

于是得到以下代码：

#include <bits/stdc++.h>
#define int long long
using namespace std;
inline int read(){
    int x=0,f=1;
    char ch=getchar();
    while (ch<'0'||ch>'9'){
        if (ch=='-') f=-1;
        ch=getchar();
    }
    while (ch>='0'&&ch<='9'){
        x=(x<<1)+(x<<3)+(ch^'0');
        ch=getchar();
    }
    return x*f;
}
int F;// ans
int G;// a^i
int a,n,p;
inline void ksm(int n){
    if (n==1) return F=G=a,void();
    ksm(n>>1ll);
    F=(F+F*G)%p;
    G=(G*G)%p;
    if (n&1){
        G=(G*a)%p;
        F=(F+G)%p;
    }
}
#undef int
int main(){
#define int long long
    a=read(),n=read(),p=read();
    ksm(n);
    printf("%lld\n",F);
}

回到本题，也能得到类似的结论，设 $F=\sum _{i=1} ^n f^i,G=f^n$ ，我们有 $F'=F+F \times G,G'=G \times G$

于是就得到了 $O(m \log m \log n)$ 的解法：

#include <bits/stdc++.h>
#define MAXN 50005
using namespace std;
const double PI=acos(-1.0);
inline int read(){
	int x=0,f=1;
	char ch=getchar();
	while (ch<'0'||ch>'9'){
		if (ch=='-') f=-1;
		ch=getchar();
	}
	while (ch>='0'&&ch<='9'){
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}
struct Complex{
	double x,y;
};
inline Complex operator + (const Complex &A,const Complex &B){
	return Complex{A.x+B.x,A.y+B.y};
}
inline Complex operator - (const Complex &A,const Complex &B){
	return Complex{A.x-B.x,A.y-B.y};
}
inline Complex operator * (const Complex &A,const Complex &B){
	return Complex{A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x};
}
int r[MAXN];
inline void FFT(Complex *A,int n,int type){
	for (register int i=0;i<n;++i) if (i<r[i]) swap(A[i],A[r[i]]);
	for (register int i=1;i<n;i<<=1){
		int R=(i<<1);
		Complex Wn=Complex{cos(2*PI/R),type*sin(2*PI/R)};
		for (register int j=0;j<n;j+=R){
			Complex w=Complex{1,0};
			for (register int k=0;k<i;++k,w=w*Wn){
				Complex x=A[j+k],y=w*A[i+j+k];
				A[j+k]=x+y;
				A[i+j+k]=x-y;
			}
		}
	}
}
struct Poly{
	int num[MAXN];
};
Poly F,G,org,temp;
int sz,m,n,o,s,u,L,mod;
inline int f(int x){
	return (o*x*x%mod+s*x%mod+u)%mod;
}
Complex a[MAXN],b[MAXN];
inline void Mul(Poly &des,const Poly &A,const Poly &B){
	for (register int i=0;i<=sz;++i) a[i]=Complex{(double)A.num[i],0},b[i]=Complex{(double)B.num[i],0};
	FFT(a,sz,1),FFT(b,sz,1);
	for (register int i=0;i<=sz;++i) a[i]=a[i]*b[i];
	FFT(a,sz,-1);
	for (register int i=1;i<=m;++i) des.num[i]=((int)(a[i].x/sz+0.5))%mod;
}
inline void Add(Poly &A,const Poly &B){
	for (register int i=1;i<=m;++i) A.num[i]=(A.num[i]+B.num[i])%mod;
}
inline void Init(int m){
	sz=1,L=0;
	while (sz<=2*m) sz<<=1,L++;
	for (register int i=0;i<=sz;++i){
		r[i]=(r[i>>1]>>1|((i&1)<<(L-1)));
	}
}
inline void ksm(int n){
	if (n==1){
		F=org,G=org;
		return ;
	}
	ksm(n>>1);
	Mul(temp,F,G);
	Add(F,temp);
	Mul(G,G,G);
	if (n&1){
		Mul(G,G,org);
		Add(F,G);
	}
}
int main(){
	m=read(),mod=read();
	Init(m);
	n=read(),o=read(),s=read(),u=read();
	for (register int i=1;i<=m;++i) org.num[i]=f(i);
	ksm(n);
	printf("%d\n",F.num[m]);
}