矩阵快速幂求递推式

75 阅读 0 评论 50 点赞

我是靠谱客的博主阔达音响，最近开发中收集的这篇文章主要介绍矩阵快速幂求递推式，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

用矩阵将快速幂可以以logn级别的时间复杂度求出递推式

典型题：求斐波那契数列第n项 n<=2^31-1 ；显然，一步一步递推o（n）算法效率不够

0 1 0　　　　f[n-3]　　f[n-2]

0 0 1　　*　 f[n-2] =　f[n-1]

1 1 3　　　　f[n-1]　　f[n-3]+f[n-2]+3*f[n-1]

在方阵右上方 n-1 * n-1 矩阵的主对角线上填1 最下面一行从右到左填递推式的系数

用矩阵实现递推式

由矩阵的基本性质知 A^n *P =A^（n-1） *A*P

所以求递推式第n项可以通过* 转移矩阵 n次实现

那么问题来了 A^n 如何快速求？

可以用二分法递归求，也可以用二进制累乘法

下面介绍二进制快速幂:

2^11=2^1 * 2^2 * 2^8

也就是按照11的二进制1011 加权算

代码如下：

while(p)
        {
            if(p&1) ans=multi(ans,b); //p是要求的幂次   矩阵ans 初始状态是单位矩阵（相当于数1）
            b=multi(b,b);　　　　　　　　//矩阵b是转移矩阵
            p>>=1;　　　　　　　　　　　　//右移迭代
        }
//  A^n    最终结果就是矩阵     ans

对于转移矩阵 * 递推首项可以把递推首项设置为方阵第一列有值，其余列为0

这样只用方阵乘法就够了

matrix multi(matrix a,matrix b)
{
    matrix ans;
    memset(ans.a,0,sizeof(ans.a));
    for(int i=0;i<maxn;i++)
        for(int j=0;j<maxn;j++)
            for(int k=0;k<maxn;k++){
                ans.a[i][j]+=a.a[i][k]*b.a[k][j];
                ans.a[i][j]%=m;
            }
    return ans;
}

以下是矩阵快速幂求递推式的模板代码 HDU1757

 1 #include <bits/stdc++.h>
 2 using namespace std;
 3 const int maxn =10;
 4 //矩阵结构体
 5 struct matrix
 6 {
 7     int a[maxn][maxn];
 8 };
 9 matrix a,b;
10 int a0[15],k,m;
11 void init()
12 {
13     memset(b.a,0,sizeof(b.a));
14     for(int i=0;i<9;i++)
15         b.a[i][i+1]=1;
16     for(int i=0;i<10;i++)
17         b.a[9][i]=a0[9-i];
18     for(int i=0;i<10;i++)
19         a.a[i][0]=i;
20 }
21 matrix multi(matrix a,matrix b)
22 {
23     matrix ans;
24     memset(ans.a,0,sizeof(ans.a));
25     for(int i=0;i<maxn;i++)
26         for(int j=0;j<maxn;j++)
27             for(int k=0;k<maxn;k++){
28                 ans.a[i][j]+=a.a[i][k]*b.a[k][j];
29                 ans.a[i][j]%=m;
30             }
31     return ans;
32 }
33 int main()
34 {
35     while(scanf("%d%d",&k,&m)!=EOF)
36     {
37         for(int i=0;i<10;i++)
38             scanf("%d",&a0[i]);
39         init();
40         int p=k-9;
41         matrix ans;
42         memset(ans.a,0,sizeof(ans.a));
43         for(int i=0;i<10;i++)
44             ans.a[i][i]=1;
45         while(p)
46         {
47             if(p&1) ans=multi(ans,b);
48             b=multi(b,b);
49             p>>=1;
50         }
51         ans=multi(ans,a);
52         cout<<ans.a[9][0]<<endl;
53     }
54     return 0;
55 }