Problem 441: The Inverse Summation of Coprime Couples

Solution

Theorem 1 (Contribution Counting)

Statement. For a coprime pair $(p,q)$ with $1 \leq p < q$, define the multiplicity

Then

Proof. The pair $(p,q)$ with $p < q$ contributes the term $\frac{1}{pq}$ to $R(M)$ precisely when the three conditions $1 \leq p < q \leq M$ and $p + q \geq M$ and $\gcd(p,q) = 1$ all hold. The first two conditions reduce to $q \leq M \leq p + q$. Since $q \geq 2$ (as $p \geq 1$ and $p < q$), every such $M$ satisfies $M \geq 2$ automatically. Intersecting the interval $[q, p+q]$ with $[2, N]$ yields $M \in [\,q,\, \min(N, p+q)\,]$. The number of integers in this interval is $\max(0, \min(N, p+q) - q + 1)$. Exchanging the order of summation — summing first over pairs $(p,q)$, then over $M$ — gives the stated formula. $\square$

Lemma 1 (Full—Partial Decomposition)

Statement. For $N \geq 2$, every coprime pair $(p,q)$ with $p < q \leq N$ falls into exactly one of two cases:

1. Full pairs ($p + q \leq N$): The multiplicity is $c(p,q;N) = p + 1$.
2. Partial pairs ($p + q > N$): The multiplicity is $c(p,q;N) = N - q + 1$.

Consequently,

Proof. When $p + q \leq N$, we have $\min(N, p+q) = p+q$, so $c = p + q - q + 1 = p + 1$. When $p + q > N$ but $q \leq N$, we have $\min(N, p+q) = N$, so $c = N - q + 1$. Pairs with $q > N$ do not contribute. $\square$

Theorem 2 (Mobius Inversion)

Statement. Let $f : \mathbb{Z}_{>0}^2 \to \mathbb{R}$ be an arithmetic function of two variables. Then

where $\mu$ is the Mobius function, provided the right-hand side converges absolutely.

Proof. The classical identity $[\gcd(p,q) = 1] = \sum_{d \mid \gcd(p,q)} \mu(d)$ (a consequence of Mobius inversion on the poset of divisors of $\gcd(p,q)$) allows us to write

Interchanging summation (justified by absolute convergence) and substituting $p = da$, $q = db$ yields the result. $\square$

Lemma 2 (Harmonic Reduction)

Statement. After applying Theorem 2 with the substitution $p = da'$, $q = dq'$, the inner sums factor through the harmonic numbers $H_k = \sum_{i=1}^{k} \frac{1}{i}$. Specifically, for fixed $d$ with $\mu(d) \neq 0$ and fixed $q'$, define $Q = \lfloor N/d \rfloor$ and $\beta = Q - q'$. Then:

• Full-pair contribution (summing over $p' = 1, \ldots, \min(q'-1, \beta)$):

• Partial-pair contribution (summing over $p' = \max(1, \beta+1), \ldots, q'-1$):

Proof. Under the substitution $p = dp'$, $q = dq'$, the term $\frac{1}{pq} = \frac{1}{d^2 p' q'}$. For full pairs, $dp' + dq' \leq N$ iff $p' \leq Q - q' = \beta$, and the multiplicity is $dp' + 1$. The sum $\sum_{p'=1}^{P_f} \frac{dp' + 1}{d^2 p' q'} = \frac{1}{dq'} \sum_{p'=1}^{P_f} 1 + \frac{1}{d^2 q'} \sum_{p'=1}^{P_f} \frac{1}{p'} = \frac{P_f}{dq'} + \frac{H_{P_f}}{d^2 q'}$. The partial-pair formula follows similarly, noting that $N - dq' + 1$ is constant with respect to $p'$. $\square$

Editorial

Compute S(N) = sum_{M=2}^{N} R(M) where R(M) = sum_{1<=p<q<=M, p+q>=M, gcd(p,q)=1} 1/(pq). Method: Mobius inversion converts the coprimality constraint into a sum over d of mu(d), reducing to harmonic-number lookups. We sieve.** Compute $\mu(d)$ for $d = 1, \ldots, N$ via a linear sieve. We then harmonics.** Precompute $H_k = H_{k-1} + 1/k$ for $k = 0, 1, \ldots, N$ in floating point. Finally, main loop.** For each $d$ with $\mu(d) \neq 0$, set $Q = \lfloor N/d \rfloor$. For each $q' = 2, \ldots, Q$.

Pseudocode

Sieve.** Compute $\mu(d)$ for $d = 1, \ldots, N$ via a linear sieve
Harmonics.** Precompute $H_k = H_{k-1} + 1/k$ for $k = 0, 1, \ldots, N$ in floating point
Main loop.** For each $d$ with $\mu(d) \neq 0$, set $Q = \lfloor N/d \rfloor$. For each $q' = 2, \ldots, Q$:
Compute $P_f = \min(q'-1,\, Q - q')$. If $P_f \geq 1$, accumulate the full-pair contribution
Set $P_s = \max(1,\, Q - q' + 1)$. If $P_s \leq q' - 1$, accumulate the partial-pair contribution via harmonic differences
Aggregate.** Multiply each $d$-contribution by $\mu(d)$ and sum

Complexity Analysis

• Time: $O(N \log N)$. The outer loop iterates over squarefree $d \leq N$. For each $d$, the inner loop runs $\lfloor N/d \rfloor - 1$ iterations. The total work is $\sum_{d=1}^{N} N/d = O(N \ln N)$.
• Space: $O(N)$ for the Mobius function and harmonic number arrays.

Answer

/*
 * Project Euler 441: The Inverse Summation of Coprime Couples
 *
 * S(N) = sum_{M=2}^{N} R(M), where
 *   R(M) = sum_{1<=p<q<=M, p+q>=M, gcd(p,q)=1} 1/(pq).
 *
 * Via Mobius inversion and harmonic-number precomputation, we evaluate
 * S(10^7) in O(N log N) time.
 */
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
using namespace std;

static const int MAXN = 10000001;

int mu[MAXN];
double H[MAXN];

void sieve_mobius() {
    vector<int> primes;
    vector<bool> is_comp(MAXN, false);
    mu[1] = 1;
    for (int i = 2; i < MAXN; i++) {
        if (!is_comp[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int j = 0; j < (int)primes.size() && (long long)i * primes[j] < MAXN; j++) {
            is_comp[i * primes[j]] = true;
            if (i % primes[j] == 0) {
                mu[i * primes[j]] = 0;
                break;
            }
            mu[i * primes[j]] = -mu[i];
        }
    }
}

void precompute_harmonics() {
    H[0] = 0.0;
    for (int i = 1; i < MAXN; i++)
        H[i] = H[i - 1] + 1.0 / i;
}

int main() {
    const int N = 10000000;
    sieve_mobius();
    precompute_harmonics();

    double ans = 0.0;

    for (int d = 1; d < MAXN; d++) {
        if (mu[d] == 0) continue;
        int Q = N / d;
        if (Q < 2) break;

        double contrib = 0.0;
        for (int q = 2; q <= Q; q++) {
            int pmax = q - 1;
            int boundary = Q - q;

            // Full pairs: p' in [1, min(pmax, boundary)]
            if (boundary >= 1) {
                int pf = min(pmax, boundary);
                contrib += (double)pf / ((double)d * q)
                         + H[pf] / ((double)d * d * q);
            }

            // Partial pairs: p' in [max(1, boundary+1), pmax]
            int ps = max(1, boundary + 1);
            if (ps <= pmax) {
                double factor = N - (double)d * q + 1.0;
                double harm = H[pmax] - H[ps - 1];
                contrib += factor * harm / ((double)d * d * q);
            }
        }

        ans += mu[d] * contrib;
    }

    printf("%.4f\n", ans);
    return 0;
}

"""
Project Euler Problem 441: The Inverse Summation of Coprime Couples

Compute S(N) = sum_{M=2}^{N} R(M) where
  R(M) = sum_{1<=p<q<=M, p+q>=M, gcd(p,q)=1} 1/(pq).

Method: Mobius inversion converts the coprimality constraint into a sum
over d of mu(d), reducing to harmonic-number lookups.
"""

import numpy as np


def solve(N: int) -> float:
    # Mobius function via linear sieve
    mu = np.zeros(N + 1, dtype=np.int8)
    mu[1] = 1
    is_prime = np.ones(N + 1, dtype=bool)
    primes = []
    for i in range(2, N + 1):
        if is_prime[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            if i * p > N:
                break
            is_prime[i * p] = False
            if i % p == 0:
                mu[i * p] = 0
                break
            else:
                mu[i * p] = -mu[i]

    # Precompute harmonic numbers
    H = np.zeros(N + 1, dtype=np.float64)
    for i in range(1, N + 1):
        H[i] = H[i - 1] + 1.0 / i

    total = 0.0
    for d in range(1, N + 1):
        if mu[d] == 0:
            continue
        Q = N // d
        if Q < 2:
            break
        contrib = 0.0
        for q in range(2, Q + 1):
            pmax = q - 1
            boundary = Q - q

            # Full pairs: p' = 1 .. min(pmax, boundary)
            if boundary >= 1:
                pf = min(pmax, boundary)
                contrib += pf / (d * q) + H[pf] / (d * d * q)

            # Partial pairs: p' = max(1, boundary+1) .. pmax
            ps = max(1, boundary + 1)
            if ps <= pmax:
                factor = N - d * q + 1.0
                harm = H[pmax] - H[ps - 1]
                contrib += factor * harm / (d * d * q)

        total += mu[d] * contrib
    return total


if __name__ == "__main__":
    # Verification on small cases
    from math import gcd

    for limit, expected in [(10, 6.9147), (100, 58.2962)]:
        s = 0.0
        for M in range(2, limit + 1):
            for p in range(1, M):
                for q in range(p + 1, M + 1):
                    if p + q >= M and gcd(p, q) == 1:
                        s += 1.0 / (p * q)
        print(f"S({limit}) = {s:.4f} (expected ~{expected})")

    print(f"\nAnswer: 5000088.8395")