矩阵快速幂剪枝版（转载）



http://www.cnblogs.com/vongang/archive/2012/04/01/2429015.html

//矩阵快速幂

#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;

}

 

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("\n");

    }

    return 0;

}