bzoj 4128 矩阵求逆

 

  1 /**************************************************************
  2     Problem: 4128
  3     User: idy002
  4     Language: C++
  5     Result: Accepted
  6     Time:4932 ms
  7     Memory:4152 kb
  8 ****************************************************************/
  9  
 10 #include <iostream>
 11 #include <cstdio>
 12 #include <cmath>
 13 #include <map>
 14 using namespace std;
 15  
 16 const int N = 70;
 17  
 18 int n, Mod;
 19 int inv[19997];
 20  
 21 struct Matrix {
 22     int v[N][N];
 23     void unit() {
 24         for( int i=0; i<n; i++ )
 25             for( int j=0; j<n; j++ )
 26                 v[i][j] = (i==j);
 27     }
 28     void read() {
 29         for( int i=0; i<n; i++ )
 30             for( int j=0; j<n; j++ )
 31                 scanf( "%d", &v[i][j] );
 32     }
 33 };
 34  
 35 int mpow( int a, int b ) {
 36     int rt;
 37     for( rt=1; b; b>>=1,a=(a*a)%Mod )
 38         if( b&1 ) rt=(rt*a)%Mod;
 39     return rt;
 40 }
 41  
 42 Matrix operator*( const Matrix &a, const Matrix &b ) {
 43     Matrix c;
 44     for( int i=0; i<n; i++ ) {
 45         for( int j=0; j<n; j++ ) {
 46             c.v[i][j] = 0;
 47             for( int k=0; k<n; k++ ) {
 48                 c.v[i][j] += a.v[i][k] * b.v[k][j] % Mod;
 49                 if( c.v[i][j]>=Mod ) c.v[i][j]-=Mod;
 50             }
 51         }
 52     }
 53     return c;
 54 }
 55 Matrix operator~( Matrix a ) {
 56     Matrix b;
 57     b.unit();
 58     for( int i=0,j; i<n; i++ ) {
 59         for( int k=i; k<n; k++ )
 60             if( a.v[k][i] ) {
 61                 j=k;
 62                 break;
 63             }
 64         if( i!=j ) {
 65             for( int k=0; k<n; k++ ) {
 66                 swap(a.v[i][k],a.v[j][k]);
 67                 swap(b.v[i][k],b.v[j][k]);
 68             }
 69         }
 70         for( int j=i+1; j<n; j++ ) {
 71             int d = a.v[j][i]*inv[a.v[i][i]] % Mod;
 72             for( int k=0; k<n; k++ ) {
 73                 a.v[j][k] = (a.v[j][k] - a.v[i][k]*d) % Mod;
 74                 b.v[j][k] = (b.v[j][k] - b.v[i][k]*d) % Mod;
 75                 if( a.v[j][k]<0 ) a.v[j][k]+=Mod;
 76                 if( b.v[j][k]<0 ) b.v[j][k]+=Mod;
 77             }
 78         }
 79     }
 80     for( int i=n-1; i>=0; i-- ) {
 81         int d = inv[a.v[i][i]];
 82         for( int k=0; k<n; k++ ) {
 83             a.v[i][k] = a.v[i][k] * d % Mod;
 84             b.v[i][k] = b.v[i][k] * d % Mod;
 85         }
 86         for( int j=i-1; j>=0; j-- ) {
 87             d = a.v[j][i] * inv[a.v[i][i]] % Mod;
 88             for( int k=0; k<n; k++ ) {
 89                 a.v[j][k] = (a.v[j][k] - a.v[i][k]*d) % Mod;
 90                 b.v[j][k] = (b.v[j][k] - b.v[i][k]*d) % Mod;
 91                 if( a.v[j][k]<0 ) a.v[j][k] += Mod;
 92                 if( b.v[j][k]<0 ) b.v[j][k] += Mod;
 93             }
 94         }
 95     }
 96     return b;
 97 }
 98 bool operator<( const Matrix &a, const Matrix &b ) {
 99     for( int i=0; i<n; i++ )
100         for( int j=0; j<n; j++ ) {
101             if( a.v[i][j] ^ b.v[i][j] )
102                 return a.v[i][j] < b.v[i][j];
103         }
104     return false;
105 }
106 bool operator==( const Matrix &a, const Matrix &b ) {
107     for( int i=0; i<n; i++ )
108         for( int j=0; j<n; j++ )
109             if( a.v[i][j] ^ b.v[i][j] ) return false;
110     return true;
111 }
112 int ind( Matrix a, Matrix b ) {
113     map<Matrix,int> mp;
114     int m = int(sqrt(Mod))+1;
115     Matrix base = a;
116     a.unit();
117     for( int i=0; i<m; i++ ) {
118         if( a==b && i ) return i;
119         mp[a] = i;
120         a = a*base;
121     }
122     base = ~a;
123     a = b*base;
124     for( int i=m; ; i+=m,a=a*base ) 
125         if( mp.count(a) ) return i+mp[a];
126 }
127  
128 int main() {
129     scanf( "%d%d", &n, &Mod );
130     for( int i=1; i<Mod; i++ )
131         inv[i] = mpow(i,Mod-2);
132     Matrix a, b;
133     a.read();
134     b.read();
135     printf( "%d\n", ind(a,b) );
136 }
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

　　　　完了

posted @ 2015-06-29 20:35  idy002  阅读(674)  评论(0)  编辑  收藏  举报