Project Euler Problem 87

Solution

Answer: 1097343

We want the number of integers below $50{,}000{,}000$ that can be written as

$$n = p_2^2 + p_3^3 + p_4^4$$

where $p_2,p_3,p_4$ are primes.

The key challenge is counting distinct values, because multiple prime triples may generate the same sum.

1. Mathematical analysis

Let the limit be

$$L = 50{,}000{,}000.$$

We need all numbers

$$p^2 + q^3 + r^4 < L$$

with $p,q,r$ prime.

Step 1: Bound the search

Since all terms are positive, each prime only needs to be generated up to a certain maximum.

For fourth powers:

$$r^4 < L \Rightarrow r < L^{1/4}.$$

Since

$$50{,}000{,}000^{1/4} \approx 84,$$

we only need primes up to 83.

For cubes:

$$q^3 < L \Rightarrow q < L^{1/3}.$$

Since

$$50{,}000{,}000^{1/3} \approx 368,$$

we only need primes up to 367.

For squares:

$$p^2 < L \Rightarrow p < \sqrt{L}.$$

Since

$$\sqrt{50{,}000{,}000}\approx 7071,$$

we only need primes up to 7071.

So the search space is actually tiny:

prime squares: about 900 primes,
prime cubes: about 70 primes,
prime fourth powers: about 20 primes.

Total combinations are roughly

$$900 \times 70 \times 20 \approx 1.26\times10^6,$$

which is trivial computationally.

Step 2: Avoid duplicates

Different triples may produce the same sum.

Example:

$$p_1^2 + q_1^3 + r_1^4 = p_2^2 + q_2^3 + r_2^4.$$

We only count the number once.

The cleanest approach is to store every valid sum in a set, since sets automatically remove duplicates.

Step 3: Efficient generation

Precompute:

all prime squares below $L$,
all prime cubes below $L$,
all prime fourth powers below $L$.

Then loop:

$$s = a+b+c$$

and insert into a set whenever

$$s < L.$$

We can prune early because the lists are sorted:

if $a+b \ge L$, stop inner loops,
if $a+b+c \ge L$, stop deepest loop.

2. Python implementation

from math import isqrt

LIMIT = 50_000_000

def sieve(n):
    """Return list of primes <= n using Sieve of Eratosthenes."""
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False

    for p in range(2, isqrt(n) + 1):
        if is_prime[p]:
            for multiple in range(p * p, n + 1, p):
                is_prime[multiple] = False

    return [i for i, prime in enumerate(is_prime) if prime]

# Generate primes only up to the needed maximum
max_prime = isqrt(LIMIT)
primes = sieve(max_prime)

# Precompute prime powers under LIMIT
prime_squares = [p**2 for p in primes if p**2 < LIMIT]
prime_cubes = [p**3 for p in primes if p**3 < LIMIT]
prime_fourths = [p**4 for p in primes if p**4 < LIMIT]

# Store distinct representable numbers
values = set()

for square in prime_squares:
    for cube in prime_cubes:
        partial = square + cube

        if partial >= LIMIT:
            break

        for fourth in prime_fourths:
            total = partial + fourth

            if total >= LIMIT:
                break

            values.add(total)

print(len(values))

3. Code walkthrough

Sieve of Eratosthenes

def sieve(n):

Builds all primes up to n.

We only need primes up to

$$\sqrt{50{,}000{,}000}$$

because square terms dominate the largest prime requirement.

Prime powers

prime_squares = [p**2 for p in primes if p**2 < LIMIT]
prime_cubes = [p**3 for p in primes if p**3 < LIMIT]
prime_fourths = [p**4 for p in primes if p**4 < LIMIT]

We precompute powers once instead of repeatedly recomputing them in loops.

Triple nested loop

for square in prime_squares:
    for cube in prime_cubes:

Compute

$$square + cube$$

first.

If already too large:

if partial >= LIMIT:
    break

we stop because cubes are increasing.

Then:

for fourth in prime_fourths:

Try adding fourth powers.

Again:

if total >= LIMIT:
    break

because fourth powers are sorted.

Deduplication

values.add(total)

A set guarantees distinct sums.

Finally:

len(values)

returns the number of valid integers.

4. Small-example verification

For numbers below 50, the program finds exactly:

$$28,33,47,49$$

which matches the problem statement (4 numbers).

That strongly validates the logic.

The overall count should be in the millions range at most, since the total number of combinations is only around $1.2$ million and duplicates exist, so a six-digit answer is plausible.

Running the computation gives:

Answer: 1097343