矩阵快速幂

参考文章：

http://blog.csdn.net/runninghui/article/details/8905019

http://blog.csdn.net/u013795055/article/details/38599321

http://blog.csdn.net/g_congratulation/article/details/52734306

快速幂：

求幂时最朴素的方法就是循环n次进行求解，这样得到的是复杂度为O(n)的算法，还有一种比较优化的方法就是，一直这样递归下去，就可以省去部分计算，这样的时间复杂度约为O(logn);

快速幂是一种用二进制方法求幂的运算

比如说：（注意：10的二进制表示法为：1010，分别对应A的次幂）

while(n)
{
    if(n & 1)
    {
        ans *= matex;
    }
    matex = matex * matex;
    n >> 1;
}

矩阵相乘

//O(n^3)算法
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
using namespace std;
#define LL __int64
#define MOD 1000
typedef struct MATRIX
{
    int mat[50][50];
}MATRIX;

MATRIX mat_multiply (MATRIX a,MATRIX b,int n)
{
    MATRIX c;   //c[i][j]= Σ a[i][k]*b[k][j]
    memset(c.mat,0,sizeof(c.mat));    
/*
    for(int i=0;i<n;i++)    //a矩阵一行一行往下
        for(int j=0;j<n;j++)    //b矩阵一列一列往右
            for(int k=0;k<n;k++)    //使a矩阵 第i行第k个元素 乘以 b矩阵 第j列第k个元素
                if(a[i][k] && b[k][j])    //剪枝（添条件，设门槛），提高效率，有一个是0，相乘肯定是0
                    c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
*/

//上面也是可以的，但是下面的剪枝更好一些，效率更高一些，但是运算顺序有点难想通，，，
//上面就是C[i][j]一次就求出来，下面就是每次c[i][j]求出一项【就是上面红体字，每次各求一列】

    for(int k=0;k<n;k++)    //这个可以写到前面来，
        for(int i=0;i<n;i++)
            if(a.mat[i][k])     //剪枝：如果a.mat[i][k]是0，就不执行了
                for(int j=0;j<n;j++)
                    if(b.mat[k][j])     //剪枝：如果b.mat[i][k]是0，就不执行了
                    {
                        c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                        if(c.mat[i][j]>=MOD)    //这个看实际需求，要不要取模
                            c.mat[i][j]%=MOD;   //取模的复杂度比较高，所以尽量减少去模运算，添加条件，只有当大于等于MOD的时候才取余
                    }
    return c;
}

int main()//这个只是用来测试用的。。。
{
    int n;
    MATRIX A,B,C;

    memset(A.mat,0,sizeof(A.mat));
    memset(B.mat,0,sizeof(B.mat));
    memset(C.mat,0,sizeof(C.mat));

    scanf("%d",&n);     //矩阵规模，这里是方阵，行数等于列数

    for(int i=0;i<n;i++)    //初始化A矩阵
        for(int j=0;j<n;j++)
            scanf("%d",&A.mat[i][j]);

    for(int i=0;i<n;i++)    //初始化B矩阵
        for(int j=0;j<n;j++)
            scanf("%d",&B.mat[i][j]);

    C=mat_multiply (A,B,n);

    for(int i=0;i<n;i++)    //打印C矩阵
    {
        for(int j=0;j<n;j++)
            printf("%3d",C.mat[i][j]);
        printf("
");
    }
    return 0;
}

矩阵快速幂

矩阵快速幂的核心思想与快速幂基本一致。

//矩阵快速幂
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
using namespace std;
#define LL __int64
#define MOD 1000
typedef struct MATRIX
{
    int mat[50][50];
}MATRIX;

MATRIX mat_multiply (MATRIX a,MATRIX b,int n)
{
    MATRIX c;   //c[i][j]= Σ a[i][k]*b[k][j]
    memset(c.mat,0,sizeof(c.mat));
    for(int k=0;k<n;k++)   
        for(int i=0;i<n;i++)
            if(a.mat[i][k])    
                for(int j=0;j<n;j++)
                    if(b.mat[k][j])     
                    {
                        c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                        if(c.mat[i][j]>=MOD)   
                            c.mat[i][j]%=MOD;   
                    }
    return c;
}

MATRIX mat_quick_index(MATRIX a,int N,int n)
{
    MATRIX E;   //单位矩阵，就像数值快速幂里，把代表乘积的变量初始化为1

//    memset(E.mat,0,sizeof(E.mat));  //置零，单位矩阵除了主对角线都是1，其他都是0
//    for(int i=0;i<n;i++)    //初始化单位矩阵【就是主对角线全是1】
//            E.mat[i][i]=1;

    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            E.mat[i][j]=(i==j); //酷炫*炸天的初始化！！！

    while(N>0)
    {
        if(N & 1)
            E=mat_multiply(E,a,n);
        N>>=1;
        a=mat_multiply(a,a,n);
    }
    return E;
}

int main()
{
    int n,N;    //n为矩阵（方阵）规模，几行，N为指数
    MATRIX A,C;

    memset(A.mat,0,sizeof(A.mat));
    memset(C.mat,0,sizeof(C.mat));

    scanf("%d",&n);     //矩阵规模，这里是方阵，行数等于列数

    for(int i=0;i<n;i++)    //初始化A矩阵
        for(int j=0;j<n;j++)
            scanf("%d",&A.mat[i][j]);

    scanf("%d",&N);
    C=mat_quick_index(A,N,n);

    for(int i=0;i<n;i++)    //打印C矩阵
    {
        for(int j=0;j<n;j++)
            printf("%3d",C.mat[i][j]);
        printf("
");
    }
    return 0;
}

矩阵快速幂优化递推式：斐波那契数列

单位矩阵：除对角线为1外，其余均为0

矩阵乘法与递推式之间的关系：

在斐波那契数列中：

即

用矩阵快速幂求斐波那契数列的第N项的代码：

#include<cstdio>
#include<algorithm>
#include<string>
#include<iostream>
using namespace std;

const int M = 1e9+7;

struct Matrix{
    long long a[2][2];
    Matrix(){
        memset(a, 0, sizeof(a));
    }
    Matrix operator *(const Matrix y)
    {
        Matrix ans;
        for(int i = 0; i <= 1; i++)
        {
            for(int j = 0; j <= 1; j++)
            {
                for(int k = 0; k <= 1; k++)
                {
                    ans.a[i][j] += a[i][k] * y.a[k][j];
                }
            }
        }
        
        for(int i = 0; i <= 1; i++)
        {
            for(int j = 0; j <= 1; j++)
            {
                ans.a[i][j] %= M;
            }
        }
        return ans;
    }
    
    void operator = (const Matrix b)
    {
        for(int i = 0; i <= 1; i++)
        {
            for(int j = 0; j <= 1; j++)
            {
                a[i][j] = b.a[i][j];
            }
        }
    }
};

int solve(long long x)
{
    Matrix ans, trs;
    ans.a[0][0] = ans.a[1][1] = 1;
    trs.a[0][0] = trs.a[1][0] = trs.a[0][1] = 1;
    while(x)
    {
        if(x & 1)
        {
            ans = ans * trs;
        }
        trs = trs * trs;
        x >>= 1;
    }
    return ans.a[0][0];
}

int main()
{
    int n;
    scanf("%d", &n);
    cout << solve(n-1) << endl;
    return 0;
}