Problem 123: Prime Square Remainders

Project Euler Problem 123 Solution

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

Problem Summary Let $p_n$ be the $n$-th prime, and define $r_n$ as the least nonnegative remainder of $$\left((p_n-1)^n+(p_n+1)^n\right)\bmod p_n^2.$$ The task is to find the smallest index $n$ for which $r_n\gt 10^{10}$. At first sight this looks like a large modular-exponentiation problem, but modulo $p_n^2$ almost every binomial term disappears. That collapse is the whole reason the problem becomes easy. Mathematical Approach The remainder is governed entirely by the first two terms of the binomial expansion. Once that is written out carefully, the problem stops being about huge powers and becomes a search over prime indices. Binomial expansion modulo $p_n^2$ Expand both powers around $p_n$: $$ (p_n\pm 1)^n=\sum_{k=0}^{n}\binom{n}{k}p_n^k(\pm 1)^{\,n-k}. $$ Adding the two expansions gives $$ (p_n-1)^n+(p_n+1)^n = \sum_{k=0}^{n}\binom{n}{k}p_n^k\left((-1)^{n-k}+1\right). $$ Any term with $k\ge 2$ is already divisible by $p_n^2$, so modulo $p_n^2$ only the cases $k=0$ and $k=1$ can survive. Which of those survives depends only on the parity of $n$. Even indices and odd indices behave differently If $n$ is even, then the constant terms add and the linear terms cancel, so $$ r_n\equiv 2 \pmod{p_n^2}. $$ Because $0\le 2\lt p_n^2$, the actual remainder is simply $r_n=2$. Even indices can therefore never exceed $10^{10}$....

Detailed mathematical approach

Problem Summary

Let $p_n$ be the $n$-th prime, and define $r_n$ as the least nonnegative remainder of

$$\left((p_n-1)^n+(p_n+1)^n\right)\bmod p_n^2.$$

The task is to find the smallest index $n$ for which $r_n\gt 10^{10}$. At first sight this looks like a large modular-exponentiation problem, but modulo $p_n^2$ almost every binomial term disappears. That collapse is the whole reason the problem becomes easy.

Mathematical Approach

The remainder is governed entirely by the first two terms of the binomial expansion. Once that is written out carefully, the problem stops being about huge powers and becomes a search over prime indices.

Binomial expansion modulo $p_n^2$

Expand both powers around $p_n$:

$$ (p_n\pm 1)^n=\sum_{k=0}^{n}\binom{n}{k}p_n^k(\pm 1)^{\,n-k}. $$

Adding the two expansions gives

$$ (p_n-1)^n+(p_n+1)^n = \sum_{k=0}^{n}\binom{n}{k}p_n^k\left((-1)^{n-k}+1\right). $$

Any term with $k\ge 2$ is already divisible by $p_n^2$, so modulo $p_n^2$ only the cases $k=0$ and $k=1$ can survive. Which of those survives depends only on the parity of $n$.

Even indices and odd indices behave differently

If $n$ is even, then the constant terms add and the linear terms cancel, so

$$ r_n\equiv 2 \pmod{p_n^2}. $$

Because $0\le 2\lt p_n^2$, the actual remainder is simply $r_n=2$. Even indices can therefore never exceed $10^{10}$.

If $n$ is odd, then the constant terms cancel and the linear terms add, giving

$$ r_n\equiv 2np_n \pmod{p_n^2}. $$

So for odd indices the remainder is controlled by the single quantity $2np_n$, but at this stage it is still a congruence rather than an exact equality.

Why the congruence becomes the exact remainder

The implementations compare $2np_n$ directly against the limit, so we need to know when $2np_n$ is already smaller than $p_n^2$. A simple counting argument is enough. For $n\ge 5$, the odd numbers

$$ 1,3,5,\dots,2n-1 $$

form a list of exactly $n$ numbers. Among them, both $1$ and $9$ are non-prime, so there are at most $n-2$ odd primes below $2n$. Including the prime $2$, there are at most $n-1$ primes not exceeding $2n$. Therefore the $n$-th prime must satisfy

$$ p_n\gt 2n \qquad (n\ge 5). $$

For odd $n\ge 5$ this implies

$$ 0\lt 2np_n\lt p_n^2, $$

so the least nonnegative residue of $2np_n \pmod{p_n^2}$ is exactly $2np_n$. Hence

$$ r_n=2np_n. $$

The first solution occurs far beyond $n=5$, so this is the formula the code can safely search with.

Worked example: $n=3$ versus $n=5$

The distinction between congruence and exact remainder matters for tiny odd indices. For $n=3$, we have $p_3=5$, so

$$ (5-1)^3+(5+1)^3=4^3+6^3=280. $$

Modulo $25$, this gives $280\bmod 25=5$. The linear formula gives $2np_n=2\cdot 3\cdot 5=30$, and indeed $30\equiv 5\pmod{25}$, but the exact remainder is $5$, not $30$.

For $n=5$, however, $p_5=11$ and

$$ (11-1)^5+(11+1)^5=10^5+12^5=348832. $$

Now $348832\bmod 121=110$, and that equals

$$ 2\cdot 5\cdot 11=110. $$

Here no wrap-around occurs because $2n\lt p_n$. The target threshold is so large that the first qualifying index lies in this stable regime.

Estimating the search range

The prime number theorem gives the rough estimate $p_n\sim n\log n$. Substituting that into the odd-index formula yields

$$ 2n^2\log n \gtrsim 10^{10}. $$

This places the crossing a little above $2\times 10^4$. So the program only needs to generate primes up to roughly the corresponding $p_n$, which is on the order of a few hundred thousand. That makes a straightforward sequential search entirely adequate.

How the Code Works

Generate primes in order

The C++, Python, and Java implementations all walk through the prime sequence in increasing order so that the current prime and its index are known simultaneously. The C++ and Java versions keep a growing list of earlier primes and test each new integer candidate only by those primes up to its square root. The Python version advances from one prime to the next through a prime-generation routine, but conceptually it performs the same ordered scan.

Check only odd indices

Each time a new prime $p_n$ is reached, the implementation first checks whether $n$ is odd. Even indices are skipped immediately because their remainder is always $2$. For odd indices, it evaluates the simple expression $2np_n$ and stops as soon as this exceeds the limit. None of the implementations ever compute $(p_n-1)^n$ or $(p_n+1)^n$ directly. One implementation also verifies the logic against the smaller threshold $10^9$ before running the full search.

Complexity Analysis

Let $P$ be the prime where the search stops. After a prime is known, the mathematical test is constant time: parity check plus one multiplication and one comparison.

In the C++ and Java implementations, the real cost is prime generation by incremental trial division. A candidate $c$ is tested only by previously found primes up to $\sqrt c$, so a standard worst-case estimate is

$$ O\!\left(\sum_{c=2}^{P}\pi(\sqrt c)\right)=O\!\left(\frac{P^{3/2}}{\log P}\right), $$

with memory usage $O(\pi(P))$ for the stored primes. For this problem $P$ is only a few $10^5$, so that cost is small in practice. The Python implementation delegates prime stepping to a library routine, but from the problem's point of view it still performs one constant-time remainder test per prime index.

Footnotes and References

Problem page: Project Euler 123 - Prime square remainders
Binomial theorem: Wikipedia - Binomial theorem
Modular arithmetic: Wikipedia - Modular arithmetic
Prime number theorem: Wikipedia - Prime number theorem
Prime numbers: Wikipedia - Prime number

Problem 123 source code

C++

#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = std::int64_t;

struct Options {
    i64 limit = 10000000000LL;
    bool run_checkpoints = true;
};

bool parse_i64_after_prefix(const std::string& arg, const std::string& prefix, i64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    i64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<i64>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_i64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.limit >= 1;
}

int solve(const i64 limit) {
    std::vector<int> primes;
    primes.reserve(300000);

    for (int candidate = 2;; ++candidate) {
        bool prime = true;
        for (int p : primes) {
            if (static_cast<i64>(p) * p > candidate) {
                break;
            }
            if (candidate % p == 0) {
                prime = false;
                break;
            }
        }
        if (!prime) {
            continue;
        }

        primes.push_back(candidate);
        const int n = static_cast<int>(primes.size());
        if ((n & 1) == 0) {
            continue;
        }

        const i64 p = candidate;
        const i64 remainder = 2LL * static_cast<i64>(n) * p;
        if (remainder > limit) {
            return n;
        }
    }
}

bool run_checkpoints() {
    if (solve(1000000000LL) != 7037) {
        std::cerr << "Checkpoint failed for limit=1e9" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

# Problem 123: Prime square remainders
# Find n where r_n = ((p_n-1)^n + (p_n+1)^n) mod p_n^2 > 10^10.
# For odd n: r_n = 2*n*p_n. Find smallest odd n where 2*n*p_n > 10^10.

from sympy import nextprime

def solve():
    limit = 10**10
    p = 2
    for n in range(1, 1000000):
        if n % 2 == 1:
            if 2 * n * p > limit:
                print(n)
                return
        p = nextprime(p)

solve()

Java

import java.util.*;

public class Euler123 {
    public static void main(String[] args) {
        long limit = 10000000000L;
        List<Integer> primes = new ArrayList<>();
        for (int c = 2;; c++) {
            boolean prime = true;
            for (int p : primes) {
                if ((long) p * p > c)
                    break;
                if (c % p == 0) {
                    prime = false;
                    break;
                }
            }
            if (!prime)
                continue;
            primes.add(c);
            int n = primes.size();
            if ((n & 1) == 0)
                continue;
            if (2L * n * c > limit) {
                System.out.println(n);
                return;
            }
        }
    }
}