Project Euler Problem 123

Solution

Answer: 21035

Let

$$r_n = (p_n-1)^n + (p_n+1)^n \pmod{p_n^2},$$

where $p_n$ is the $n$-th prime.

We seek the least $n$ such that

$$r_n > 10^{10}.$$

---

Mathematical analysis

The key observation is that the expression can be simplified dramatically using the binomial theorem.

Let $p=p_n$.

We expand:

$$(p+1)^n = \sum_{k=0}^{n} \binom{n}{k} p^k$$

and

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

We only care modulo $p^2$, so every term containing $p^2$ or higher vanishes.

Thus modulo $p^2$,

$$(p+1)^n \equiv 1 + np$$

and

$$(p-1)^n \equiv (-1)^n + n(-1)^{n-1}p.$$

Adding:

$$r_n \equiv (1+(-1)^n) + np\left(1+(-1)^{n-1}\right)\pmod{p^2}.$$

Now split into parity cases.

Case 1: $n$ even

Then

$$(-1)^n=1,\qquad (-1)^{n-1}=-1.$$

So

$$r_n \equiv 2 + np(1-1)=2 \pmod{p^2}.$$

Hence for even $n$,

$$r_n = 2.$$

This can never exceed $10^{10}$.

---

Case 2: $n$ odd

Then

$$(-1)^n=-1,\qquad (-1)^{n-1}=1.$$

So

$$r_n \equiv 0 + np(1+1)=2np \pmod{p^2}.$$

Therefore, for odd $n$,

$$r_n = 2p_n n$$

(as long as $2np_n < p_n^2$, which is true in the relevant range).

So the problem reduces to finding the smallest odd $n$ such that

$$2n p_n > 10^{10}.$$

---

Estimating the range

Using the Prime Number Theorem,

$$p_n \sim n\log n.$$

Thus

$$2n p_n \sim 2n^2\log n.$$

We solve approximately:

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

This suggests $n$ is around $2\times 10^4$.

We now compute exactly.

---

Python implementation

from math import isqrt

LIMIT = 10**10

def generate_primes(limit_count):
    """
    Generate the first 'limit_count' primes using
    a standard incremental primality test.
    """
    primes = [2]
    candidate = 3

    while len(primes) < limit_count:
        is_prime = True
        root = isqrt(candidate)

        for p in primes:
            if p > root:
                break
            if candidate % p == 0:
                is_prime = False
                break

        if is_prime:
            primes.append(candidate)

        candidate += 2

    return primes

# We know the answer is around 2e4, so 30000 primes is plenty.
primes = generate_primes(30000)

answer = None

for n in range(1, len(primes) + 1, 2):   # only odd n matter
    p = primes[n - 1]                    # p_n (1-indexed)

    remainder = 2 * n * p

    if remainder > LIMIT:
        answer = n
        break

print(answer)

---

Code walkthrough

Prime generation

def generate_primes(limit_count):

Creates the first limit_count primes.

---

primes = [2]
candidate = 3

Start with the first prime and test only odd numbers afterward.

---

root = isqrt(candidate)

We only test divisibility up to $\sqrt{\text{candidate}}$.

---

if candidate % p == 0:

If divisible by any earlier prime, it is composite.

---

if is_prime:
    primes.append(candidate)

Store newly found primes.

---

Main search

for n in range(1, len(primes) + 1, 2):

Only odd $n$ can produce large remainders.

---

p = primes[n - 1]

Python lists are zero-indexed, but $p_n$ is one-indexed.

---

remainder = 2 * n * p

Using the derived formula

$$r_n = 2np_n.$$

---

if remainder > LIMIT:

Stop at the first $n$ exceeding $10^{10}$.

---

Verification on the example

For $n=3$:

• $p_3 = 5$
• Formula gives

$$r = 2 \cdot 3 \cdot 5 = 30.$$

Modulo $25$,

$$30 \equiv 5 \pmod{25},$$

matching the problem statement:

$$4^3 + 6^3 = 280 \equiv 5 \pmod{25}.$$

Correct.

Also, the problem says the least $n$ for $10^9$ is $7037$. Plugging into the formula:

$$2 \cdot 7037 \cdot p_{7037}$$

indeed first exceeds $10^9$, confirming consistency.

---

Final verification

The solution must be:

• odd (even $n$ give remainder $2$),
• around $2\times10^4$ from asymptotics,
• and the first threshold crossing.

Running the computation gives:

$$n = 21035.$$

---

Answer: 21035