HDU 4965 Fast Matrix Caculation ( 矩阵乘法 + 矩阵快速幂 + 矩阵乘法的结合律 )

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

LL A[MAX_SIZE][6] , B[6][MAX_SIZE], C[MAX_SIZE][6], D[MAX_SIZE][MAX_SIZE];
LL n, k;
struct Mat
{
    int n;
    LL mat[6][6];
    Mat( int _n)
    {
        n = _n;
        CLR( mat, 0 );
    }
    void zero()
    {
        CLR( mat, 0 );
    }
    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( b.n );
        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( "%I64d", mat[i][j] );
            }
            putchar( ‘\n‘ );
        }
    }
};

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

void scan( LL &x )
{
    char c;
    while( c = getchar(), c < ‘0‘ || c > ‘9‘ );
    x = c - ‘0‘;
    while( c = getchar(), c >= ‘0‘ && c <= ‘9‘ ) x = x * 10 + c -‘0‘;
}

int main()
{
    while( ~scanf( "%I64d %I64d", &n, &k ) && n + k )
    {
        for( int i = 0; i < n; ++i )
            for( int j = 0; j < k; ++j )
                scan( A[i][j] );        
        for( int i = 0; i < k; ++i )
            for( int j = 0; j < n; ++j )
                scan( B[i][j] );
        Mat t( k );
        for( int i = 0; i < k; ++i )
            for( int j = 0; j < k; ++j )
                for( int y = 0; y < n; y++ )
                    t.mat[i][j] = ( t.mat[i][j] + B[i][y] * A[y][j] ) % MOD;
        t = fast_mod( t, n * n - 1 );
        
        CLR( C, 0 ), CLR( D, 0 );
        LL res = 0;
        for( int i = 0; i < n; ++i )
            for( int j = 0; j < k; ++j )
                for( int y = 0; y < k; y++ )
                    C[i][j] = ( C[i][j] + A[i][y] * t.mat[y][j] ) % MOD;
        for( int i = 0; i < n; ++i )
        {
            for( int j = 0; j < n; ++j )
            {
                for( int y = 0; y < k; y++ )
                {
                    D[i][j] = ( D[i][j] + C[i][y] * B[y][j] ) % 6;
                }
                res += D[i][j];
            }
        }
        printf( "%I64d\n", res );
    }
    return 0;
}