矩阵运算快速幂来快速计算线性递推式

92 阅读 0 评论 61 点赞

我是靠谱客的博主粗犷导师，这篇文章主要介绍矩阵运算快速幂来快速计算线性递推式，现在分享给大家，希望可以做个参考。

斐波那契数列
f(0)=0; f(1)=1; f(n)=f(n-1)+f(n-2),n>1

从上面这个方程中我们可以看到很明显的递推关系，当n>1的时候很明显发现会有一个关系式，但是实际上我们在做运算的时候，如果一步一步的按照递推式计算，将会消耗大量的时间（最短也是O(n)的时间复杂度），于是我们这个时候就需要引入矩阵乘法和快速幂来减少时间复杂度

矩阵乘法：
设A为mp的矩阵，B为pn的矩阵，那么称m*n的矩阵C为矩阵A与B的乘积，其中矩阵C中的第i行第j列元素可以表示为A的第i行与B的第j列对应元素乘积和

矩阵相乘：矩阵相乘最重要的方法是一般矩阵乘积。它只有在第一个矩阵的列数（column）和第二个矩阵的行数（row）相同时才有意义。一般单指矩阵乘积时，指的便是一般矩阵乘积。一个m×n的矩阵就是m×n个数排成m行n列的一个数阵。由于它把许多数据紧凑的集中到了一起，所以有时候可以简便地表示一些复杂的模型。
设A为mp的矩阵，B为pn的矩阵，那么称m*n的矩阵C为矩阵A与B的乘积，记作C=AB：

首先：快速幂取模模板
例如：求x的n次幂并取模

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
	int quickmod(long long n,long long x,long long r) //整型的
	{
		int ans=1;
		while(n)
		{
			if(n&1)
				ans=ans*x%r;
			n>>=1;
			x=x*x%r;
		}
		return ans;
	}

矩阵 (struct)来定义

复制代码

1
2
3
4
struct   Matrix{
	int m[2][2];
}

矩阵的乘法运算

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
    Matrix mul(Matrix k,Matrix n)
    {
        matrix ans;
        for(int i=0;i<2;i++)
            for(int j=0;j<2;j++)
            {
                ans.m[i][j]=0;
                for(int k=0;k<2;k++)
                //ans.m[i][j]+=A.m[i][k]*B.m[k][j]%mod;
        		//修正
        		ans.m[i][j] +=A.m[i][k]*B.m[k][k];
        		ans.m[i][j] = ams.m[i][j]%mod;
            }
        return ans;
    }

两者结合一下即可
最终代码

复制代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include <iostream>
using namespace std;
const long long mod=1e9+7;
struct  Matrix{
	int m[2][2];
};
Matrix mul(Matrix A,Matrix B)
{
    Matrix ans;
    for(int i=0;i<2;i++)
        for(int j=0;j<2;j++)
        {
            ans.m[i][j]=0;
            for(int k=0;k<2;k++)
                //ans.m[i][j]+=A.m[i][k]*B.m[k][j]%mod;
        		//修正
        		ans.m[i][j] +=A.m[i][k]*B.m[k][k];
        		ans.m[i][j] = ams.m[i][j]%mod;
        }
    return ans;
}
Matrix quickmod(long long n,Matrix x)
{
	Matrix ans;
	ans.m[0][0]=1;
	ans.m[0][1]=0;
	ans.m[1][0]=0;
	ans.m[1][1]=1;
	while(n)
	{
		if(n&1)
			ans=mul(ans,x);
		n>>=1;
		x=mul(x,x);
	}
	return ans;
}
int main()
{
	 Matrix A,ans;
	 int n;
	 cin>>n;
	 ans.m[0][0]=1;
	 ans.m[0][1]=0;
  	 A.m[0][0]=1;
 	 A.m[0][1]=1;
	 A.m[1][0]=1;
	 A.m[1][1]=0;
	 ans=mul(ans,quickmod(n-1,A));
	 cout<<ans.m[0][0]<<endl;
	 return 0;
}

然后如何将线性递推式转化成矩阵的形式来求解呢
先看一道hdu的题A Simple Math Problem hdu 1757

Problem Description
Lele now is thinking about a simple function f(x).
If x < 10 f(x) = x.
If x >= 10 f(x) = a0 * f(x-1) + a1 * f(x-2) + a2 * f(x-3) + …… + a9 * f(x-10);
And ai(0<=i<=9) can only be 0 or 1 .
Now, I will give a0 ~ a9 and two positive integers k and m ,and could you help Lele to caculate f(k)%m.
Input
The problem contains mutiple test cases.Please process to the end of file.
In each case, there will be two lines.
In the first line , there are two positive integers k and m. ( k<2*10^9 , m < 10^5 )
In the second line , there are ten integers represent a0 ~ a9.
Output
For each case, output f(k) % m in one line.
Sample Input
10 9999
1 1 1 1 1 1 1 1 1 1
20 500
1 0 1 0 1 0 1 0 1 0
Sample Output
45
104

这里当x<10的时候 f(x) = x;
当x>=10的时候 f(x)=a0 * f(x-1) + a1 * f(x-2) + a2 * f(x-3) + …… + a9 * f(x-10)这也就是一个线性递推式了
我们就可以把他变成矩阵的形式来求解。

复制代码

#include <iostream>
#include <cstring>
using namespace std;
typedef long long ll;
const int MAX =10; //最多也就一个10*10的矩阵乘法
struct matrix
{
    ll m[MAX][MAX];
};
matrix mul(matrix A,matrix B,ll mo)
{
    matrix ans;
    for(int i=0;i<MAX;i++){
        for(int j=0;j<MAX;j++){
            ans.m[i][j] = 0;
            for(int k=0;k<MAX;k++){
                ans.m[i][j] += A.m[i][k]*B.m[k][j];
                ans.m[i][j] = ans.m[i][j]%mo;
            }
        }
    }
    return ans;
}
matrix quickmod(matrix A,ll k,ll mo)
{
    matrix ans;
    memset(ans.m,0,sizeof(ans.m));
    for(int i=0;i<MAX;i++){
        ans.m[i][i] = 1;
    }
    while(k){
        if(k&1)
            ans = mul(ans,A,mo);
        k>>=1;
        A = mul(A,A,mo);
    }
    return ans;
}
int main()
{
    ll k,m;
    matrix Init;
    for(int i=0;i<MAX;i++){
        Init.m[i][0] = MAX-i-1;
    }
    while(cin>>k>>m){
        matrix T;
        if(k<10){
            cout<<k<<endl;
            continue;
        }
        memset(T.m,0,sizeof(T));
        for(int i=0;i<MAX;i++){cin>>T.m[0][i];} //记录下a0,a1,a2,a3,a4……
        for(int i=1;i<MAX;i++){ // 初始化后面的根据矩阵
            T.m[i][i-1] = 1;
        }

cout << mul(quickmod(T,k-9,m),Init,m).m[0][0] <<endl;
    }

return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include <iostream>
#include <cstring>
using namespace std;
typedef long long ll;
const int MAX =10; //最多也就一个10*10的矩阵乘法
struct matrix
{
    ll m[MAX][MAX];
};
matrix mul(matrix A,matrix B,ll mo)
{
    matrix ans;
    for(int i=0;i<MAX;i++){
        for(int j=0;j<MAX;j++){
            ans.m[i][j] = 0;
            for(int k=0;k<MAX;k++){
                ans.m[i][j] += A.m[i][k]*B.m[k][j];
                ans.m[i][j] = ans.m[i][j]%mo;
            }
        }
    }
    return ans;
}
matrix quickmod(matrix A,ll k,ll mo)
{
    matrix ans;
    memset(ans.m,0,sizeof(ans.m));
    for(int i=0;i<MAX;i++){
        ans.m[i][i] = 1;
    }
    while(k){
        if(k&1)
            ans = mul(ans,A,mo);
        k>>=1;
        A = mul(A,A,mo);
    }
    return ans;
}
int main()
{
    ll k,m;
    matrix Init;
    for(int i=0;i<MAX;i++){
        Init.m[i][0] = MAX-i-1;
    }
    while(cin>>k>>m){
        matrix T;
        if(k<10){
            cout<<k<<endl;
            continue;
        }
        memset(T.m,0,sizeof(T));
        for(int i=0;i<MAX;i++){cin>>T.m[0][i];} //记录下a0,a1,a2,a3,a4……
        for(int i=1;i<MAX;i++){ // 初始化后面的根据矩阵
            T.m[i][i-1] = 1;
        }

        cout << mul(quickmod(T,k-9,m),Init,m).m[0][0] <<endl;
    }

    return 0;
}