Day 169-Q5 — 2026-10-04 | Pohlig–Hellman（p−1 が小さい素因数だけなら、10^18 の素数でも離散対数が解ける。g が原始根でないときの最小解）

📋 コピーimport sys
from math import isqrt


def sieve(n):
    s = bytearray([1]) * (n + 1)
    s[0] = s[1] = 0
    for i in range(2, isqrt(n) + 1):
        if s[i]:
            s[i * i::i] = bytearray(len(range(i * i, n + 1, i)))
    return [i for i in range(n + 1) if s[i]]


PRIMES = sieve(100000)


def bsgs(a, b, q, p):
    # 位数 q の元 a について a^x = b (0 <= x < q) を解く。必ず解がある前提
    m = isqrt(q) + 1
    table = {}
    e = 1
    for j in range(m):
        table.setdefault(e, j)
        e = e * a % p
    step = pow(a, -m, p)       # a^{-m}
    y = b
    for i in range(m):
        j = table.get(y)
        if j is not None:
            return i * m + j
        y = y * step % p
    raise ValueError("no log in subgroup")


def solve(p, g, h):
    g %= p
    h %= p
    if p == 2:
        return 0                # g = h = 1
    # p - 1 を素因数分解（素因数はすべて 1e5 以下の保証）
    n = p - 1
    fac = []
    for q in PRIMES:
        if q * q > n and n > 1 and n <= 100000:
            fac.append((n, 1))
            n = 1
        if n == 1:
            break
        if n % q == 0:
            e = 0
            while n % q == 0:
                n //= q
                e += 1
            fac.append((q, e))
    # g の位数 ord を p-1 から素因数を削って求める
    order = p - 1
    ofac = []
    for q, e in fac:
        for _ in range(e):
            if pow(g, order // q, p) == 1:
                order //= q
            else:
                break
        cnt = 0
        t = order
        while t % q == 0:
            t //= q
            cnt += 1
        if cnt:
            ofac.append((q, cnt))
    # 巡回群なので h が <g> に入る ⇔ h^ord = 1
    if pow(h, order, p) != 1:
        return -1
    # Pohlig–Hellman: 各 q^e について x mod q^e を q 進の桁ごとに求め、CRT でまとめる
    x_all, mod_all = 0, 1
    for q, e in ofac:
        qe = q ** e
        g0 = pow(g, order // qe, p)          # 位数 q^e
        h0 = pow(h, order // qe, p)
        gamma = pow(g0, q ** (e - 1), p)     # 位数 q
        ginv = pow(g0, -1, p)
        x = 0
        qk = 1
        for k in range(e):
            hk = pow(h0 * pow(ginv, x, p) % p, q ** (e - 1 - k), p)
            d = bsgs(gamma, hk, q, p)
            x += d * qk
            qk *= q
        # CRT: x_all mod mod_all と x mod qe を合わせる
        t = (x - x_all) * pow(mod_all, -1, qe) % qe
        x_all += mod_all * t
        mod_all *= qe
    return x_all % order if order > 1 else 0


def main():
    data = sys.stdin.buffer.read().split()
    t = int(data[0])
    out = []
    for i in range(t):
        p, g, h = int(data[1 + 3 * i]), int(data[2 + 3 * i]), int(data[3 + 3 * i])
        out.append(str(solve(p, g, h)))
    print("\n".join(out))


main()