Problem 518: Prime Triples and Geometric Sequences

Project Euler Problem 518 Solution

EulerSolve provides an optimized solution for Project Euler Problem 518, Prime Triples and Geometric Sequences, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must sum $a+b+c$ over all prime triples $(a,b,c)$ with $a<b<c<n$ such that $a+1$, $b+1$, and $c+1$ form a geometric progression. If we write $$x=a+1,\qquad y=b+1,\qquad z=c+1,$$ then the defining condition is $y^2=xz$. A naive search over all prime triples below $n$ is infeasible for the real input, so the implementations enumerate the hidden arithmetic structure instead of scanning primes directly. Mathematical Approach The key idea is to classify every valid triple by a common scale factor and a coprime pair that describes the geometric progression exactly. Step 1: Convert the geometric progression into one Diophantine equation Three positive numbers are in geometric progression exactly when the square of the middle term equals the product of the outer terms. Therefore $$\frac{b+1}{a+1}=\frac{c+1}{b+1}\iff (b+1)^2=(a+1)(c+1).$$ With $x=a+1$, $y=b+1$, and $z=c+1$, the problem becomes $$y^2=xz,\qquad x<y<z.$$ So we no longer have to reason about geometric progressions directly; we only need a complete description of the positive integer solutions to $y^2=xz$. Step 2: Parameterize all integer solutions Let $$d=\gcd(x,z),\qquad x=d r,\qquad z=d s,\qquad \gcd(r,s)=1.$$ Substituting into $y^2=xz$ gives $$y^2=d^2rs.$$ Hence $rs$ must be a perfect square....

Detailed mathematical approach

Problem Summary

We must sum $a+b+c$ over all prime triples $(a,b,c)$ with $a<b<c<n$ such that $a+1$, $b+1$, and $c+1$ form a geometric progression. If we write

$$x=a+1,\qquad y=b+1,\qquad z=c+1,$$

then the defining condition is $y^2=xz$. A naive search over all prime triples below $n$ is infeasible for the real input, so the implementations enumerate the hidden arithmetic structure instead of scanning primes directly.

Mathematical Approach

The key idea is to classify every valid triple by a common scale factor and a coprime pair that describes the geometric progression exactly.

Step 1: Convert the geometric progression into one Diophantine equation

Three positive numbers are in geometric progression exactly when the square of the middle term equals the product of the outer terms. Therefore

$$\frac{b+1}{a+1}=\frac{c+1}{b+1}\iff (b+1)^2=(a+1)(c+1).$$

With $x=a+1$, $y=b+1$, and $z=c+1$, the problem becomes

$$y^2=xz,\qquad x<y<z.$$

So we no longer have to reason about geometric progressions directly; we only need a complete description of the positive integer solutions to $y^2=xz$.

Step 2: Parameterize all integer solutions

Let

$$d=\gcd(x,z),\qquad x=d r,\qquad z=d s,\qquad \gcd(r,s)=1.$$

Substituting into $y^2=xz$ gives

$$y^2=d^2rs.$$

Hence $rs$ must be a perfect square. Since $r$ and $s$ are coprime, each of them must already be a perfect square, so we may write

$$r=u^2,\qquad s=v^2,\qquad \gcd(u,v)=1.$$

Therefore every solution has the form

$$x=du^2,\qquad y=duv,\qquad z=dv^2,$$

and every positive choice of $d,u,v$ with $\gcd(u,v)=1$ satisfies $y^2=xz$. Because $x<z$, we only need $u<v$.

Step 3: Translate back to the prime triple

Undoing the substitution yields

$$a=du^2-1,\qquad b=duv-1,\qquad c=dv^2-1.$$

Since $u<v$, we automatically get

$$du^2<duv<dv^2,$$

so the ordering $a<b<c$ is built into the parameterization. The coprimality condition $\gcd(u,v)=1$ removes duplicate descriptions, which means each valid prime triple corresponds to one unique search state $(d,u,v)$.

Step 4: Use parity to isolate the only odd-$d$ case

If $a$, $b$, and $c$ are all odd primes, then $a+1$, $b+1$, and $c+1$ are all even, so $x$, $y$, and $z$ are even. In the parameterization above that forces $d$ to be even, because with odd $d$ the coprime pair $u,v$ cannot make all three numbers $du^2$, $duv$, and $dv^2$ even.

Now suppose $d$ is odd. Since $b$ is a prime greater than $2$, $b+1=y$ is even, so $uv$ must be even. Because $\gcd(u,v)=1$, exactly one of $u$ and $v$ is even.

The largest prime $c$ is also greater than $2$, hence $c+1=z$ must be even. Therefore $v$ must be even and $u$ must be odd. But then $x=du^2$ is odd, so $a=x-1$ is even. The only even prime is $2$, so

$$a=2,\qquad x=a+1=3.$$

Thus

$$du^2=3,$$

which forces

$$u=1,\qquad d=3.$$

So the odd-$d$ part of the search is not a general family at all: it is the single special branch

$$a=2,\qquad b=3v-1,\qquad c=3v^2-1,$$

with $v$ even.

Step 5: Bound the finite search region

From the largest value we get

$$c=dv^2-1<n,$$

so, because everything is integral,

$$dv^2\le n,\qquad d\le \left\lfloor\frac{n}{v^2}\right\rfloor.$$

Since $d\ge 1$, this also implies

$$v^2\le n,\qquad v\le \lfloor\sqrt{n}\rfloor.$$

Therefore the search can be organized as: loop over $v$, loop over coprime $u<v$, and then loop over all admissible scale factors $d$. The main branch uses even $d$; the only odd branch is $d=3$, $u=1$, $v$ even.

Worked Example: $n=100$

The special odd branch already gives one valid triple. Taking

$$u=1,\qquad v=2,\qquad d=3$$

produces

$$a=2,\qquad b=5,\qquad c=11.$$

From the even branch, for example,

$$u=1,\qquad v=2,\qquad d=6\Rightarrow (a,b,c)=(5,11,23),$$

and

$$u=2,\qquad v=3,\qquad d=2\Rightarrow (a,b,c)=(7,11,17).$$

Continuing the enumeration below $100$ gives exactly 11 triples:

$(2,5,11)$, $(2,11,47)$, $(5,11,23)$, $(5,17,53)$, $(7,11,17)$, $(7,23,71)$, $(11,23,47)$, $(17,23,31)$, $(17,41,97)$, $(31,47,71)$, $(71,83,97)$.

Their total is

$$18+60+39+75+35+101+81+71+155+149+251=1035,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same plan. First they build an odd-only prime table up to $n$, so every primality test becomes a constant-time table lookup for odd numbers, plus one explicit check for $2$.

Next they loop over $v=2,3,\dots,\lfloor\sqrt{n}\rfloor$. For each $v$, they compute $v^2$ and the largest allowable scale factor $\left\lfloor n/v^2 \right\rfloor$. Then they scan all $u$ with $1\le u<v$, keeping only the coprime pairs.

For each such pair, the implementation first checks the unique odd-$d$ branch $d=3$, $u=1$, $v$ even. After that it runs through all even values of $d$ up to the bound, forms

$$a=du^2-1,\qquad b=duv-1,\qquad c=dv^2-1,$$

and adds $a+b+c$ whenever all three numbers are prime.

The C++ implementation also performs two short internal validations before the full run: the known value $S(100)=1035$, and a fast-versus-direct comparison at $n=500$. The Python and Java implementations keep the same optimized search structure without the extra direct checker.

Complexity Analysis

Building the odd-only sieve up to $n$ costs $O(n\log\log n)$ time and $O(n)$ memory.

The search phase examines

$$\sum_{v\le \sqrt{n}}\sum_{\substack{1\le u<v\\ \gcd(u,v)=1}} O\!\left(\frac{n}{v^2}\right)$$

candidate scale factors, which is $O(n\log n)$ overall. The gcd checks fit within the same order, so the enumeration is dominated by arithmetic generation rather than by primality testing. The full method therefore runs in $O(n\log n)$ time and uses $O(n)$ memory.

Footnotes and References

Problem page: https://projecteuler.net/problem=518
Geometric progression: Wikipedia — Geometric progression
Coprime integers: Wikipedia — Coprime integers
Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes

Problem 518 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = std::int64_t;

struct OddPrimeTable {
    // index i stores primality of number (2*i + 1)
    std::vector<std::uint8_t> odd_is_prime;

    explicit OddPrimeTable(const int limit) : odd_is_prime(static_cast<std::size_t>((limit + 1) / 2), 1U) {
        if (!odd_is_prime.empty()) {
            odd_is_prime[0] = 0U;  // 1 is not prime
        }
        const int max_i = static_cast<int>(odd_is_prime.size()) - 1;
        for (int i = 1; (2 * i + 1) * (2 * i + 1) <= max_i * 2 + 1; ++i) {
            if (!odd_is_prime[static_cast<std::size_t>(i)]) {
                continue;
            }
            const int p = 2 * i + 1;
            int start = (p * p) >> 1;
            for (int j = start; j <= max_i; j += p) {
                odd_is_prime[static_cast<std::size_t>(j)] = 0U;
            }
        }
    }

    bool is_prime(const int x) const {
        if (x == 2) {
            return true;
        }
        if (x < 2 || (x & 1) == 0) {
            return false;
        }
        return odd_is_prime[static_cast<std::size_t>(x >> 1)] != 0U;
    }
};

u64 solve(const int n) {
    OddPrimeTable prime_table(n);
    const int limit = static_cast<int>(std::sqrt(static_cast<long double>(n)));

    u64 sum = 0ULL;
    for (int v = 2; v <= limit; ++v) {
        const int v2 = v * v;
        const int max_d = n / v2;
        for (int u = 1; u < v; ++u) {
            if (std::gcd(u, v) != 1) {
                continue;
            }

            const int u2 = u * u;
            const int uv = u * v;

            // Odd d can only work when a=2, i.e. d*u^2=3 -> u=1, d=3.
            if (u == 1 && (v % 2 == 0) && max_d >= 3) {
                const int d = 3;
                const int a = u2 * d - 1;
                const int b = uv * d - 1;
                const int c = v2 * d - 1;
                if (prime_table.is_prime(a) && prime_table.is_prime(b) && prime_table.is_prime(c)) {
                    sum += static_cast<u64>(a + b + c);
                }
            }

            for (int d = 2; d <= max_d; d += 2) {
                const int a = u2 * d - 1;
                const int b = uv * d - 1;
                const int c = v2 * d - 1;
                if (prime_table.is_prime(a) && prime_table.is_prime(b) && prime_table.is_prime(c)) {
                    sum += static_cast<u64>(a + b + c);
                }
            }
        }
    }
    return sum;
}

u64 brute(const int n) {
    OddPrimeTable prime_table(n);
    std::vector<int> primes;
    for (int x = 2; x < n; ++x) {
        if (prime_table.is_prime(x)) {
            primes.push_back(x);
        }
    }

    u64 sum = 0ULL;
    for (std::size_t i = 0; i < primes.size(); ++i) {
        const int a = primes[i];
        for (std::size_t j = i + 1; j < primes.size(); ++j) {
            const int b = primes[j];
            const i64 lhs = static_cast<i64>(b + 1) * static_cast<i64>(b + 1);
            for (std::size_t k = j + 1; k < primes.size(); ++k) {
                const int c = primes[k];
                const i64 rhs = static_cast<i64>(a + 1) * static_cast<i64>(c + 1);
                if (lhs == rhs) {
                    sum += static_cast<u64>(a + b + c);
                }
            }
        }
    }
    return sum;
}

bool run_checkpoints() {
    if (solve(100) != 1'035ULL) {
        std::cerr << "Checkpoint failed: S(100)\n";
        return false;
    }
    if (solve(500) != brute(500)) {
        std::cerr << "Checkpoint failed: fast/brute mismatch at 500\n";
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr int n = 100'000'000;
    std::cout << solve(n) << '\n';
    return 0;
}

Python

import math

def solve():
    N = 100_000_000

    # Odd sieve for primality
    half = (N + 1) // 2
    sieve = bytearray([1]) * half
    sieve[0] = 0  # 1 is not prime
    for i in range(1, half):
        p = 2 * i + 1
        if p * p > N: break
        if sieve[i]:
            start = (p * p) >> 1
            for j in range(start, half, p):
                sieve[j] = 0

    def is_prime(x):
        if x == 2: return True
        if x < 2 or x % 2 == 0: return False
        return sieve[x >> 1] != 0

    limit = int(N ** 0.5)
    total = 0
    for v in range(2, limit + 1):
        v2 = v * v
        max_d = N // v2
        for u in range(1, v):
            if math.gcd(u, v) != 1: continue
            u2 = u * u
            uv = u * v

            # Odd d=3 special case
            if u == 1 and v % 2 == 0 and max_d >= 3:
                d = 3
                a = u2 * d - 1
                b = uv * d - 1
                c = v2 * d - 1
                if is_prime(a) and is_prime(b) and is_prime(c):
                    total += a + b + c

            # Even d
            for d in range(2, max_d + 1, 2):
                a = u2 * d - 1
                b = uv * d - 1
                c = v2 * d - 1
                if is_prime(a) and is_prime(b) and is_prime(c):
                    total += a + b + c

    return str(total)

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

Java

public class Euler518 {

    static class OddPrimeTable {
        byte[] oddIsPrime;

        OddPrimeTable(int limit) {
            int size = (limit + 1) / 2;
            oddIsPrime = new byte[size];
            for (int i = 0; i < size; i++)
                oddIsPrime[i] = 1;
            if (size > 0)
                oddIsPrime[0] = 0;

            int maxI = size - 1;
            for (int i = 1; (long) (2 * i + 1) * (2 * i + 1) <= maxI * 2 + 1; i++) {
                if (oddIsPrime[i] == 1) {
                    int p = 2 * i + 1;
                    int start = (p * p) >> 1;
                    for (int j = start; j <= maxI; j += p) {
                        oddIsPrime[j] = 0;
                    }
                }
            }
        }

        boolean isPrime(int x) {
            if (x == 2)
                return true;
            if (x < 2 || (x & 1) == 0)
                return false;
            return oddIsPrime[x >> 1] == 1;
        }
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int temp = a % b;
            a = b;
            b = temp;
        }
        return a;
    }

    static long solve(int n) {
        OddPrimeTable primeTable = new OddPrimeTable(n);
        int limit = (int) Math.sqrt(n);

        long sum = 0;
        for (int v = 2; v <= limit; v++) {
            int v2 = v * v;
            int maxD = n / v2;
            boolean vEven = (v % 2 == 0);

            for (int u = 1; u < v; u++) {
                if (gcd(u, v) != 1)
                    continue;

                int u2 = u * u;
                int uv = u * v;

                if (u == 1 && vEven && maxD >= 3) {
                    int d = 3;
                    int a = u2 * d - 1;
                    int b = uv * d - 1;
                    int c = v2 * d - 1;
                    if (primeTable.isPrime(a) && primeTable.isPrime(b) && primeTable.isPrime(c)) {
                        sum += (long) (a + b + c);
                    }
                }

                for (int d = 2; d <= maxD; d += 2) {
                    int a = u2 * d - 1;
                    int b = uv * d - 1;
                    int c = v2 * d - 1;
                    if (primeTable.isPrime(a) && primeTable.isPrime(b) && primeTable.isPrime(c)) {
                        sum += (long) (a + b + c);
                    }
                }
            }
        }
        return sum;
    }

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