Power of Matrix（uva11149+矩阵快速幂）

Power of Matrix

Time Limit:3000MS     Memory Limit:0KB     64bit IO Format:%lld & %llu

Submit Status Practice UVA 11149

Appoint description: System Crawler  (2015-03-15)

Description

 

Problem B : Power of Matrix

Time limit: 10 seconds

Consider an n-by-n matrix A. 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 = . Then A² =  = , thus:

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

 

Input

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

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 their last digit. Print a blank line after each case.

 

Sample Input

3 20 2 00 0 20 0 00 0

 

Sample Output

0 2 40 0 20 0 0

 

 

 

首先我们来想一下计算A+A^2+A^3...+A^k。

如果A=2,k=6。那你怎么算                        

2+2²+2³+2⁴+2⁵+2⁶ = ？= (2+2²+2³)*(1+2³)

 

如果A=2,k=7。那你怎么算                        

2+2²+2³+2⁴+2⁵+2⁶+2⁷ = ？= (2+2²+2³)*(1+2³)+2⁷

 

^{so....同理：}

当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。

 

转载请注明出处：寻找&星空の孩子

题目链接：UVA 11149

 

#include<cstdio>#include<cstring>#include<iostream>#include<algorithm>using namespace std;#define LL __int64#define mmax 45struct matrix{    int mat[mmax][mmax];};int N;matrix multiply(matrix a,matrix b){    matrix c;    memset(c.mat,0,sizeof(c.mat));    for(int i=0; i<N; i++)    {        for(int j=0; j<N; j++)        {            if(a.mat[i][j]==0)continue;            for(int k=0; k<N; k++)            {                if(b.mat[j][k]==0)continue;                c.mat[i][k]=(c.mat[i][k]+a.mat[i][j]*b.mat[j][k])%10;            }        }    }    return c;}matrix quickmod(matrix a,int n){    matrix res;    for(int i=0; i<N; i++) //单位阵        for(int j=0; j<N; j++)            res.mat[i][j]=(i==j);    while(n)    {        if(n&1)            res=multiply(a,res);        a=multiply(a,a);        n>>=1;    }    return res;}matrix add (matrix a,matrix b){    matrix 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;}matrix solve(matrix a,int k){    if(k==1) return a;    matrix ans;    for(int i=0; i<N; i++)        for(int j=0; j<N; j++)            ans.mat[i][j]=(i==j);    if(k==0) return ans;    ans=multiply((add(quickmod(a,(k>>1)),ans)),solve(a,(k>>1)));    if(k%2) ans=add(quickmod(a,k),ans);    return ans;}int main(){    int k;    while(scanf("%d%d",&N,&k)!=EOF)    {        if(!N)break;        matrix ans;        for(int i=0;i<N;i++)        {            for(int j=0;j<N;j++)            {                int temp;                scanf("%d",&temp);                ans.mat[i][j]=temp%10;            }        }        ans=solve(ans,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]);        }        printf("\n");    }    return 0;}