uva11149 - Power of Matrix 矩阵倍增

最新推荐文章于 2019-04-04 18:18:19 发布

原创最新推荐文章于 2019-04-04 18:18:19 发布 · 673 阅读

0 ·

CC 4.0 BY-SA版权

数学专栏收录该内容

46 篇文章

订阅专栏

本文介绍了一种利用矩阵快速幂高效计算矩阵A+A^2+...+A^k的方法，其中k可达百万级别。文章提供了完整的C++实现代码，并通过实例解释了算法的基本原理及其在图论中路径计数问题的应用。

Consider an n-by-n matrixA. We define A^k = A *A * ... * A (k times). Here, * denotes the usual matrix multiplication.

You are to write a program that computes the matrix A + A² + A³ + ... + A^k.

Example

Suppose A = . ThenA² = =, thus:

Such computation has various applications. For instance, the above example actually counts all the paths in the following graph:

Input consists of no more than 20 test cases. The first line for each case contains two positive integersn (≤ 40) and k (≤ 1000000). This is followed by n lines, each containing n non-negative integers, giving the matrixA.

Input is terminated by a case where n = 0. This case need NOT be processed.

Output

For each case, your program should compute the matrix A + A² + A³ + ... + A^k. Since the values may be very large, you only need to print theirlast digit. Print a blank line after each case.

Sample Input

Sample Output

0 2 4
0 0 2
0 0 0

计算A+A^2+A^3...+A^k。

当k是偶数，A+A^2+A^3...+A^k=（E+A^(k/2)）*(A+A^2...+A^(k/2))。当k是奇数，A+A^2+A^3...+A^k=（E+A^(k/2)）*(A+A^2...+A^(k/2))+A^k。

类似二分的思想。话说被重载运算符优先级坑了，看来以后重载运算符要打括号啊。。

#include<iostream>
#include<algorithm>
#include<queue>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<cstring>
#define INF 0x3f3f3f3f
#define MAXN 45
#define MAXM 110
#define MOD 1000000007
#define MAXNODE 4*MAXN
#define eps 1e-9
using namespace std;
typedef long long LL;
int N,K;
struct Mat{
    int mat[MAXN][MAXN];
    Mat(){
        memset(mat,0,sizeof(mat));
    }
};
Mat operator * (Mat a,Mat b){
    int i,j,k;
    Mat ret;
    for(k=0;k<N;k++){
        for(i=0;i<N;i++){
            if(a.mat[i][k]==0) continue;
            for(j=0;j<N;j++){
                if(b.mat[k][j]==0) continue;
                ret.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                ret.mat[i][j]%=10;
            }
        }
    }
    return ret;
}
Mat operator + (Mat a,Mat b){
    Mat ret;
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++) ret.mat[i][j]=(a.mat[i][j]+b.mat[i][j])%10;
    return ret;
}
Mat operator ^ (Mat a,int n){
    Mat ret,t=a;
    int i,j;
    for(i=0;i<N;i++)
        for(j=0;j<N;j++) ret.mat[i][j]=(i==j);
    while(n){
        if(n&1) ret=ret*t;
        t=t*t;
        n>>=1;
    }
    return ret;
}
Mat solve(Mat a,int k){
    if(k==1) return a;
    Mat ans;
    for(int i=0;i<N;i++) ans.mat[i][i]=1;
    if(k==0) return ans;
    ans=((a^(k>>1))+ans)*solve(a,k>>1);
    if(k%2) ans=(a^k)+ans;
    return ans;
}
Mat a,ans;
int main(){
    freopen("in.txt","r",stdin);
    while(scanf("%d%d",&N,&K)!=EOF&&N){
        for(int i=0;i<N;i++)
            for(int j=0;j<N;j++){
                int t;
                scanf("%d",&t);
                a.mat[i][j]=t%10;
            }
        ans=solve(a,K);
        for(int i=0;i<N;i++){
            for(int j=0;j<N-1;j++) printf("%d ",ans.mat[i][j]);
            printf("%d\n",ans.mat[i][N-1]);
        }
        puts("");
    }
    return 0;
}