poj 3233 Matrix Power Series （矩阵快速幂）

这里介绍两种做法：

一：思路：
做这题时用的是类似二分的方法做，仔细观察可以发现
Sk = A + A2 + A3 + … + Ak   
    =(1+Ak/2)*(A + A2 + A3 + … + Ak/2  )+{Ak}
    =(1+Ak/2)*(Sk/2 )+{Ak}
当k为偶数时不要大括号里面的数；
代码：

//800ms
#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
int n , m ;
struct matrix {//矩阵
    int a[30][30];
    void init(  )
    {//初始化
        int i ;
        memset ( a, 0 ,sizeof ( a ) ) ;
        for ( i = 0 ; i < 30 ; i ++ )
            a[i][i] = 1 ;
    }
};
void print ( matrix s )
{//打印一个矩阵
    int i , j ;
    for ( i = 0 ; i < n ; i ++ )
    {
        for ( j = 0 ; j < n ; j ++ )
        {
            if ( j )
                cout<<' ';
            cout<<s.a[i][j]%m;
        }
        cout<<endl;
    }
}
matrix matrixadd ( matrix a , matrix b )
{//矩阵加法，计算时记得要取模，避免超int
    int i , j ;
    matrix c ;
    for ( i = 0 ; i < n ; i ++ )
        for ( j = 0  ; j < n ; j ++ )
            c.a[i][j]=((a.a[i][j]+b.a[i][j])%m);
    return c ;
}
matrix matrixmul ( matrix a , matrix b )
{//矩阵乘法，计算时记得要取模，避免超int
    int i , j , k ;
    matrix c ;
    for ( i = 0 ; i < n ; i ++ )
    {
        for ( j = 0 ; j < n ; j ++ )
        {
            c.a[i][j]=0;
            for ( k = 0 ; k < n ; k ++ )
                c.a[i][j] +=((a.a[i][k]*b.a[k][j])%m) ;
            c.a[i][j] %= m ;
        }
    }
    return c ;
}
matrix mul ( matrix s , int k  )
{//矩阵的k次方，快速幂
    matrix ans ;
    ans .init () ;
    while ( k >= 1 )
    {
        if ( k & 1 )
            ans = matrixmul ( ans , s ) ;
        k = k >> 1 ;
        s = matrixmul ( s , s ) ;
    }
    return ans ;
}
matrix sum ( matrix s , int k )
{//举证前k项求和
    if ( k == 1 )
        return s ;
    matrix tmp ;//用来保存答案
    tmp.init();//初始化
    tmp = matrixadd ( tmp , mul ( s , k >> 1 ) );    //计算1+A^(k/2)
    tmp = matrixmul ( tmp , sum ( s , k >> 1 ) ) ;    //计算(1+A^(k/2))*(A + A^2 + A^3 + … + A^(k/2)  )
    if ( k&1 )//考虑是否要+A^k
        tmp = matrixadd ( tmp , mul ( s , k ) ) ;
    return tmp ;//返回前n项的值
}

int main()
{
    int k ;
    while ( cin >> n >> k >> m )
    {
        int i , j ;
        matrix s ;
        for ( i = 0 ; i < n ; i ++ )
            for ( j = 0 ; j < n ; j ++ )
                cin >> s.a[i][j] ;
        s = sum ( s , k ) ;
        print(s);    
    }
    return 0;
}

二：

构造一个矩阵T

T={A  I}  

      0   I

所以T*T={A*A   A*I+I*I}
                   0           I

那么容易得到T^(K+1)={A^(K+1)           I+A+A^2+A^3+...+A^K}
                                         0                                I

这样我们只需要算T的k+1次幂就可以了；

代码：

#include <cstdio>
#include <cstring>
#include <iostream>
#include <cmath>
#include <algorithm>
const int maxn = 66;
using namespace std;

struct mat{
   int v[maxn][maxn];
   mat() {
       memset(v, 0, sizeof(v));
   }
};
int n, m;

mat matrix_mul(mat a, mat b) {
    mat c;
    int i, j, k;
    for(i = 0; i < 2*n; i++)
        for(j = 0; j < 2*n; j++)
                for(k = 0; k < 2*n; k++) {
                    c.v[i][k] += (a.v[i][j]*b.v[j][k])%m;
                    c.v[i][k] %= m;
                }
    return c;
}

mat matrix_mi(mat p, int k) {
    mat t;
    for(int i = 0; i < 2*n; i++)
        t.v[i][i] = 1;
    while(k) {
        if (k & 1)
            t = matrix_mul(t, p);
        k >>= 1;
        p = matrix_mul(p, p);
    }
    return t;
}

int main() {
    int k;
    cin>>n>>k>>m;
    mat p, t;
    for(int i = 0; i < n; i++) {
        p.v[i][i+n] = p.v[i+n][i+n] = 1;
        for(int j = 0; j < n; j++)
            cin>>p.v[i][j];
        }
        t = matrix_mi(p, k + 1);
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++)
            if (i != j)
                printf("%d ", t.v[i][j+n]);
            else printf("%d ", (t.v[i][j+n] + m - 1)%m);
        printf("\n");
    }
	return 0;
}