ProPro

你可以把一个数mm分成不超过nn个数aia_i,每一种方案对答案的贡献是f(ai)\prod f(a_i),求最后的贡献值。

SolSol

很容易列出O(nm2)O(nm^2)dpdp方程

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

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

f[x]f[x]代表多项式ff的第xx项代表的值,f(x)f(x)代表题目所说的函数)

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

于是发现dpi=fidp_i=f^i

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


先看一道简化版的题目:

给你a,n,pa,n,p,求i=1naimodp\sum _{i=1}^n a^i \mod p的值。(n1018n \le 10^{18}

不要用数学方法乱搞。

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

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

发现F=F+F×G,G=G×GF'=F+F \times G,G'=G \times G

nn为奇数的时候还有处理一下边角,G=G×a,F=F+GG=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=i=1nfi,G=fnF=\sum _{i=1} ^n f^i,G=f^n,我们有F=F+F×G,G=G×GF'=F+F \times G,G'=G \times G

于是就得到了O(mlogmlogn)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]);
}