计算小练习（3）Treasure 另解
NaCly_Fish
·
·
个人记录
先讲个笑话：这篇题解我花了超过 5h 才写出来。
此题的组合推导比较简单，如果需要的话，可以先去题解区看看。
我们希望求出所有 a_{i,2i-m} 的值，不过我们并不需要像原题解中那样，找出这一条线上的递推式。如果能以 \Theta(1) 从 a_{n,k} 附近的值推出 a_{n,k+1} 和 a_{n+1,k}，那么问题就解决了。
原递推式为
a_{n,k}=a_{n-2,k-1}+a_{n-1,k-1}+a_{n-1,k}
对列建立生成函数，得到
C_k(x)=(x^2+x)C_{k-1}(x)+xC_k(x)
那么对于 k \ge 1，可以直接解得
C_k(x) = \frac{2x}{1-x^2}\left( \frac{x(1+x)}{1-x}\right)^k
$$(n-k)a_{n+1,k}=2ka_{n,k}+(n-k)a_{n-1,k} $$
****
现在考虑另一个方向上的递推，设第 $n$ 行的生成函数：
$$R_n(x)=(1+x)R_{n-1}(x)+xR_{n-2}(x)$$
设
$$u^2-(1+x)u-x=0 \quad \texttt{A}$$
取常数项非零的那个根，则第 $n$ 行就可以用 $u$ 来表示：
$$R_n(x) \equiv \frac{2u^{n+1}}{u^2+2u-1} \pmod{x^{n+1}}$$
这里稍微展开介绍一下 $R_n(x)$ 的整式递推是怎么来的，设
$$f=\frac{2x^{n+1}}{x^2+2x-1}$$
它满足 ODE
$$((n-1)x^2+2nx-(n+1))f=x(x^2+2x-1)f'$$
要求出 $F=f \circ u$ 的 ODE，按机械化方法首先需要把 $u'$ 用 $u$ 的整式来表示。把式 $\texttt{A}$ 求导，整理一下得到
$$u'= \frac{u+1}{2u-(1+x)} $$
可以证明对于任意代数方程，$u'$ 都可以表示为关于 $u$ 的整式。关键在于求 $2u-(1+x)$ 的逆元（在满足 $\texttt{A}$ 的情况下），具体求法应用多项式 exgcd 即可：
$$u'=\frac{1-x-(3+x)u}{1+6x+x^2} \quad \texttt{B}$$
对 $f$ 的 ODE 两边同时右复合 $u$ 得到
$$((n-1)u^2+2nu-(n+1))(f \circ u)=u(u^2+2u-1)(f'\circ u)$$
这里我们逆用链式法则 $(f\circ u)'=(f'\circ u)u'$，就可以得到
$$((n-1)u^2+2nu-(n+1))u'F=u(u^2+2u-1)F'$$
接下来的操作我愿意称之为「约化」（Reduce），首先利用式 $\texttt{B}$ 把 $u'$ 都用 $u$ 来替换；然后利用式 $\texttt{A}$ 来降次（也就是多项式取模），把 $u$ 的次数降到不超过 $1$。于是就有
$$\begin{aligned}((-(n+1)+(4n-2)x+(n-1)x^2)+((5n-3)+(6n-4)x+(n-1)x^2)u) & F \\ -(1+6x+x^2)(x + (2 + x) u) & F' = 0 \end{aligned}$$
当然，$x+(2+x)u$ 也可以求逆，这样就能让 $F'$ 表示为 $F$ 乘上关于 $u$ 的整式。由此得出
$$(-((n+1)+(3n+1)x)+((n+1)-(n-1)x)u)F +x(1+6x+x^2) F'=0 $$
按照套路继续求导，然后约化。注意其中产生的 $F'$ 还可以用 $F$ 来表示，这样得到的结果中，$F''$ 可以由 $F'$ 乘上 $u$ 的整式来表示：
$$\begin{aligned} (((n^2+n)+(7n^2-n-8)x+(7n^2-17n-6)x^2+(n^2-7n-2)x^3) & - \\((n^2+n)+(5n^2-3n-8)x-(5n^2-7n+6)x^2-(n^2-3n+2)x^3) u) & F \\ -x^2(1+6x+x^2)^2 & F'' = 0 \end{aligned}$$
终于到了激动人心的导出 ODE 时间，把上面两个式子线性组合起来，把 $F$ 前面的 $u$ 消掉就好了。最终得到的是
$$\begin{aligned} 2n(1-n)(1+n)(1-x) & F \\ +((n + n^2) + (5 n^2 - 3 n - 8) x - (5 n^2 - 7 n + 6) x^2 - (n^2 -
3 n + 2) x^3) & F' \\ -x(1+6x+x^2)((n+1)-(n-1)x)& F'' = 0 \end{aligned}$$
由此可以得到整式递推式为
$$\begin{aligned} (n+1)(k+1)(n-k) & a_{n,k+1} \\ +(2n(1-n^2)+(5n^2-3n-8)k-(7+5n)k(k-1)) & a_{n,k} \\ -(2n(1-n^2)+(5n^2-7n+6)(k-1)+(7-5n)(k-1)(k-2)) & a_{n,k-1} \\ -(k-2)(n-1)(n-k+1) & a_{n,k-2} = 0 \end{aligned}$$
****
根据递推式的阶数，我们维护一个 $2$ 行 $3$ 列的矩形范围内 $a$ 的值即可。有一点常数优化的地方，就是根据 $b_{n,k}=a_{n+1,k}-a_{n,k}$，再利用竖直方向的递推式，就可以直接根据矩形内的值求出 $b$。
注意到递推边界处有些不便于操作，那直接 $\mathcal O(n)$ 求出那些位置的值即可，总时间复杂度还是不变的。
update：[飞雨烟雁](https://www.luogu.com.cn/user/375984)同学指出：其实不一定总是要纵向和横向递推，也可以纵向+斜向。比如此题中做代换 $a'_{n,k} = a_{n+k,k}$ 的话，$a'_{n,k}$ 的横向递推为简单的二阶整式递推，也就对应了 $a_{n,k}$ 的斜方向递推。这样就可以最大幅度地减少人类的手算工作。
可以发现此做法有以下优势：跑得比高斯消元更快、便于处理取 $a$ 的位置更复杂的情况、~~有一定手算 ODE 的可行性~~。参考代码如下：
```cpp
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<vector>
#define N 6000003
#define ll long long
#define p 998244353
using namespace std;
inline int power(int a,int t){
int res = 1;
while(t){
if(t&1) res = (ll)res*a%p;
a = (ll)a*a%p;
t >>= 1;
}
return res;
}
int f[2][3],inv[N],fac[N],ifac[N]; // f[i][j] = a_{n-i,k-j}
int tn,tk;
void init(int n){
inv[1] = fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
for(int i=2;i<=n;++i) fac[i] = (ll)fac[i-1]*i%p;
ifac[n] = power(fac[n],p-2);
for(int i=n-1;i;--i) ifac[i] = (ll)ifac[i+1]*(i+1)%p;
for(int i=2;i<=n;++i) inv[i] = (ll)fac[i-1]*ifac[i]%p;
}
void down(){
int a[3]; // a[j] = a_{n+1,k-j}
for(int j=0;j<3;++j) a[j] = (2ll*(tk-j)*inv[tn-(tk-j)]%p*f[0][j]+f[1][j])%p;
for(int j=0;j<3;++j) f[1][j] = f[0][j],f[0][j] = a[j];
++tn;
}
void right(){
int a[2]; // a[i] = a_{n-i,k+1};
for(int i=0;i<2;++i){
a[i] = (-(2*tn*(1-(ll)tn*tn%p)%p + (5ll*tn*tn-3ll*tn-8)%p*tk%p - (7+5ll*tn)*tk%p*(tk-1)%p)*f[i][0]%p +
(2*tn*(1-(ll)tn*tn%p)%p + (5ll*tn*tn-7ll*tn+6)%p*(tk-1)%p + (7-5ll*tn)*(tk-1)%p*(tk-2)%p)*f[i][1]%p -
(ll)(tk-2)*(1-tn)%p*(tn-tk+1)%p*f[i][2])%p * inv[tn+1]%p*inv[tk+1]%p*inv[tn-tk]%p;
tn--;
}
tn += 2,++tk;
for(int i=2;i>0;--i) f[0][i] = f[0][i-1],f[1][i] = f[1][i-1];
f[0][0] = a[0],f[1][0] = a[1];
}
inline int binom(int n,int m){
if(n<m||m<0) return 0;
return (ll)fac[n]*ifac[m]%p*ifac[n-m]%p;
}
inline int solve(int n,int k){ // [x^{n-k-1}] 2 (1+x)^(k-1)/(1-x)^(k+1)
if(k==0) return 1;
n -= k+1;
if(n<0) return 0;
ll res = 0;
for(int i=0;i<=n;++i) res += (ll)binom(k-1,i)*binom(n-i+k,k)%p;
return res*2%p;
}
int n,k,m;
int a[N],b[N]; // a[i] = a_{i,2i-m}
ll ans;
int main(){
scanf("%d%d%d",&n,&k,&m);
init(n<<1);
tn = (m>>1)+2;
tk = 2*tn-m;
for(int i=0;i<2;++i)
for(int j=0;j<3;++j)
f[i][j] = solve(tn-i,tk-j);
a[tn-1] = f[1][2],b[tn-1] = f[0][2]-f[1][2];
if(tk==4) a[tn-2] = 1;
while(tn<n-1&&tk<tn){
a[tn] = f[0][0],b[tn] = ((3ll*tk-tn)*inv[tn-tk]%p*f[0][0]+f[1][0])%p;
right(),right();
down();
}
for(int i=tn-1;i<=n;++i){
a[i] = solve(i,2*i-m);
b[i] = solve(i+1,2*i-m)-a[i];
if(a[i]==0&&b[i]==0) break;
}
for(int i=(m>>1);i<n;++i)
if(2*i>=m)
ans += ((ll)a[i]*binom(2*(n-i-1),k-2*i+m-2)+(ll)b[i]*binom(2*(n-i)-1,k-2*i+m-1))%p;
if(k+m==2*n) ans += a[n]+b[n];
printf("%lld",(ans%p+p)%p);
return 0;
}
```