这篇博客介绍了如何解决LOJ #6672的数学问题，涉及到生成函数和整式递推的概念。作者通过建立数学模型，解析出整式递推公式，并利用生成函数进行求解。文章还提到了计算过程中遇到的挑战，如计算大量项的问题，并分享了AC代码以及整式递推求逆元的离线解决方案。
题目
没有latex就没有推式子的动力怎么破？
设
f
i
,
j
f_{i,j}
fi,j表示在前
j
j
j个回合里，血量为
0
0
0了
i
i
i个回合且第
j
j
j个回合血量为
0
0
0的方案数。
f
i
,
j
=
∑
k
=
0
j
−
1
f
i
−
1
,
k
g
j
−
k
f_{i,j} = \sum_{k=0}^{j-1} f_{i-1,k}g_{j-k}
fi,j=∑k=0j−1fi−1,kgj−k
其中
g
i
g_{i}
gi表示经过
i
+
1
i+1
i+1个回合只有第一个和最后一个回合血量为
0
0
0的方案数。
g
i
=
h
i
−
2
g_{i} = h_{i-2}
gi=hi−2
其中
h
i
h_{i}
hi表示经过
i
+
1
i+1
i+1个回合（血量不为负数）且第一个和最后一个回合血量为
0
0
0的方案数，那么上式的意义就是我第一步往上，最后一步往下，中间血量不
<
1
<1
<1。
计算
h
h
h时，通过枚举第一次除起点外血量为
0
0
0的时机来转移，
N
a
m
e
l
y
:
Namely:
Namely:
h
i
=
∑
j
=
1
i
g
j
h
i
−
j
=
∑
j
=
1
i
h
j
−
2
h
i
−
j
h_{i} = \sum_{j=1}^i g_jh_{i-j} = \sum_{j=1}^i h_{j-2}h_{i-j}
hi=∑j=1igjhi−j=∑j=1ihj−2hi−j
但是因为
h
−
1
h_{-1}
h−1不存在，所以上式应写成：
h
i
=
g
1
h
i
−
1
+
∑
j
=
2
i
h
j
−
2
h
i
−
j
=
h
i
−
1
+
∑
j
=
2
i
h
j
−
2
h
i
−
j
h_{i} = g_1h_{i-1}+\sum_{j=2}^i h_{j-2} h_{i-j} = h_{i-1}+\sum_{j=2}^i h_{j-2}h_{i-j}
hi=g1hi−1+∑j=2ihj−2hi−j=hi−1+∑j=2ihj−2hi−j
这个显然是卡特兰数的亲戚，我们马上想到用整式递推。
开始逆着推生成函数回去。
H
(
x
)
=
∑
i
=
0
h
i
x
i
=
H
(
x
)
x
+
x
2
H
(
x
)
2
+
1
H(x) = \sum_{i=0} h_ix^i = H(x)x + x^2H(x)^2 + 1
H(x)=∑i=0hixi=H(x)x+x2H(x)2+1
说明
H
(
x
)
H(x)
H(x)在有理分式域是代数的，存在整式递推式。
x
2
H
(
x
)
2
+
(
x
−
1
)
H
(
x
)
+
1
=
0
x^2H(x)^2 + (x-1)H(x) + 1 = 0
x2H(x)2+(x−1)H(x)+1=0
H
(
x
)
=
1
−
x
±
−
3
x
2
−
2
x
+
1
2
x
2
H(x) = \frac {1-x \pm \sqrt{-3x^2-2x+1}}{2x^2}
H(x)=2x21−x±−3x2−2x+1
H
(
x
)
H(x)
H(x)应为形式幂级数，故此处应取负号。
G
(
x
)
=
∑
i
=
0
g
i
x
i
=
H
(
x
)
x
2
+
x
G(x) = \sum_{i=0}g_ix^i =H(x)x^2+x
G(x)=∑i=0gixi=H(x)x2+x这里
+
x
+x
+x指的是
g
1
=
1
g_1=1
g1=1
F
i
(
x
)
=
∑
j
=
0
g
j
x
j
=
F
i
−
1
(
x
)
G
(
x
)
=
F
i
−
1
(
x
)
(
x
+
x
2
H
(
x
)
)
F_i(x) = \sum_{j=0}g_jx^j=F_{i-1}(x)G(x) = F_{i-1}(x)(x+x^2H(x))
Fi(x)=∑j=0gjxj=Fi−1(x)G(x)=Fi−1(x)(x+x2H(x))
F
1
(
x
)
=
1
F_1(x) = 1
F1(x)=1
所以
F
k
(
x
)
=
(
x
+
x
2
H
(
x
)
)
k
−
1
=
(
1
+
x
−
−
3
x
2
−
2
x
+
1
2
)
k
−
1
F_k(x) = (x+x^2H(x))^{k-1} = (\frac {1+x-\sqrt{-3x^2-2x+1}}2)^{k-1}
Fk(x)=(x+x2H(x))k−1=(21+x−−3x2−2x+1)k−1
然后我们发现我们的所有过程都是符合整式递推数列的生成函数变换要求的。
（
v
(
x
)
=
1
+
x
−
−
3
x
2
−
2
x
+
1
2
v(x)=\frac {1+x-\sqrt{-3x^2-2x+1}}2
v(x)=21+x−−3x2−2x+1为代数形式幂级数，而
u
(
x
)
=
x
k
−
1
u(x)=x^{k-1}
u(x)=xk−1也显然微分有限，所以
F
k
(
x
)
=
u
(
v
(
x
)
)
F_k(x)=u(v(x))
Fk(x)=u(v(x))也微分有限，为
P
−
r
e
c
u
r
s
i
v
e
P-recursive
P−recursive）。
但是我们发现一个很难的事情，那就是
F
k
(
x
)
F_k(x)
Fk(x)的前
k
−
1
k-1
k−1项都为
0
0
0，我们显然不能用这些项求出一个较短的整式递推式，所以我们就用第
k
k
k项之后的项来求即可。
那么问题来了，怎么算
k
<
=
1
e
7
k<=1e7
k<=1e7之后的几项？
简单，
F
k
(
x
)
=
(
1
+
x
−
−
3
x
2
−
2
x
+
1
2
x
)
k
−
1
x
k
−
1
F_k(x) = (\frac {1+x-\sqrt{-3x^2-2x+1}}{2x})^{k-1}x^{k-1}
Fk(x)=(2x1+x−−3x2−2x+1)k−1xk−1
那么算
(
1
+
x
−
−
3
x
2
−
2
x
+
1
2
x
)
k
−
1
(\frac {1+x-\sqrt{-3x^2-2x+1}}{2x})^{k-1}
(2x1+x−−3x2−2x+1)k−1的前几项即可。
所以我们可以打一个多项式开根
那么我们把
H
(
x
)
H(x)
H(x)用
n
2
n^2
n2递推式算出来，然后算出
1
+
x
H
(
x
)
1+xH(x)
1+xH(x)
再快速幂暴力卷
k
−
1
k-1
k−1次即可。
剩下的就是整式递推的事了.
PS:整式递推求逆元可以离线
O
(
n
)
O(n)
O(n)求。
一开始快速幂打错了疯狂assertfail，拿zzq的原版板子也assertfail的调了半天
充分学习了assert对于调试端来说是多么的重要
A C C o d e \mathrm{AC \ Code} AC Code(内附输出整式递推式器)
#include<bits/stdc++.h>
#define rep(i,j,k) for(int i=(j),LIM=(k);i<=LIM;i++)
#define per(i,j,k) for(int i=(j),LIM=(k);i>=LIM;i--)
#define vc vector
#define vi vc<int>
#define Ct const
#define mod 998244353
#define LL long long
#define maxn 10000007
#define St 13
using namespace std;
int n,K,h[maxn],f[maxn];
void MUL(int *A,int *B,int *C){
static int t[St+1];
rep(i,0,St){ t[i]=0;
rep(j,0,i) t[i]=(t[i]+(LL)A[j]*B[i-j])%mod;}
rep(i,0,St) C[i]=t[i];
}
int Pow(int b,int k){ int r=1;for(;k;k>>=1,b=1ll*b*b%mod) if(k&1) r=1ll*r*b%mod;return r; }
vc<vi >PR(Ct vi&a,int D){
int N=a.size(),B=(N+2)/(D+2),C=(D+1)*B,R=N-(B-1),c=C,p;
vc<vi >b(R,vi(C));
rep(i,0,R-1) rep(j,0,B-1) for(int k=0,x=a[i+j];k<=D;k++,x=(LL)x*i%mod) b[i][j*(D+1)+k]=x;
rep(i,0,C-1){ rep(j,(p=-1,i),R-1) if(b[j][i]){p=j;break;}
if(p==-1){c=i;break;}swap(b[i],b[p]);
p=Pow(b[i][i],mod-2);rep(j,i,C-1) b[i][j]=(LL)b[i][j]*p%mod;
rep(j,i+1,R-1) if(p=b[j][i]) rep(k,i,C-1) b[j][k]=(b[j][k]-(LL)b[i][k]*p)%mod;
}assert(c^C);
per(i,c-1,0) if(b[i][c]) rep(j,0,i-1) b[j][c]=(b[j][c]-(LL)b[j][i]*b[i][c])%mod;
int M=c/(D+1)+1,S;
vc<vi >r(M,vi(D+1));
rep(i,(r[0][c%(D+1)]=1,0),c-1) r[M-i/(D+1)-1][i%(D+1)]=-b[i][c];
for(int i=0;i<M;i++) rep(j,0,D) rep(k,(S=(j+1ll)*(-M+1)%mod,j+1),D)
r[i][j]=(r[i][j]+(LL)S*r[i][k])%mod,S=S*(-M+1ll)%mod*(k+1)%mod*Pow(k-j+1,mod-2)%mod;
return r;
}
inline int cal(Ct vi&a,int x){ int s=1,r=0;rep(i,0,a.size()-1) r=(r+(LL)s*a[i])%mod,s=(LL)s*x%mod; return r; }
int val[maxn],pr[maxn],sc,inv[maxn];
int pan(int a){ a=(a+mod)%mod;if(a>(mod/2)) a-=mod;return a; }
int main(){
scanf("%d%d",&n,&K);h[0]=h[1]=1;
rep(i,2,St){ h[i]=h[i-1];
rep(j,2,i) h[i]=(h[i]+(LL)h[j-2]*h[i-j])%mod;}
per(i,St,1) h[i]=h[i-1];
h[1]=f[0]=1;
for(int t=K-1;t;t>>=1,MUL(h,h,h)) if(t&1) MUL(f,h,f);
vc<vi >r=PR(vi(f,f+St+1),2);
/* printf("@%d\n",r.size());
rep(i,0,r.size()-1){
rep(j,0,r[i].size()-1)
printf("%dx^%d%c",pan(r[i][j]),j,"+\n"[j==r[i].size()-1]);
puts("");
}*/
rep(i,St+1,n-K+1) val[i]=cal(r[0],i),pr[i]=(LL)(i==St+1?1:pr[i-1])*val[i]%mod;
sc = Pow(pr[n-K+1],mod-2);
per(i,n-K+1,St+1) inv[i]=(LL)sc*(i==St+1?1:pr[i-1])%mod,sc=(LL)sc*val[i]%mod;
rep(i,St+1,n-K+1){
rep(j,1,r.size()-1) f[i]=(f[i]-(LL)f[i-j]*cal(r[j],i))%mod;
f[i]=1ll*f[i]*inv[i]%mod;
}
printf("%d\n",(f[n-K+1]+mod)%mod);
}
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
