一年多前初看《挑战程序设计竞赛》，看到如下这样的话
事实上，要求m项递推式的第n项的值可以不使用矩阵，而是使用初项的线性表示，通过快速幂在O(m^2 log n)的时间内求出答案。有兴趣的读者可以试着思考看看。
当时就百度了一下，发现各种解法都是首先要学一堆定理，然后还要学习各种多项式的知识。于是搁置了很久。
接近一年之后，又因为某些原因尝试再次学习这个递推，发现找到了一份很短的代码。琢磨了一下，发现可以从另一个角度去理解这个
k
2
∗
l
o
g
(
n
)
k^2*log(n)
k2∗log(n)的递推。
我们分两步走：
第一步
现有长度是4的递推式
a
n
=
x
1
a
n
−
1
+
x
2
a
n
−
2
+
x
3
a
n
−
3
+
x
4
a
n
−
4
a_n=x_1 a_{n-1}+x_2 a_{n-2}+x_3 a_{n-3}+x_4 a_{n-4}
an=x1an−1+x2an−2+x3an−3+x4an−4
我们这样存入rec数组里面
int rec[]={x1,x2,x3,x4}
类似矩阵快速幂，我们倍增地去考虑计算一个递推式。
比如我们已经计算有了
a
n
=
p
0
a
0
+
p
1
a
1
+
p
2
a
2
+
p
3
a
3
a_n=p_0 a_0+p_1 a_1+p_2 a_2+p_3 a_3
an=p0a0+p1a1+p2a2+p3a3
注意到相当于我们计算出了
a
2
n
=
p
0
a
n
+
p
1
a
n
+
1
+
p
2
a
n
+
2
+
p
3
a
n
+
3
a_{2n}=p_0 a_{n}+p_1 a_{n+1}+p_2 a_{n+2}+p_3 a_{n+3}
a2n=p0an+p1an+1+p2an+2+p3an+3
回代式子，得到：
a
2
n
=
p
0
(
p
0
a
0
+
p
1
a
1
+
p
2
a
2
+
p
3
a
3
)
+
p
1
a
n
+
2
+
p
2
a
n
+
3
+
p
3
a
n
+
4
a_{2n}=p_0 (p_0 a_0+p_1 a_1+p_2 a_2+p_3 a_3)+p_1 a_{n+2}+p_2 a_{n+3}+p_3 a_{n+4}
a2n=p0(p0a0+p1a1+p2a2+p3a3)+p1an+2+p2an+3+p3an+4
回代所有的项后，得到：
a
2
n
=
p
0
p
0
a
0
+
(
p
1
p
0
+
p
0
p
1
)
a
1
+
.
.
.
+
(
p
3
p
2
+
p
2
p
3
)
a
5
+
p
3
p
3
a
6
a_{2n}=p_0p_0 a_0+(p_1p_0+p_0p_1) a_1+...+(p_3 p_2+p_2 p_3)a_5+p_3p_3a_6
a2n=p0p0a0+(p1p0+p0p1)a1+...+(p3p2+p2p3)a5+p3p3a6
其系数就是多项式
(
p
0
+
p
1
+
p
2
+
p
3
)
(
p
0
+
p
1
+
p
2
+
p
3
)
(p_0+p_1+p_2+p_3)(p_0+p_1+p_2+p_3)
(p0+p1+p2+p3)(p0+p1+p2+p3)
展开之后第i项的系数
第二步
这个时候我们完成了第一步，计算出
a
2
n
a_{2n}
a2n，但是有个问题，这个递推式的长度也倍增了，看起来似乎它会和我们的项数一样长。这个时候原来的递推式上场。
因为
a
n
=
x
1
a
n
−
1
+
x
2
a
n
−
2
+
x
3
a
n
−
3
+
x
4
a
n
−
4
a_n=x_1 a_{n-1}+x_2 a_{n-2}+x_3 a_{n-3}+x_4 a_{n-4}
an=x1an−1+x2an−2+x3an−3+x4an−4
所以
a
6
=
p
3
p
3
(
x
1
a
5
+
x
2
a
4
+
x
3
a
3
+
x
4
a
2
)
a_6=p_3p_3(x_1 a_5+x_2 a_4+x_3 a_3+x_4 a_2)
a6=p3p3(x1a5+x2a4+x3a3+x4a2)
然后我们干掉了
a
6
a_6
a6,递推式长度减少了1
然后我们再继续对
a
5
a_5
a5,
a
4
a_4
a4做同样的操作，我们就得到了
a
2
n
a_{2n}
a2n关于
a
0
,
a
1
,
a
2
,
a
3
a_0,a_1,a_2,a_3
a0,a1,a2,a3的递推式了。
上代码。
#define rep(i,a,b) for(int i=a,_b=b;i<_b;i++)
typedef long long ll;
void mul(vector<ll>&rec,ll a[],ll b[],int k){
ll c[2*k]={};
rep(i,0,k) rep(j,0,k)
c[i+j]=(c[i+j]+a[i]*b[j])%mod2;
for(int i=k*2-2;i>=k;i--)
for(int j=1;j<=k;j++)//use recursion to go back
c[i-j]=(c[i-j]+rec[j]*c[i])%mod2;
rep(i,0,k) a[i]=c[i];
}
// recursion initial value nth item
ll linear(vector<ll>&a,ll b[],ll n){
int k=a.size()-1;ll res[2*k]={},c[2*k]={};
c[1]=res[0]=1;
for(;n;n/=2,mul(a,c,c,k))
if(n&1) mul(a,res,c,k);
ll ret=0;
rep(i,0,k) ret=(ret+b[i]*res[i])%mod2;
return ret;
}
linear函数的意思是，给定a这个递推式，给定初值b[]，求第n项。
mul函数就是给定递推式rec,a[]表示第x项用初项表示的系数,b[]表示第y项用初项表示的系数，现计算x+y项，结果存入a[],k是递推式的长度。
综上，我们实现了k^2log(n)计算递推式。
扩展：
今年年初我写了篇文章：
Number of Battlefields UVA - 11885 矩阵快速幂（高斯消元解线性递推）
原文章中提到了BM算法，其作用是在递推式长度的平方时间复杂度内计算出递推式。
后面我利用了矩阵快速幂来求解这个递推式的递推。
以下设递推式的长度为
k
k
k
所以整体时间复杂度是高斯消元的
k
3
k^3
k3+ 矩阵快速幂的
k
3
l
o
g
(
n
)
k^3log(n)
k3log(n)
如果我们即使用了BM算法在
k
2
k^2
k2时间复杂度内求出递推式，瓶颈仍然在快速幂。
但现在如果递推式优化到了
k
2
l
o
g
(
n
)
k^2log(n)
k2log(n),那么BM突然间又有用起来了。
当然知道往后学仍然有 k l o g ( k ) l o g ( n ) klog(k)log(n) klog(k)log(n)求递推式的方法，但似乎需要大量多项式的知识。
需要快速求解递推式的题还挺少。也许因为这个科技还比较少知道吧。
这个是求递推式的，但是递推式挺短的：
Covering
或者和上篇文章一样：
Number of Battlefields UVA - 11885
参考代码：代码
后来看到cf上TLE写了篇BM算法的文章，其中也提到了类似的递推式计算。
添加链接描述
update 5：
今天(6.9)重新去oi wiki上看了 OI Wiki常系数齐次线性递推的内容，发现看懂了。其实本文的内容大概从“尽量避免解释多项式”的角度去解释了一个多项式乘法和多项式取模的过程，其内容和OIwiki上本质是一样的。
347
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
