POJ 3233 - Matrix Power Series ( 矩阵快速幂 + 二分)

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long LL;
#define MAX_SIZE 30
#define CLR( a, b ) memset( a, b, sizeof(a) )

int MOD = 0;
int n, k;

struct Mat
{
    LL mat[MAX_SIZE][MAX_SIZE];
    Mat() 
    {
        CLR( mat, 0 );
    }
    void zero()
    {
        CLR( mat, 0 );
    }
    void setv( int v )
    {
        for( int i = 0; i < n; ++i )
            for( int j = 0; j < n; ++j )
                mat[i][j] = v;
    }
    void init()
    {
        for( int i = 0; i < n; ++i )
            for( int j = 0; j < n; ++j )
                mat[i][j] = ( i == j );    
    }
    Mat operator + ( const Mat &b ) const
    {
        Mat c;
        for( int i = 0; i < n; ++i )
            for( int j = 0; j < n; ++j )
                c.mat[i][j] = ( mat[i][j] + b.mat[i][j] ) % MOD;
        return c;
    }    
    Mat operator * ( const Mat &b ) const
    {
        Mat c;
        for( int k = 0; k < n; ++k )
            for( int i = 0; i < n; ++i ) if( mat[i][k] )
                for( int j = 0; j < n; ++j )
                    c.mat[i][j] = ( c.mat[i][j] + mat[i][k] * b.mat[k][j] ) % MOD;
        return c;
    }
    void debug()
    {
        for( int i = 0; i < n; ++i )
        {
            for( int j = 0; j < n; ++j )
            {
                if( j != 0 )    putchar( ‘ ‘ );
                printf( "%lld", mat[i][j] );
            }
            putchar( ‘\n‘ );
        }
    }
};


Mat fast_mod( Mat a, int b )
{
    Mat res;
    res.init();
    while( b )
    {
        if( b & 1 )    res = res * a;
        a = a * a;
        b >>= 1; 
    }
    return res;
}


Mat solve( Mat a, int b )
{
    if( b == 1 )
        return a;
    Mat res;
    res.init();
    res = res + fast_mod( a, b >> 1);
    res = res * solve( a, b >> 1 );
    if( b & 1 ) 
        res = res + fast_mod( a, b );
    return res;
}

void Orz()
{
    while( ~scanf( "%d %d %d", &n, &k, &MOD ) )
    {
        Mat c;
        for( int i = 0; i < n; ++i )
            for( int j = 0; j < n; ++j )
                scanf( "%d", &c.mat[i][j] );
        Mat res = solve( c, k );
        res.debug();
    }
}

int main()
{
    Orz();
    return 0;
}