2019杭电多校第五场 I. discrete logarithm problem

Aug 5, 2019 由 小羊

题目链接：HDU 6632

已知数字 $a,b,p$，并且 $p=2^a3^b+1$。想求解离散对数 $a^x \equiv b \mod p$。

需要一个东西叫做 Pohlig-Hellman 算法。以下是来自维基百科内容的翻译，将使用群论语言。

p^e 阶循环群

这个算法用于计算 $p^e$ 阶循环群 $G=\langle g \rangle$ 的离散对数 $g^x = h$。

1. 初始化 $x_0 := 0$

2. 计算 $\gamma := g^{p^{e-1}}$，由拉格朗日配集分解计数定理知，$ord~\gamma = p$。

3. 对于 $k \in \lbrace 0,\cdots,e-1 \rbrace$，分别做

   1. 计算 $h_k := (g^{-x_k}h)^{p^{e-1-k}}$。由此构造可知，此元素的阶一定整除 $p$，即阶为 $1$ 或 $p$，因此 $h_k \in \langle \gamma \rangle$。

   2. 使用 BSGS 算法，计算 $d_k \in \lbrace 0, \cdots, p-1 \rbrace$，其中 $\gamma^{d_k} = h_k$。这个操作的复杂度为 $O(\sqrt{p})$

   3. 令 $x_{k+1} := x_k + p^k d_k$

4. 返回 $x_e$

当 $e \ll p$ 时，这个算法优于 BSGS，复杂度为 $O(e\sqrt{p})$。

一般合数阶循环群

首先我们将 $G = \langle g \rangle$ 的阶数刻画为 $n = \prod_{i=1}^r {p_i}^{e_i}$。现在我们解方程 $g^x = h$。

1. 对循环群阶数做因数分解

2. 对于 $i \in \lbrace 1, \cdots, r \rbrace$，分别做

   1. 计算 $g_i := g^{n / {p_i}^{e_i}}$，由拉格朗日定理知此元素的阶为 ${p_i}^{e_i}$。

   2. 计算 $h_i := h^{n/{p_i}^{e_i}}$，由此构造知 $h_i \in \langle g_i \rangle$

   3. 利用上面算法计算 $x_i$ 使得 ${g_i}^{x_i} = h_i$

3. 利用中国剩余定理解方程组 $x \equiv x_i \mod {p_i}^{e_i}$

4. 返回 $x$

因此总的时间复杂度为

$T(n) = O\left( \sum_{i=1}^r e_i (\log n + \sqrt{p_i}) \right)$

回到本题

显然因为 $p_i = 2, 3$，我们不用 BSGS……好的。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;

lld fpow(lld _a, lld k, lld MOD) {
    __int128 b = 1, a = _a % MOD;
    while (k) {
        if (k & 1) b = b * a % MOD;
        a = a * a % MOD;
        k >>= 1;
    }
    return b;
}

__int128 extgcd(__int128 a, __int128 b, __int128 &x, __int128 &y) {
    if (!b) return x = 1, y = 0, a;
    __int128 d = extgcd(b, a % b, y, x);
    return y -= (a / b) * x, d;
}

int T, e2, e3, G;
lld MOD, n, p2, p3;

lld solvePE(lld g, lld h, lld p, lld pe, int e)
{
    lld x = 0, pp = 1;
    lld y = fpow(g, pe/p, MOD);

    for (int k = 0; k < e; k++)
    {
        __int128 hk = fpow(g, pe-x, MOD);
        hk = hk * h % MOD;
        hk = fpow(hk, pe/p/pp, MOD);
        __int128 d, yy;
        for (d = 0, yy = 1; d < p; d++, yy = yy * y % MOD)
            if (yy == hk) break;
        x = (x + pp * d) % MOD, pp *= p;
    }

    return x;
}

__int128 solveN(lld h)
{
    if (e3 == 0) return solvePE(G, h, 2, p2, e2);
    else if (e2 == 0) return solvePE(G, h, 3, p3, e3);
    __int128 x2 = solvePE(fpow(G, p3, MOD), fpow(h, p3, MOD), 2, p2, e2);
    __int128 x3 = solvePE(fpow(G, p2, MOD), fpow(h, p2, MOD), 3, p3, e3);
    __int128 invp2, invp3; assert(extgcd(p2, p3, invp2, invp3) == 1);
    return (x2 * p3 * invp3 + x3 * p2 * invp2) % (MOD-1);
}

int main()
{
    cin >> T;

    while (T--)
    {
        long long a, b;
        cin >> MOD >> a >> b;
        n = MOD-1, e2 = e3 = 0, p2 = p3 = 1;
        while (n % 2 == 0) n /= 2, p2 *= 2, e2++;
        while (n % 3 == 0) n /= 3, p3 *= 3, e3++;
        for (G = 2; (e2 > 0 && fpow(G, (MOD-1)/2, MOD) == 1)
                 || (e3 > 0 && fpow(G, (MOD-1)/3, MOD) == 1); G++);
        __int128 inda = solveN(a), indb = solveN(b), x, y;
        __int128 d = extgcd(inda, MOD-1, x, y);
        lld ans = indb % d ? -1
            : lld((indb / d * x % (MOD-1) + (MOD-1)) % ((MOD-1)/d));
        cout << ans << endl;
    }

    return 0;
}

我竟然都忘了 extgcd 是怎么求多解的，我好菜啊QAQ