矩阵快速幂——入门

题目链接POJ307

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

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

先聊聊快速幂

首先，如果要你自己写一个pow()函数，需要求2²⁰，怎么求

#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时，这个复杂度的时间差别就很大了

怎么得到这个快速幂的代码的呢？我们假设需要求2⁷⁴

74转化为二进制为1001010，根据二进制对应位的1，可以把74拆分为64+8+2 ，所以
74 = 2⁶ + 2³ + 2¹
2⁷⁴ = 2²⁶ + 2²³ + 2²¹

 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 = 2ⁿ | 2¹| 2² | 2⁴ | 2⁸ | 2¹⁶| 2³² | 2⁶⁴ | 2¹²⁸
ans | 2⁰| 2⁰ | 2⁰×2² | 2⁰×2² | 2⁰×2²×2⁸| 2⁰×2²×2⁸ | 2⁰×2²×2⁸ | 2⁰×2²×2⁸×2⁶⁴

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

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

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

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

• 宏观上看矩阵

上图中，乘号右边是一系列的矩阵。我们称为B₁、B₂、B₃、B₄ …… B_n
该序列中的同行为一个斐波那契数列。乘号左边是一个转移矩阵A，B_n = A×B_n-1

• 微观上看矩阵

  
  

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

• 总结一波
  现在需要求的是B_n%10000，B_n = B₁ × A^n-1，B₁是[1,1]，所以 A^n-1就是需要的答案。那把求B_n改成了求A^n-1有什么好处呢？见下文

矩阵快速幂的实际运用

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

如果斐波那契第一项是A，第二项是A²，第三项是A³，第4项是A⁴ ……
那第8项的值是，A⁸ = A⁴×A⁴

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

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

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

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

代码拆分

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

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]; //赋值返回数据
}

A¹，也就是转移矩阵，题中已经给出了，用数组表示是[1,1,1,0]。这里还需要一个A⁰，A⁰其实就是单位矩阵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数组一开始是A⁰，temp数组一开始是A¹。我们来模拟下求第13项的斐波那契值是怎么求的。

第一步:n = 13，执行Matrix(cot,temp)和Matrix(temp,temp)，得到cot数组为A¹，t数组为A²，n = 6

第二步:n = 6，执行Matrix(temp,temp)，得到cot数组为A¹，t数组为A⁴，n = 3

第三步:n = 3，执行Matrix(cot,temp)和Matrix(temp,temp)，得到cot数组为A⁵，t数组为A⁸，n = 1

第四步:n = 1，执行Matrix(cot,temp)和Matrix(temp,temp)，得到cot数组为A¹³，t数组为A¹⁶，n = 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;
}

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

后续学习的传送门

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