BSGS 求离散对数

发表于 分类于 ，

概述:

BSGS (Baby-Step-Giant-Step) 算法解决的是这样一个问题：当 \((a,p)=1\) 时，求最小的 \(x\), 使得 \(a^x \equiv b(mod\ p)\). 这个 \(x\) 即为离散对数 \(log_ab\).

先考虑最简单的情况，即 \(p\) 为素数.

根据费马小定理，这时只需检查 \(0, 1, ..., p-1\) 是否为解即可。因为 \(a^{p-1} \equiv 1(mod\ p)\), 当 \(x\) 超过 \(p-1\) 时 \(a^x\) 就循环了.

直接检查复杂度为 \(O(p)\), 显然不太够.

实现:

可以利用分块的思想。设置一个块的大小 \(m\). 对于前 \(m\) 项，即 \(a^0, a^1, ..., a^{m-1}\), 直接判断是否为 b. 若相等则直接返回， 否则将 \(a^i\ mod\ p\) 存入 \(e_i\) 之中。最后求出 \(a^m\) 的逆 \(a^{-m}\).

然后考虑下一个块 \(a^m, a^{m+1}, ..., a^{2m-1}\), 若它们之中有解，则说明存在 \(e_i \times a^m \equiv b(mod\ p)\), 左乘 \(a^{-m}\) 即有 \(e_i \equiv b \times a^{-m}(mod\ p)\), 即我们只需要检查 \(e_i\) 中是否存在 \(b \times a^{-m}\) 即可.

以此类推，在第一个块求得 \(e_i\) 之后， 后面的块只需要检查 \(e_i\) 之中是否存在 \(b \times a^{-im}\) 即可.

若使用 std::map 存储 \(e_i\), 时间复杂度即为 \(O((m+p/m)log\ p)\), 当 \(m = p^{1/2}\) 时时间最优， 此时时间复杂度为 \(O(p^{1/2}log\ p)\), 若采用哈希表存储 \(e_i\) 的话平均时间复杂度即为 \(O(\sqrt{p})\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
int BSGS(int a,int p,int b) {
    int m = ceil(sqrt(p));
    unordered_map<int,int> e;
    e[1] = 0;   //a^0 = 1
    int x = a%p;
    for(int i(1);i<m;++i) {
        if(x==b) {
            return i;
        }
        if(!e.count(x)) {
            e[x] = i;   //a^i = x
        }
        x = x*a%p;
    }
    int anm = inv(x,p); //inv为求逆元
    x = (b*anm)%p;
    for(int i(1);i<=m;++i) {
        if(e.count(x)) {
            return i*m + e[x];
        }
        x = (x*anm)%p;
    }
    return -1;
}

当 \(p\) 不是素数时解法大致相同， 只是证明部分的费马小定理换成更为普适的欧拉定理，且上界为 \(\phi(p)\). 同时求逆元时使用扩展欧几里得算法.

(实际实现时上界仍然为 p, 因为求 \(\phi(p)\) 的时间复杂度较高)

扩展:

当 \((a,p) \ne 1\) 时， 上面的普通 BSGS 将会失效，因为此时无法求得逆元.

此时需要将其转化为 \((a,p) = 1\) 的形式.

首先给出一个定理:

\[ \frac{a}{d} \equiv \frac{b}{d} (mod\ \frac{m}{d})\ (1)\]

其中 \(d | (a,b,m)\)(证明见最下).

若 \(d \nmid b\), 除了 \(b=1\) 的情况下 \(x=0\), 否则根据裴蜀定理，原式无解，返回 - 1. (\(a \cdot a^{x-1} +py = b\)) 否则可以利用这个公式将原式:

\[ a^x \equiv b (mod\ p) \]

化为:

\[ a^{x-1} \cdot \frac{a}{d} \equiv \frac{b}{d} (mod\ \frac{p}{d}) \]

其中 \(d = (a,p)\).

继续变换，可得:

\[ a^{x-k} \cdot \frac{a^k}{\prod^{k}_{i}d_i} \equiv \frac{b}{\prod^{k}_{i}d_i}(mod\ \frac{p}{\prod^{k}_{i}d_i}) \]

其中 \(d_i\) 为每次变换时的 \((a,\frac{p}{\prod^{i-1}_{j}d_j})\)

此时 \((a,\frac{p}{\prod^{k}_{i}d_i})=1\). 令 \(c = \frac{a^k}{\prod^{k}_{i}d_i}\), \(b' = \frac{b}{\prod^{k}_{i}d_i}\), \(p' = \frac{p}{\prod^{k}_{i}d_i}\), 上式可写为:

\[ \begin{aligned} & a^{x-k} \cdot c \equiv b'(mod\ p') \\ \Rightarrow & a^{x-k} \equiv b' \cdot c^{-1}(mod\ p') \end{aligned} \]

这样，我们就把 \((a,p) \neq 1\) 的情况转化为了 \((a,p) = 1\) 的情况，然后用普通的 BSGS 解决就行了.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
int exbsgs(int a,int p,int b) {
    if(b==1) {
        return 0;
    }
    int c = 1,g,k(0);
    while((g=gcd(a,p))!=1) {    //将(a,p) != 1的情况化归
        if(b%g!=0) {
            return -1;
        }
        c = c*(a/g)%p;
        b /= g;
        p /= g;
        ++k;
        if(b == 1) {
            return k;
        }
    }
    b = b*inv(c,p)%p;
    /* 普通BSGS流程 */
    int m = ceil(sqrt(p));
    unordered_map<int,int> e;
    e[1] = 0;
    int x = a%p;
    for(int i(1);i<m;++i) {
        if(x==b) {
            return i+k; //记得答案要+k, 因为普通BSGS求得的是x-k
        }
        if(!e.count(x)) {
            e[x] = i;
        }
        x = x*a%p;
    }
    int anm = inv(x,p);
    x = (b*anm)%p;
    for(int i(1);i<=m;++i) {
        if(e.count(x)) {
            return i*m + e[x] + k;  //同上
        }
        x = (x*anm)%p;
    }
    return -1;
}

证明:

首先是公式 \((1)\):

\[ \begin{aligned} & a \equiv b(mod\ m) \\ \Rightarrow & a-b = km (k \in \mathbb{Z}) \\ \Rightarrow & \frac{a}{d} - \frac{b}{d} = k \cdot \frac{m}{d} \\ \Rightarrow & \frac{a}{d} \equiv \frac{b}{d} (mod\ \frac{m}{d}) \end{aligned} \]

接着是文字证明为何 \(c\) 在 \(mod\ p'\) 意义下有逆元:

通过化归的过程，有: \((a,p') = 1\)

这意味着对 \(a,p'\) 进行质因数分解后， 它们两个没有相同的质因数.

同时 \(c = \frac{a^k}{\prod^{k}_{i}d_i}\), 不会从 \(a\) 中产生新的质因数，所以 \(c,p'\) 也没有相同的质因数，\((c,p') = 1\), \(c\) 在 \(mod\ p'\) 意义下逆元存在.