Solution

Answer: 326227335

Let

$$A_n \equiv i \pmod{p_i}\qquad (1\le i\le n),$$

where $p_i$ is the $i$-th prime, and $A_n$ is the smallest positive solution.

We must compute

$$S(300000)=\sum_{\substack{p\le 300000\ p\mid A_n\text{ for some }n}} p.$$

1. Mathematical analysis

1.1 Constructing $A_n$

Define

$$P_n=\prod_{i=1}^n p_i .$$

Suppose we already know $A_n$.

Then every solution to the first $n$ congruences is

$$A_n+kP_n.$$

To obtain $A_{n+1}$, we need

$$A_n+kP_n \equiv n+1 \pmod{p_{n+1}}.$$

Since $P_n$ is coprime to $p_{n+1}$, there is a unique solution:

$$k \equiv (n+1-A_n),P_n^{-1}\pmod{p_{n+1}}.$$

Choosing the smallest nonnegative $k$ gives $A_{n+1}$.

Thus:

$$A_{n+1}=A_n+t_nP_n,$$

where

$$t_n=((n+1)-A_n),P_n^{-1}\pmod{p_{n+1}}.$$

1.2 Key observation about divisibility

Fix a prime $q=p_k$.

For every $n\ge k$,

$$A_n \equiv k \pmod{p_k},$$

because the defining system already includes the congruence modulo $p_k$.

Therefore:

if $p_k\mid A_n$, then necessarily $n<k$.

So for each prime $p_k$, we only need to track the residues

$$A_1,A_2,\dots,A_{k-1}\pmod{p_k}.$$

1.3 Updating residues efficiently

We never need the gigantic integers $A_n$ themselves.

For each future prime $p_k$, maintain:

$r_k = A_n \bmod p_k$,
$c_k = P_n \bmod p_k$.

Initially:

$$A_1=1,\qquad P_1=2.$$

Hence for all $k>1$:

$$r_k=1,\qquad c_k=2.$$

Now suppose we move from $A_n$ to $A_{n+1}$.

We already computed:

$$A_{n+1}=A_n+t_nP_n.$$

Therefore modulo any future prime $p_k$,

$$r_k \leftarrow r_k+t_n c_k \pmod{p_k},$$

$$c_k \leftarrow c_k,p_{n+1}\pmod{p_k}.$$

Whenever $r_k=0$, we have found a prime dividing some $A_n$.

1.4 Complexity

There are

$$\pi(300000)=25997$$

primes up to $300000$.

The algorithm performs a triangular set of updates:

$$\sum_{i=1}^{m} (m-i)=O(m^2),$$

which is about $3.4\times 10^8$ modular updates.

Using vectorized NumPy operations, this runs quickly in practice.

2. Python implementation

import numpy as np
from math import isqrt

LIMIT = 300000

# ------------------------------------------------------------
# Sieve of Eratosthenes
# ------------------------------------------------------------
is_prime = bytearray(b"\x01") * (LIMIT + 1)
is_prime[0:2] = b"\x00\x00"

for p in range(2, isqrt(LIMIT) + 1):
    if is_prime[p]:
        is_prime[p * p : LIMIT + 1 : p] = (
            b"\x00" * (((LIMIT - p * p) // p) + 1)
        )

primes_list = [i for i, v in enumerate(is_prime) if v]
primes = np.array(primes_list, dtype=np.int64)

m = len(primes)

# ------------------------------------------------------------
# r[k] = current A_n mod p_k
# c[k] = current P_n mod p_k
# ------------------------------------------------------------
r = np.ones(m, dtype=np.int64)

c = np.empty(m, dtype=np.int64)
c[0] = 0
c[1:] = 2

# hit[k] becomes True if p_k divides some A_n
hit = np.zeros(m, dtype=bool)

# ------------------------------------------------------------
# Main recurrence
# ------------------------------------------------------------
for n in range(1, m - 1):

    q = int(primes[n])  # p_(n+1)

    # t_n = ((n+1) - A_n) * P_n^{-1} mod p_(n+1)
    t = ((n + 1) - int(r[n])) * pow(int(c[n]), -1, q) % q

    # Update all future residues simultaneously
    sl = slice(n + 1, m)

    r[sl] = (r[sl] + t * c[sl]) % primes[sl]
    c[sl] = (c[sl] * q) % primes[sl]

    # Any zero residue means divisibility
    hit |= (r == 0)

# ------------------------------------------------------------
# Final answer
# ------------------------------------------------------------
answer = int(primes[hit].sum())

print(answer)

3. Code walkthrough

Prime generation

The sieve builds all primes up to $300000$.

primes_list = [i for i, v in enumerate(is_prime) if v]

This gives:

$$2,3,5,7,11,\dots$$

Residue arrays

r[k]

stores

$$A_n \bmod p_k.$$

c[k]

stores

$$P_n \bmod p_k.$$

These are enough to update all future residues without ever constructing $A_n$.

Computing $t_n$

t = ((n + 1) - int(r[n])) * pow(int(c[n]), -1, q) % q

implements

$$t_n=((n+1)-A_n),P_n^{-1}\pmod{p_{n+1}}.$$

Vectorized updates

r[sl] = (r[sl] + t * c[sl]) % primes[sl]

implements

$$A_{n+1}\equiv A_n+t_nP_n.$$

Similarly:

c[sl] = (c[sl] * q) % primes[sl]

updates

$$P_{n+1}=P_n\cdot p_{n+1}.$$

4. Verification on the example

For $n\le 50$, the program finds exactly:

$$5,\ 23,\ 41$$

as the primes dividing some $A_n$.

Thus

$$S(50)=5+23+41=69,$$

matching the statement.

5. Final verification

The computation produces:

$$S(300000)=326227335.$$

The result is plausible:

only a relatively small fraction of primes divide some $A_n$,
the sum is comfortably below the total sum of all primes up to $300000$,
the sample case matches perfectly.

Answer: 326227335

Project Euler Problem 552