bzoj 4128 矩阵求逆

/**************************************************************
    Problem: 4128
    User: idy002
    Language: C++
    Result: Accepted
    Time:4932 ms
    Memory:4152 kb
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <cmath>
#include <map>
using namespace std;
 
const int N = 70;
 
int n, Mod;
int inv[19997];
 
struct Matrix {
    int v[N][N];
    void unit() {
        for( int i=0; i<n; i++ )
            for( int j=0; j<n; j++ )
                v[i][j] = (i==j);
    }
    void read() {
        for( int i=0; i<n; i++ )
            for( int j=0; j<n; j++ )
                scanf( "%d", &v[i][j] );
    }
};
 
int mpow( int a, int b ) {
    int rt;
    for( rt=1; b; b>>=1,a=(a*a)%Mod )
        if( b&1 ) rt=(rt*a)%Mod;
    return rt;
}
 
Matrix operator*( const Matrix &a, const Matrix &b ) {
    Matrix c;
    for( int i=0; i<n; i++ ) {
        for( int j=0; j<n; j++ ) {
            c.v[i][j] = 0;
            for( int k=0; k<n; k++ ) {
                c.v[i][j] += a.v[i][k] * b.v[k][j] % Mod;
                if( c.v[i][j]>=Mod ) c.v[i][j]-=Mod;
            }
        }
    }
    return c;
}
Matrix operator~( Matrix a ) {
    Matrix b;
    b.unit();
    for( int i=0,j; i<n; i++ ) {
        for( int k=i; k<n; k++ )
            if( a.v[k][i] ) {
                j=k;
                break;
            }
        if( i!=j ) {
            for( int k=0; k<n; k++ ) {
                swap(a.v[i][k],a.v[j][k]);
                swap(b.v[i][k],b.v[j][k]);
            }
        }
        for( int j=i+1; j<n; j++ ) {
            int d = a.v[j][i]*inv[a.v[i][i]] % Mod;
            for( int k=0; k<n; k++ ) {
                a.v[j][k] = (a.v[j][k] - a.v[i][k]*d) % Mod;
                b.v[j][k] = (b.v[j][k] - b.v[i][k]*d) % Mod;
                if( a.v[j][k]<0 ) a.v[j][k]+=Mod;
                if( b.v[j][k]<0 ) b.v[j][k]+=Mod;
            }
        }
    }
    for( int i=n-1; i>=0; i-- ) {
        int d = inv[a.v[i][i]];
        for( int k=0; k<n; k++ ) {
            a.v[i][k] = a.v[i][k] * d % Mod;
            b.v[i][k] = b.v[i][k] * d % Mod;
        }
        for( int j=i-1; j>=0; j-- ) {
            d = a.v[j][i] * inv[a.v[i][i]] % Mod;
            for( int k=0; k<n; k++ ) {
                a.v[j][k] = (a.v[j][k] - a.v[i][k]*d) % Mod;
                b.v[j][k] = (b.v[j][k] - b.v[i][k]*d) % Mod;
                if( a.v[j][k]<0 ) a.v[j][k] += Mod;
                if( b.v[j][k]<0 ) b.v[j][k] += Mod;
            }
        }
    }
    return b;
}
bool operator<( const Matrix &a, const Matrix &b ) {
    for( int i=0; i<n; i++ )
        for( int j=0; j<n; j++ ) {
            if( a.v[i][j] ^ b.v[i][j] )
                return a.v[i][j] < b.v[i][j];
        }
    return false;
}
bool operator==( const Matrix &a, const Matrix &b ) {
    for( int i=0; i<n; i++ )
        for( int j=0; j<n; j++ )
            if( a.v[i][j] ^ b.v[i][j] ) return false;
    return true;
}
int ind( Matrix a, Matrix b ) {
    map<Matrix,int> mp;
    int m = int(sqrt(Mod))+1;
    Matrix base = a;
    a.unit();
    for( int i=0; i<m; i++ ) {
        if( a==b && i ) return i;
        mp[a] = i;
        a = a*base;
    }
    base = ~a;
    a = b*base;
    for( int i=m; ; i+=m,a=a*base ) 
        if( mp.count(a) ) return i+mp[a];
}
 
int main() {
    scanf( "%d%d", &n, &Mod );
    for( int i=1; i<Mod; i++ )
        inv[i] = mpow(i,Mod-2);
    Matrix a, b;
    a.read();
    b.read();
    printf( "%d
", ind(a,b) );
}

收获:

　　1. 矩阵进行如下操作可以相当于用一个矩阵乘以它:

　　　　将一行上的所有数乘以k

　　　　将一行加到另一行上

　　　　交换两行

　　2. 求逆的过程

　　　　如果要求矩阵A的逆矩阵A^-1,先得到一个单位矩阵B,

　　　　然后用上面1中的三种操作将A变成单位矩阵(不能变成单位矩阵则说明该矩阵行列式为0,即该矩阵不存在逆)

　　　　将对A的所有操作同样地应用于B,最终B就是A^-1

　　3. 求逆的正确性

　　　　我们对A进行了一系列变换,等同于用一个矩阵C乘以A使得 C*A = I

　　　　即C是A的逆矩阵, 将同样的操作作用于B,得到的矩阵为 C*B = C*I = C

　　　　即最终B的结果就是我们要求的逆

　　4. 高斯消元的另一种理解

　　　　A*X = B

　　　　C*A*X = C*B

　　　　完了