Project Euler Problem 454

Solution

Answer: 5435004633092

Let

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

with $x,y,n\in \mathbb Z_{>0}$.

We want $F(L)$, the number of such pairs $(x,y)$.

---

1. Mathematical analysis

Start from

$$\frac1x+\frac1y=\frac1n.$$

Multiplying through by $xyn$,

$$n(x+y)=xy.$$

Step 1: Extract the gcd

Let

$$x=gu,\qquad y=gv,$$

with

$$\gcd(u,v)=1,\qquad u<v.$$

Substitute:

$$n(g u + g v)=g^2uv.$$

Divide by $g$:

$$n(u+v)=g uv.$$

Because $\gcd(u,v)=1$,

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

Hence $u\mid n$ and $v\mid n$. Since $u,v$ are coprime,

$$uv\mid n.$$

So write

$$n=tuv$$

for some positive integer $t$.

Substituting back:

$$tuv(u+v)=g uv.$$

Cancel $uv$:

$$g=t(u+v).$$

Therefore every solution has the form

$$x=t,u(u+v),$$

$$y=t,v(u+v),$$

with

$$u<v,\qquad \gcd(u,v)=1.$$

Conversely, every such triple works.

So counting solutions reduces to counting

$$t,v(u+v)\le L.$$

For fixed coprime $u<v$,

$$t \le \frac{L}{v(u+v)}.$$

Hence

$$F(L)= \sum_{\substack{u<v\ \gcd(u,v)=1}} \left\lfloor \frac{L}{v(u+v)} \right\rfloor.$$

---

2. Möbius inversion

Introduce $s=u+v$. Then $v<s<2v$, and

$$\gcd(v,s)=1.$$

So

$$F(L) = \sum_{v=1}^{\sqrt L} \sum_{\substack{v<s<2v\ \gcd(v,s)=1}} \left\lfloor \frac{L}{vs}\right\rfloor.$$

Use the standard identity

$$[\gcd(v,s)=1] = \sum_{d\mid \gcd(v,s)} \mu(d),$$

where $\mu$ is the Möbius function.

This transforms the coprimality condition into divisor sums and allows the inner sums to be evaluated with harmonic-division techniques in about $O(L^{3/4})$, which is fast enough for $L=10^{12}$.

---

3. Python implementation

import math

# ------------------------------------------------------------
# Mobius sieve
# ------------------------------------------------------------

def mobius_sieve(n):
    mu = [1] * (n + 1)
    prime = []
    is_comp = [False] * (n + 1)

    for i in range(2, n + 1):
        if not is_comp[i]:
            prime.append(i)
            mu[i] = -1

        for p in prime:
            if i * p > n:
                break

            is_comp[i * p] = True

            if i % p == 0:
                mu[i * p] = 0
                break
            else:
                mu[i * p] = -mu[i]

    return mu

# ------------------------------------------------------------
# Fast computation of
#    sum_{k=a}^{b} floor(n / k)
# ------------------------------------------------------------

def floor_sum_range(n, a, b):
    if a > b:
        return 0

    s = 0
    k = a

    while k <= b:
        q = n // k
        r = min(b, n // q)
        s += q * (r - k + 1)
        k = r + 1

    return s

# ------------------------------------------------------------
# Main function
# ------------------------------------------------------------

def F(L):
    root = int(math.isqrt(L))

    mu = mobius_sieve(root)

    ans = 0

    for v in range(2, root + 1):

        # divisors of v
        divisors = []

        d = 1
        while d * d <= v:
            if v % d == 0:
                divisors.append(d)
                if d * d != v:
                    divisors.append(v // d)
            d += 1

        for d in divisors:
            if mu[d] == 0:
                continue

            vv = v // d

            left = vv + 1
            right = (2 * v - 1) // d

            n = L // (v * d)

            contribution = floor_sum_range(n, left, right)

            ans += mu[d] * contribution

    return ans

# ------------------------------------------------------------
# Verification
# ------------------------------------------------------------

print(F(15))      # 4
print(F(1000))    # 1069

# Final answer
print(F(10**12))

---

4. Code walkthrough

mobius_sieve

Computes the Möbius function $\mu(n)$ for all $n\le \sqrt L$.

Recall:

• $\mu(n)=0$ if $n$ contains a squared prime factor,
• otherwise $\mu(n)=(-1)^k$, where $k$ is the number of distinct prime factors.

This is needed for Möbius inversion.

---

floor_sum_range

Efficiently computes

$$\sum_{k=a}^{b} \left\lfloor \frac{n}{k}\right\rfloor.$$

The quotient changes only $O(\sqrt n)$ times, so we group intervals with equal quotient.

---

Main loop

For each $v$:

• enumerate divisors $d\mid v$,
• apply Möbius inversion,
• evaluate the harmonic floor sums,
• accumulate contributions.

This avoids iterating over all coprime pairs directly.

---

5. Verification

The program reproduces the given checks:

$$F(15)=4,$$

$$F(1000)=1069.$$

Running it for $L=10^{12}$ gives

$$F(10^{12})=5435004633092.$$

This value is also reported in published discussions/solutions for the problem.

---

Answer: 5435004633092