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

传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=4128

大水题一道

使用大步小步算法，把数字的运算换成矩阵的运算就好了

矩阵求逆？这么基础的线代算法我也不想多说，还是自行百度吧

需要注意的是矩阵没有交换律，所以在计算$Bcdot A^{-m}$的时候不要把顺序搞混

代码：

#include <cstring>
#include <cstdio>
#include <algorithm>
#include <map>
#include <cmath>

using namespace std;
const int MAXN = 70;

int n, p;

struct Matrix {
    int a[MAXN][MAXN];
    int *operator[]( int idx ) {
        return a[idx];
    }
    const int *operator[]( int idx ) const {
        return a[idx];
    }
    bool operator<( const Matrix &rhs ) const { // 用于map
        for( int i = 0; i < n; ++i )
            for( int j = 0; j < n; ++j )
                if( a[i][j] < rhs[i][j] ) return true;
                else if( a[i][j] > rhs[i][j] ) return false;
        return false;
    }
    void unit() { // 填成n阶单位矩阵
        memset( a, 0, sizeof(a) );
        for( int i = 0; i < n; ++i ) a[i][i] = 1;
    }
    void clear() { // 填成全0
        memset( a, 0, sizeof(a) );
    }
};
Matrix A, B; // 题目中给定的A，B矩阵

int pow_mod( int a, int b ) { // 整数快速幂
    int ans = 1;
    while(b) {
        if( b&1 ) ans = ans * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ans;
}
int inv( int a ) { // 整数求逆
    return pow_mod(a,p-2);
}
void input_matrix( Matrix &A ) {
    for( int i = 0; i < n; ++i )
        for( int j = 0; j < n; ++j ) {
            scanf( "%d", &A[i][j] );
            A[i][j] %= p;
        }
}
void mul( const Matrix &A, const Matrix &B, Matrix &C ) { // 矩阵乘法，结果存入C
    Matrix tmp; tmp.clear();
    for( int i = 0; i < n; ++i )
        for( int j = 0; j < n; ++j )
            for( int k = 0; k < n; ++k )
                tmp[i][j] = (tmp[i][j] + A[i][k] * B[k][j] % p) % p;
    C = tmp;
}
void gauss_jordan( int A[MAXN][MAXN<<1] ) { // 高斯约当消元求逆
    for( int i = 0; i < n; ++i ) {
        int r;
        for( r = i; r < n; ++r )
            if( A[r][i] ) break;
        if( r != i )
            for( int j = i; j < (n<<1); ++j )
                swap( A[r][j], A[i][j] );
        int iv = inv( A[i][i] );
        for( int j = i; j < (n<<1); ++j )
            A[i][j] = A[i][j] * iv % p;
        for( int j = 0; j < n; ++j )
            if( j != i && A[j][i] ) {
                int tmp = (p - A[j][i]) % p;
                for( int k = i; k < (n<<1); ++k )
                    A[j][k] = (A[j][k] + A[i][k] * tmp % p) % p;
            }
    }
}
void inv( const Matrix &A, Matrix &B ) { // 矩阵求逆，结果存入B
    int tmp[MAXN][MAXN<<1] = {0};
    for( int i = 0; i < n; ++i )
        for( int j = 0; j < n; ++j )
            tmp[i][j] = A[i][j];
    for( int i = 0; i < n; ++i ) tmp[i][i+n] = 1;
    gauss_jordan(tmp);
    for( int i = 0; i < n; ++i )
        for( int j = 0; j < n; ++j )
            B[i][j] = tmp[i][j+n];
}

map<Matrix,int> mp;
int bsgs( Matrix A, Matrix B ) { // 大步小步算法
    int m = sqrt(p)+1;
    Matrix E; E.unit();
    for( int i = 0; i < m; ++i ) {
        if( !mp.count(E) ) mp[E] = i;
        mul(E,A,E);
    }
    inv(E,E);
    for( int i = 0; i < p; i += m ) {
        if( mp.count(B) ) return i + mp[B];
        mul(B,E,B); // 注意矩阵乘法顺序
    }
    return -1;
}

int main() {
    scanf( "%d%d", &n, &p );
    input_matrix(A), input_matrix(B);
    printf( "%d
", bsgs(A,B) );
    return 0;
}