poj 3233 矩阵快速幂

题意：给定一个n*n矩阵A，一个常数k和另一个常数m，求S = A + A^2 + A^3 + … + A^k。

思路：构造新矩阵T={A I }，则T^k={A^K   I+A+A^2+A^3+...+A^(K-1)}
                                   I  I                   0                     I

显然只要求出T^k即可，方法可以使用快速幂。即将k分解成若干2的幂相加的形式，那位为1则加上相应的T的幂。（http://blog.youkuaiyun.com/wangjian8006/article/details/7868864）

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#define clr(s,t) memset(s,t,sizeof(s))
#define N 66
int n,k,m,b;
struct matrix{
    int ar[N][N];
}res,s;
struct matrix mul(struct matrix s,struct matrix t){//矩阵相乘
    int i,j,p;
    struct matrix ans;
    clr(ans.ar, 0);
    for(i = 0;i<b;i++)
        for(j = 0;j<b;j++){
            for(p = 0;p<b;p++)
                ans.ar[i][j] += (s.ar[i][p] * t.ar[p][j] % m);
            ans.ar[i][j] %= m;
        }
    return ans;
}
void test(){
    while(k){
        if(k&1)
            res = mul(res,s);
        s = mul(s,s);
        k >>= 1;
    }
}
int main(){
    int i,j;
    scanf("%d %d %d",&n,&k,&m);
    clr(res.ar, 0);
    clr(s.ar, 0);
    for(i = 0;i<n;i++)
        for(j = 0;j<n;j++)
            scanf("%d",&s.ar[i][j]);
    for(i = 0;i<n;i++){
        s.ar[i][i+n] = 1;
        res.ar[i][i] = 1;
        s.ar[i+n][i+n] = res.ar[i+n][i+n] = 1;
    }
    b = n<<1;
    test();
    for(i = 0;i<n;i++){
        for(j = 0;j<n;j++){
            if(i==j)
                printf("%d ",(res.ar[i][j]+res.ar[i][j+n]-1)%m);
            else
                printf("%d ",(res.ar[i][j]+res.ar[i][j+n])%m);
        }
        printf("\n");
    }
    return 0;
}