矩阵快速幂——入门

题目链接POJ307

题意:求第n位斐波那契数mod 10000的大小。其中n的大小高达1000000000

由于我对矩阵快速幂也不是很了解,所以只能浅谈,不敢深谈,希望对新手有点帮助。

先聊聊快速幂

首先,如果要你自己写一个pow()函数,需要求220,怎么求

#include <iostream>
using namespace std;
int quickPow(int a,int b){   //快速幂
    int ans = 1;
    while(b){
        if(b&1) ans*=a;
        a *= a;
        b /= 2;
    }
    return ans;
}
int pow(int a,int b){   //普通乘
    int ans = 1;
    for(int i=1; i<=b; i++)
        ans *= a;
    return ans;
}
int main(int argc, const char * argv[]) {
    cout<<quickPow(2, 20)<<endl;   //第一种方式
    cout<<pow(2,20)<<endl;    //第二种方式
    return 0;
}

上面的两种方式,一种复杂度是O(logn),一种复杂度是O(n)。当你数量级很大的时候,比如10000000000时,这个复杂度的时间差别就很大了

怎么得到这个快速幂的代码的呢?我们假设需要求274

74转化为二进制为1001010,根据二进制对应位的1,可以把74拆分为64+8+2 ,所以
74 = 26 + 23 + 21
274 = 226 + 223 + 221

 while(b){    //快速幂核心代码
      if(b&1) ans*=a;
      a *= a;
      b /= 2;
 }

下面看一下上面这段快速幂程序中的ans和a的用法 (74二进制1001010为从右到左下面表格从1-7)

初始值 0 1(ans变了) 0 1(ans变了) 0 0 1(ans变了)
a = 2n 21 22 24 28 216 232 264
ans 20 20 20×22 20×22 20×22×28 20×22×28 20×22×28 ×264

这里有两个变量,第一个变量a默认用来存2n,,变量ans存储的是最后的结果,当遇到第k位的二进制位是1的时候,则ans ×= 2k

举个例子,74的二进制表示是1001010,因为是除2运算,所以从右往左遍历,到右边第二个时,ans ×= 22,到右边第四个时,ans ×= 28,到右边第7个时,ans ×= 264,最后ans = 274,至于ans乘的过程中的那些22、28、264怎么来?a变量不是一直在不断增加吗,ans乘的就是a变量对应位的值,对应上面的表格可能看的更清楚一点

不知道,你看懂没……如果还是没懂,就给个面子装懂好嘛……

矩阵?怎么过渡到斐波那契

  • 宏观上看矩阵

上图中,乘号右边是一系列的矩阵。我们称为B1、B2、B3、B4 …… Bn
该序列中的同行为一个斐波那契数列。乘号左边是一个转移矩阵A,Bn = A×Bn-1

  • 微观上看矩阵


原始的序列是 [A,B],需要得到的序列是[A+B,A],根据矩阵相乘定律,可以得到转移矩阵A

  • 总结一波
    现在需要求的是Bn%10000,Bn = B1 × An-1,B1[1,1],所以 An-1就是需要的答案。那把求Bn改成了求An-1有什么好处呢?见下文

矩阵快速幂的实际运用

上述矩阵代替斐波那契数列的方式。已经把斐波那契数列从原先的递加变成了现在的递乘。那变成递乘有什么用?

如果斐波那契第一项是A,第二项是A2,第三项是A3,第4项是A4 ……
那第8项的值是,A8 = A4×A4

换句话说,第4项的值 × 第4项的值 = 第8项的值

求第16项的值
原来的方式是:1->2->3->4……->15->16需要经过15次运算
现在的方式是:1->2->4->8->16只需要经过4次运算

为什么只需要经过4次?从A->A2得到第2项,我要求第4项不需要经过第3项,直接A2 × A2 -> A4就可以得到第4项。这又回到了刚开始的快速幂

具体怎么得到,下面请看代码解析

代码拆分

代码不长,先看下面这个函数。这里使用了引用传递,即a数组发生变化后,传入的实参也会改变。例如传入数组a数组b,得到的数组a = 数组a × 数组b。即Matrix(cot,temp)执行的结果是cot *= temp,且cottemp都表示的是矩阵数组

void Matrix(int (&a)[2][2],int b[2][2]){   //矩阵乘法
    int tmp[2][2] = {0};
    for(int i = 0; i < 2; ++i)
        for(int j = 0; j < 2; ++j)
            for(int k = 0; k < 2; ++k)
                tmp[i][j] = (tmp[i][j] + a[i][k] * b[k][j]) % N;  //矩阵相乘,注意mod
    for(int i = 0; i < 2; ++i)
        for(int j = 0; j < 2; ++j)  
            a[i][j] = tmp[i][j]; //赋值返回数据
}

A1,也就是转移矩阵,题中已经给出了,用数组表示是[1,1,1,0]。这里还需要一个A0,A0其实就是单位矩阵E(E × 任意矩阵A = A)。单位矩阵用数组表示就是[1,0,0,1](对角线都是1的矩阵)

while(n){
       if(n & 1) Matrix(cot,temp);   //如果是奇数
       Matrix(temp,temp);   
       n /= 2;  //不断除2
}

上述5行代码是矩阵快速幂关键(至于怎么得到这5行代码的,在上述快速幂的解释中应该已经可以得到了)。cot数组一开始是A0temp数组一开始是A1。我们来模拟下求第13项的斐波那契值是怎么求的。

第一步:n = 13,执行Matrix(cot,temp)Matrix(temp,temp),得到cot数组为A1t数组为A2n = 6

第二步:n = 6,执行Matrix(temp,temp),得到cot数组为A1t数组为A4n = 3

第三步:n = 3,执行Matrix(cot,temp)Matrix(temp,temp),得到cot数组为A5t数组为A8n = 1

第四步:n = 1,执行Matrix(cot,temp)Matrix(temp,temp),得到cot数组为A13t数组为A16n = 0

退出while循环并输出cot数组对应的值。一共进行了4步

完整代码

#include <iostream>
using namespace std;
int N = 10000,n;
void Matrix(int (&a)[2][2],int b[2][2]){
    int tmp[2][2] = {0};
    for(int i = 0; i < 2; ++i)
        for(int j = 0; j < 2; ++j)
            for(int k = 0; k < 2; ++k)
                tmp[i][j] = (tmp[i][j] + a[i][k] * b[k][j]) % N;
    for(int i = 0; i < 2; ++i)
        for(int j = 0; j < 2; ++j)
            a[i][j] = tmp[i][j];
}
int main(){
    while(cin>>n && n!=-1){
        int temp[2][2] = {1, 1, 1, 0},cot[2][2] = {1, 0, 0, 1};
        while(n){
            if(n & 1) Matrix(cot,temp);
            Matrix(temp,temp);
            n /= 2;
        }
        cout<<cot[0][1]<<endl;
    }
    return 0;
}

这篇文章主要利用斐波那契数列来浅谈矩阵快速幂的用法,矩阵快速幂还有更大的使用空间,后续有空会补上矩阵快速幂的其他用法,另外送上一篇文章 十个利用矩阵乘法解决的经典题目

后续学习的传送门

矩阵快速幂——进阶
矩阵快速幂——实战

推荐阅读更多精彩内容