Möbius Inversion for Coprime Counting

Intuition, Vocabulary and Foundations

A gcd condition couples several values, which makes direct counting awkward. Möbius inversion replaces that global condition with a signed sum of simpler divisibility conditions. Instead of asking whether two values have gcd $1$, count how many are divisible by each $d$, then combine those counts with the Möbius function.

The central identity is

$$ \mathbf 1_{\gcd(x_1,\ldots,x_k)=1} = \sum_{d\mid\gcd(x_1,\ldots,x_k)}\mu(d). $$

The divisor sum equals $1$ when the gcd is $1$ and $0$ otherwise because $\sum_{d\mid n}\mu(d)=\mathbf1_{n=1}$. This cancellation resembles inclusion–exclusion over prime divisors. Terms with $\mu(d)=1$ add, terms with $\mu(d)=-1$ subtract, and terms with $\mu(d)=0$ ignore divisors containing a squared prime.

For a multiset of positive integers, let $C_d=\#\{i:d\mid a_i\}$. There are $\binom{C_d}{2}$ unordered index pairs whose values are both divisible by $d$. Thus the number of unordered coprime index pairs is $\sum_{d\ge1}\mu(d)\binom{C_d}{2}$. For ordered $k$-tuples from $\{1,\ldots,n\}$ with repetition, the corresponding count is $\sum_{d=1}^{n}\mu(d)\lfloor n/d\rfloor^k$. The combinatorial model determines the inner count.

Step-by-Step Derivations and Algorithms

For an unordered index pair $(i,j)$, insert the indicator identity:

$$ \mathbf1_{\gcd(a_i,a_j)=1} = \sum_{d\mid a_i,\ d\mid a_j}\mu(d). $$

Summing over pairs and exchanging sums gives

$$ \sum_{i\lt j}\mathbf1_{\gcd(a_i,a_j)=1} = \sum_{d\ge1}\mu(d) \sum_{\substack{i\lt j\\d\mid a_i,\ d\mid a_j}}1 = \sum_{d\ge1}\mu(d)\binom{C_d}{2}. $$

Every pair with gcd $g$ receives total weight $\sum_{d\mid g}\mu(d)$, which is one only for $g=1$. To compute all $C_d$ for values bounded by $V$, build an occurrence array and sum over multiples. The number of visits is $\sum_{d=1}^{V}\lfloor V/d\rfloor=O(V\log V)$.

For $\sum_{d\le n}\mu(d)\lfloor n/d\rfloor^k$, precompute the Mertens prefix. If $q=\lfloor n/l\rfloor$, the quotient remains constant through $r=\lfloor n/q\rfloor$, and the block contributes $(M(r)-M(l-1))q^k$. Only $O(\sqrt n)$ distinct quotients occur.

To count gcd exactly $g$, divide eligible values by $g$ and count coprime objects in the scaled domain, or recover all exact-gcd counts by descending divisor inversion. The best approach depends on whether there is one target or an entire gcd distribution.

flowchart TD
    A["Define ordered or unordered objects"] --> B["Replace gcd condition by divisor indicator"]
    B --> C{"Explicit multiset?"}
    C -->|Yes| D["Build value frequencies"]
    D --> E["Accumulate counts over multiples"]
    E --> F["Apply mu[d] times combinatorial count"]
    C -->|No, dense interval| G["Precompute mu and Mertens prefix"]
    G --> H{"Repeated floor quotients?"}
    H -->|Yes| I["Group quotient intervals"]
    H -->|No| J["Sum directly"]
    F --> K["Verify against gcd oracle"]
    I --> K
    J --> K

Implementation, Complexity and Engineering

The reference implementation computes $\mu$ with a linear sieve, counts unordered coprime index pairs in a positive-integer multiset, and counts ordered coprime tuples from a dense interval. Independent brute-force functions verify both counting models.

Python integers prevent count overflow. In fixed-width languages, $\binom{C_d}{2}$ may require 64 bits even when values use 32 bits, and $\lfloor n/d\rfloor^k$ grows rapidly. Use checked wide arithmetic or reduce modulo a specified output modulus.

from itertools import product
from math import gcd


def mobius_linear(limit):
    if limit < 1:
        raise ValueError("limit must be at least 1")
    primes = []
    smallest = [0] * (limit + 1)
    mu = [0] * (limit + 1)
    mu[1] = 1
    for value in range(2, limit + 1):
        if smallest[value] == 0:
            smallest[value] = value
            primes.append(value)
            mu[value] = -1
        for prime in primes:
            composite = value * prime
            if composite > limit:
                break
            smallest[composite] = prime
            if prime == smallest[value]:
                mu[composite] = 0
                break
            mu[composite] = -mu[value]
    return mu


def coprime_unordered_pairs(values):
    if any(value <= 0 for value in values):
        raise ValueError("all values must be positive")
    if len(values) < 2:
        return 0
    maximum = max(values)
    frequency = [0] * (maximum + 1)
    for value in values:
        frequency[value] += 1
    mu = mobius_linear(maximum)
    answer = 0
    for divisor in range(1, maximum + 1):
        divisible = sum(
            frequency[multiple]
            for multiple in range(divisor, maximum + 1, divisor)
        )
        answer += mu[divisor] * divisible * (divisible - 1) // 2
    return answer


def mertens_prefix(limit):
    mu = mobius_linear(limit)
    prefix = [0] * (limit + 1)
    for value in range(1, limit + 1):
        prefix[value] = prefix[value - 1] + mu[value]
    return mu, prefix


def ordered_coprime_tuples(n, arity):
    if n < 1 or arity < 1:
        raise ValueError("n and arity must be positive")
    _, mertens = mertens_prefix(n)
    answer = 0
    left = 1
    while left <= n:
        quotient = n // left
        right = n // quotient
        answer += (
            mertens[right] - mertens[left - 1]
        ) * quotient**arity
        left = right + 1
    return answer


def brute_unordered_pairs(values):
    return sum(
        gcd(values[i], values[j]) == 1
        for i in range(len(values))
        for j in range(i + 1, len(values))
    )


def tuple_gcd(items):
    result = 0
    for item in items:
        result = gcd(result, item)
    return result


def brute_ordered_tuples(n, arity):
    return sum(
        tuple_gcd(items) == 1
        for items in product(range(1, n + 1), repeat=arity)
    )


sample = [2, 3, 4, 6]
assert coprime_unordered_pairs(sample) == 2
assert coprime_unordered_pairs(sample) == brute_unordered_pairs(sample)

with_duplicates = [2, 2, 3, 3, 6, 7]
assert coprime_unordered_pairs(with_duplicates) == brute_unordered_pairs(with_duplicates)

for maximum in range(1, 11):
    values = list(range(1, maximum + 1))
    assert coprime_unordered_pairs(values) == brute_unordered_pairs(values)

for n in range(1, 9):
    for arity in (1, 2, 3):
        assert ordered_coprime_tuples(n, arity) == brute_ordered_tuples(n, arity)

mu = mobius_linear(12)
assert mu[1:13] == [1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0]

try:
    coprime_unordered_pairs([0, 1, 2])
    raise AssertionError("zero input was accepted")
except ValueError:
    pass

print("all tests passed")

Applications, Case Studies and Exercises

Möbius inversion appears whenever common-divisor counts are easy but exact gcd constraints are hard. Typical uses include coprime array pairs, primitive lattice points, visible grid directions, tuples with gcd one, and full gcd histograms.