矩阵快速幂

矩阵快速幂其实并不是一个很难的算法，相对来说其实还算是比较简单，但是初看可能会难以理解他到底在干什么，其实有些我们只要动手模拟模拟就好，还有就是公式比较重要，看到有些博客拿斐波那契数列举例子，但是又不给出公式，很是让人头疼。

 

1.快速幂

在了解矩阵快速幂之前我们需要先了解快速幂。何为快速幂？即加快求幂运算的一种方式。

例如让我们求解2的7次方，我们之前的算法是将2重复的乘7次，也就是循环7次，但是有了矩阵快速幂之后，可以将循环次数减少到三次。我们将2的7次方进行分解，分解为2*(2*2)*(2*2*2*2)。可以划分为4个2* 2个2 * 1个2，这刚好与我们7的二进制的每一位数字对应，也就是111,我们知道了两个2相乘的结果，那么利用两个2相乘的结果去求解4个2相乘就只需要乘1次。先说到这里，剩下的看代码之后应该就能理解其中的含义了。

// n代表几次幂
int pow(int a,int n){
    int b = 1;
    while(n){
        if(n&1){
            b *= a;
        }
        // a进行翻倍，这样就能少乘几次了
        a *= a;
        n >>= 1;
    }
    return b;
}

上面的7次可能会有点特殊，因为其二进制的每一位都是1。可以测试一下8次吧，8的二进制只有第4位才为1。

 

2.矩阵快速幂

这里的矩阵就是我们线性代数中的矩阵，看这个之前必须要先了解线性代数中的矩阵乘法，这个就自行百度吧。

我们还是以斐波那契数列数列来讲解，首先要给出公式，斐波那契数列的递推公式为，我们如何用矩阵乘法表示它呢？

我们知道了矩阵的乘法，然后我们就知道上面f[i] = 1 * f[i-1] + 1 * f[i-2]，这与我们的斐波那契数列公式是不是一模一样，再往下看

有了这个式子，我们只要算出即可直接由和求出。那么怎么快速计算出呢？快速幂即可！但有个需要注意的地方：当指数为零时，如果是数字，，但矩阵应该怎么办呢？我们需要用到神奇的单位矩阵，它的主对角线上的数字为，其它都为，它有个有用的性质：任何矩阵乘以它都会等于本身。

其实你可能会有所疑问，为什么算出就可以知道f(n)的结果，因为f(n)就是前面无数个f(2)和f(1)的和，算矩阵乘法的过程中其实就是一直去加f(2)+f(1)，快速幂矩阵就是能一次就很多个f(2)和f(1)，最终得到f(n)，可以先简单这样理解吧，剩下的还是要看代码。看代码的过程记得带上面的公式

代码实现：

#define N 2

int mod=10000;
int res[N][N];
int temp[N][N];
int a[N][N];
// 线性代数矩阵乘法
void Mul(int a[][N],int b[][N])
{
    memset(temp,0,sizeof(temp));
    for(int i=0; i<N; i++)
        for(int j=0; j<N; j++)
            for(int k=0; k<N; k++)
                temp[i][j]=(temp[i][j]+a[i][k]*b[k][j])%mod;
    for(int i=0; i<N; i++)
        for(int j=0; j<N; j++)
            a[i][j]=temp[i][j];
}
void fun(int a[][N],int n)
{
    memset(res,0,sizeof(res));
    // 单位矩阵初始化，单位矩阵与我们的数字1是一种概念，任意一个矩阵乘单位矩阵都不会变化
    for(int i=0; i<N; i++)
        res[i][i]=1;
    // 快速幂
    while(n)
    {
        if(n&1)
        {
            Mul(res,a);
        }
        Mul(a,a);
        n>>=1;
    }
}

int main()
{
    printf("%d",pow(2,31)-1);
    int n;
    while(~scanf("%d",&n)&&n!=-1)
    {
        // 与我们上面给定的公式是一样的
        // 这里的值只是普通的斐波那契是这个样子
        // 如果矩阵的对应关系不一样，矩阵中的值会有所变化
        // 其实知道原理之后，有些值还是很好推
        a[0][0]=1;
        a[1][0]=1;
        a[0][1]=1;
        a[1][1]=0;
        if(n==0||n==1)
            printf("%d\n",n);
        else{
            // n代表的是斐波那契数列的第几个
            // 注意下标开始的位置
            fun(a,n);
            printf("%d\n",res[1][0]);
        }
    }
}

3.小应用

上面的斐波那契算是一个经典的应用，下面这个题也是一个经典的应用。

HDU 2157:http://acm.hdu.edu.cn/showproblem.php?pid=2157

一看是个图论问题，刚开始看这个题真没发现与矩阵快速幂有多大的关系。

简述：给你一个n个点，给你然后告诉你哪个点能到哪个点，注意是有向的，最后给你两个点，问你从点a到点b经过k个点有多少种方案，注意：可走重复边，即一个边可以经过多次。

首先将给定的有向图转换为邻接矩阵，这就得到了我们快速幂矩阵中的矩阵

假设给定了4个点，让我们求0到3这个点有多少方案，那么点i到点j经过两个点的方案就有f(i,j) = f (i,k) + f(k,j) k是中间点，i，j是起点和终点，我们就需要去枚举这个中间点。枚举的过程就是矩阵的乘法的过程。三个点的方案也就用之前两个点去乘，因为路径可以重复走，如果不懂的话，就将邻接矩阵画出来模拟一下吧。

代码实现：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define N 21

int mod = 1000;
int map[N][N];
int temp[N][N];
int ts[N][N];
int res[N][N];

void matrix(int a[][N],int b[][N],int n){
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            temp[i][j] = 0;
            for(int k=0;k<n;k++){
                temp[i][j] = (temp[i][j] + a[i][k] * b[k][j])%mod;
            }
        }
    }

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            a[i][j] = temp[i][j];
        }
    }
}

void quickPow(int a[][N],int n,int size){
    memset(res,0,sizeof(res));
    for(int i=0;i<size;i++){
        res[i][i] = 1;
    }
    while(n){
        if(n&1){
            matrix(res,a,size);
        }
        matrix(a,a,size);
        n >>= 1;
    }
}

int main()
{
    int n,m,x,y,t,k;
    while(scanf("%d%d",&n,&m) == 2 && (n+m)){
        memset(map,0,sizeof(map));
        for(int i=0;i<m;i++){
            scanf("%d%d",&x,&y);
            map[x][y] = 1;
        }
        scanf("%d",&t);
        for(int i=0;i<t;i++){
            scanf("%d%d%d",&x,&y,&k);
            memcpy(ts,map,sizeof(map));
            quickPow(ts,k,n);
            printf("%d\n",res[x][y]);
        }
    }
    return 0;
}