Day 176-Q5 — 2026-10-12 | 合成数を法とする平方根（Hensel の持ち上げと中国剰余定理で、解の個数と解の 1 つを求める）

📋 コピーimport sys


def factorize(m):
    fs = []
    d = 2
    while d * d <= m:
        if m % d == 0:
            k = 0
            while m % d == 0:
                m //= d
                k += 1
            fs.append((d, k))
        d += 1 if d == 2 else 2
    if m > 1:
        fs.append((m, 1))
    return fs


_TS = {}


def sqrt_mod_p(a, p):
    """奇素数 p、a は p と互いに素な平方剰余。p-1 = q·2^s の s が小さければ Tonelli–Shanks、
    大きければ Cipolla（Tonelli は最悪 s^2/2 回の 2 乗が要り、s=30 で遅い）"""
    if p % 4 == 3:
        return pow(a, (p + 1) // 4, p)
    if p not in _TS:
        q, s = p - 1, 0
        while q % 2 == 0:
            q //= 2
            s += 1
        z = 2
        while pow(z, (p - 1) // 2, p) != p - 1:
            z += 1
        _TS[p] = (q, s, pow(z, q, p))
    q, s, c = _TS[p]
    if s > 24:
        return cipolla(a, p)
    m, t, r = s, pow(a, q, p), pow(a, (q + 1) // 2, p)
    while t != 1:
        i, t2 = 0, t
        while t2 != 1:
            t2 = t2 * t2 % p
            i += 1
        b = pow(c, 1 << (m - i - 1), p)
        m, c, t, r = i, b * b % p, t * b * b % p, r * b % p
    return r


def cipolla(a, p):
    """w = t^2 - a が非剰余になる t を探し、F_p[√w] で (t + √w)^((p+1)/2) を計算する"""
    t = 1
    while True:
        w = (t * t - a) % p
        if pow(w, (p - 1) // 2, p) == p - 1:
            break
        t += 1
    e = (p + 1) // 2
    x0, x1 = 1, 0
    b0, b1 = t, 1
    while e:
        if e & 1:
            x0, x1 = (x0 * b0 + x1 * b1 % p * w) % p, (x0 * b1 + x1 * b0) % p
        b0, b1 = (b0 * b0 + b1 * b1 % p * w) % p, 2 * b0 * b1 % p
        e >>= 1
    return x0


def unit_root(a, p, j):
    """p と互いに素な a について y^2 ≡ a (mod p^j) の (解の個数, 解の 1 つ)。j >= 1"""
    if p == 2:
        if j == 1:
            return 1, 1
        if j == 2:
            return (2, 1) if a % 4 == 1 else (0, None)
        if a % 8 != 1:
            return 0, None
        y = 1                                   # mod 8 の解。1 ビットずつ持ち上げる
        for i in range(3, j):                   # y は mod 2^i の解 → mod 2^(i+1) の解に
            if (y * y - a) % (1 << (i + 1)):
                y += 1 << (i - 1)
        return 4, y % (1 << j)
    if pow(a % p, (p - 1) // 2, p) != 1:
        return 0, None
    y = sqrt_mod_p(a % p, p)
    pj = p ** j
    mod = p
    while mod < pj:                             # Hensel: 精度を 2 倍ずつ（ニュートン法）
        mod = min(mod * mod, pj)
        y = (y - (y * y - a) * pow(2 * y, -1, mod)) % mod
    return 2, y


def solve_prime_power(a, p, k):
    pk = p ** k
    a %= pk
    if a == 0:
        return p ** (k // 2), 0
    v = 0
    while a % p == 0:
        a //= p
        v += 1
    if v % 2:
        return 0, None
    cnt, y = unit_root(a, p, k - v)
    if cnt == 0:
        return 0, None
    h = v // 2
    return cnt * p ** h, (p ** h * y) % pk


def main():
    data = sys.stdin.buffer.read().split()
    m, q = int(data[0]), int(data[1])
    fs = factorize(m)
    mods = [p ** k for p, k in fs]
    # CRT の係数（M/pk の pk を法とする逆元）を前計算
    coef = [(m // pk) * pow(m // pk, -1, pk) % m for pk in mods]
    out = []
    for i in range(q):
        a = int(data[2 + i])
        total, x = 1, 0
        for (p, k), c in zip(fs, coef):
            cnt, r = solve_prime_power(a, p, k)
            if cnt == 0:
                total = 0
                break
            total *= cnt
            x += r * c
        if total == 0:
            out.append("0 -1")
        else:
            out.append(f"{total} {x % m}")
    print("\n".join(out))


main()