Day 088-Q1 — Pohlig-Hellman法（離散対数問題・$p-1$ の滑らかさを利用）

問題

素数 $p$ と、$\bmod p$ における原始根 $g$、および整数 $y$（$1 \le y \le p-1$）が与えられる。$g^x \equiv y \pmod p$ を満たす最小の非負整数 $x$（$0 \le x \le p-2$）を求めよ。$p-1$ の素因数はすべて $10^6$ 以下であることが保証される。

入出力例

23 5 17

7

41 6 1

0

例1: $5^7 \bmod 23 = 17$、$p-1=22=2\times11$。例2: $g^0=1$ が最小解。

概念図: $p-1$ を素因数分解して部分群に分割し CRT で統合

ヒント

def pohlig_hellman_prime_power(g, y, q, e, p):
    gamma = pow(g, q ** (e - 1), p)  # 位数 q の元
    x = 0
    g_inv = pow(g, -1, p)
    for k in range(e):
        exp = q ** (e - 1 - k)
        yk = pow((y * pow(g_inv, x, p)) % p, exp, p)
        xk = bsgs_small(gamma, yk, p, q)
        x += xk * (q ** k)
    return x

模範解答

import sys
import math

def factorize(n):
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors

def bsgs_small(g, h, p, order):
    m = math.isqrt(order) + 1
    table = {}
    e = 1
    for j in range(m):
        table.setdefault(e, j)
        e = e * g % p
    ginv_m = pow(g, -m, p)
    cur = h % p
    for i in range(m):
        if cur in table:
            x = i * m + table[cur]
            if x < order:
                return x
        cur = cur * ginv_m % p
    return -1

def pohlig_hellman_prime_power(g, y, q, e, p):
    gamma = pow(g, q ** (e - 1), p)
    x = 0
    g_inv = pow(g, -1, p)
    for k in range(e):
        exp = q ** (e - 1 - k)
        yk = pow((y * pow(g_inv, x, p)) % p, exp, p)
        xk = bsgs_small(gamma, yk, p, q)
        x += xk * (q ** k)
    return x

def crt_combine(residues, moduli):
    x, m = 0, 1
    for r, mi in zip(residues, moduli):
        t = ((r - x) * pow(m, -1, mi)) % mi
        x = x + m * t
        m *= mi
    return x % m

def solve():
    p, g, y = map(int, sys.stdin.readline().split())
    n = p - 1
    factors = factorize(n)
    residues, moduli = [], []
    for q, e in factors.items():
        qe = q ** e
        gi = pow(g, n // qe, p)
        yi = pow(y, n // qe, p)
        xi = pohlig_hellman_prime_power(gi, yi, q, e, p)
        residues.append(xi)
        moduli.append(qe)
    x = crt_combine(residues, moduli)
    print(x)

solve()

計算量: 因数分解 $O(B)$（$B=10^6$）。各成分 $O(e_i\sqrt{q_i})$。全体で単純 BSGS の $O(\sqrt p)$ より大幅高速。

Step-by-Step 解説

Step 1: 群の位数を素因数分解する

$(\mathbb{Z}/p\mathbb{Z})^*$ の位数 $p-1$ を $q_1^{e_1}\cdots q_k^{e_k}$ に分解。$q_i\le10^6$ なので試し割りで十分高速。

Step 2: 各素数冪の部分群に帰着

$g_i=g^{(p-1)/q_i^{e_i}}$（位数 $q_i^{e_i}$）、$y_i=y^{(p-1)/q_i^{e_i}}$ とし $g_i^x=y_i$ の $x\bmod q_i^{e_i}$ を求める。

Step 3: 桁ごとの Pohlig-Hellman

Step 4: CRT で統合

各成分の $x\bmod q_i^{e_i}$ が求まったら中国剰余定理で $x\bmod(p-1)$ を復元する。

よくあるミス

次のステップ

• 発展問題: $p-1$ が滑らかでない場合の Index Calculus 法との併用
• 高速化: 因数分解を事前計算した素数篩で行う

自己評価