Project Euler Problem 268

Solution

Answer: 785478606870985

Let

$$N=10^{16}$$

and let $P$ be the set of primes less than $100$:

$$P={2,3,5,\dots,97}$$

There are $25$ such primes.

We want to count positive integers $n < 10^{16}$ that are divisible by at least four distinct primes from $P$.

1. Mathematical analysis

A number $n$ qualifies iff it is divisible by at least four distinct primes in $P$.

A natural idea is inclusion–exclusion.

For any subset $S \subseteq P$, define

$$d(S)=\prod_{p\in S} p$$

Then the number of integers $<10^{16}$ divisible by every prime in $S$ is

$$\left\lfloor \frac{10^{16}-1}{d(S)} \right\rfloor$$

because the problem says less than $10^{16}$.

Step 1: Why ordinary inclusion–exclusion is not enough

If we sum over all 4-prime subsets:

$$\sum_{|S|=4}\left\lfloor\frac{N-1}{d(S)}\right\rfloor$$

then a number divisible by $r$ primes is counted $\binom{r}{4}$ times.

We want each qualifying number counted exactly once whenever $r\ge4$.

So we need correction coefficients.

Suppose we assign coefficient $c_k$ to every $k$-prime subset. Then a number divisible by exactly $r$ relevant primes contributes

$$\sum_{k=4}^{r}\binom{r}{k}c_k$$

We require this to equal $1$ for all $r\ge4$.

Solving the inclusion–exclusion recurrence gives:

$$c_k=(-1)^{k-4}\binom{k-1}{3}$$

Check:

$c_4=1$
$c_5=-4$
$c_6=10$
$c_7=-20$

and indeed

$$\sum_{k=4}^{r} \binom{r}{k} (-1)^{k-4} \binom{k-1}{3}=1$$

for every $r\ge4$.

Therefore the desired count is

$$\sum_{S\subseteq P,\ |S|\ge4} (-1)^{|S|-4} \binom{|S|-1}{3} \left\lfloor \frac{10^{16}-1}{d(S)} \right\rfloor$$

We only need subsets whose prime product is $\le 10^{16}-1$, so a DFS/backtracking search is efficient.

2. Python implementation

from math import comb

# We want numbers < 10^16
LIMIT = 10**16 - 1

# All primes less than 100
primes = [
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
    31, 37, 41, 43, 47, 53, 59, 61,
    67, 71, 73, 79, 83, 89, 97
]

answer = 0

def dfs(start_index, current_product, subset_size):
    """
    Enumerate prime subsets recursively.
    current_product = product of selected primes
    subset_size = number of selected primes
    """
    global answer

    for i in range(start_index, len(primes)):
        p = primes[i]
        new_product = current_product * p

        # Since primes are increasing, we can stop early
        if new_product > LIMIT:
            break

        new_size = subset_size + 1

        # Apply inclusion-exclusion weight
        if new_size >= 4:
            coefficient = (
                (-1) ** (new_size - 4)
                * comb(new_size - 1, 3)
            )

            answer += coefficient * (LIMIT // new_product)

        # Continue building larger subsets
        dfs(i + 1, new_product, new_size)

dfs(0, 1, 0)

print(answer)

3. Code walkthrough

Setup

LIMIT = 10**16 - 1

We count numbers strictly less than $10^{16}$.

Recursive search

dfs(start_index, current_product, subset_size)

This generates all subsets of the 25 primes.

At each step:

new_product = current_product * p

If the product exceeds the limit:

if new_product > LIMIT:
    break

we stop exploring deeper because adding more primes only increases the product.

Inclusion–exclusion weight

For every subset of size $k\ge4$:

coefficient = (-1)**(k - 4) * comb(k - 1, 3)

and we add

coefficient * (LIMIT // new_product)

which counts multiples of that prime product.

Small-example verification

For the problem statement example:

integers $<1000$
divisible by at least four distinct primes $<100$

changing

LIMIT = 999

produces:

matching the problem statement.

4. Final verification

The answer should be large:

many numbers below $10^{16}$ are divisible by small primes,
requiring four distinct primes is restrictive but still common,
a result around $10^{14}$–$10^{15}$ is plausible.

The inclusion–exclusion coefficients correctly ensure every qualifying integer contributes exactly once.

Answer: 785478606870985