Steven_MengのBlog

$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$

于是得到以下代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#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)$的解法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#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]);
}

 评论