Problem 808: Reversible Prime Squares

Project Euler Problem 808 Solution

EulerSolve provides an optimized solution for Project Euler Problem 808, Reversible Prime Squares, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We seek numbers $s$ with three simultaneous properties: $s$ must be the square of a prime, $s$ must not be palindromic in base 10, and the decimal reversal of $s$ must again be the square of a prime. The task is to list these reversible prime squares in increasing order and sum the first $50$ of them. If $s$ qualifies and $t=\operatorname{rev}(s)$, then $t$ is also a reversible prime square unless $s=t$, which is exactly the palindromic case that must be excluded. The search therefore focuses on a very special subset of prime squares rather than on all integers. Mathematical Approach Let $$\mathcal{R}=\left\{p^2 : p\in\mathbb{P},\ p^2\ne \operatorname{rev}(p^2),\ \operatorname{rev}(p^2)=q^2,\ q\in\mathbb{P}\right\}.$$ The implementations search this set directly. The mathematics is simple but precise: restrict the search to prime roots, reverse the decimal digits, and certify that the reversed value is also a prime square. Step 1: Restrict the Search Space to Prime Squares Every admissible value has the form $$s=p^2,\qquad p\in\mathbb{P}.$$ So there is no reason to test arbitrary integers. It is enough to enumerate primes $p$ and examine their squares. Because the map $p\mapsto p^2$ is strictly increasing for positive integers, scanning prime roots in increasing order automatically scans candidate squares in increasing order as well....

Detailed mathematical approach

Problem Summary

We seek numbers $s$ with three simultaneous properties: $s$ must be the square of a prime, $s$ must not be palindromic in base 10, and the decimal reversal of $s$ must again be the square of a prime. The task is to list these reversible prime squares in increasing order and sum the first $50$ of them.

If $s$ qualifies and $t=\operatorname{rev}(s)$, then $t$ is also a reversible prime square unless $s=t$, which is exactly the palindromic case that must be excluded. The search therefore focuses on a very special subset of prime squares rather than on all integers.

Mathematical Approach

Let

$$\mathcal{R}=\left\{p^2 : p\in\mathbb{P},\ p^2\ne \operatorname{rev}(p^2),\ \operatorname{rev}(p^2)=q^2,\ q\in\mathbb{P}\right\}.$$

The implementations search this set directly. The mathematics is simple but precise: restrict the search to prime roots, reverse the decimal digits, and certify that the reversed value is also a prime square.

Step 1: Restrict the Search Space to Prime Squares

Every admissible value has the form

$$s=p^2,\qquad p\in\mathbb{P}.$$

So there is no reason to test arbitrary integers. It is enough to enumerate primes $p$ and examine their squares. Because the map $p\mapsto p^2$ is strictly increasing for positive integers, scanning prime roots in increasing order automatically scans candidate squares in increasing order as well.

Step 2: Reverse the Decimal Expansion and Exclude Fixed Points

For a candidate square $s$, define its decimal reversal by

$$r=\operatorname{rev}(s).$$

The problem requires a genuinely different reversed value, so palindromes are rejected immediately:

$$s\ne r.$$

This is more than a cosmetic filter. If $s=r$, then the reversal gives back the same number instead of a distinct reversible partner. Also, for any prime $p>5$, the square $p^2$ ends in $1$ or $9$, so its reversal begins in $1$ or $9$; there is no ambiguity from leading zeros.

Step 3: Characterize When the Reversal Is Also a Prime Square

A positive integer is a prime square exactly when it is a perfect square and its positive square root is prime. Therefore $s$ belongs to $\mathcal{R}$ if and only if

$$s=p^2,\qquad s\ne \operatorname{rev}(s),\qquad \operatorname{rev}(s)=q^2,\qquad p,q\in\mathbb{P}.$$

So the test for a reversed value $r$ is exact:

$$a=\left\lfloor\sqrt{r}\right\rfloor,\qquad a^2=r,\qquad a\in\mathbb{P}.$$

If either the square test fails or the primality test fails, the original square is discarded.

Step 4: Reversible Prime Squares Come in Distinct Pairs

Suppose $s=p^2\in\mathcal{R}$ and let $t=\operatorname{rev}(s)=q^2$. Reversing again gives

$$\operatorname{rev}(t)=s.$$

Because palindromes are excluded, we have $s\ne t$. Hence reversible prime squares typically appear as two-way pairs $(s,t)$. The algorithm does not need any special pairing logic: once it tests every prime square in increasing order, both members of a valid pair will be encountered naturally.

Step 5: Worked Example

The smallest example comes from

$$13^2=169,\qquad \operatorname{rev}(169)=961=31^2.$$

Since both $13$ and $31$ are prime, $169$ is valid. The reverse direction also works:

$$31^2=961,\qquad \operatorname{rev}(961)=169=13^2,$$

so $961$ is valid as well.

Two quick counterexamples illustrate the other filters. First,

$$11^2=121,\qquad \operatorname{rev}(121)=121,$$

so $121$ is rejected for being palindromic. Second,

$$17^2=289,\qquad \operatorname{rev}(289)=982,$$

and $982$ is not a perfect square, so $289$ is rejected.

Step 6: Why an Expanding Prime Bound Is Enough

The implementations generate primes up to a current bound $B$, test every square $p^2$ with $p\le B$, and keep all valid values found in that pass. If fewer than $50$ values have been collected, the bound is enlarged and the search is repeated. This works because the candidate squares are examined in strictly increasing order. Once a pass finds at least $50$ reversible prime squares, the first $50$ found are exactly the first $50$ elements of $\mathcal{R}$.

How the Code Works

The C++, Python, and Java implementations all follow the same logical pipeline. They begin by constructing a prime sieve up to a current bound and then iterate through the primes in ascending order. Each prime is squared, and palindromic squares are discarded before any more expensive work is done.

For every remaining square, the implementation reverses its decimal digits, computes an exact integer square root of the reversed value, and checks whether the reversal is a perfect square. If it is, the square root is then tested for primality. The primality check uses deterministic Miller-Rabin after an initial screen by a short list of small primes, so the result is exact for the integer range under consideration.

Whenever the reversed value is confirmed to be a prime square, the original square is appended to the answer list. If the list reaches $50$ entries, the implementations sum those values and return the total. Otherwise they increase the sieve bound and restart with a larger search range. The arithmetic details differ slightly by language, but the mathematical test is identical in all three implementations.

Complexity Analysis

Let $B$ be the current sieve bound on the prime root. Building the sieve costs $O(B\log\log B)$ time and $O(B)$ memory. The number of candidates tested is $\pi(B)$, one for each prime up to $B$. For each candidate, digit reversal and integer square root use $O(\log B)$ digit work, while deterministic Miller-Rabin uses a fixed witness set, so in this setting the primality certification is effectively constant-time per candidate and more formally requires $O(\log B)$ modular multiplications. In practice the sieve dominates the growth, with small extra work per prime root.

Footnotes and References

Problem page: https://projecteuler.net/problem=808
Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes
Miller-Rabin primality test: Wikipedia — Miller-Rabin primality test
Integer square root: Wikipedia — Integer square root
Palindromic number: Wikipedia — Palindromic number

Problem 808 source code

C++

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

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

static inline u64 mul_mod(u64 a, u64 b, u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) % static_cast<u128>(mod));
}

static u64 pow_mod(u64 a, u64 e, u64 mod) {
    u64 r = 1 % mod;
    a %= mod;
    while (e > 0) {
        if (e & 1ULL) {
            r = mul_mod(r, a, mod);
        }
        a = mul_mod(a, a, mod);
        e >>= 1ULL;
    }
    return r;
}

static bool is_prime(u64 n) {
    if (n < 2) {
        return false;
    }
    for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (n == p) {
            return true;
        }
        if (n % p == 0) {
            return false;
        }
    }

    u64 d = n - 1;
    int s = 0;
    while ((d & 1ULL) == 0ULL) {
        d >>= 1ULL;
        ++s;
    }

    static constexpr u64 WITNESSES[] = {2ULL, 325ULL, 9'375ULL, 28'178ULL, 450'775ULL, 9'780'504ULL, 1'795'265'022ULL};
    for (u64 a : WITNESSES) {
        if (a % n == 0) {
            continue;
        }
        u64 x = pow_mod(a, d, n);
        if (x == 1 || x == n - 1) {
            continue;
        }
        bool comp = true;
        for (int r = 1; r < s; ++r) {
            x = mul_mod(x, x, n);
            if (x == n - 1) {
                comp = false;
                break;
            }
        }
        if (comp) {
            return false;
        }
    }
    return true;
}

static std::vector<u32> sieve_primes(const u32 limit) {
    std::vector<bool> is_prime_vec(static_cast<std::size_t>(limit + 1), true);
    if (limit >= 0) {
        is_prime_vec[0] = false;
    }
    if (limit >= 1) {
        is_prime_vec[1] = false;
    }
    for (u32 p = 2; static_cast<u64>(p) * p <= limit; ++p) {
        if (!is_prime_vec[p]) {
            continue;
        }
        for (u32 q = p * p; q <= limit; q += p) {
            is_prime_vec[q] = false;
        }
    }
    std::vector<u32> primes;
    primes.reserve(static_cast<std::size_t>(limit / 10));
    for (u32 i = 2; i <= limit; ++i) {
        if (is_prime_vec[i]) {
            primes.push_back(i);
        }
    }
    return primes;
}

static u64 reverse_digits(u64 x) {
    u64 r = 0;
    while (x > 0) {
        r = r * 10 + (x % 10);
        x /= 10;
    }
    return r;
}

static bool is_palindrome(u64 x) {
    return x == reverse_digits(x);
}

static u64 isqrt_u64(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;
}

static std::vector<u64> first_reversible_prime_squares(int need) {
    for (u32 limit = 1U << 16U;; limit <<= 1U) {
        const auto primes = sieve_primes(limit);
        std::vector<u64> vals;
        vals.reserve(static_cast<std::size_t>(need + 16));

        for (u32 p : primes) {
            const u64 sq = static_cast<u64>(p) * static_cast<u64>(p);
            if (is_palindrome(sq)) {
                continue;
            }
            const u64 rev = reverse_digits(sq);
            const u64 r = isqrt_u64(rev);
            if (r * r == rev && is_prime(r)) {
                vals.push_back(sq);
                if (static_cast<int>(vals.size()) >= need) {
                    vals.resize(static_cast<std::size_t>(need));
                    return vals;
                }
            }
        }
        if (limit >= (1U << 30U)) {
            break;
        }
    }

    return {};
}

int main() {
    const auto vals = first_reversible_prime_squares(50);

    assert(vals[0] == 169ULL);
    assert(vals[1] == 961ULL);

    u64 sum = 0;
    for (u64 v : vals) {
        sum += v;
    }

    std::cout << sum << '\n';
    return 0;
}

Python

import math

def solve():
    def mod_pow(base, exp, mod):
        r = 1 % mod; a = base % mod
        while exp > 0:
            if exp & 1: r = r * a % mod
            a = a * a % mod
            exp >>= 1
        return r

    def is_prime(n):
        if n < 2: return False
        for p in [2,3,5,7,11,13,17,19,23,29,31,37]:
            if n == p: return True
            if n % p == 0: return False
        d = n - 1; s = 0
        while d % 2 == 0: d >>= 1; s += 1
        for a in [2,325,9375,28178,450775,9780504,1795265022]:
            if a % n == 0: continue
            x = mod_pow(a, d, n)
            if x == 1 or x == n - 1: continue
            comp = True
            for _ in range(1, s):
                x = x * x % n
                if x == n - 1: comp = False; break
            if comp: return False
        return True

    def rev_digits(x):
        r = 0
        while x > 0:
            r = r * 10 + x % 10
            x //= 10
        return r

    def is_palindrome(x): return x == rev_digits(x)

    need = 50
    limit = 1 << 16
    while True:
        sieve = bytearray([1]) * (limit + 1)
        sieve[0] = sieve[1] = 0
        for p in range(2, int(limit**0.5) + 1):
            if sieve[p]:
                for q in range(p*p, limit+1, p):
                    sieve[q] = 0
        primes = [i for i in range(2, limit+1) if sieve[i]]
        vals = []
        for p in primes:
            sq = p * p
            if is_palindrome(sq): continue
            rev = rev_digits(sq)
            r = math.isqrt(rev)
            if r * r == rev and is_prime(r):
                vals.append(sq)
                if len(vals) >= need: break
        if len(vals) >= need:
            return str(sum(vals[:need]))
        limit <<= 1

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

Java

import java.util.ArrayList;

public class Euler808 {

    static long mulMod(long a, long b, long mod) {
        long x0 = a & 0xFFFFFFL;
        long x1 = (a >>> 24) & 0xFFFFFFL;
        long y0 = b & 0xFFFFFFL;
        long y1 = (b >>> 24) & 0xFFFFFFL;

        long p0 = x0 * y0;
        long p1 = x1 * y0 + x0 * y1;

        return (p0 + (p1 << 24)) % mod;
    }

    static long powMod(long a, long e, long mod) {
        long r = 1 % mod;
        a %= mod;
        while (e > 0) {
            if ((e & 1L) == 1L) {
                // simple modular product since mod is up to ~ 10^16, we can use
                // Math.multiplyHigh, but since values are small, BigInteger or custom is better
                // if overflow happens.
                // Wait, java 9+ has Math.multiplyHigh, but here, mod is at most sqrt(10^16)^2?
                // Wait.
                // For Miller-Rabin, n can be up to 10^16.
                // We can just use BigInteger for powMod to be safe if mod is large.
                r = java.math.BigInteger.valueOf(r).multiply(java.math.BigInteger.valueOf(a))
                        .mod(java.math.BigInteger.valueOf(mod)).longValue();
            }
            a = java.math.BigInteger.valueOf(a).multiply(java.math.BigInteger.valueOf(a))
                    .mod(java.math.BigInteger.valueOf(mod)).longValue();
            e >>= 1L;
        }
        return r;
    }

    static boolean isPrime(long n) {
        if (n < 2)
            return false;
        long[] smallPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
        for (long p : smallPrimes) {
            if (n == p)
                return true;
            if (n % p == 0)
                return false;
        }

        long d = n - 1;
        int s = 0;
        while ((d & 1L) == 0L) {
            d >>= 1L;
            s++;
        }

        long[] witnesses = { 2L, 325L, 9375L, 28178L, 450775L, 9780504L, 1795265022L };
        for (long a : witnesses) {
            if (a % n == 0)
                continue;
            long x = powMod(a, d, n);
            if (x == 1 || x == n - 1)
                continue;
            boolean comp = true;
            for (int r = 1; r < s; r++) {
                x = java.math.BigInteger.valueOf(x).multiply(java.math.BigInteger.valueOf(x))
                        .mod(java.math.BigInteger.valueOf(n)).longValue();
                if (x == n - 1) {
                    comp = false;
                    break;
                }
            }
            if (comp)
                return false;
        }
        return true;
    }

    static ArrayList<Long> sievePrimes(int limit) {
        boolean[] isPrimeVec = new boolean[limit + 1];
        for (int i = 2; i <= limit; i++)
            isPrimeVec[i] = true;

        int r = (int) Math.sqrt(limit);
        for (int p = 2; p <= r; p++) {
            if (isPrimeVec[p]) {
                for (int q = p * p; q <= limit; q += p) {
                    isPrimeVec[q] = false;
                }
            }
        }

        ArrayList<Long> primes = new ArrayList<>();
        for (int i = 2; i <= limit; i++) {
            if (isPrimeVec[i]) {
                primes.add((long) i);
            }
        }
        return primes;
    }

    static long reverseDigits(long x) {
        long r = 0;
        while (x > 0) {
            r = r * 10 + (x % 10);
            x /= 10;
        }
        return r;
    }

    static boolean isPalindrome(long x) {
        return x == reverseDigits(x);
    }

    static long isqrt(long n) {
        if (n < 0)
            return 0;
        long r = (long) Math.sqrt(n);
        while ((r + 1) * (r + 1) <= n && (r + 1) > r) {
            r++;
        }
        while (r * r > n) {
            r--;
        }
        return r;
    }

    static ArrayList<Long> firstReversiblePrimeSquares(int need) {
        for (int limit = 1 << 16;; limit <<= 1) {
            ArrayList<Long> primes = sievePrimes(limit);
            ArrayList<Long> vals = new ArrayList<>();

            for (long p : primes) {
                long sq = p * p;
                if (isPalindrome(sq))
                    continue;
                long rev = reverseDigits(sq);
                long r = isqrt(rev);
                if (r * r == rev && isPrime(r)) {
                    vals.add(sq);
                    if (vals.size() >= need) {
                        return new ArrayList<>(vals.subList(0, need));
                    }
                }
            }
            if (limit >= (1 << 30))
                break;
        }
        return new ArrayList<>();
    }

    public static String solve() {
        ArrayList<Long> vals = firstReversiblePrimeSquares(50);
        long sum = 0;
        for (long v : vals) {
            sum += v;
        }
        return Long.toString(sum);
    }

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