在讲解矩阵快速幂之前我先来讲解一下快速幂
主要是为了了解快速幂的思想。
快速幂：
————————————————————————————————————————————
快速幂解决的是求 n ^ k 的值的问题；
需要注意的是这里的k非常大，往往是1e7之上的。(这里假设最后求得的值不会溢出)；
那么用朴素的求法来求 n^ k的过程是：
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.快速幂的结果是一个实数，矩阵快速幂的结果是一个矩阵。
为了更好地讲解矩阵快速幂，我们先了解矩阵和矩阵乘法：
矩阵的定义：
矩阵乘法：
简单来说就是:乘积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;
}———————————————————————————————————————————————
京公网安备 11010502036488号