常系数齐次线性递推学习小记

适用范围：

若遇到形如此类的递推式：$f_i={\sum }_{j=1}^k{f}_{i-j}\ast a_j$。
我们称其为线性递推。遇到此类问题，我们可以用矩阵快速幂求第$n$项，复杂度是$O(k^3\log n)。$
这个复杂度在$k$比较大的时候很吃力，我们就有了一种优化方案：
也就是常系数齐次线性递推。

具体推导(不含定理证明)：

我们设其转移矩阵为$A$。
初始向量为：$ST=${$f_0,f_1......,f_{k-1}$}。
我们最终要求的就是：$ST\ast A^n$。
我们考虑把$A^n$表示成如下形式$$A^n=\sum _{i=0}^{k-1}{c}_i{A}^i$$
那么我们要求的就是${\sum }_{i=0}^{k-1}{c}_i\ast ST\ast A^i$。(乘法满足分配律)。
考虑到$ST\ast A^i=S{T}_i$。相当于初始向量转移$i$次。
那么也就是说我们只要求出$ci$，那么只要将初始向量对应乘上去就可以得到第 $n$项了。
这个怎么求呢？
我们可以把一个多项式表示成如下形式：$A^n=C(x)=P(x)Q(x)+R(x)$。(在这里，$x$是一个矩阵。)
若当$x=A$时，$Q(x)$是零向量，那么问题转化为：$C(x)\equiv R(x)(modQ(x))$。
我们只需要快速幂加多项式取模即可求出$C(x)$即$c_i$。
这样的$Q(x)$一定存在，但我好像不太会证。
结论就是：$Q_k=1,Q_i=-a_{k-i}(0\le i\le k-1)$。这里的$a_i$是递推系数。
暴力多项式取模复杂度是$O(k^2\log n)$
用$NTT$可以优化到$O(k\log k\log n)$。

例题 NOI2017 泳池：

题意：有一个长度为$N$，高度为$1001$的网格图。每一个格子有$q$的概率是安全的。
我们定义一个合法的子矩形为：矩形中格子都是安全的，且矩形下边界贴着网格图下边界。
求最大子矩形面积为$K$的概率。

数据范围：$n\le 10^9,k\le 10^3$

题解：
为了方便$DP$，我们用最大面积至多为$K$的概率减去面积至多为$K-1$的概率。
考虑逐列递推，设$fi$表示到第$i$列的答案。考虑添加一个最大子矩形，枚举矩形的长。
需要其他$DP$辅助转移。
设一个$g_{i,j},h_{i,j}$分别表示长度为$j$，高度$i$以下全为安全格，以上任意，但最大子矩形面积不超过$k$的概率，长度为$j$，高度$i$以下全为安全格，且第$i+1$行最右边的格子不安全的，最大子矩形面积不超过$k$的概率。
枚举子矩形在哪里被一个不合法格子隔开，可以得到转移。
初始化：$g_{k,1}=q^k\ast (1-q),g_{i,0}=h_{i,0}=1$。
$$g_{i,j}=\sum _{l=0}^j{h}_{i,l}\ast g_{i+1,j-l}$$ $$h_{i,j}=\sum _{l=0}^{j-1}{h}_{i,l}\ast g_{i+1,j-l-1}\ast (1-q)\ast q^i$$
此时我们就可以递推$f$了。 $f_i={\sum }_{j=1}^k{f}_{i-j-1}\ast g_{1,j}\ast (1-q)$。
我们设$a_i=g_{1,i-1}\ast (1-q)$。就有：$f_i={\sum }_{j=1}^{k+1}{f}_{i-j}\ast a_j$。
这是一个线性递推式，暴力多项式取模快速幂优化即可。
复杂度：$O(k^2\log n)$。

Code：

# include<cstdio>
# include<cstring>
# include<algorithm>
using namespace std;
const int N = 1e3 + 5;
const int mo = 998244353;
typedef long long ll;
int g[N][N],h[N][N],a[N],f[N],q[N];
int t[N << 1],c[N << 1],d[N << 1],tmp[N << 1];
inline int pow(int x,int p)
{
	int ret = 1;
	for (; p ; p >>= 1,x = (ll)x * x % mo)
	if (p & 1) ret = (ll)ret * x % mo;
	return ret;
}
inline void mul(int *a,int *b,int len)
{
	memset(tmp,0,sizeof(tmp));
	for (int i = 0 ; i <= len ; ++i)
		for (int j = 0 ; j <= len ; ++j) tmp[i + j] = (tmp[i + j] + (ll)a[i] * b[j] % mo) % mo;
	for (int i = 0 ; i <= (len << 1) ; ++i) a[i] = tmp[i];
}
inline void mod(int *a,int fir,int sec)
{
	for (int i = fir ; i > sec ; --i)
	{
		int x = a[i];
		for (int j = i - sec - 1 ; j <= i ; ++j) a[j] = (a[j] - (ll)x * t[j - i + sec + 1] % mo + mo) % mo;
	}
}
inline int solve(int n,int K)
{
	if (!K) return pow(mo + 1 - q[0],n);
	memset(g,0,sizeof(g)),memset(h,0,sizeof(h));
	g[K][1] = (ll)q[K - 1] * (mo + 1 - q[0]) % mo,g[K][0] = 1;
	for (int i = K - 1 ; i ; --i)
	{
		int lim = K / i; g[i][0] = h[i][0] = 1;
		for (int j = 1 ; j <= lim ; ++j)
			for (int k = 0 ; k < j ; ++k) h[i][j] = (h[i][j] + (ll)h[i][k] * g[i + 1][j - 1 - k] % mo * q[i - 1] % mo * (mo + 1 - q[0]) % mo) % mo;
		for (int j = 1 ; j <= lim ; ++j)
			for (int k = 0 ; k <= j ; ++k) g[i][j] = (g[i][j] + (ll)h[i][k] * g[i + 1][j - k] % mo) % mo;
	}
	for (int i = 1 ; i <= K + 1 ; ++i) a[i] = (ll)g[1][i - 1] * (mo + 1 - q[0]) % mo;
	f[0] = t[K + 1] = 1;
	for (int i = K ; ~i ; --i) t[i] = mo - a[K + 1 - i];
	memset(c,0,sizeof(c)),memset(d,0,sizeof(d));
	int p = n; c[d[0] = 1] = 1;
	for (; p ; p >>= 1,mul(c,c,K),mod(c,K << 1,K))
	if (p & 1) mul(d,c,K),mod(d,K << 1,K);
	for (int i = 1 ; i <= K ; ++i)
	{
		f[i] = g[1][i];
		for (int j = 1 ; j <= i ; ++j) f[i] = (f[i] + (ll)f[i - j] * a[j] % mo) % mo;
	}
	int ret = 0; for (int i = 0 ; i <= K ; ++i) ret = (ret + (ll)f[i] * d[i] % mo) % mo;
	return ret;
}
int main()
{
	int n,K,x,y; scanf("%d%d%d%d",&n,&K,&x,&y); q[0] = (ll)x * pow(y,mo - 2) % mo;
	for (int i = 1 ; i <= K - 1 ; ++i) q[i] = (ll)q[i - 1] * q[0] % mo;
	printf("%d\n",(solve(n,K) - solve(n,K - 1) + mo) % mo);
	return 0;
}