Day 179-Q5 — 2026-10-15 | 法 $m\le10^{18}$ の平方根を全部求める（Pollard ρ で素因数分解し、奇素数べきは Tonelli–Shanks＋Hensel、2 べきは 1 ビットずつ、最後に CRT で全組み合わせ）

📋 コピーimport sys
from math import gcd
from random import randrange


def is_prime(n):
    if n < 2:
        return False
    for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
        if n % p == 0:
            return n == p
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2; s += 1
    for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):   # n < 3.3e24 で決定的
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                break
        else:
            return False
    return True


def rho(n):
    """Brent 版 Pollard ρ。gcd を 128 歩ごとにまとめて取る"""
    if n % 2 == 0:
        return 2
    while True:
        c = randrange(1, n)
        y, m, g, r, q = randrange(n), 128, 1, 1, 1
        while g == 1:
            x = y
            for _ in range(r):
                y = (y * y + c) % n
            k = 0
            while k < r and g == 1:
                ys = y
                for _ in range(min(m, r - k)):
                    y = (y * y + c) % n
                    q = q * abs(x - y) % n
                g = gcd(q, n)
                k += m
            r *= 2
        if g == n:                      # まとめすぎて n になったら 1 歩ずつ戻す
            g = 1
            while g == 1:
                ys = (ys * ys + c) % n
                g = gcd(abs(x - ys), n)
        if g != n:
            return g


def factor(n, out):
    for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47):
        while n % p == 0:
            out[p] = out.get(p, 0) + 1
            n //= p
    stack = [n] if n > 1 else []
    while stack:
        x = stack.pop()
        if x == 1:
            continue
        if is_prime(x):
            out[x] = out.get(x, 0) + 1
            continue
        d = rho(x)
        stack.append(d); stack.append(x // d)
    return out


def sqrt_mod_prime(a, p):
    """Tonelli–Shanks。a は p の平方剰余（p 奇素数）"""
    if p % 4 == 3:
        return pow(a, (p + 1) // 4, p)
    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
    m, c, t, r = s, pow(z, q, p), pow(a, q, p), pow(a, (q + 1) // 2, p)
    while t != 1:
        i, tt = 0, t
        while tt != 1:
            tt = tt * tt % 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 roots_prime_power(a, p, e):
    """x^2 ≡ a (mod p^e) の解（gcd(a,p)=1）"""
    pe = p ** e
    a %= pe
    if p == 2:
        if e == 1:
            return [1]
        if e == 2:
            return [1, 3] if a % 4 == 1 else []
        if a % 8 != 1:
            return []
        r = 1                              # 法 8 で 1。k ビット目ずつ直す
        for k in range(3, e):
            if (r * r - a) % (1 << (k + 1)):
                r += 1 << (k - 1)
        h = 1 << (e - 1)
        return sorted({r % pe, (-r) % pe, (r + h) % pe, (-r + h) % pe})
    if pow(a % p, (p - 1) // 2, p) != 1:
        return []
    r = sqrt_mod_prime(a % p, p)
    k = 1
    while k < e:                           # Hensel（ニュートン）で桁を倍々に
        k = min(2 * k, e)
        mod = p ** k
        r = (r - (r * r - a) * pow(2 * r, -1, mod)) % mod
    return [r, pe - r] if r != pe - r else [r]


def main():
    data = sys.stdin.buffer.read().split()
    q = int(data[0])
    out = []
    for i in range(q):
        a, m = int(data[1 + 2 * i]), int(data[2 + 2 * i])
        if m == 1:
            out.append("1 0")
            continue
        fs = factor(m, {})
        sols, mod = [0], 1
        for p, e in sorted(fs.items()):
            rs = roots_prime_power(a, p, e)
            if not rs:
                sols = []
                break
            pe = p ** e
            inv = pow(mod, -1, pe)          # CRT: x = s + mod * ((r - s) * inv mod pe)
            sols = [s + mod * ((r - s) * inv % pe) for s in sols for r in rs]
            mod *= pe
        if sols:
            out.append(f"{len(sols)} {min(sols)}")
        else:
            out.append("0 -1")
    print("\n".join(out))


main()