Möbius Function + Inversion, in depth: the divisor-sum identity, linear sieve, coprime and squarefree counting, and Mertens via Du's sieve

Definition and first values

For a positive integer n, μ(n) is defined by its prime factorisation. μ(1) = 1. If n has a squared prime factor, μ(n) = 0. Otherwise n is a product of k distinct primes and μ(n) = (-1)^k. So primes get -1, products of two distinct primes get +1, and anything divisible by 4, 9, 25 and so on gets 0.

μ is multiplicative: μ(ab) = μ(a)μ(b) whenever gcd(a, b) = 1. That property is what lets a sieve compute it for every n up to a limit in linear time, and it is why μ behaves well inside sums over divisors.

The identity everything rests on

The identity that drives everything is that the sum of μ(d) over the divisors d of n equals 1 when n = 1 and 0 otherwise. Write it as Σ_{d|n} μ(d) = [n = 1], where the bracket is 1 if the condition holds and 0 if not.

Proof: for n = 1 the only divisor is 1 and μ(1) = 1. For n greater than 1, let n have k ≥ 1 distinct prime factors. Divisors with a squared prime contribute 0, so only squarefree divisors count, and those are products of a subset of the k primes. A subset of size j contributes (-1)^j and there are C(k, j) such subsets, so the sum is Σ_j C(k, j)(-1)^j = (1 - 1)^k = 0. Check with n = 12: μ(1) + μ(2) + μ(3) + μ(4) + μ(6) + μ(12) = 1 - 1 - 1 + 0 + 1 + 0 = 0.

In practice the identity is a filter. Whenever you need [gcd(a, b) = 1] inside a sum, replace it with Σ_{d | gcd(a,b)} μ(d), which equals Σ over d dividing both a and b of μ(d). Then swap the order of summation so d is on the outside. The inner sums become counts of multiples of d, which are just floor divisions. That one move solves most competitive programming problems tagged Möbius.

Möbius inversion

Möbius inversion states: if g(n) = Σ_{d|n} f(d) for all n, then f(n) = Σ_{d|n} μ(d) g(n/d), and the converse also holds. In words, if g accumulates f over divisors, μ undoes the accumulation.

Proof: substitute the definition of g. Σ_{d|n} μ(d) g(n/d) = Σ_{d|n} μ(d) Σ_{e | n/d} f(e). The pairs (d, e) range over all d, e with d·e dividing n, so regroup by e: Σ_{e|n} f(e) Σ_{d | n/e} μ(d). By the key identity the inner sum is 1 only when n/e = 1, so the whole expression is f(n).

The figure uses a classic pair. Every integer n equals the sum of φ(d) over its divisors, so with g(n) = n and f = φ, inversion gives φ(n) = Σ_{d|n} μ(d)·n/d. For n = 12 that is 12 - 6 - 4 + 2 = 4, and the four numbers up to 12 coprime to it are 1, 5, 7 and 11. The totient article covers φ itself; here φ is only the test case.

Inversion is inclusion-exclusion on the divisor lattice. Ordinary inclusion-exclusion over sets is the same theorem on the lattice of subsets, where the Möbius function of a pair of sets S ⊆ T is (-1)^{|T| - |S|}. Recognising that both are one idea helps when a problem mixes them.

Computing μ: the linear sieve

For a single n, compute μ by trial division in O(√n): divide out each prime factor, return 0 if any prime divides twice, and flip the sign for each prime found. For every value up to a limit, use the linear sieve, which visits each composite exactly once through its smallest prime factor.

def mobius_sieve(n: int) -> list[int]:
    """mu[0..n] in O(n). mu[0] is unused and left as 0."""
    mu = [0] * (n + 1)
    if n >= 1:
        mu[1] = 1
    is_comp = bytearray(n + 1)
    primes = []
    for i in range(2, n + 1):
        if not is_comp[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            ip = i * p
            if ip > n:
                break
            is_comp[ip] = 1
            if i % p == 0:          # p already divides i, so p*p divides ip
                mu[ip] = 0
                break               # ip will not be reached again: linear time
            mu[ip] = -mu[i]         # one more distinct prime flips the sign
    return mu

assert mobius_sieve(12)[1:] == [1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0]

Memory is the real limit. In Python a list of 10^7 values costs about 80 MB in pointers alone, and the loop is slow; use an array of signed bytes or numpy for large limits. In C++ an int8_t array of 10^8 entries is 100 MB, the primes list adds about 23 MB (5.76 million primes as 32-bit ints), and the sieve runs in roughly a second.

Worked problem: coprime pairs

Count ordered pairs (a, b) with 1 ≤ a, b ≤ N and gcd(a, b) = 1. Apply the filter: the count is Σ_a Σ_b Σ_{d | gcd(a,b)} μ(d). Move d outside. For a fixed d, a and b must both be multiples of d, and there are ⌊N/d⌋ choices for each, so the answer is Σ_{d=1}^{N} μ(d)·⌊N/d⌋².

Worked example, N = 10. Only squarefree d contribute: d = 1 gives +100, d = 2 gives -25, d = 3 gives -9, d = 5 gives -4, d = 6 gives +1, d = 7 gives -1, d = 10 gives +1, and d = 4, 8 and 9 give 0. The total is 63. Cross-check another way: the pairs with a ≤ b number Σ φ(b) for b from 1 to 10, which is 32, and doubling minus the single diagonal pair (1, 1) gives 2·32 - 1 = 63.

The same derivation gives k-tuples with gcd 1 as Σ μ(d)·⌊N/d⌋^k, pairs from two different ranges as Σ μ(d)·⌊N/d⌋·⌊M/d⌋, and pairs with gcd exactly g by running the formula on ⌊N/g⌋. Raise to the power with fast modular exponentiation when the answer is wanted modulo a prime.

Floor-division blocks

Summing over every d up to N costs O(N). The quotient ⌊N/d⌋ takes only about 2√N distinct values, and it is constant on blocks of consecutive d. With a prefix sum of μ, each block costs O(1).

from itertools import accumulate

def coprime_pairs(N: int, M_prefix: list[int]) -> int:
    """M_prefix[x] = mu(1) + ... + mu(x), valid for x <= N."""
    total, d = 0, 1
    while d <= N:
        q = N // d
        last = N // q                       # largest d' with N // d' == q
        total += (M_prefix[last] - M_prefix[d - 1]) * q * q
        d = last + 1
    return total

mu = mobius_sieve(10)
assert coprime_pairs(10, list(accumulate(mu))) == 63

With many queries, sieve once to the largest N and answer each query in O(√N). For two ranges N and M, the block ends at min(N // (N // d), M // (M // d)).

Counting squarefree numbers

A number is squarefree when no square greater than 1 divides it. The count of squarefree numbers up to N is Q(N) = Σ_{d=1}^{⌊√N⌋} μ(d)·⌊N/d²⌋. The reason: the indicator of squarefree n is Σ over d with d² dividing n of μ(d), because the largest square dividing n is s², and the sum over d dividing s of μ(d) is [s = 1]. Swap the sums and count multiples of d².

Check N = 10: d = 1 gives 10, d = 2 gives -⌊10/4⌋ = -2, d = 3 gives -⌊10/9⌋ = -1, so Q(10) = 7. The squarefree numbers are 1, 2, 3, 5, 6, 7 and 10. Because the sieve only needs to reach √N, N around 10^14 needs a sieve to 10^7, which is cheap.

The Mertens function at large N

Some problems need the Mertens function M(n) = Σ_{k≤n} μ(k) at n far beyond what you can sieve, say 10^11. Summing the key identity over all n up to N gives Σ_{k=1}^{N} M(⌊N/k⌋) = 1, so M(N) = 1 - Σ_{k=2}^{N} M(⌊N/k⌋). The right side only needs M at the O(√N) distinct quotients, and grouping k into blocks makes each evaluation O(√n). Sieving the small values up to about N^{2/3} and memoising the large ones gives the O(N^{2/3}) method usually called Du's sieve.

def make_mertens(limit: int):
    small = list(accumulate(mobius_sieve(limit)))    # small[x] = M(x) for x <= limit
    memo: dict[int, int] = {}

    def M(n: int) -> int:
        if n <= limit:
            return small[n]
        if n in memo:
            return memo[n]
        res, k = 1, 2
        while k <= n:
            q = n // k
            last = n // q
            res -= (last - k + 1) * M(q)
            k = last + 1
        memo[n] = res
        return res
    return M

M = make_mertens(1000)
assert M(10) == -1
assert make_mertens(100)(10**5) == make_mertens(10**5)(10**5)   # recursion agrees with a plain sieve

Choose limit near N^{2/3} times a small constant; for N = 10^10 that is a few million. The recursion depth stays small because each call divides n by at least 2, but in Python convert it to an explicit loop over the quotients in increasing order if you push N past 10^11.

A sieve-free alternative

You do not always need μ. To count pairs with gcd exactly g for every g at once, start from cnt(g) = ⌊N/g⌋², the pairs whose gcd is a multiple of g, and subtract exact counts of the larger multiples, working from g = N down to 1.

def pairs_with_exact_gcd(N: int) -> list[int]:
    f = [0] * (N + 1)
    for g in range(N, 0, -1):
        f[g] = (N // g) ** 2
        for m in range(2 * g, N + 1, g):
            f[g] -= f[m]
    return f

assert pairs_with_exact_gcd(10)[1] == 63

This runs in O(N log N) by the harmonic sum, works when the per-multiple counts come from data rather than a formula (frequencies of array values, for example), and has no sign bookkeeping. The μ formula wins when you need one value fast, when N is too large to iterate, or when the sum must be split into floor blocks.

Beyond gcd: periodic objects

Inversion also counts objects that are periodic. The number of monic irreducible polynomials of degree n over a field with q elements is (1/n)·Σ_{d|n} μ(d)·q^{n/d}. For q = 2 and n = 4 that is (16 - 4 + 0)/4 = 3, matching x⁴+x+1, x⁴+x³+1 and x⁴+x³+x²+x+1. The same pattern counts aperiodic necklaces (Lyndon words), which shows up in de Bruijn sequence construction and in combinatorics problems about rotations.

Failure modes

• Negative values under a modulus. μ is -1 for many d; in C++ or Java, add the modulus before taking the remainder.
• Overflow. ⌊N/d⌋² overflows a signed 64-bit integer once ⌊N/d⌋ passes about 3·10^9, which happens at small d for N = 10^10; use 128-bit arithmetic or reduce modulo early.
• Sieving to the wrong bound. Squarefree counting needs μ to √N, coprime pairs need it to N, Du's sieve needs it to about N^{2/3}.
• Forgetting μ(1) = 1. A sieve that only sets values at primes and composites leaves the most important term at 0.
• Off-by-one in blocks. The block end is N // (N // d), inclusive; use prefix[last] - prefix[d - 1].
• Counting ordered versus unordered pairs. The formula counts ordered pairs; convert explicitly.

What to do next

1. Implement mobius_sieve and check it against trial division for every n up to 10^5.
2. Solve the coprime pairs problem for N up to 10^7 with the block sum, then for two different ranges.
3. Count squarefree numbers up to 10^14 and compare with the 6/π² density estimate.
4. Implement the Mertens recursion and test M(10^9) against a sieve run on a machine with enough memory.
5. Rewrite one solution with the exact-gcd-by-multiples method and compare speed and simplicity.
6. Read the advanced number theory page next for primitive roots and discrete logs.