Day 173-Q5 — 2026-10-09 | 安全素数でない p の離散対数（大きい素因数 q は指数計算法、残りは Pohlig–Hellman、CRT で合わせる）

📋 コピーimport sys
import random
from math import isqrt


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


def factor_small(n, ps, prod):
    # n が ps だけで割り切れる（B-smooth）なら指数の辞書、そうでなければ None
    if pow(prod, 40, n) != 0:                  # n < 2^40 なので、smooth ⇔ n | prod^40
        return None
    e = {}
    for i, q in enumerate(ps):
        if n % q == 0:
            c = 0
            while n % q == 0:
                n //= q
                c += 1
            e[i] = c
            if n == 1:
                return e
    return e if n == 1 else None


def bsgs(g, h, order, p):
    # g の位数が order の部分群で g^x = h を解く（x in [0, order)）
    m = isqrt(order) + 1
    table = {}
    cur = 1
    for j in range(m):
        table.setdefault(cur, j)
        cur = cur * g % p
    step = pow(g, (p - 2) * m, p)              # g^{-m}
    cur = h
    for i in range(m + 1):
        if cur in table:
            return (i * m + table[cur]) % order
        cur = cur * step % p
    return None


def solve_mod_q(rows, rhs, k, q):
    # rows[r][i] x_i = rhs[r] (mod q) を掃き出し、決まった変数だけ返す
    a = [row[:] + [r] for row, r in zip(rows, rhs)]
    where = [-1] * k
    r = 0
    for c in range(k):
        piv = next((i for i in range(r, len(a)) if a[i][c] % q), None)
        if piv is None:
            continue
        a[r], a[piv] = a[piv], a[r]
        inv = pow(a[r][c], q - 2, q)
        a[r] = [x * inv % q for x in a[r]]
        for i in range(len(a)):
            if i != r and a[i][c] % q:
                f = a[i][c]
                ai, ar = a[i], a[r]
                a[i] = [(x - f * y) % q for x, y in zip(ai, ar)]
        where[c] = r
        r += 1
    sol = {}
    for c in range(k):
        if where[c] >= 0:
            row = a[where[c]]
            if all(row[j] % q == 0 for j in range(k) if j != c):   # 自由変数に依存しない行だけ使う
                sol[c] = row[k] % q
    return sol


def main():
    data = sys.stdin.read().split()
    p, g, t = int(data[0]), int(data[1]), int(data[2])
    hs = [int(x) for x in data[3:3 + t]]
    rng = random.Random(12345)
    n = p - 1
    # p-1 = (小さい素因数の積) × q（q は 10^4 より大きい素数、高々 1 個）
    fac = {}
    m = n
    for r in small_primes(10 ** 4):
        while m % r == 0:
            fac[r] = fac.get(r, 0) + 1
            m //= r
    q = m if m > 1 else None

    # ---- 大きい素数 q の成分: 指数計算法で log mod q ----
    if q is not None:
        B = 1000
        ps = small_primes(B)
        k = len(ps)
        prod = 1
        for r in ps:
            prod *= r
        rows, rhs = [], []
        while len(rows) < k + 20:
            e = rng.randrange(1, n)
            ex = factor_small(pow(g, e, p), ps, prod)
            if ex is not None:
                row = [0] * k
                for i, c in ex.items():
                    row[i] = c
                rows.append(row)
                rhs.append(e % q)
        logq = solve_mod_q(rows, rhs, k, q)

    out = []
    for h in hs:
        res, mods = [], []
        # 小さい素数冪 r^a の成分: Pohlig–Hellman（桁ごとに BSGS）
        for r, a in fac.items():
            gr = pow(g, n // r, p)              # 位数 r の元
            x = 0
            ra = 1
            for j in range(a):
                hj = pow(h * pow(g, n - x, p) % p, n // (ra * r), p)
                d = bsgs(gr, hj, r, p)
                x += d * ra
                ra *= r
            res.append(x)
            mods.append(r ** a)
        if q is not None:
            while True:                         # h * g^s が B-smooth で、出てくる素数の log が全部分かるまで
                s = rng.randrange(n)
                ex = factor_small(h * pow(g, s, p) % p, ps, prod)
                if ex is not None and all(i in logq for i in ex):
                    res.append((sum(c * logq[i] for i, c in ex.items()) - s) % q)
                    mods.append(q)
                    break
        x, md = 0, 1                            # CRT
        for r_, m_ in zip(res, mods):
            t_ = (r_ - x) * pow(md, -1, m_) % m_
            x += md * t_
            md *= m_
        out.append(x % n)
    print('\n'.join(map(str, out)))


if __name__ == "__main__":
    main()