矩阵快速幂求线性常系数递推矩阵快速幂求线性常系数递推

75 阅读 0 评论 50 点赞

我是靠谱客的博主生动毛衣，最近开发中收集的这篇文章主要介绍矩阵快速幂求线性常系数递推矩阵快速幂求线性常系数递推，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

矩阵快速幂求线性常系数递推

这个东西十分好用，但是搞懂这个需要先知道两个要点

——矩阵乘法

——常系数递推

矩阵乘法

我们有两个矩阵

X = ⎡ ⎣ ⎢ 213101343 ⎤ ⎦ ⎥ Y = ⎡ ⎣ ⎢ 112 1 - 1 0 ⎤ ⎦ ⎥

$begin {matrix} X=begin{bmatrix} 2 & 1 & 3\ 1 & 0 & 4\ 3 & 1 & 3 end{bmatrix} & Y = begin{bmatrix} 1 & 1 \ 1 & -1 \ 2 & 0 end{bmatrix} end{matrix}$
然后我们定义一个新矩阵 Z ，

Z=X×Y Z = X × Y $Z=X times Y$

其中X为a行b列，Y必须为b行c列。

那么对于矩阵中每一个元素 $Z_{i ;j} = sum {_{k=1}^b} {;} X_{i;k}times Y_{k;j}$

而Z本身就是一个a*c的矩阵。

Z = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 2 * 1 + 1 * 1 + 3 * 2 1 * 1 + 0 * 1 + 4 * 1 3 * 1 + 1 * 1 + 3 * 2 2 * 1 + 1 * (- 1) + 3 * 0 1 * 1 + 0 * (- 1) + 4 * 0 3 * 1 + 1 * (- 1) + 3 * 0 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$begin{matrix} Z=begin{bmatrix} 2*1+1*1+3*2 qquad & 2*1+1*(-1)+3*0\ \ 1*1+0*1+4*1 qquad & 1*1+0*(-1)+4*0\ \ 3*1+1*1+3*2 qquad & 3*1+1*(-1)+3*0\ end{bmatrix} end{matrix}$
于是我们就做完了一次矩阵乘法。

目前貌似算矩阵乘法只能枚举i，j，k，时间复杂度为 $O({n^3})$ ，是很高的，但是解决常系数递推问题是足够的。

线性常系数递推

设 $f(n)=a_1f(n-1)+a_2f(n-2)+a_3f(n-3)+cdots+a_nf(1)+k$

则这样的可以由这样的式子转移过来的f(n)，我们就称为它满足线性常系数递推关系。

这种式子很常见，比如斐波那契数列就是这样的式子。

斐波那契数列： $f(n)=f(n-1)+f(n-2)$

如果给出一个题，求斐波那契的第N项，那么我们非常自然的想到 $O(N)$ 递推过去。但如果数据范围 $N<=10^{18}$ ,显然 $O(N)$ 的方法是不行的。此时，利用矩阵乘法，我们可以做到 $O{(log_2(N))}$ 。

如何解决？

我们考虑将斐波那契数列转为1个1*2的矩阵。

我们命名为fib

F i b = [f [i - 1] f [i - 2]]

$begin{matrix} Fib=begin{bmatrix} quad f[i-1] quad \ quad f[i-2] quad end{bmatrix} end{matrix}$

然后我们再定义一个2*2的矩阵

T r a n s f e r = [1110]

$begin{matrix} Transfer=begin{bmatrix} 1 & 1\ 1 & 0\ end{bmatrix} end{matrix}$
然后我们定义一个ans矩阵

对fib和transfer做一下矩阵乘法，就会发现

a n s = [f [i] f [i - 1]]

$begin {matrix} ans = begin {bmatrix} f[i]\ f[i-1] end {bmatrix} end {matrix}$
于是

ans11 a n s 1 1 $ans_{ 1;1 }$ 就是答案。

有人可能会说，这比线性递推慢很多啊，这样算一次得进行4次运算，但是线性递推1次就能出结果

那么我们再来想，如果要再求一遍，是不是只需要让ans再乘一遍transfer就可以了，那么我们求n项，就乘n次transfer，根据乘法结合律，我们可以先把n个transfer相乘算出来，这个时候用快速幂，就能优化到 $O(log_2(N))$ 。这就是用处。

那么我们先来看一段整数快速幂运算

int ksm(int x,int y){//求x的y次方
    int res=1;
    while(y){
        if(y&1){
            res*=x;
        }
        x*=x;
        y>>=1;
    }
    return res;
}

我们可以看到，返回的答案res，在一开始的是时候赋值为1，因为任何数乘以1都等于它本身，这样就不会修改值了。

那么矩阵快速幂，返回的答案为ans，它在一开始的时候怎么赋值呢？显然不能全是1，参考上面的矩乘定义。

其实这里是有通法的，我们定义一个矩阵unit，如果任意矩阵乘以unit都等于它本身，则我们将unit称之为单位矩阵

u n i t = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 1000 \dots 0 0100000100 0 \dots 0 0 \dots 0 0 \dots 0 1 \dots 0 0 \dots 0 00001 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$begin{matrix} unit=begin{bmatrix} 1&0&0&0cdots0&0\ 0&1&0&0cdots0&0\ 0&0&1&0cdots0&0\ 0&0&0&1cdots0&0\ cdots\ 0&0&0&0cdots0&1 end{bmatrix} end{matrix}$
简单来说，就是一个除主对角线是1以外都是0的矩阵。

至于为什么是这样大家可以自己手推一下。

因为如果有一个矩阵与这个矩阵相乘，肯定是单位矩阵的第一行乘以这个矩阵的一列，无论这个矩阵这一列是谁，都只会保留第一个数，因为第一行第二个往后的数都是0，这个矩阵的第一列第二个往后的数乘上它们都得0，对答案是没有影响的。接下来同理。

所以我们就得到了一个单位矩阵。

接下来就是代码部分：

#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cstring>
using namespace std;
const long long mod=1e9+7;
struct mat{
    int m[2][2];
}; 
mat zhuanyi={
    1,1,
    1,0
};
mat danwei={
    1,0,
    0,1
};
mat mat_mul(mat a,mat b){
    mat c;
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            c.m[i][j]=0;
            for(int k=0;k<2;k++){
                c.m[i][j]+=(a.m[i][k]*b.m[k][j])%mod;   
            }
            c.m[i][j]%=mod;
        }
    }
    return c;
}
mat mat_pow(int k){
    mat ans,p;
    ans=danwei;
    p=zhuanyi;
    while(k){
        if(k&1){
            ans=mat_mul(ans,p);
        }
        p=mat_mul(p,p);
        k>>=1;
    }
    return ans;
}
int main(){
    int n;
    cin>>n;
    mat ans=mat_pow(n-1);
    cout<<ans.m[0][0]<<endl;
    return 0;
}

我们将矩阵直接定义为结构体，这样在返回和调用的时候就不需要指针操作，非常便捷。

另外这一题由于递推式的项数只有2，所以矩乘的复杂度为8，是个常数级别，就可以直接忽略了。但是对于项数超多的递推式，就要斟酌一下线性还是矩乘了。

mat_pow是矩阵快速幂，可以看到，基本上和一般快速幂没有区别，只是乘法改为调用上面的mat_mul矩乘函数。

那么最后还有一个问题，就是如何推出一个线性常系数递推关系式的转移矩阵？

也是有通法的。对于一个线性常系数齐次递推 $f(n)=a_1times f(n-1)+a_2times f(n-2)+a_3times f(n-3)+cdots+a_ktimes f(n-k)$ ，有如下转移矩阵:

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ a 1 100 \dots 0 a 2 0100 a 3 0010 a 4 \dots a k - 1 0 \dots 0 0 \dots 0 0 \dots 0 0 \dots 1 a k 0000 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$begin{matrix} begin{bmatrix} a_1&a_2&a_3&quad a_4 cdots a_{k-1}&a_k\ 1&0&0&0cdots ;0&0\ 0&1&0&0cdots ;0&0\ 0&0&1&0cdots ;0&0\ cdots\ 0&0&0&0cdots ;1&0 end{bmatrix} end{matrix}$
至于为何是这样，还是可以手玩一下自己算一算，用中文很难讲的比较明白。

而另一个数列也就显而易见了

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ f (k - 1) f (k - 2) f (k - 3) \dots f (1) f (0) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$begin{matrix} begin{bmatrix} f(k-1)\ f(k-2)\ f(k-3)\ cdots\ f(1)\ f(0)\ end{bmatrix} end{matrix}$

对于一个线性常系数非齐次递推

$f(n)=a_1times f(n-1)+a_2times f(n-2)+a_3times f(n-3)+cdots+a_ktimes f(n-k)+b$ ，有如下转移矩阵：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ a 1 100 \dots 00 a 2 01000 a 3 00100 a 4 \dots a k - 1 0 \dots 0 0 \dots 0 0 \dots 0 0 \dots 1 0 \dots 0 a k 00000 b 00001 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$begin{matrix} begin{bmatrix} a_1&a_2&a_3&quad a_4 cdots a_{k-1}&a_k&b\ 1&0&0&0cdots ;0&0&0\ 0&1&0&0cdots ;0&0&0\ 0&0&1&0cdots ;0&0&0\ cdots\ 0&0&0&0cdots ;1&0&0\ 0&0&0&0cdots ;0&0&1 end{bmatrix} end{matrix}$
另一个数列为