Day 172-Q5 — 2026-10-08 | 二平方和（Pollard の rho で分解し、4k+1 型の素数を Cornacchia で a²+b² に、ガウス整数の積で組み立てる。表し方の数は 4∏(e+1)）

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


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):   # この 12 個の底で n < 3×10^23 まで決定的（記憶）
        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):
    # Pollard の rho（Brent 版、gcd を 128 回分まとめて取る）
    if n % 2 == 0:
        return 2
    while True:
        y, c, m = random.randrange(1, n), random.randrange(1, n), 128
        g = r = q = 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:
            g = 1
            while g == 1:
                ys = (ys * ys + c) % n
                g = gcd(abs(x - ys), n)
        if g != n:
            return g


def factor(n, res):
    if n == 1:
        return
    if is_prime(n):
        res[n] = res.get(n, 0) + 1
        return
    r = isqrt(n)
    if r * r == n:                             # 平方数は先に開く（p^2 の rho は実測で遅かった）
        factor(r, res)
        factor(r, res)
        return
    d = rho(n)
    factor(d, res)
    factor(n // d, res)


def two_squares_prime(p):
    # p ≡ 1 (mod 4) を a^2 + b^2 に。x^2 ≡ -1 (mod p) を作り、Cornacchia（ユークリッドの互除法を √p 未満で止める）
    c = 2
    while pow(c, (p - 1) // 2, p) != p - 1:   # 平方非剰余 c を探す（半分の c が当たる）
        c += 1
    x = pow(c, (p - 1) // 4, p)               # x^2 = c^((p-1)/2) = -1
    a, b = p, x
    r = isqrt(p)
    while b > r:
        a, b = b, a % b
    return b, isqrt(p - b * b)


def gmul(x, y):
    # ガウス整数の積 (a+bi)(c+di)
    return (x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0])


def gpow(x, e):
    r = (1, 0)
    while e:
        if e & 1:
            r = gmul(r, x)
        x = gmul(x, x)
        e >>= 1
    return r


def solve(n):
    f = {}
    small = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47)
    for p in small:                            # 小さい素因数は割り算で先に取る（rho の負担を減らす）
        while n % p == 0:
            f[p] = f.get(p, 0) + 1
            n //= p
    factor(n, f)
    count = 4
    z = (1, 0)
    for p, e in f.items():
        if p % 4 == 3:
            if e % 2:
                return 0, None                 # 4k+3 型の素数が奇数乗 → 表せない
            z = (z[0] * p ** (e // 2), z[1] * p ** (e // 2))
        elif p == 2:
            z = gmul(z, gpow((1, 1), e))       # 2 = (1+i)(1-i)
        else:
            count *= e + 1                     # r2(n) = 4 × Π(e+1)（p ≡ 1 mod 4）
            z = gmul(z, gpow(two_squares_prime(p), e))
    a, b = sorted((abs(z[0]), abs(z[1])))
    return count, (a, b)


def main():
    data = sys.stdin.buffer.read().split()
    t = int(data[0])
    out = []
    for i in range(1, t + 1):
        n = int(data[i])
        cnt, ab = solve(n)
        out.append(f"{cnt} {ab[0]} {ab[1]}" if cnt else "0")
    print("\n".join(out))


if __name__ == "__main__":
    main()