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

;
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\n", ind(a,b) );
}

View Code

收获:

　　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

　　　　完了