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, so $p_1=2$, $p_2=3$, $p_3=5$, $p_4=7$, and let $r_n$ be the remainder when $(p_n-1)^n+(p_n+1)^n$ is divided by $p_n^2$. For $n=3$ the statement computes $4^3+6^3=280$, and $280=11\cdot25+5$, so $r_3=5$. We are told that the least $n$ for which the remainder exceeds $10^9$ is $n=7037$, and we must find the least $n$ for which it exceeds $10^{10}$. ^[1]

The whole solution in one sentence: when the two powers are multiplied out, every term except the two lowest is a multiple of $p_n^2$, so the remainder is $2$ for even $n$ and exactly $2n\,p_n$ for odd $n\ge5$, and we only have to walk through the primes in order until $2n\,p_n$ first exceeds $10^{10}$ at an odd index.

Mathematical Approach

We compute the first remainders by hand, explain their pattern with the binomial theorem, prove when the simple formula is exact, and check it on the value $7037$ from the statement.

1. Small cases

Here are the first remainders, each sum written as a multiple of $p_n^2$ plus the remainder:

$n=1$, $p_1=2$: $1+3=4=1\cdot4+0$, so $r_1=0$;
$n=2$, $p_2=3$: $2^2+4^2=20=2\cdot9+2$, so $r_2=2$;
$n=3$, $p_3=5$: $4^3+6^3=280=11\cdot25+5$, so $r_3=5$;
$n=4$, $p_4=7$: $6^4+8^4=5392=110\cdot49+2$, so $r_4=2$;
$n=5$, $p_5=11$: $10^5+12^5=348832$${}=2882\cdot121+110$, so $r_5=110$;
$n=6$, $p_6=13$: $12^6+14^6=10515520$${}=62222\cdot169+2$, so $r_6=2$;
$n=7$, $p_7=17$: $16^7+18^7=880655488$${}=3047250\cdot289+238$, so $r_7=238$.

The sums grow very fast, but the remainders stay small. Every even $n$ gives $2$, and from $n=5$ on every odd $n$ gives $2n\,p_n$: $110=2\cdot5\cdot11$, $238=2\cdot7\cdot17$, and the next odd index gives $r_9=414=2\cdot9\cdot23$. The two smallest odd indices are different: there $2n\,p_n$ would be $4$ and $30$, but $r_1=0$ and $r_3=5$. These remainders form the sequence A379544 of the OEIS. ^[2]

2. Remainders and congruences

For an integer $x$ and a divisor $m\ge1$ there is exactly one way to write $x=qm+s$ with integers $q$ and $0\le s\lt m$; the remainder $s$ is written $x\bmod m$. Two integers are congruent modulo $m$, written $x\equiv y\pmod m$, when $m$ divides $x-y$. For example $280\equiv5\pmod{25}$, because $280-5=275=11\cdot25$. Congruences can be added and multiplied like equations, so in a calculation modulo $m$ every multiple of $m$ may be dropped at any moment. ^[3]

3. Only two terms of each power survive

Write $p$ for $p_n$. The binomial theorem multiplies out a power of a sum: ^[4]

$$\begin{aligned}(p+1)^n&=\sum_{k=0}^{n}\binom nk p^k,\\(p-1)^n&=\sum_{k=0}^{n}\binom nk p^k(-1)^{n-k}.\end{aligned}$$

Here $\binom nk=\frac{n!}{k!\,(n-k)!}$, with $n!=1\cdot2\cdots n$ and $0!=1$, counts the ways to choose the $k$ factors, out of $n$, that contribute $p$; for $n=3$ these numbers are $1,3,3,1$. Every term with $k\ge2$ contains the factor $p^2$, so modulo $p^2$ only the terms $k=0$ and $k=1$ remain. In the second sum they are $(-1)^n$ and $np\,(-1)^{n-1}=-(-1)^n\,np$, hence

$$\begin{aligned}(p+1)^n&\equiv1+np,\\(p-1)^n&\equiv(-1)^n(1-np)\pmod{p^2}.\end{aligned}$$

For $n=3$ and $p=5$: $6^3=1+15+75+125=216$ and $4^3=-1+15-75+125=64$. Dropping the multiples $75$ and $125$ of $25$ leaves $16=1+3\cdot5$ and $14=-(1-3\cdot5)$, and indeed $216=8\cdot25+16$ and $64=2\cdot25+14$.

4. Even and odd indices

When the two congruences are added, the terms $\pm np$ cancel for even $n$, and the terms $\pm1$ cancel for odd $n$. So the sum is congruent to $2$ for even $n$ and to $2np$ for odd $n$, modulo $p^2$.

For even $n$ we have $p\ge3$, so $2$ already lies between $0$ and $p^2-1$, and the remainder is exactly $r_n=2$. For odd $n$, divide $2n$ by $p$, say $2n=qp+s$ with $0\le s\lt p$. Then $2np=q\,p^2+s\,p$ with $0\le sp\lt p^2$, so for odd $n$

$$r_n=p_n\cdot\bigl(2n\bmod p_n\bigr).$$

For $n=1$ this gives $2\cdot(2\bmod2)=0$, for $n=3$ it gives $5\cdot(6\bmod5)=5$, and for $n=5$ it gives $11\cdot10=110$, the values of Section 1.

5. The prime outruns $2n$

By Section 4, $r_n=2n\,p_n$ for odd $n$ exactly when $2n\lt p_n$. This fails for the first four primes, since $2\le2$, $3\le4$, $5\le6$ and $7\le8$, but it holds from then on.

Claim. $p_n\gt2n$ for every $n\ge5$.

Proof. The primes up to $2n$ are $2$ and the odd primes among the $n$ odd numbers $1,3,5,\dots,2n-1$. Two of these odd numbers are not prime: $1$, and $9$, which is on the list because $2n-1\ge9$. So at most $n-2$ odd primes, and at most $n-1$ primes in all, are smaller than $2n$, and the $n$-th prime must be larger than $2n$.

Hence $2n\bmod p_n=2n$ for every odd $n\ge5$, and Section 4 gives the exact value $r_n=2n\,p_n$.

6. The least index

Let $L\ge5$ be a limit. The indices $n=1,2,3,4$ give the remainders $0,2,5,2$, and every even $n$ gives $2$, so none of them exceeds $L$. For odd $n\ge5$ the remainder $2n\,p_n$ grows strictly with $n$, because both factors grow. Hence the least $n$ with $r_n\gt L$ is the first odd $n\ge5$ with $2n\,p_n\gt L$, and from there on every odd index exceeds $L$ as well. The answer is therefore always odd; for $5\le L\le109$, for instance, it is $n=5$, because $r_5=110$.

7. Worked example: the checkpoint

For $L=10^9$ the primes at the end of the search are $p_{7035}=71023$, $p_{7036}=71039$ and $p_{7037}=71059$. By Section 5,

$$\begin{aligned}r_{7035}&=2\cdot7035\cdot71023=999\,293\,610,\\r_{7037}&=2\cdot7037\cdot71059=1\,000\,084\,366,\end{aligned}$$

while $r_{7036}=2$. The first value is below $10^9$, and so are the remainders of all earlier odd indices, because they grow with $n$; the second value is above. So $7037$ is the least index, as the statement says. Without Section 5 we would have to reduce two powers of $34141$ digits each.

8. How far the search runs

Rosser's theorem says that $p_n\gt n\ln n$ for every $n\ge1$, where $\ln$ is the natural logarithm. ^[5] Hence $2n\,p_n\gt2n^2\ln n$, and $2n^2\ln n$ already exceeds $10^{10}$ at the odd index $n=22345$. The search for $10^{10}$ therefore stops at an odd $n\le22345$, and it needs only primes up to $p_{22345}=253811$. At its end the compared number $2n\,p_n$ exceeds $10^{10}$, far more than the largest signed 32-bit integer $2^{31}-1=2147483647$; the first odd index with $2n\,p_n\gt2^{31}-1$ is $n=10117$, where $2n\,p_n=2148102142$. So the product must be formed in 64-bit arithmetic.

How the Code Works

None of the three programs computes a power: they walk through the primes in increasing order and apply Section 6. The C++ program ^[6] keeps its settings in the structure Options: limit $=10^{10}$ and run_checkpoints = true. parse_arguments accepts --skip-checkpoints and --limit= followed by digits, which parse_i64_after_prefix accumulates in an i64 without an overflow check, so a value above $2^{63}-1$ has undefined behaviour; any other argument, or a limit below $1$, ends the program with exit code $1$.

solve(limit) builds the primes itself by trial division. ^[7] It stores the primes found so far in a std::vector<int> with room reserved for $300000$ of them. Each candidate $2,3,4,\dots$ is divided by the stored primes $p$ in increasing order while $p^2\le$ candidate; a zero remainder marks a composite, and a candidate that survives is prime, because a composite number has a prime factor no larger than its square root. A new prime is appended, and its index $n$ is the length of the list. Even $n$ are skipped by the test (n & 1) == 0, since their remainder is $2$ (Section 4); for odd $n$ the program forms remainder = 2LL * static_cast<i64>(n) * p in the 64-bit type i64 and returns $n$ as soon as it exceeds the limit (Sections 5 and 6). run_checkpoints first requires solve(1000000000) $=7037$, the value from the statement, and otherwise ends the program with exit code $2$; then main prints solve(limit).

The program uses $2n\,p_n$ also for $n=1$ and $n=3$, where the true remainders are $0$ and $5$ (Section 5). This is harmless for the default limit and the checkpoint, but a small --limit gives a wrong index: the limits $5$ to $29$ make it print $3$ instead of $5$, and the limits $1$, $2$, $3$ make it print $1$ instead of $2$, $3$, $3$; the limit $4$ and every limit from $30$ up to about $4.5\cdot10^{17}$ are handled correctly. Beyond that the needed primes would pass $2^{31}-1$, and the int counter candidate would overflow first.

The Python version has the limit $10^{10}$ built in, without options or a checkpoint. Its function solve starts from p = 2 and runs n through range(1, 1000000), a bound that Section 8 shows is never reached. At an odd n it prints n and returns as soon as 2 * n * p > limit; otherwise it moves on with SymPy's nextprime, which, apart from a short cached list of small primes, tests the candidates of the form $6j\pm1$ with SymPy's primality test isprime. ^[8] Python integers have unlimited precision, so no overflow can occur. The Java version repeats the loop of solve in main, with an ArrayList<Integer> of primes and the same fixed limit; the cast in (long) p * p and the long literal 2L in 2L * n * c make Java evaluate both products in 64-bit arithmetic, whereas a 32-bit product would overflow silently from $n=10117$ on. ^[9]

Complexity and Verification

For the default limit let $P$ be the prime at which the search stops; by Section 8, $P\le253811$. In C++ and Java an even index costs one parity check, and an odd index two multiplications and a comparison more; the real work is the trial division. A candidate $c$ needs at most $\pi(\sqrt c)$ divisions, where $\pi(x)$ counts the primes up to $x$, so all candidates together need at most $P\,\pi(\sqrt P)$, which is less than $2.6\,P^{3/2}/\ln P$ because $\pi(x)\lt1.25506\,x/\ln x$ for $x\gt1$. ^[5] For $P=253811$ the sum of all $\pi(\sqrt c)$ is about $1.7\cdot10^7$, but the real count is far smaller, because almost three quarters of the composites are rejected by $2$ or $3$ at once: the search for $10^{10}$ needs about $2.2$ million divisions, and the checkpoint search for $10^9$, which the C++ program runs first, needs $484988$. Both programs store the at most $22345$ primes found; the C++ reservation of $300000$ four-byte int entries takes $1.2$ MB. The Python version keeps only the current prime and calls nextprime once per index, so its cost is that of SymPy's primality tests. A sieve of Eratosthenes up to $253811$ would need $O(P\log\log P)$ operations, ^[10] but at this size each program finishes in well under a second on a current computer.

The C++ program checks itself with the checkpoint from the statement, and all three programs print the same index. Separate computations for this article, not part of the published solutions, confirmed the following. We computed the true remainders for every $n\le33860$, that is, for every prime below $4\cdot10^5$, by modular exponentiation, which reduces modulo $p_n^2$ after every multiplication; ^[11] they equal $2$ for even $n$ and $2n\,p_n$ for odd $n$, except for $n=1$ and $n=3$, and the first $57$ agree with the data of A379544. Searching these true remainders gives $7037$ for $10^9$ and, for $10^{10}$, the index that the three programs print. Finally, the C++ program with --limit agrees with the true remainders for the limits $4$, $30$ to $399$ and $10^3,10^4,\dots,10^9$, and differs exactly for the small limits listed above.

Footnotes and References

The references below supply the background facts used above. The formula for the remainder, the bound $p_n\gt2n$ and the stopping rule are derived explicitly in this article.

Project Euler 123 — Prime Square Remainders. The official statement: pₙ is the n-th prime and r is the remainder when (pₙ − 1)ⁿ + (pₙ + 1)ⁿ is divided by pₙ²; for n = 3, 4³ + 6³ = 280 ≡ 5 (mod 25). The least n for which the remainder first exceeds 10⁹ is 7037, and the task asks for the least n for which it first exceeds 10¹⁰.
OEIS A379544. The sequence a(n) = ((p − 1)ⁿ + (p + 1)ⁿ) mod p², with p the n-th prime: 0, 2, 5, 2, 110, 2, 238, …, given up to n = 57, with the example n = 3 and a link to Problem 123.
Wikipedia — Modular arithmetic. Two integers are congruent modulo m when their difference is a multiple of m; congruence is compatible with addition, subtraction, multiplication and exponentiation, and a mod m is the unique r with 0 ≤ r < m and r ≡ a (mod m).
Wikipedia — Binomial theorem. (x + y)ⁿ is the sum of C(n, k) xⁿ⁻ᵏ yᵏ for k = 0, …, n, where C(n, k) = n!/(k!(n − k)!).
Wikipedia — Rosser's theorem. J. Barkley Rosser proved in 1939 that the n-th prime satisfies pₙ > n log n for every n ≥ 1, with the natural logarithm. Wikipedia — Prime-counting function: π(x) < 1.25506 x / log x for every x > 1, and π(x) is asymptotic to x / log x by the prime number theorem.
C++ — Euler123.cpp. The linked immutable C++ revision contains parse_arguments, solve and run_checkpoints with the checkpoint for 10⁹. It is the implementation record for this article.
Wikipedia — Trial division. Prime divisors up to the square root of a number suffice, because a composite number has a factor no larger than its square root.
SymPy — nextprime and isprime. nextprime(n, ith=1) returns the ith prime greater than n, and its search uses that potential primes are located at 6j ± 1; isprime is definitive for n < 2⁶⁴.
Java Language Specification — Types, Values, and Variables: int and long are 32-bit and 64-bit two's-complement integers, a multiplication with a long operand is carried out in 64-bit precision, and integer operators do not indicate overflow. Python — numeric types: integers have unlimited precision.
Wikipedia — Sieve of Eratosthenes. The sieve marks the multiples of each prime as composite and finds all primes up to n with O(n log log n) operations and O(n) memory.
Wikipedia — Modular exponentiation. Reducing modulo m after every multiplication keeps the numbers small, and with exponentiation by squaring bᵉ mod m takes O(log e) steps.

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;
            }
        }
    }
}