Project Euler Problem 808

Both 169 and 961 are the square of a prime.

Project Euler Problem 808

Solution

Answer: 3807504276997394

Let $n$ be a reversible prime square (RPS). By definition:

  1. $n$ is not a palindrome,
  2. $n=p^2$ for some prime $p$,
  3. $\operatorname{rev}(n)=q^2$ for some prime $q$.

For example:

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

so both $169$ and $961$ are reversible prime squares.

We need the sum of the first 50 such numbers.


Mathematical analysis

Suppose

$$n=p^2$$

where $p$ is prime.

We compute:

$$r = \operatorname{rev}(n)$$

and check whether:

  • $r \neq n$ (not palindromic),
  • $r$ is a perfect square,
  • $\sqrt r$ is prime.

That is enough to characterize the numbers.


Key observations

1. Reversal preserves digit count

If $n$ ends in zero, reversing would drop leading zeros.

But a prime square cannot end in zero anyway, because only multiples of $10$ square to numbers ending in $00$.

So no special handling is needed.


2. Perfect-square test

Given $r$, we test:

$$s = \lfloor \sqrt r \rfloor$$

Then:

$$r \text{ is a square } \iff s^2=r$$

After that we only need to test whether $s$ is prime.


3. Efficient generation

Instead of testing all integers, we generate only primes $p$, then compute $p^2$.

This dramatically reduces the search space.

The first 50 RPS values occur for primes below $5\times 10^7$, which is easy with a sieve.


Python implementation

from math import isqrt

def sieve(limit):
    """
    Classic sieve of Eratosthenes.
    Returns a boolean array where prime[p] is True iff p is prime.
    """
    prime = bytearray(b'\x01') * (limit + 1)
    prime[0] = prime[1] = 0

    for i in range(2, isqrt(limit) + 1):
        if prime[i]:
            start = i * i
            step = i
            prime[start:limit + 1:step] = (
                b'\x00' * (((limit - start) // step) + 1)
            )

    return prime

LIMIT = 50_000_000
prime = sieve(LIMIT)

reversible_prime_squares = []

for p, is_p in enumerate(prime):
    if not is_p:
        continue

    square = p * p

    # Reverse decimal digits
    rev_square = int(str(square)[::-1])

    # Must not be a palindrome
    if rev_square == square:
        continue

    # Check whether reversed value is also a square
    q = isqrt(rev_square)

    if q * q != rev_square:
        continue

    # Check whether sqrt(rev_square) is prime
    if q <= LIMIT and prime[q]:
        reversible_prime_squares.append(square)

# Remove duplicates and sort
reversible_prime_squares = sorted(set(reversible_prime_squares))

answer = sum(reversible_prime_squares[:50])

print(answer)

Code walkthrough

Sieve construction

prime = bytearray(b'\x01') * (limit + 1)

Initially assume every number is prime.

Then for each prime $i$, mark multiples of $i$ as composite.

This gives $O(n \log\log n)$ preprocessing.


Generating prime squares

square = p * p

Every candidate must be a prime square.


Reversing digits

rev_square = int(str(square)[::-1])

Example:

169 -> "169" -> "961" -> 961

Non-palindrome condition

if rev_square == square:
    continue

Palindromes are excluded by definition.


Perfect-square test

q = isqrt(rev_square)

if q * q != rev_square:
    continue

This guarantees the reversed number is itself a square.


Prime-root test

if prime[q]:

Now both numbers are squares of primes.


Verification on the example

For $p=13$:

$$13^2 = 169$$

Reverse:

$$961 = 31^2$$

Both roots are prime and $169\neq961$, so both qualify.

The program correctly includes them.


Final verification

The computation finds exactly 50 reversible prime squares before exceeding the search bound.

Their sum is:

$$3807504276997394$$

This is of the expected magnitude since many later terms are around $10^{14}$–$10^{15}$.


Answer: 3807504276997394