Day 170-Q5 — 2026-10-05 | 指数計算法（index calculus）で多数の離散対数（安全素数 p=2q+1。小さい素数の対数を連立方程式で一度だけ求め、各クエリは「滑らかな数」を引くまで乱択）

📋 コピーimport sys
import random


def main():
    data = sys.stdin.buffer.read().split()
    p, g, qn = int(data[0]), int(data[1]), int(data[2])
    hs = [int(x) for x in data[3:3 + qn]]
    q = (p - 1) // 2                     # p = 2q + 1、q は素数
    rng = random.Random(12345)

    # 1. 因子基底: B 以下の素数
    B = 1000
    sieve = bytearray([1]) * (B + 1)
    sieve[0:2] = b"\x00\x00"
    for i in range(2, int(B ** 0.5) + 1):
        if sieve[i]:
            sieve[i * i::i] = bytearray(len(sieve[i * i::i]))
    fb = [i for i in range(2, B + 1) if sieve[i]]
    K = len(fb)
    P = 1
    for x in fb:
        P *= x
    E = p.bit_length()                   # r < p < 2^E なので、素因数の重複は E 回まで

    def factor(r):
        # r が B-smooth なら指数ベクトル（疎: {番号: 指数}）、そうでなければ None
        if pow(P % r, E, r) != 0:
            return None
        vec = {}
        for idx, f in enumerate(fb):
            if r % f == 0:
                c = 0
                while r % f == 0:
                    r //= f
                    c += 1
                vec[idx] = c
                if r == 1:
                    break
        return vec

    # 2. 関係式を集め、F_q 上で消去する（行を 1 本ずつ入れて、ピボット行を保つ）
    piv = [None] * K                     # piv[c]: 列 c がピボットの行（係数リスト, 右辺）
    rank = 0
    k = rng.randrange(1, p - 1)
    r = pow(g, k, p)
    while rank < K:
        k += 1
        r = r * g % p                    # g^k を 1 回の掛け算で更新する
        vec = factor(r)
        if vec is None:
            continue
        row = [0] * K
        for idx, c in vec.items():
            row[idx] = c
        rhs = k % q
        for c in range(K):
            if row[c] == 0:
                continue
            if piv[c] is None:
                inv = pow(row[c], q - 2, q)
                row = [x * inv % q for x in row]
                rhs = rhs * inv % q
                piv[c] = (row, rhs)
                rank += 1
                break
            prow, prhs = piv[c]
            t = row[c]
            row = [(x - t * y) % q for x, y in zip(row, prow)]
            rhs = (rhs - t * prhs) % q
    # 後退代入で log(fb[c]) mod q を求める
    lg = [0] * K
    for c in range(K - 1, -1, -1):
        row, rhs = piv[c]
        s = rhs
        for d in range(c + 1, K):
            if row[d]:
                s -= row[d] * lg[d]
        lg[c] = s % q

    # 3. 各クエリ: 偶奇はオイラーの規準、mod q は h·g^k を smooth にして足し合わせる
    out = []
    for h in hs:
        par = 0 if pow(h, q, p) == 1 else 1
        k = rng.randrange(0, p - 1)
        r = h * pow(g, k, p) % p
        while True:
            vec = factor(r)
            if vec is not None:
                break
            k += 1
            r = r * g % p
        xq = (sum(c * lg[idx] for idx, c in vec.items()) - k) % q
        x = xq if xq % 2 == par else xq + q       # CRT: x ≡ xq (mod q), x ≡ par (mod 2)
        out.append(x)
    print("\n".join(map(str, out)))


main()