矩阵快速幂

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/u013377068/article/details/79960366

在讲解矩阵快速幂之前我先来讲解一下快速幂

主要是为了了解快速幂的思想。


快速幂:

————————————————————————————————————————————

快速幂解决的是求 n ^ k 的值的问题;

需要注意的是这里的k非常大,往往是1e7之上的。(这里假设最后求得的值不会溢出);

那么用朴素的求法来求 n^ k的过程是:

n * n * n * n*n ..... *n;

这样的求法往往会超时。

那么我们就要换一种求法了。

我们来个例子:

在现实生活中,如果让你求2^20次方,你会怎么求?

一般的求法是 1. 让 2 * 2 = 4,这样就能算出2^2的值了。

                     2. 让 4 * 4 =16,这样就能算出2^4的值了。

                     3. 让16 * 16 = 256,这样就能算出2 ^ 8的值了

                     4. 让256 * 256 = 65536,这样就能算出2 ^ 16的值了

                     5. 让65536 * (2 ^ 4) = 1048576,这样就能算出2 ^ 20的值了。


这样我们只需要5步就能算出2^20的值了,如果是按朴素的求法需要20步。

而且随着k的增大,他们的差距也会越来越大。(O(n) 和 O( log2(n) )的差距。 

但是第二种求法的话,我们需计算出最后一步的a * b中的b (也就是举例中的65536 * (2 ^ 4) 中的(2 ^ 4)) ,

这样一来,当k比较的话,最后一步也很麻烦。

在这里我们采用一种和第二种方法思想上类似的方法(没错,就是二分) 

我们以a ^ 19来举例:如图


我们从最后一层开始计算:

    1.我们可以发现最后一层从左边开始a 和 a 是一样的(听上去好像是废话),我们设res为上一层的第一部分

        的值。显然此时res = a *  a

    2.倒数第二层的第一部分和第二部分是一样的,所以res = res * res。

    3.倒数第三层的第一部分和第二部分是一样的,所以res = res * res。

    4.倒数第四层的第一部分和第二部分是一样的,所以res = res * res。到了这一步我们就求处倒数第五层(

        也就是正数第一层)的res。如果按照一般情况,这时候res = a ^19, 但是我为让每一层的第一部分

        和第二部分相等,我把每一层a的数量当成偶数来分了,所以如图,我们还有两个a没乘上。我们可以在

        这个过程中用ans来存这部分忽略的值,所以a ^ 19 = res * ans;

快速幂就是这样一个过程;

代码如下:

  

#include<stdio.h>
int FastPower(int a , int k);	//这里假设a ^ k不会溢出 
int main()
{
	int a ,k;
	
	scanf("%d %d",&a, &k);
	
	int result = FastPower( a,  k);
	
	printf("%d\n", result);
	
	
	return 0;
}
int FastPower(int a , int k)
{
	int res = a ; //最初res是等于a的。
	 
	int ans = 1;   
	
	while(k)	//当k ==0 时,结束计算 
	{
		if(k & 1) //相当于 k % 2 . 是用来判断奇偶的,用位运算会快点 
		ans = ans * res;
		
		res = res * res;
		
		k = k >> 1;		//相当于k = k /2,用位运算会快点 
	} 
	
	return ans;
} 

———————————————————————————————————————————————


接下来就是矩阵快速幂了:

———————————————————————————————————————————————

矩阵快速幂基本上和快速幂是一样的,所以这里我用快速幂来类比矩阵快速幂

它们有3点不同:

1.快速幂中的基本元素是一个正实数,矩阵快速幂的基本元素是一个矩阵。

2.快速幂中的基本操作是乘法,矩阵快速幂的基本操作是矩阵乘法。

3.快速幂的结果是一个实数,矩阵快速幂的结果是一个矩阵。


为了更好地讲解矩阵快速幂,我们先了解矩阵和矩阵乘法:

矩阵的定义:

由 m × n 个数aij排成的m行n列的数表称为m行n列的矩阵,简称m × n矩阵。记作:


矩阵乘法:


简单来说就是:乘积C的第m行第n列的元素等于矩阵A的第m行的元素与矩阵B的第n列对应元素乘积之和。

需要注意的是:矩阵相乘的条件是:矩阵A的行数等于矩阵B的列数

                                                      矩阵A的列数等于矩阵B的行数


代码如下:

#include<stdio.h>
struct Mat
{
	int n,m;
	int data[10][10];
};

Mat MatMul(Mat a, Mat b);		//矩阵乘法 

Mat ans; 		 

Mat MatFastPower(Mat a, int k);


int main()
{
	
	int n,m;
	
	scanf("%d %d",&n, &m);		//输入矩阵的长宽 
	Mat a;
	
	a.n =n, a.m =m;
	
	for(int i = 1; i <= n; i++)		
	{
		for(int j =1; j <= m; j++)
		scanf("%d",&a.data[i][j]);
	}
	
	ans.n =m, ans.m =n;
	
	for(int i =1 ;i <= n; i++)		//使得ans="1"; 
	{
		for(int j =1; j <= m; j++)
		{
			if( i== j)
			ans.data[i][j] = 1;
		}
	}
	
	int k;
	
	scanf("%d", &k);			//输入幂次 
	

	Mat result = MatFastPower( a,  k);		//矩阵快速幂 
	
	for(int i = 1; i <= result.n; i++)		//输出矩阵a的k次方 
	{
		for(int j = 1; j <= result.m; j++)
		printf("%d ",result.data[i][j]) ;
		
		printf("%\n");
	} 
	
	return 0;
}

Mat MatFastPower(Mat a, int k)
{
	Mat res = a;
	
	while(k)
	{
		if(k & 1)
		ans = MatMul(ans,res);
		
		res = MatMul(res, res);
		k = k>> 1;
	}
	return ans;
	
}


Mat MatMul(Mat a, Mat b)		//矩阵乘法 
{
	Mat temp;
	
	temp.n = a.n;
	
	temp.m = b.m;
	
	for(int i = 1; i <= temp.n; i++)
	{
		
		for(int j = 1; j <= temp.m; j++)
		{
			temp.data[i][j]=0;
			
			for(int k =1 ; k <= a.m; k++)
			
			temp.data[i][j] += (a.data[i][k] * b.data[k][j]); 
			
		}
	}
	return temp;
	
}

可以看出除了乘法的过程不一样,其他基本一样。

———————————————————————————————————————————————


例题

———————————————————————————————————————————————

51Nod - 1242    

点击打开链接



斐波那契数列的定义如下:


F(0) = 0
F(1) = 1
F(n) = F(n - 1) + F(n - 2) (n >= 2)


(1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, ...)
给出n,求F(n),由于结果很大,输出F(n) % 1000000009的结果即可。
Input
输入1个数n(1 <= n <= 10^18)。
Output
输出F(n) % 1000000009的结果。
Sample Input
11
Sample Output

89


可以看出n是一个很大的值,我们也不能求出F(n);

我们可以把F(n)转换成a * b ^ n 的形式;

根据递推式F(n) = F(n-1) + F(n-2)

我们可以把它转化成F(n) = F(n-1) * 1 + F(n-2)  * 1;

以下是我推导的过程


这样我们就得到F(n)的一个公式,用矩阵快速幂来做,求出F(N)的复杂度log2(N) (不包括矩阵相乘的过程)

这道题的代码如下

#include<stdio.h>
#define ll long long
const ll Mod =1000000009;
struct Mat
{
	int n, m;
	
	ll data[4][4];
};

Mat ans;

Mat MatMul(Mat a, Mat b);

Mat MatMod(Mat a, ll mod);		//矩阵取模 

Mat MatFastPoswer(Mat a, ll k);

int main()
{
	
	ll n;
	scanf("%lld", &n);
	
	Mat a;			//构建矩阵{1,1,1,0} ;
	a.n = 2, a.m = 2; 
	a.data[1][1] = 1, a.data[1][2]=1, a.data[2][1]=1, a.data[2][2]=0;
	
	Mat result = MatFastPoswer(a, n-1);
	
	ll res = 1 * result.data[1][1];
	
	printf("%lld\n",res % Mod);
	
	
	
	
	return 0;
}

Mat MatFastPoswer(Mat a, ll k)
{
	ans.m=2, ans.n = 2;
	
	ans.data[1][1] = 1, ans.data[1][2] = 0;
	
	ans.data[2][1] = 0, ans.data[2][2] = 1;
	
	Mat res = a;
	
	while(k)
	{
		if(k & 1)
		ans = MatMod(MatMul(ans,res), Mod);
		
		res = MatMod(MatMul(res, res), Mod);
		
		k = k >> 1;
	}
	
	
	
	return ans;
	
}

Mat MatMul(Mat a, Mat b)
{
	Mat temp ;
	
	temp.n = a.n, temp.m = b.m;
	
	for(int i = 1; i <= temp.n; i++)
	{
		for(int j =1; j <= temp.m; j++)
		{
			temp.data[i][j] = 0;
			
			for(int k = 1; k <= temp.m; k++)
			{
				temp.data[i][j] += (a.data[i][k] * b.data[k][j]);
			}
		}
	}
	return temp;
}
Mat MatMod(Mat a, ll mod)
{
	Mat temp = a; 
	
	for(int i = 1; i <= temp.n; i++)
	{
		for(int j = 1; j <= temp.m; j++)
		{
			temp.data[i][j] %= mod;
		}
	}
	
	return temp;
}

———————————————————————————————————————————————

阅读更多

扫码向博主提问

糖宋元明清

非学,无以致疑;非问,无以广识
去开通我的Chat快问
想对作者说点什么?

博主推荐

换一批

没有更多推荐了,返回首页