Project Euler Problem 355

Solution

Answer: 1726545007

Let

$$\operatorname{Co}(n)$$

be the maximum possible sum of a subset of ${1,2,\dots,n}$ whose elements are pairwise coprime.

The key observation is that pairwise coprime means:

No prime factor may appear in more than one selected number.

So the problem is really about assigning primes to selected integers.

Step 1: A natural baseline

For every prime $p \le n$, if we use $p$ at all, the best standalone contribution is the largest power of $p$ not exceeding $n$:

$$M(p)=\max{p^k\le n}.$$

Since prime powers of distinct primes are mutually coprime, a valid baseline set is

$${1}\cup {M(p): p\le n},$$

with sum

$$B = 1 + \sum_{p\le n} M(p).$$

Example for $n=10$:

$$M(2)=8,\quad M(3)=9,\quad M(5)=5,\quad M(7)=7$$

giving

$$1+8+9+5+7=30,$$

which is optimal.

However, sometimes combining primes into one composite beats separate prime powers.

Example ($n=30$):

separate: $16+25=41$ for primes $2,5$,
combined: $27?$ impossible (contains only $3$),
but $2\cdot 13=26$ may beat $16+13=29$? no.

In general, for any squarefree prime set $S$, define:

$$V(S)=\max{x\le n:\text{ prime factors of }x\text{ are exactly }S}.$$

Replacing the separate prime powers by $V(S)$ changes the score by

$$\Delta(S) = V(S)-\sum_{p\in S} M(p).$$

Only sets with $\Delta(S)>0$ matter.

The optimization becomes:

Choose disjoint prime sets $S$ maximizing total improvement $\sum \Delta(S)$.

This is a weighted set-packing / binary integer optimization problem:

variable $x_S\in{0,1}$,
each prime may appear in at most one chosen set,
maximize total gain.

Since $n=200000$, the number of profitable sets is modest (~25k), and modern MILP solves it instantly.

Python implementation

from sympy import primerange
from scipy.optimize import milp, LinearConstraint, Bounds
import scipy.sparse as sp
import numpy as np

def Co(n):
    primes = list(primerange(2, n + 1))

    # Largest prime power <= n for each prime
    M = {}
    for p in primes:
        x = p
        while x * p <= n:
            x *= p
        M[p] = x

    baseline = 1 + sum(M.values())

    # Smallest prime factor sieve
    spf = list(range(n + 1))
    for i in range(2, int(n**0.5) + 1):
        if spf[i] == i:
            for j in range(i * i, n + 1, i):
                if spf[j] == j:
                    spf[j] = i

    # For each prime set, keep the largest number <= n using exactly those primes
    best = {}

    for x in range(2, n + 1):
        y = x
        prime_set = []

        while y > 1:
            p = spf[y]
            prime_set.append(p)
            while y % p == 0:
                y //= p

        prime_set = tuple(prime_set)

        if x > best.get(prime_set, 0):
            best[prime_set] = x

    # Keep only profitable replacements
    candidates = []

    for prime_set, x in best.items():
        gain = x - sum(M[p] for p in prime_set)

        if gain > 0 and len(prime_set) >= 2:
            candidates.append((gain, prime_set))

    # Build sparse constraint matrix:
    # each prime may be used at most once
    p_index = {p: i for i, p in enumerate(primes)}

    rows, cols, data = [], [], []

    for j, (_, prime_set) in enumerate(candidates):
        for p in prime_set:
            rows.append(p_index[p])
            cols.append(j)
            data.append(1)

    A = sp.coo_matrix(
        (data, (rows, cols)),
        shape=(len(primes), len(candidates))
    ).tocsr()

    # MILP:
    # maximize gains == minimize negative gains
    c = -np.array([gain for gain, _ in candidates], dtype=float)

    constraints = LinearConstraint(
        A,
        -np.inf * np.ones(len(primes)),
        np.ones(len(primes))
    )

    integrality = np.ones(len(candidates), dtype=int)
    bounds = Bounds(
        np.zeros(len(candidates)),
        np.ones(len(candidates))
    )

    result = milp(
        c=c,
        constraints=constraints,
        integrality=integrality,
        bounds=bounds
    )

    return int(round(baseline - result.fun))

# Verification
print(Co(30))   # 193
print(Co(100))  # 1356

# Final answer
print(Co(200000))

Code walkthrough

Compute $M(p)$

For each prime, find its largest power $\le n$. These form the baseline solution. 2. Factor every integer $x\le n$

Using a smallest-prime-factor sieve, determine the distinct prime factors of $x$. 3. Record best representative for each prime set

For a prime subset $S$, keep only the largest number using exactly $S$. 4. Compute gains

If replacing separate prime powers with a composite increases the sum, record the positive gain. 5. MILP optimization

Binary variable per candidate replacement. Constraints enforce that each prime appears at most once. 6. Add gains to baseline

Final optimum:

$$\operatorname{Co}(n) = B+\max \sum \Delta(S).$$

Verification

The program reproduces the given checks:

$\operatorname{Co}(30)=193$
$\operatorname{Co}(100)=1356$

For $n=200000$, it computes:

Answer: 1726545007