Project Euler Problem 381

Solution

Answer: 139602943319822

For primes $p \ge 5$, define

$$S(p)=\sum_{k=1}^{5}(p-k)! \pmod p.$$

We must compute

$$\sum_{5\le p<10^8} S(p),$$

where the sum runs over primes.

---

1. Mathematical analysis

The key idea is to use Wilson’s theorem:

$$(p-1)! \equiv -1 \pmod p$$

for every prime $p$.

We repeatedly “divide backward” to express the other factorials in terms of $(p-1)!$.

---

Step 1: Express each factorial modulo $p$

Since

$$(p-1)!=(p-1)(p-2)!,$$

and $p-1\equiv -1 \pmod p$,

$$- (p-2)! \equiv -1 \pmod p,$$

so

$$(p-2)! \equiv 1 \pmod p.$$

Next,

$$(p-2)!=(p-2)(p-3)! \equiv (-2)(p-3)!,$$

hence

$$(p-3)! \equiv -\frac12 \pmod p.$$

Similarly,

$$(p-4)! \equiv \frac16 \pmod p,$$

and

$$(p-5)! \equiv -\frac1{24} \pmod p.$$

Therefore,

$$S(p)\equiv -1+1-\frac12+\frac16-\frac1{24}.$$

The first two terms cancel:

$$S(p)\equiv -\frac12+\frac16-\frac1{24}.$$

Using denominator $24$,

$$S(p)\equiv \frac{-12+4-1}{24} = -\frac9{24} = -\frac38 \pmod p.$$

So the entire problem reduces to

$$\boxed{S(p)\equiv -3\cdot 8^{-1}\pmod p}.$$

This is the crucial simplification.

---

Step 2: Efficient computation

Instead of factorials, we only need a modular inverse:

S(p) = (-3 * pow(8, -1, p)) % p

for each prime $p$.

The complexity is then dominated by prime generation up to $10^8$.

There are about $5.76\times10^6$ primes below $10^8$, so a sieve of Eratosthenes is entirely feasible.

---

2. Python implementation

from math import isqrt

LIMIT = 10**8

# Sieve of Eratosthenes (odd numbers only)
n = LIMIT
sieve = bytearray(b"\x01") * (n // 2)
sieve[0] = 0  # 1 is not prime

for i in range(3, isqrt(n) + 1, 2):
    if sieve[i // 2]:
        start = i * i
        step = 2 * i
        sieve[start // 2 :: i] = b"\x00" * (((n - start - 1) // step) + 1)

total = 0

# Iterate over odd primes
for p in range(5, n, 2):
    if sieve[p // 2]:
        total += (-3 * pow(8, -1, p)) % p

print(total)

---

3. Code walkthrough

Sieve construction

sieve = bytearray(b"\x01") * (n // 2)

We store only odd numbers to halve memory usage.

• index i represents number 2*i + 1.

---

Crossing out composites

for i in range(3, isqrt(n) + 1, 2):

We test only odd candidate primes.

If i is prime:

sieve[start // 2 :: i] = ...

marks all odd multiples of i starting from $i^2$.

---

Computing $S(p)$

Using the derived formula:

(-3 * pow(8, -1, p)) % p

where

pow(8, -1, p)

computes the modular inverse of $8$ modulo $p$.

---

Verification on the example

For $p=7$:

$$8^{-1}\equiv 1 \pmod 7,$$

so

$$S(7)\equiv -3\equiv 4 \pmod 7,$$

matching the statement.

Also,

sum(S(p) for 5 <= p < 100) = 480

which agrees with the problem’s check value.

---

4. Final verification

• The derivation from Wilson’s theorem is exact.
• The reduction to a modular inverse eliminates all large factorials.
• The implementation uses a standard prime sieve and modular arithmetic only.
• The known check value $480$ is reproduced exactly.

Running the program yields:

$$139602943319822$$

---

Answer: 139602943319822