算法分类整理+模板①：矩阵快速幂

一直有一个想法，感觉自己很多基础算法不是很扎实，想要找个机会写一些算法的整理，顺便自己总结一些实用的模板。

最近偶然在训练赛中连续做了2道思维+矩阵快速幂的题目，碰巧有时间，就以矩阵快速幂作为这个系列博客的开始吧。

 

如果想要了解矩阵快速幂，首先要了解什么叫做快速幂。

举个例子，如果让你求2^10的值，一个for循环可以轻松地解决问题，但是如果是2^10000000000000呢？

且不管这个值能否表示出来，单单说for循环的时间复杂度O(n)就注定不能直接暴力求解。

当然，为了求得这个解，我们一般要求答案对于某个数取模，常用的MOD值有10007,1000000007。

由此，我们可以看出，当问题存在超时+取模的限制时，我们需要一种新的算法，即快速幂。

快速幂是基于二分思想的一种时间复杂度为O(lgn)的算法。

我们可以考虑一个例子，如果要求2^10的值，我们能否这样算：

首先把2^10分解成(2^5)*(2^5)，其次把2^5分解成2*(2^4)，然后将2^4分成(2^2)*(2^2)，最后把2^2变成2*2。

这样，我们就将2^10变成了：(2*(2*2)^2)^2。这样我们只需要计算4次乘法就可以得到2^10的值，而线性的算法需要10次，快速幂进行了极大地优化。

一般地，对于a^b来说，当b为偶数时，我们可以写成(a^(b/2))^2；当b为奇数时，可以写成a*(a^(b-1))。

所以，经过快速幂算法优化后的quick_pow计算只需要log(b)次！b越大，这个优化就越明显！

模板代码如下：

 

#include <stdio.h>
#include <iostream>
#define MOD 1000000007
typedef long long int LL
using namespace std;
LL quick_pow(int a,int b){
    LL res=1;
    while(b){
        if(b&1)
            res=((a%MOD)*(res%MOD))%MOD;
        res=(res*res)%MOD;
        b>>=1;
    }
    return res;
}

 

 

了解了快速幂的基本思想与代码实现，我们就要来看看矩阵快速幂。其实矩阵快速幂的基本思想和普通快速幂是一样的。

对于矩阵，我相信学过线性代数的同学应该对其深恶痛绝……不要怕，我们的矩阵快速幂只涉及到矩阵的乘法，是不是很简单啊~(对于不会矩阵乘法的同学请自行百度)

对于矩阵乘法，模板如下：

 

#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <string.h>
#define MOD 1000000007
typedef long long int LL;
using namespace std;
struct Matrix{
    int a[2][2];
    Matrix(){
        memset(a,0,sizeof(a));
    }
    Matrix operator* (const Matrix &p){
        Matrix res;
        for(int i=0;i<2;i++){
            for(int j=0;j<2;j++){
                for(int k=0;k<2;k++){
                    res.a[i][j]+=(a[i][k]*p.a[k][j]%MOD);
                }
                res.a[i][j]%=MOD;
            }
        }
        return res;
    }
}init,unit;

 

首先来看一个例子：如果问题是求解斐波那契数列的第1000000000项对于MOD取模的值，我们应该如何去做？

从快速幂我们知道对于这个问题线性的计算是肯定超时的，所以我们依旧采用快速幂的思想。

对于斐波那契数列我们有如下规律(以下为矩阵)：

f(n+1) f(n)    乘以    1     1     结果是     f(n+1)+f(n)  f(n+1)　　即：f(n+2)  f(n+1)

f(n)  f(n-1)             1     0　　　　　　  f(n)+f(n-1)  f(n)　　　　     f(n+1)  f(n)

这是我们会惊奇的发现 当我们用斐波那契数列组成的矩阵乘以一个特定矩阵时会得到下一个斐波那契数的值。

所以不难想象，我们只要知道数列的前3项，用他们组成的矩阵乘以这个特定矩阵的k次幂就能得到任意项的斐波那契数，并且时间复杂度是O(lgn)的！

所以，到这里我们就要知道如何去找这个特定矩阵。

一般地，如果有通项公式：f(n)=a*f(n-1)+b*f(n-2)+c*f(n-3)……(这里我们以3个为例)，若f(1)==p,f(2)==q,f(3)==r

我们设定一个init矩阵表示初始值：

r  0  0

q  0  0

p  0  0

一个unit矩阵表示那个特定矩阵：

a  b  c

1  0  0

0  1  0

这样，unit矩阵左乘init矩阵等于：

ap+bq+cs  0  0

            p  0  0

            q  0  0

这样我们就构造出了一个矩阵，表示出了整个数列的递推关系：

 a  b  c　　　　　　　　　　　f(3)  0  0　　　　　  f(n+2)  0  0

(1  0  0)     的n次幂   乘以     f(2)  0  0　　等于  　f(n+1)  0  0

 0  1  0　　　　　　　　　　  f(1)  0  0　　　　　      f(n)  0  0

 

当然，构造这样的矩阵的方法不一，此方法只是较为通用的方法，对于某些通项公式可以找到更简便的矩阵使得矩阵快速幂成立。

 

矩阵快速幂模板如下：

#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <string.h>
#define MOD 1000000007
typedef long long int LL;
using namespace std;
struct Matrix{
    int a[2][2];
    Matrix(){
        memset(a,0,sizeof(a));
    }
    Matrix operator* (const Matrix &p){
        Matrix res;
        for(int i=0;i<2;i++){
            for(int j=0;j<2;j++){
                for(int k=0;k<2;k++){
                    res.a[i][j]+=(a[i][k]*p.a[k][j]%MOD);
                }
                res.a[i][j]%=MOD;
            }
        }
        return res;
    }
}init,unit;
Matrix quick_pow(Matrix unit,Matrix init,int k){
    while(k){
        if(k&1)
            init=init*unit;
        unit=unit*unit;
        k>>=1;
    }
    return init;
}
void init_Matrix(){
    init.a[0][0]=1;
    init.a[0][1]=0;
    init.a[1][0]=0;
    init.a[1][1]=1;
    unit.a[0][0]=1;
    unit.a[0][1]=1;
    unit.a[1][0]=1;
    unit.a[1][1]=0;
}

 

主函数中首先对init和unit矩阵进行初始化，然后再调用quick_pow()。

 

 

小Tips：关于如何识别矩阵快速幂的问题。

一般来说，题目中如果有”答案对于MOD取模“这句话时，并且操作次数巨大，我们就可以考虑使用快速幂或矩阵快速幂。

关于题目中有”取模“的说法时，一般来说有几种可能。一是快速幂，二是dp，三是组合数学。

另外推荐两道矩阵快速幂的题目和题解：

http://www.cnblogs.com/Torrance/p/5412802.html

http://www.cnblogs.com/Torrance/p/5410755.html