Problem 454: Diophantine Reciprocals III

Project Euler Problem 454 Solution

EulerSolve provides an optimized solution for Project Euler Problem 454, Diophantine Reciprocals III, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must count the number $F(L)$ of positive integer triples $(x,y,n)$ satisfying $$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n.$$ The symmetry between $x$ and $y$ means that the strict inequality $x \lt y$ removes the mirror-image duplicate of every solution. The published checkpoints are $F(15)=4$ and $F(1000)=1069$, while the real target is $F(10^{12})$. Mathematical Approach Step 1: Extract the gcd and obtain a coprime parametrization Let $$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$$ with $u \lt v$ and $\gcd(u,v)=1$. Substituting into the reciprocal equation gives $$n(u+v)=g\,u\,v.$$ Because $\gcd(u,u+v)=1$ and $\gcd(v,u+v)=1$, we can read divisibility information directly from this identity. Since $u\mid n(u+v)$ and $u$ is coprime to $u+v$, it follows that $u\mid n$. The same argument gives $v\mid n$. As $\gcd(u,v)=1$, we conclude that $uv\mid n$, so we may write $$n=tuv,\qquad t\in\mathbb Z_{>0}.$$ Then the previous equation simplifies to $g=t(u+v)$. Therefore every admissible solution has the form $$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v,$$ with $u,v,t \in \mathbb Z_{>0}$, $u \lt v$, and $\gcd(u,v)=1$. Conversely, every triple of this form satisfies $\frac1x+\frac1y=\frac1n$, so this is a bijection and not just a one-way construction....

Detailed mathematical approach

Problem Summary

We must count the number $F(L)$ of positive integer triples $(x,y,n)$ satisfying

$$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n.$$

The symmetry between $x$ and $y$ means that the strict inequality $x \lt y$ removes the mirror-image duplicate of every solution. The published checkpoints are $F(15)=4$ and $F(1000)=1069$, while the real target is $F(10^{12})$.

Mathematical Approach

Step 1: Extract the gcd and obtain a coprime parametrization

Let

$$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$$

with $u \lt v$ and $\gcd(u,v)=1$. Substituting into the reciprocal equation gives

$$n(u+v)=g\,u\,v.$$

Because $\gcd(u,u+v)=1$ and $\gcd(v,u+v)=1$, we can read divisibility information directly from this identity. Since $u\mid n(u+v)$ and $u$ is coprime to $u+v$, it follows that $u\mid n$. The same argument gives $v\mid n$. As $\gcd(u,v)=1$, we conclude that $uv\mid n$, so we may write

$$n=tuv,\qquad t\in\mathbb Z_{>0}.$$

Then the previous equation simplifies to $g=t(u+v)$. Therefore every admissible solution has the form

$$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v,$$

with $u,v,t \in \mathbb Z_{>0}$, $u \lt v$, and $\gcd(u,v)=1$. Conversely, every triple of this form satisfies $\frac1x+\frac1y=\frac1n$, so this is a bijection and not just a one-way construction.

Step 2: Rewrite the search space in the form used by the algorithm

Introduce the new variables

$$a=v,\qquad s=u+v.$$

Then $u=s-a$. The condition $1\le u \lt v$ becomes

$$a+1\le s\le 2a-1.$$

The coprimality condition is preserved because

$$\gcd(a,s)=\gcd(v,u+v)=\gcd(u,v)=1.$$

The upper bound on $y$ is now

$$y=t\,a\,s\le L.$$

So for a fixed coprime pair $(a,s)$, the number of admissible scale factors $t$ is exactly

$$1\le t\le \left\lfloor\frac{L}{as}\right\rfloor,$$

which contributes $\left\lfloor L/(as)\right\rfloor$ solutions. Hence

$$F(L)=\sum_{a=1}^{A_{\max}}\ \sum_{s=a+1}^{U(a)} \mathbf{1}_{\gcd(a,s)=1}\left\lfloor\frac{L}{as}\right\rfloor,$$

where

$$U(a)=\min\left(2a-1,\left\lfloor\frac{L}{a}\right\rfloor\right).$$

The outer loop can stop once no valid $s$ exists. Since the smallest allowed $s$ is $a+1$, we need $a(a+1)\le L$. Therefore

$$A_{\max}=\left\lfloor\frac{\sqrt{1+4L}-1}{2}\right\rfloor.$$

Step 3: Count coprime values in an interval with Möbius inversion

For a fixed outer parameter $a$, define

$$C_a(\alpha,\beta)=\#\left\{s\in\mathbb Z:\alpha\le s\le \beta,\ \gcd(a,s)=1\right\}.$$

The standard Möbius identity is

$$\mathbf{1}_{\gcd(a,s)=1}=\sum_{d\mid \gcd(a,s)}\mu(d).$$

Summing this over an interval gives

$$C_a(\alpha,\beta)=\sum_{d\mid a}\mu(d)\left(\left\lfloor\frac{\beta}{d}\right\rfloor-\left\lfloor\frac{\alpha-1}{d}\right\rfloor\right).$$

Only squarefree divisors contribute, because $\mu(d)=0$ whenever $d$ contains a repeated prime factor. The implementation therefore precomputes, for each outer value, exactly those divisors with nonzero Möbius sign and reuses them for every interval query.

Step 4: Group equal quotients into harmonic blocks

Fix $a$ and set

$$M=\left\lfloor\frac{L}{a}\right\rfloor,\qquad q(s)=\left\lfloor\frac{M}{s}\right\rfloor.$$

As $s$ increases, $q(s)$ remains constant across contiguous ranges. If the current block starts at $\ell$ and

$$q=\left\lfloor\frac{M}{\ell}\right\rfloor,$$

then the same quotient remains valid up to

$$r=\min\left(U(a),\left\lfloor\frac{M}{q}\right\rfloor\right).$$

Every $s\in[\ell,r]$ contributes the same multiplier $q$, so the whole block contributes

$$q\cdot C_a(\ell,r).$$

This is the main speedup: instead of scanning each inner value one by one, the algorithm jumps directly between maximal ranges on which the floor quotient is constant.

Worked Example: $L=15$

Here

$$A_{\max}=\left\lfloor\frac{\sqrt{61}-1}{2}\right\rfloor=3.$$

For $a=1$, the interval $a+1\le s\le 2a-1$ is empty, so there is no contribution.

For $a=2$, the only possible value is $s=3$. It is coprime to $2$, and it contributes

$$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2.$$

For $a=3$, the admissible values are $s=4$ and $s=5$. Both are coprime to $3$, giving

$$\left\lfloor\frac{15}{3\cdot 4}\right\rfloor=1,\qquad \left\lfloor\frac{15}{3\cdot 5}\right\rfloor=1.$$

Therefore

$$F(15)=2+1+1=4,$$

which matches the checkpoint exactly.

How the Code Works

The C++, Python, and Java implementations all follow this derivation. They first compute $A_{\max}$ using an integer square root so the outer range is exact. Then they build the Möbius function up to that limit with a sieve, and from it they assemble, for every outer value, the signed list of squarefree divisors needed by the interval formula above.

During the main summation, the implementation determines the valid inner interval, advances through maximal quotient blocks, and for each block evaluates the coprime count using the precomputed signed divisors. The C++ and Java versions store this data in compact flattened form to reduce allocation overhead, while the Python version keeps direct per-value lists. The Java implementation also parallelizes independent outer ranges across available CPU cores.

Complexity Analysis

Let $A=A_{\max}=O(\sqrt L)$. Computing the Möbius function is sieve-like and near-linear in $A$. Building the stored squarefree-divisor data costs

$$\sum_{\substack{d\le A\\ \mu(d)\ne 0}}\left\lfloor\frac{A}{d}\right\rfloor=O(A\log A),$$

and the memory needed for that data has the same order.

The summation stage is much smaller than a naive search over all $x$ and $y$. It processes quotient blocks rather than individual inner values, and each block performs one Möbius interval count over the stored divisors of the current outer value. In practice this behaves close to $O(\sqrt L\log L)$, which is why the method remains feasible for $L=10^{12}$, whereas direct enumeration would be hopelessly large.

Footnotes and References

Problem page: https://projecteuler.net/problem=454
Möbius function: Wikipedia — Möbius function
Möbius inversion formula: Wikipedia — Möbius inversion formula
Unit-fraction background: Wikipedia — Egyptian fraction

Problem 454 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using i64 = std::int64_t;
using u128 = __uint128_t;

struct Options {
    u64 l = 1'000'000'000'000ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        out = static_cast<u64>(std::stoull(tail));
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--l=", options.l)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.l >= 1ULL;
}

u64 integer_sqrt(const u64 n) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((r + 1ULL) <= n / (r + 1ULL)) {
        ++r;
    }
    while (r > n / r) {
        --r;
    }
    return r;
}

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }
    std::string out;
    while (value > 0) {
        const unsigned digit = static_cast<unsigned>(value % 10);
        out.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(out.begin(), out.end());
    return out;
}

u64 solve(const u64 L) {
    const u64 root = integer_sqrt(1ULL + 4ULL * L);
    const u32 max_b = static_cast<u32>((root - 1ULL) / 2ULL);

    std::vector<int> mu(static_cast<std::size_t>(max_b) + 1U, 0);
    std::vector<u32> primes;
    primes.reserve(static_cast<std::size_t>(max_b / 10U));
    std::vector<bool> is_composite(static_cast<std::size_t>(max_b) + 1U, false);
    mu[1] = 1;

    for (u32 i = 2U; i <= max_b; ++i) {
        if (!is_composite[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (const u32 p : primes) {
            const u64 v = static_cast<u64>(i) * static_cast<u64>(p);
            if (v > max_b) {
                break;
            }
            is_composite[static_cast<std::size_t>(v)] = true;
            if (i % p == 0U) {
                mu[static_cast<std::size_t>(v)] = 0;
                break;
            }
            mu[static_cast<std::size_t>(v)] = -mu[i];
        }
    }

    std::vector<u32> divisor_count(static_cast<std::size_t>(max_b) + 1U, 0U);
    for (u32 d = 1U; d <= max_b; ++d) {
        if (mu[d] == 0) {
            continue;
        }
        for (u32 b = d; b <= max_b; b += d) {
            ++divisor_count[b];
        }
    }

    std::vector<u32> offset(static_cast<std::size_t>(max_b) + 2U, 0U);
    for (u32 b = 1U; b <= max_b; ++b) {
        offset[static_cast<std::size_t>(b) + 1U] =
            offset[static_cast<std::size_t>(b)] + divisor_count[b];
    }

    std::vector<int> signed_divisor(offset[static_cast<std::size_t>(max_b) + 1U], 0);
    std::vector<u32> fill = offset;

    for (u32 d = 1U; d <= max_b; ++d) {
        const int mud = mu[d];
        if (mud == 0) {
            continue;
        }
        const int sd = mud * static_cast<int>(d);
        for (u32 b = d; b <= max_b; b += d) {
            signed_divisor[fill[b]++] = sd;
        }
    }

    u128 answer = 0;
    for (u32 b = 1U; b <= max_b; ++b) {
        const u64 M = L / b;
        const u64 upper = std::min<u64>(2ULL * b - 1ULL, M);
        if (upper <= b) {
            continue;
        }

        u64 c = static_cast<u64>(b) + 1ULL;
        while (c <= upper) {
            const u64 q = M / c;
            const u64 c2 = std::min<u64>(upper, M / q);

            const u64 left_minus_1 = c - 1ULL;
            i64 coprime_count = 0;
            for (u32 idx = offset[b]; idx < offset[static_cast<std::size_t>(b) + 1U]; ++idx) {
                const int sd = signed_divisor[idx];
                const u32 d = static_cast<u32>(sd > 0 ? sd : -sd);
                const i64 term =
                    static_cast<i64>(c2 / d) - static_cast<i64>(left_minus_1 / d);
                coprime_count += (sd > 0) ? term : -term;
            }

            answer += static_cast<u128>(q) * static_cast<u128>(coprime_count);
            c = c2 + 1ULL;
        }
    }

    return static_cast<u64>(answer);
}

bool run_checkpoints() {
    if (solve(15ULL) != 4ULL) {
        std::cerr << "Checkpoint failed: L(15)\n";
        return false;
    }
    if (solve(1000ULL) != 1069ULL) {
        std::cerr << "Checkpoint failed: L(1000)\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const u64 answer = solve(options.l);
    std::cout << answer << '\n';
    return 0;
}

Python

import math

def solve():
    L = 1000000000000
    root = math.isqrt(1 + 4 * L)
    max_b = (root - 1) // 2

    # Mobius sieve
    mu = [0] * (max_b + 1)
    mu[1] = 1
    is_comp = bytearray(max_b + 1)
    primes = []
    for i in range(2, max_b + 1):
        if not is_comp[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            v = i * p
            if v > max_b: break
            is_comp[v] = 1
            if i % p == 0:
                mu[v] = 0
                break
            mu[v] = -mu[i]

    # For each b, store signed divisors d with mu[d] != 0
    div_lists = [[] for _ in range(max_b + 1)]
    for d in range(1, max_b + 1):
        if mu[d] == 0: continue
        sd = mu[d] * d
        for b in range(d, max_b + 1, d):
            div_lists[b].append(sd)

    answer = 0
    for b in range(1, max_b + 1):
        M = L // b
        upper = min(2 * b - 1, M)
        if upper <= b: continue
        c = b + 1
        while c <= upper:
            q = M // c
            c2 = min(upper, M // q)
            left_m1 = c - 1
            coprime_count = 0
            for sd in div_lists[b]:
                d = abs(sd)
                term = c2 // d - left_m1 // d
                coprime_count += term if sd > 0 else -term
            answer += q * coprime_count
            c = c2 + 1

    return str(answer)

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;

public class Euler454 {

    static long integerSqrt(long n) {
        long r = (long) Math.sqrt(n);
        while ((r + 1) <= n / (r + 1)) {
            r++;
        }
        while (r > n / r) {
            r--;
        }
        return r;
    }

    public static String solve() {
        long L = 1000000000000L;
        long root = integerSqrt(1L + 4L * L);
        final int maxB = (int) ((root - 1L) / 2L);

        int[] mu = new int[maxB + 1];
        List<Integer> primes = new ArrayList<>(maxB / 10);
        boolean[] isComposite = new boolean[maxB + 1];
        mu[1] = 1;

        for (int i = 2; i <= maxB; i++) {
            if (!isComposite[i]) {
                primes.add(i);
                mu[i] = -1;
            }
            for (int p : primes) {
                long v = (long) i * p;
                if (v > maxB)
                    break;
                isComposite[(int) v] = true;
                if (i % p == 0) {
                    mu[(int) v] = 0;
                    break;
                }
                mu[(int) v] = -mu[i];
            }
        }

        int[] divisorCount = new int[maxB + 1];
        for (int d = 1; d <= maxB; d++) {
            if (mu[d] == 0)
                continue;
            for (int b = d; b <= maxB; b += d) {
                divisorCount[b]++;
            }
        }

        int[] offset = new int[maxB + 2];
        for (int b = 1; b <= maxB; b++) {
            offset[b + 1] = offset[b] + divisorCount[b];
        }

        int[] signedDivisor = new int[offset[maxB + 1]];
        int[] fill = new int[maxB + 2];
        System.arraycopy(offset, 0, fill, 0, maxB + 2);

        for (int d = 1; d <= maxB; d++) {
            int mud = mu[d];
            if (mud == 0)
                continue;
            int sd = mud * d;
            for (int b = d; b <= maxB; b += d) {
                signedDivisor[fill[b]++] = sd;
            }
        }

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        int initialChunkSize = maxB / threads;
        final int chunkSize = initialChunkSize == 0 ? 1 : initialChunkSize;

        int totalChunks = (maxB + chunkSize - 1) / chunkSize;

        long answer = IntStream.range(0, totalChunks).parallel().mapToLong(chunkIdx -> {
            int startB = chunkIdx * chunkSize + 1;
            int endB = Math.min(maxB, startB + chunkSize - 1);
            long localAns = 0;

            for (int b = startB; b <= endB; b++) {
                long M = L / b;
                long upper = Math.min(2L * b - 1L, M);
                if (upper <= b)
                    continue;

                long c = (long) b + 1L;
                int offStart = offset[b];
                int offEnd = offset[b + 1];

                while (c <= upper) {
                    long q = M / c;
                    long c2 = Math.min(upper, M / q);

                    long leftMinus1 = c - 1L;
                    long coprimeCount = 0;

                    for (int idx = offStart; idx < offEnd; idx++) {
                        int sd = signedDivisor[idx];
                        int d = sd > 0 ? sd : -sd;
                        long term = (c2 / d) - (leftMinus1 / d);
                        if (sd > 0)
                            coprimeCount += term;
                        else
                            coprimeCount -= term;
                    }

                    localAns += q * coprimeCount;
                    c = c2 + 1L;
                }
            }
            return localAns;
        }).sum();

        return Long.toString(answer);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}