Dirichlet Convolution + Multiplicative Functions, in depth: the convolution ring, inverses, prime-by-prime transforms and the hyperbola method

The operation and its algebra

An arithmetic function is any map from the positive integers to numbers. Define (f * g)(n) = sum over d | n of f(d) g(n/d). Each divisor pair (d, n/d) appears once, so the operation is commutative. Associativity follows because (f * g) * h at n sums f(a) g(b) h(c) over all ordered triples with abc = n, which does not depend on grouping.

The identity is ε, with ε(1) = 1 and ε(n) = 0 otherwise: (f * ε)(n) keeps only the term d = n. Pointwise addition distributes over *, so arithmetic functions form a commutative ring. The constant function 1, the identity id(n) = n and the power functions id_k(n) = n^k are the building blocks.

Inverses and where Möbius comes from

A function f has a Dirichlet inverse g, with f * g = ε, exactly when f(1) ≠ 0. Solve for g one value at a time: at n = 1 the equation is f(1) g(1) = 1, and at n > 1 it is f(1) g(n) + sum over d | n, d < n of g(d) f(n/d) = 0. Every term in the sum uses a smaller argument, so the recursion is well founded:

g(1) = 1/f(1), and g(n) = -(1/f(1)) · sum over d | n, d < n of g(d) f(n/d).

So the functions with f(1) ≠ 0 form an abelian group under *. Over the integers the inverse stays integral only when f(1) = ±1; otherwise work in rationals or modulo a prime that does not divide f(1). The Möbius function is simply the inverse of the constant 1, which is why Möbius inversion works: if F = f * 1 then F * μ = f * (1 * μ) = f * ε = f.

Multiplicative functions and closure

A function is multiplicative if f(1) = 1 and f(mn) = f(m) f(n) whenever gcd(m, n) = 1. It is then fixed by its values on prime powers. Two closure facts make the whole theory practical:

• If f and g are multiplicative, so is f * g. A divisor d of mn with gcd(m, n) = 1 splits uniquely as d = ab with a | m and b | n, and the sum factors into a product.
• If f is multiplicative, so is its inverse. The inverse is built from f by the recursion above, and an induction on n shows it inherits the factorisation.

The consequence is that on a prime power the convolution is a short sum: (f * g)(p^k) = sum over i from 0 to k of f(p^i) g(p^(k-i)). That is how τ(p^k) = k + 1 and σ(p^k) = 1 + p + ... + p^k drop out, and it is what a linear sieve exploits when it builds f * g in O(N) by tracking the exponent of the smallest prime.

One trap: the pointwise sum of multiplicative functions is not multiplicative in general, and neither is the pointwise product with a non-multiplicative function. Closure holds for the convolution, the pointwise product of two multiplicative functions, and the inverse.

Behind all this sits the Dirichlet series F(s) = sum f(n)/n^s. Multiplying two series collects n^-s terms with ab = n, so the product of series is the series of the convolution. ζ(s) is the series of 1, ζ(s)^2 is the series of τ, and 1/ζ(s) is the series of μ. A multiplicative function has an Euler product over primes, which is the analytic face of the prime-power rule above.

Computing convolutions up to N

Three routines cover most needs. The first convolves any two functions up to N by looping over multiples: for each i it touches N/i entries, so the total is about N ln N. The second inverts any f with f(1) = ±1 using the recursion, pushing each new g(d) forward to its multiples. The third is the prime-by-prime transform: F = f * 1 computed by treating each prime as one axis of a prefix sum over exponent vectors, at a cost of N times the sum of 1/p, about N ln ln N.

def convolve(f, g, N):
    """h = f * g on 1..N, O(N log N). Lists are indexed from 0; index 0 is unused."""
    h = [0] * (N + 1)
    for i in range(1, N + 1):
        if f[i]:
            fi = f[i]
            for j in range(1, N // i + 1):
                h[i * j] += fi * g[j]
    return h

def dirichlet_inverse(f, N):
    """g with f * g = eps. Requires f[1] in (1, -1) for an integer result."""
    g = [0] * (N + 1)
    g[1] = f[1]                          # 1/f(1) when f(1) is +-1
    acc = [0] * (N + 1)                  # acc[n] = sum over d|n, d<n of g(d) f(n/d)
    for d in range(1, N + 1):
        if d > 1:
            g[d] = -acc[d] * g[1]
        if g[d]:
            for k in range(2, N // d + 1):
                acc[d * k] += g[d] * f[k]
    return g

def zeta_transform(f, N, primes):
    """F = f * 1 in O(N log log N). Ascending i is required."""
    F = f[:]
    for p in primes:
        for i in range(1, N // p + 1):
            F[i * p] += F[i]
    return F

def mobius_transform(F, N, primes):
    """f = F * mu, the exact inverse of zeta_transform. Descending i is required."""
    f = F[:]
    for p in primes:
        for i in range(N // p, 0, -1):
            f[i * p] -= f[i]
    return f

The direction of the inner loop matters. Ascending i in the zeta transform lets F[i] already include contributions along the p-axis, so F[ip] accumulates every power of p. Running it descending would add only one step and produce f(n) + f(n/p). The inverse is the mirror image: descending i subtracts the unmodified neighbour, which is a difference along the axis.

Checks that passed: convolve(μ, 1) = ε; dirichlet_inverse(1) equals μ from a sieve; convolve(μ, id) matches a gcd-count totient for n < 300; zeta_transform(id) equals σ and mobius_transform undoes it; and 200 random integer functions with f(1) = ±1 each satisfied f * inverse(f) = ε up to N = 2,000.

Worked example at n = 12, and measured cost

Take n = 12 with divisors 1, 2, 3, 4, 6, 12. For τ = 1 * 1 each divisor contributes 1, giving 6. For σ = id * 1 the terms are the divisors themselves, summing to 28. For φ = μ * id the nonzero μ values are μ(1) = 1, μ(2) = -1, μ(3) = -1 and μ(6) = 1, giving 12 - 6 - 4 + 2 = 4, which matches the four residues 1, 5, 7, 11 coprime to 12.

Check multiplicativity the cheap way: 12 = 4 · 3 and τ(4) τ(3) = 3 · 2 = 6, σ(4) σ(3) = 7 · 4 = 28, φ(4) φ(3) = 2 · 2 = 4. A test that compares h(ab) with h(a) h(b) over coprime pairs confirmed τ and φ as computed by the generic convolution are multiplicative, which is a useful smoke test for any new function you build.

Measured cost at N = 1,000,000: computing σ with the harmonic double loop performed 13,970,034 inner updates and took 4.91 s; the prime-by-prime transform performed 2,853,708 updates and took 1.58 s, with identical output. When one side of the convolution is the constant 1, always use the transform.

The hyperbola method for summatory functions

Often you need the summatory function S(x) = sum over n ≤ x of (f * g)(n) for x far beyond any sieve, say 10^12. Write it as a sum over pairs: S(x) = sum over ab ≤ x of f(a) g(b). Let F and G be the prefix sums of f and g, and split at u = v = ⌊√x⌋:

S(x) = sum over a ≤ u of f(a) G(x/a) + sum over b ≤ v of g(b) F(x/b) - F(u) G(v).

Each pair with a ≤ u is counted in the first sum, each pair with b ≤ v in the second, and pairs with both are subtracted once. No pair can have a > √x and b > √x, so nothing is missed. If F and G are cheap to evaluate, the cost is O(√x).

For f = g = 1 this is the classic divisor summatory function:

from math import isqrt

def divisor_summatory(x):
    """D(x) = sum of tau(n) for n <= x, in O(sqrt x)."""
    r = isqrt(x)
    return 2 * sum(x // i for i in range(1, r + 1)) - r * r

It agrees with the naive sum of ⌊x/i⌋ for every x below 3,000. D(100) = 482, D(10^12) = 27,785,452,449,086, and D(10^14) = 3,239,062,263,181,054 took 5.52 s in pure Python for its 10^7 floor divisions. When F or G is itself expensive, such as the Mertens function, pair the hyperbola split with memoised recursion; that is the idea behind Du's sieve on the Möbius page.

Operational guidance

• Pick the cheapest tool for the shape. Generic f * g up to N: harmonic loop. f * 1 or its inverse: prime-by-prime transform. Both multiplicative: linear sieve with prime-power rules. One summatory value at huge x: hyperbola method.
• Integer width. σ up to 10^7 fits in 64 bits, but summatory functions such as D(10^14) are near 3.2 · 10^15 and sums of σ grow like x^2. In C++ or Java, use 128-bit intermediates or reduce modulo a prime throughout.
• Memory. Python lists of ints cost about 8 bytes per pointer plus objects; at N = 10^7 prefer numpy or array('q'). The prime-by-prime transform is in place and vectorises per prime with numpy slicing.
• Testing. Keep a brute-force divisor-sum oracle and compare for n ≤ 2,000 on every change, plus the identity f * inverse(f) = ε on random inputs.

Failure modes

• Confusing convolutions. The Dirichlet product sums over divisor pairs; the ordinary (Cauchy) product used for polynomials sums over a + b = n. Swapping them produces plausible numbers that are wrong.
• Inverting with f(1) not ±1 over integers. Integer division silently truncates. Use fractions or a modulus.
• Wrong loop direction in the prime-by-prime transforms, as described above; the output still looks like a divisor sum.
• Assuming closure that does not exist. Sums of multiplicative functions are not multiplicative, so a sieve that relies on f(p^k m) = f(p^k) f(m) gives wrong results for them.
• Off-by-one at the square. Forgetting to subtract F(u) G(v), or using a floating-point square root that rounds the wrong way at perfect squares. Use isqrt.

Trade-offs

What to do next

1. Implement convolve and a brute-force divisor-sum oracle; verify the six identities in the table for n ≤ 2,000.
2. Implement dirichlet_inverse and check that the inverse of 1 is μ and the inverse of id is μ(n)·n.
3. Replace every f * 1 in your code with the prime-by-prime transform and measure the difference at your N.
4. Compute D(10^12) with the hyperbola method and compare with the value above.
5. Read the Möbius page for Du's sieve, the totient page for φ in practice, and prime counting for sub-linear sums over primes; the number theory overview has the gcd and modular tools these rely on.