【题解】Matrix BZOJ 4128 矩阵求逆 离散对数 大步小步算法
传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=4128
大水题一道
使用大步小步算法，把数字的运算换成矩阵的运算就好了
矩阵求逆？这么基础的线代算法我也不想多说，还是自行百度吧
需要注意的是矩阵没有交换律，所以在计算$B\cdot A^{-m}$的时候不要把顺序搞混
代码：
1 #include <cstring> 2 #include <cstdio> 3 #include <algorithm> 4 #include <map> 5 #include <cmath> 6 7 using namespace std; 8 const int MAXN = 70; 9 10 int n, p; 11 12 struct Matrix { 13 int a[MAXN][MAXN]; 14 int *operator[]( int idx ) { 15 return a[idx]; 16 } 17 const int *operator[]( int idx ) const { 18 return a[idx]; 19 } 20 bool operator<( const Matrix &rhs ) const { // 用于map 21 for( int i = 0; i < n; ++i ) 22 for( int j = 0; j < n; ++j ) 23 if( a[i][j] < rhs[i][j] ) return true; 24 else if( a[i][j] > rhs[i][j] ) return false; 25 return false; 26 } 27 void unit() { // 填成n阶单位矩阵 28 memset( a, 0, sizeof(a) ); 29 for( int i = 0; i < n; ++i ) a[i][i] = 1; 30 } 31 void clear() { // 填成全0 32 memset( a, 0, sizeof(a) ); 33 } 34 }; 35 Matrix A, B; // 题目中给定的A，B矩阵 36 37 int pow_mod( int a, int b ) { // 整数快速幂 38 int ans = 1; 39 while(b) { 40 if( b&1 ) ans = ans * a % p; 41 a = a * a % p; 42 b >>= 1; 43 } 44 return ans; 45 } 46 int inv( int a ) { // 整数求逆 47 return pow_mod(a,p-2); 48 } 49 void input_matrix( Matrix &A ) { 50 for( int i = 0; i < n; ++i ) 51 for( int j = 0; j < n; ++j ) { 52 scanf( "%d", &A[i][j] ); 53 A[i][j] %= p; 54 } 55 } 56 void mul( const Matrix &A, const Matrix &B, Matrix &C ) { // 矩阵乘法，结果存入C 57 Matrix tmp; tmp.clear(); 58 for( int i = 0; i < n; ++i ) 59 for( int j = 0; j < n; ++j ) 60 for( int k = 0; k < n; ++k ) 61 tmp[i][j] = (tmp[i][j] + A[i][k] * B[k][j] % p) % p; 62 C = tmp; 63 } 64 void gauss_jordan( int A[MAXN][MAXN<<1] ) { // 高斯约当消元求逆 65 for( int i = 0; i < n; ++i ) { 66 int r; 67 for( r = i; r < n; ++r ) 68 if( A[r][i] ) break; 69 if( r != i ) 70 for( int j = i; j < (n<<1); ++j ) 71 swap( A[r][j], A[i][j] ); 72 int iv = inv( A[i][i] ); 73 for( int j = i; j < (n<<1); ++j ) 74 A[i][j] = A[i][j] * iv % p; 75 for( int j = 0; j < n; ++j ) 76 if( j != i && A[j][i] ) { 77 int tmp = (p - A[j][i]) % p; 78 for( int k = i; k < (n<<1); ++k ) 79 A[j][k] = (A[j][k] + A[i][k] * tmp % p) % p; 80 } 81 } 82 } 83 void inv( const Matrix &A, Matrix &B ) { // 矩阵求逆，结果存入B 84 int tmp[MAXN][MAXN<<1] = {0}; 85 for( int i = 0; i < n; ++i ) 86 for( int j = 0; j < n; ++j ) 87 tmp[i][j] = A[i][j]; 88 for( int i = 0; i < n; ++i ) tmp[i][i+n] = 1; 89 gauss_jordan(tmp); 90 for( int i = 0; i < n; ++i ) 91 for( int j = 0; j < n; ++j ) 92 B[i][j] = tmp[i][j+n]; 93 } 94 95 map<Matrix,int> mp; 96 int bsgs( Matrix A, Matrix B ) { // 大步小步算法 97 int m = sqrt(p)+1; 98 Matrix E; E.unit(); 99 for( int i = 0; i < m; ++i ) { 100 if( !mp.count(E) ) mp[E] = i; 101 mul(E,A,E); 102 } 103 inv(E,E); 104 for( int i = 0; i < p; i += m ) { 105 if( mp.count(B) ) return i + mp[B]; 106 mul(B,E,B); // 注意矩阵乘法顺序 107 } 108 return -1; 109 } 110 111 int main() { 112 scanf( "%d%d", &n, &p ); 113 input_matrix(A), input_matrix(B); 114 printf( "%d\n", bsgs(A,B) ); 115 return 0; 116 }
