题解Matrix BZOJ 4128 矩阵求逆 离散对数 大步小步算法

Posted 2020-09-08 mlystdcall

传送门：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 }