bzoj4128: Matrix | jeffswt: the ac moments

Description

给定矩阵 \(A, B\) 和模数 \(p\), 求最小的 \(x\) 满足 \(A^x = B (\bmod p)\).

Input

第一行两个整数 \(n\) 和 \(p\), 表示矩阵的阶和模数, 接下来一个 \(n \times n\) 的矩阵 \(A\). 接下来一个 \(n \times n\) 的矩阵 \(B\).

Output

输出一个正整数, 表示最小的可能的 \(x\), 数据保证在 \(p\) 内有解.

Sample Input

Sample Output

Data Range

对于 100% 的数据,\(n \leq 70, p \leq 19997, 0 \leq A_{i,j}, B_{i,j} \lt p\), 且满足 \(p\) 为质数.

同时保证 \(A\) 有逆矩阵.

Explanation

这道题用的大小步算法 (又称拔山盖世算法).

原文地址:http://blog.csdn.net/clover_hxy/article/details/50683832 (这一篇讲的是对于整数的大小步算法)

先把 x=i*m-j, 其中 m=ceil (sqrt (C)) , (ceil 是向上取整).

这样原式就变为 A^ (i*m-j)=B (mod C),

再变为 A^j×B=A (m*i) (mod C).

枚举 j (范围 0-m), 将 A^j×B 存入 hash 表

枚举 i (范围 1-m), 从 hash 表中寻找第一个满足 A^j×B=A (m*i) (mod C).

此时 x=i*m-j 即为所求.

在网上看到的其他题解大多用的是 x=im+j, 也可以做, 只是会牵扯的求逆元, 所以比较麻烦. 使 x=im-j 就可以轻松避免这个问题了.

那么肯定有人会有疑问为何只计算到 m=ceil (sqrt (C)) 就可以确定答案呢?

x=i*m-j 也就是 x 的最大值不会超过 p, 那超过 p 的怎么办?

有一个公式 a^ (k mod p)=a^k (mod p) 这个公式的推导需要用到费马小定理

k mod p 可以看做 k-mp, 原式可化成 a^{k/ (a}p)^m=ak (mod p)

根据费马小定理 a^p=1 (mod p) 其中 p 为质数, a, p 互质, 可得 a^k/1m=a^k (mod p) a^k=ak (mod p) 得证.

原文地址:http://blog.csdn.net/Lcomyn/article/details/46714569 (这篇内容和下一篇差不多, 而且没有那个详细)

原文地址:http://blog.csdn.net/xym_csdn/article/details/50675542

枚举 j, 接下来有两种方法, 一个是存哈希表, 但蒟蒻觉得太麻烦, 所以就学习了一下随机化, 随便搞一个 n×1 的矩阵 C, 在等式两边同乘, 变成 (A^m) i×C=B×A^j×C (注意顺序! 不满足交换律!) , 然后枚举 i, 和得到的所有结果比较, 碰到就可以输出了

主过程时间复杂度: 存结果时 O (2n^{2√p) 判断时 O (√p (n}2+√pn)) , 好像是这个吧

能让 n^{3 变成 n}2, 原本跑 1. 5s 的程序就可以优化到 0. 6s 左右了==, 而且比 hash 快一些

Source Code


#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>

#define rep(_var,_begin,_end) for(int _var=_begin;_var<=_end;_var++)
using namespace std;
typedef long long lli;
const int maxn = 80, maxm = 150; // maxm^2 > p

int n, m, modr;

struct matrix {
    lli val[maxn][maxn];
    matrix& operator *= (matrix b) {
        static lli c[maxn][maxn];
        rep(i, 1, n) rep(j, 1, n) {
            c[i][j] = 0;
            rep(k, 1, n) c[i][j] += val[i][k] * b.val[k][j];
            c[i][j] %= modr;
        }
        rep(i, 1, n) rep(j, 1, n)
            val[i][j] = c[i][j];
        return *this;
    }
    matrix(int v) {
        rep(i, 1, n) rep(j, 1, n)
            val[i][j] = i == j ? v : 0;
        return ;
    }
    matrix(void) {
        rep(i, 1, n) rep(j, 1, n)
            val[i][j] = 0;
        return ;
    }
    void print(void) {
        rep(i, 1, n) rep(j, 1, n)
            printf("%d ", val[i][j]);
        printf("\n");
        return ;
    }
};

matrix operator ^ (matrix a, lli powr)
{
    matrix res(1);
    for (powr; powr > 0; powr >>= 1, a *= a)
        if (powr & 1)
            res *= a;
    return res;
}

// A magic comparator that compares the hashes between hashed multiplication
// results
bool operator == (matrix a, matrix b)
{
    rep(i, 1, n) if (a.val[i][1] != b.val[i][1])
        return false;
    return true;
}

// A magic multiplication that hashes the matrix multiplication, avoiding
// excessive computation as we don't require the exact results
matrix operator * (matrix a, matrix b)
{
    matrix c;
    rep(i, 1, n) {
        rep(k, 1, n) c.val[i][1] += a.val[i][k] * b.val[k][1];
        c.val[i][1] %= modr;
    }
    return c;
}

matrix hash[maxm];

int main(int argc, char** argv)
{
    scanf("%d%d", &n, &modr);
    matrix a, b, c__, c;
    rep(i, 1, n) rep(j, 1, n)
        scanf("%lld", &a.val[i][j]);
    rep(i, 1, n) rep(j, 1, n)
        scanf("%lld", &b.val[i][j]);
    // Baby step giant step
    m = sqrt(modr);
    rep(i, 1, n)
        c__.val[i][1] = i % modr;
    c = c__;
    for (int i = 0; i <= m; i++) {
        hash[i] = b * c;
        c = a * c;
    }
    a = a ^ m;
    c = c__;
    int res = 0;
    for (int i = 1; i <= m; i++) {
        c = a * c;
        for (int j = 0; j <= m; j++)
            if (hash[j] == c) {
                res = i * m - j;
                i = m + 1;
                break;
            }
    }
    printf("%d\n", res);
    return 0;
}