Project Euler Problem 574

Solution

Answer: 5780447552057000454

Let

$$P(q)=\prod_{r<q,\ r\text{ prime}} r$$

(the primorial over primes smaller than $q$).

We seek the minimum $A$ such that for a prime $p$,

$$p=A+B \quad \text{or} \quad p=A-B,$$

with:

$\gcd(A,B)=1$,
every prime $r<q$ divides $AB$,
$p<q^2$.

Because $\gcd(A,B)=1$, every prime $r<q$ divides exactly one of $A,B$. Thus for each such prime $r$,

$$A \equiv 0 \pmod r \quad\text{or}\quad A \equiv p \pmod r,$$

since in both sum and difference cases, $B\equiv p-A$ or $A-p$, giving the same congruence condition modulo $r$.

So for fixed $q$, $A$ is determined by a system of congruences:

$$A \equiv \varepsilon_r \pmod r, \qquad \varepsilon_r\in{0,p}.$$

By the Chinese Remainder Theorem, each assignment gives a unique residue class modulo $P(q)$.

For each CRT residue $x$, the valid candidates are:

sum case $p=A+B$:

$A\in[\lceil p/2\rceil,p)$, $A\equiv x\pmod{P(q)}$,

difference case $p=A-B$:

$A>p$, $A\equiv x\pmod{P(q)}$.

We recursively build the CRT solution prime-by-prime and keep the smallest feasible $A$.

The following Python computes $V(p)$ and then

$$S(3800)=\sum_{p<3800}V(p).$$

from math import gcd, isqrt

def sieve(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False

    for i in range(2, isqrt(n) + 1):
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = False

    return [i for i in range(n + 1) if is_prime[i]

primes = sieve(100)

# primes less than q
primes_below_q = {
    q: [r for r in primes if r < q]
    for q in primes
}

def V(p):
    best = float('inf')

    for q in primes:
        if q * q <= p:
            continue

        ps = primes_below_q[q]

        def dfs(i, x, modulus):
            nonlocal best

            # pruning
            if x >= best:
                return

            # finished CRT construction
            if i == len(ps):

                # ----- sum case: p = A + B -----
                low = (p + 1) // 2
                A = x

                if A < low:
                    A += ((low - A + modulus - 1) // modulus) * modulus

                if A < p and gcd(A, p) == 1:
                    best = min(best, A)

                # ----- difference case: p = A - B -----
                A = x

                if A <= p:
                    A += ((p + 1 - A + modulus - 1) // modulus) * modulus

                if gcd(A, p) == 1:
                    best = min(best, A)

                return

            r = ps[i]
            inv = pow(modulus % r, -1, r)
            target = p % r

            # A ≡ 0 (mod r)
            k = ((-x) % r) * inv % r
            dfs(i + 1, x + modulus * k, modulus * r)

            # A ≡ p (mod r)
            k = ((target - x) % r) * inv % r
            dfs(i + 1, x + modulus * k, modulus * r)

        dfs(0, 0, 1)

    return best

def S(n):
    return sum(V(p) for p in sieve(n - 1))

print(S(10))    # 10
print(S(200))   # 7177
print(S(3800))

Verification:

$S(10)=10$ ✔️
$S(200)=7177$ ✔️
$V(37)=22$ ✔️
$V(151)=165$ ✔️

Therefore,

Answer: 302121311939