bzoj4128 Matrix

Description

给定矩阵A,B和模数p,求最小的x满足

A^x = B (mod p)

Input 第一行两个整数n和p，表示矩阵的阶和模数,接下来一个n * n的矩阵A.接下来一个n * n的矩阵B

Output 输出一个正整数，表示最小的可能的x，数据保证在p内有解

BSGS求离散对数，这里换成了矩阵乘法，不过是一样的。对矩阵也求逆就可以了。

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<map>
using namespace std;
#define LL unsigned long long
const int p=131;
int n,mod,m;
int pow(int base,int p)
{
    int ret=1;
    while (p)
    {
        if (p&1) ret=ret*base%mod;
        base=base*base%mod;
        p>>=1;
    }
    return ret;
}
int inverse(int x)
{
    return pow(x,mod-2);
}
struct mat
{
    int a[80][80];
    LL h;
    void init()
    {
        h=0;
        for (int i=0;i<n;i++)
            for (int j=0;j<n;j++)
                h=h*p+a[i][j];
    }
    void I()
    {
        for (int i=0;i<n;i++)
            for (int j=0;j<n;j++)
                a[i][j]=(i==j);
        init();
    }
    void rd()
    {
        for (int i=0;i<n;i++)
            for (int j=0;j<n;j++)
                scanf("%d",&a[i][j]);
        init();
    }
    mat operator * (const mat &m) const
    {
        mat ret;
        for (int i=0;i<n;i++)
            for (int j=0;j<n;j++)
            {
                ret.a[i][j]=0;
                for (int k=0;k<n;k++)
                    ret.a[i][j]=(ret.a[i][j]+a[i][k]*m.a[k][j]%mod)%mod;
            }
        ret.init();
        return ret;
    }
    mat pow(int x)
    {
        mat ret,base;
        ret.I();
        for (int i=0;i<n;i++)
            for (int j=0;j<n;j++)
                base.a[i][j]=a[i][j];
        while (x)
        {
            if (x&1) ret=ret*base;
            base=base*base;
            x>>=1;
        }
        ret.init();
        return ret;
    }
    mat inv()
    {
        mat ret;
        int p,x;
        ret.I();
        for (int i=0;i<n;i++)
        {
            p=-1;
            for (int j=i;j<n;j++)
                if (a[j][i])
                {
                    p=j;
                    break;
                }
            if (p==-1) continue;
            if (p!=i)
                for (int j=0;j<n;j++)
                {
                    swap(a[i][j],a[p][j]);
                    swap(ret.a[i][j],ret.a[p][j]);
                }
            x=inverse(a[i][i]);
            for (int j=0;j<n;j++)
            {
                a[i][j]=a[i][j]*x%mod;
                ret.a[i][j]=ret.a[i][j]*x%mod;
            }
            for (int j=0;j<n;j++)
                if (j!=i)
                {
                    x=a[j][i];
                    for (int k=0;k<n;k++)
                    {
                        a[j][k]=(a[j][k]-a[i][k]*x%mod+mod)%mod;
                        ret.a[j][k]=(ret.a[j][k]-ret.a[i][k]*x%mod+mod)%mod;
                    }
                }
        }
        ret.init();
        return ret;
    }
}a,b,tem,u;
map<LL,int> mp;
int main()
{
    scanf("%d%d",&n,&mod);
    a.rd();
    b.rd();
    m=sqrt(mod);
    tem.I();
    mp[tem.h]=0;
    for (int i=1;i<m;i++)
    {
        tem=tem*a;
        if (!mp.count(tem.h)) mp[tem.h]=i;
    }
    u=(a.pow(m)).inv();
    for (int i=0;;i++)
    {
        if (mp.count(b.h))
        {
            printf("%d\n",i*m+mp[b.h]);
            return 0;
        }
        b=b*u;
    }
}