Problem 268: At Least Four Distinct Prime Factors Less Than 100

Project Euler Problem 268 Solution

EulerSolve provides an optimized solution for Project Euler Problem 268, At Least Four Distinct Prime Factors Less Than 100, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must count integers in the half-open range $1 \le x < N$ with $N=10^{16}$ such that $x$ is divisible by at least four distinct primes below $100$. Multiplicity does not matter: a prime like $2$ contributes only once, whether $x$ is divisible by $2$ or by $2^{20}$. The C++ program exposes the bound as --limit=<N> , so the same method works for any $N\ge 1$. Mathematical Approach 1. The Relevant Prime Set Let $$\mathcal P=\{p \text{ prime}: p<100\}.$$ There are $25$ such primes. For a subset $S\subseteq\mathcal P$, define its squarefree product $$q_S=\prod_{p\in S}p.$$ An integer $x$ is divisible by every prime in $S$ exactly when $q_S\mid x$. Therefore the number of integers $x<N$ divisible by all primes of $S$ is $$\left\lfloor\frac{N-1}{q_S}\right\rfloor.$$ The code uses $N-1$ because the searched range is $1\le x<N$, not $1\le x\le N$. 2. Grouping by Subset Size For each $k\ge4$, define $$T_k=\sum_{|S|=k}\left\lfloor\frac{N-1}{q_S}\right\rfloor.$$ So $T_k$ counts, with multiplicity, how many numbers are divisible by every prime in some $k$-element subset of $\mathcal P$. If a number $x$ has exactly $r$ distinct prime divisors from $\mathcal P$, then: 1. it appears in $T_k$ exactly $\binom{r}{k}$ times, because we may choose any $k$ of those $r$ primes, 2....

Detailed mathematical approach

Problem Summary

We must count integers in the half-open range $1 \le x < N$ with $N=10^{16}$ such that $x$ is divisible by at least four distinct primes below $100$. Multiplicity does not matter: a prime like $2$ contributes only once, whether $x$ is divisible by $2$ or by $2^{20}$.

The C++ program exposes the bound as --limit=<N>, so the same method works for any $N\ge 1$.

Mathematical Approach

1. The Relevant Prime Set

Let

$$\mathcal P=\{p \text{ prime}: p<100\}.$$

There are $25$ such primes. For a subset $S\subseteq\mathcal P$, define its squarefree product

$$q_S=\prod_{p\in S}p.$$

An integer $x$ is divisible by every prime in $S$ exactly when $q_S\mid x$. Therefore the number of integers $x<N$ divisible by all primes of $S$ is

$$\left\lfloor\frac{N-1}{q_S}\right\rfloor.$$

The code uses $N-1$ because the searched range is $1\le x<N$, not $1\le x\le N$.

2. Grouping by Subset Size

For each $k\ge4$, define

$$T_k=\sum_{|S|=k}\left\lfloor\frac{N-1}{q_S}\right\rfloor.$$

So $T_k$ counts, with multiplicity, how many numbers are divisible by every prime in some $k$-element subset of $\mathcal P$.

If a number $x$ has exactly $r$ distinct prime divisors from $\mathcal P$, then:

1. it appears in $T_k$ exactly $\binom{r}{k}$ times, because we may choose any $k$ of those $r$ primes,

2. it contributes nothing to $T_k$ when $k>r$.

So the whole problem is to find coefficients $a_k$ such that

$$\sum_{k=4}^{r} a_k\binom{r}{k}= \begin{cases} 0,&r<4,\\ 1,&r\ge4. \end{cases}$$

3. The Correct Inclusion-Exclusion Weights

The program uses

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

Thus the final answer is

$$\sum_{k=4}^{25} a_kT_k =\sum_{k=4}^{25}(-1)^{k-4}\binom{k-1}{3}T_k.$$

The first few coefficients are

$$a_4=1,\qquad a_5=-4,\qquad a_6=10,\qquad a_7=-20,\dots$$

So the formula starts as

$$T_4-4T_5+10T_6-20T_7+\cdots.$$

4. Why Those Weights Give “At Least Four”

Take a number with exactly $r$ eligible prime factors. Its total weight becomes

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

For $r<4$, this sum is empty, so $w(r)=0$, exactly as required.

For $r\ge4$, there is a neat closed-form evaluation. First use

$$\binom{r}{k}\binom{k-1}{3}=\binom{r}{4}\binom{r-4}{k-4}\frac{4}{k}.$$

Now write $j=k-4$. Then

$$w(r)=4\binom{r}{4}\sum_{j=0}^{r-4}(-1)^j\binom{r-4}{j}\frac{1}{j+4}.$$

Using $\frac{1}{j+4}=\int_0^1 x^{j+3}\,dx$, we get

$$w(r)=4\binom{r}{4}\int_0^1 x^3(1-x)^{r-4}\,dx.$$

This is a beta integral:

$$\int_0^1 x^3(1-x)^{r-4}\,dx=\frac{3!(r-4)!}{r!}.$$

Therefore

$$w(r)=4\cdot\frac{r!}{4!(r-4)!}\cdot\frac{3!(r-4)!}{r!}=1.$$

So every number with at least four distinct eligible primes contributes exactly once, and every number with fewer than four contributes zero.

5. A Concrete Overcount Example

The number $2310=2\cdot3\cdot5\cdot7\cdot11$ has $r=5$ eligible prime factors.

It appears in $T_4$ exactly $\binom54=5$ times, once for each choice of four of its five primes. It appears in $T_5$ exactly $\binom55=1$ time. It appears in no $T_k$ with $k\ge6$.

Its final weighted contribution is therefore

$$1\cdot5+(-4)\cdot1=1.$$

This is the simplest way to see why the coefficients correct the heavy overcount from $T_4$.

6. DFS Enumeration and Pruning

The program never enumerates all $2^{25}$ subsets explicitly. Instead it performs a depth-first search over increasing prime indices.

At a DFS state, the code knows:

1. the next prime index it may use,

2. how many primes have already been taken,

3. the current squarefree product $q_S$.

If taken >= 4, the code immediately adds

$$\left\lfloor\frac{N-1}{q_S}\right\rfloor$$

to subset_sums[taken].

Before multiplying by a new prime $p$, it checks

$$q_S \le \frac{N-1}{p},$$

implemented as product > max_value / p. This avoids overflow and prunes any branch whose product already exceeds $N-1$.

7. Small Checkpoints

The first checkpoint is

$$\texttt{solve}(1000)=23.$$

Indeed, the numbers below $1000$ with at least four distinct primes below $100$ are exactly $23$ in number; examples include $210,330,390,462,\dots,990$.

The second checkpoint compares

$$\texttt{solve}(100000)$$

against a brute-force routine that explicitly counts distinct prime divisors for every $x<100000$. This confirms that the weighted subset formula and the DFS implementation match the direct definition.

How the Code Works

primes_below_100() generates the 25 primes less than $100$ with a small linear sieve.

solve(limit) sets max_value = limit - 1, runs the DFS, stores all $T_k$-type totals in subset_sums, then forms the weighted sum with coefficients $(-1)^{k-4}\binom{k-1}{3}$.

binom is only used for those small coefficients, so 64-bit signed integers are sufficient.

The command-line interface accepts --limit=<N> and --skip-checkpoints.

Complexity Analysis

The theoretical subset space has size $2^{25}$, but in practice most branches are cut immediately once the squarefree product becomes too large. The main cost is therefore the number of valid subset products $q_S\le N-1$, not the full power set.

The final weighted sum over $k=4,\dots,25$ is tiny compared with the DFS itself, and memory usage is $O(25)$ beyond the recursion stack because only the prime list and the 26-element accumulator array are stored.

Problem 268 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = long long;

struct Options {
    u64 limit = 10000000000000000ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<u64>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 1;
}

std::vector<int> primes_below_100() {
    constexpr int max_n = 100;
    std::array<bool, max_n + 1> composite{};
    std::vector<int> primes;
    for (int i = 2; i < max_n; ++i) {
        if (!composite[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
        for (int p : primes) {
            const int v = i * p;
            if (v >= max_n) {
                break;
            }
            composite[static_cast<std::size_t>(v)] = true;
            if (i % p == 0) {
                break;
            }
        }
    }
    return primes;
}

i64 binom(const int n, const int k) {
    if (k < 0 || k > n) {
        return 0;
    }
    if (k == 0 || k == n) {
        return 1;
    }
    i64 numer = 1;
    i64 denom = 1;
    for (int i = 1; i <= k; ++i) {
        numer *= (n - k + i);
        denom *= i;
    }
    return numer / denom;
}

u64 solve(const u64 limit) {
    if (limit <= 1) {
        return 0;
    }

    const std::vector<int> primes = primes_below_100();
    std::array<u64, 26> subset_sums{};

    const u64 max_value = limit - 1;

    const auto dfs = [&](const auto& self, const std::size_t index, const int taken, const u64 product) -> void {
        if (taken >= 4) {
            subset_sums[static_cast<std::size_t>(taken)] += max_value / product;
        }

        if (index >= primes.size()) {
            return;
        }

        for (std::size_t i = index; i < primes.size(); ++i) {
            const u64 p = static_cast<u64>(primes[i]);
            if (product > max_value / p) {
                continue;
            }
            self(self, i + 1, taken + 1, product * p);
        }
    };

    dfs(dfs, 0U, 0, 1ULL);

    i64 answer = 0;
    for (int k = 4; k <= static_cast<int>(primes.size()); ++k) {
        const i64 coeff = ((k - 4) % 2 == 0 ? 1 : -1) * binom(k - 1, 3);
        answer += coeff * static_cast<i64>(subset_sums[static_cast<std::size_t>(k)]);
    }

    return static_cast<u64>(answer);
}

u64 brute(const int limit) {
    const std::vector<int> primes = primes_below_100();
    u64 count = 0;
    for (int x = 1; x < limit; ++x) {
        int distinct = 0;
        int t = x;
        for (int p : primes) {
            if (t % p == 0) {
                ++distinct;
                while (t % p == 0) {
                    t /= p;
                }
            }
        }
        if (distinct >= 4) {
            ++count;
        }
    }
    return count;
}

bool run_checkpoints() {
    if (solve(1000ULL) != 23ULL) {
        std::cerr << "Checkpoint failed for limit 1000" << '\n';
        return false;
    }
    if (solve(100000ULL) != brute(100000)) {
        std::cerr << "Checkpoint failed for brute-force cross-check at 100000" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

def primes_below_100():
    max_n = 100
    composite = [False] * (max_n + 1)
    primes = []
    for i in range(2, max_n):
        if not composite[i]:
            primes.append(i)
        for p in primes:
            v = i * p
            if v >= max_n:
                break
            composite[v] = True
            if i % p == 0:
                break
    return primes

def binom(n, k):
    if k < 0 or k > n:
        return 0
    if k == 0 or k == n:
        return 1
    numer = 1
    denom = 1
    for i in range(1, k + 1):
        numer *= (n - k + i)
        denom *= i
    return numer // denom

def solve_limit(limit):
    if limit <= 1:
        return 0

    primes = primes_below_100()
    subset_sums = [0] * 26
    max_value = limit - 1

    def dfs(index, taken, product):
        if taken >= 4:
            subset_sums[taken] += max_value // product

        if index >= len(primes):
            return

        for i in range(index, len(primes)):
            p = primes[i]
            if product > max_value // p:
                continue
            dfs(i + 1, taken + 1, product * p)

    dfs(0, 0, 1)

    answer = 0
    for k in range(4, len(primes) + 1):
        coeff = (1 if (k - 4) % 2 == 0 else -1) * binom(k - 1, 3)
        answer += coeff * subset_sums[k]

    return answer

def solve():
    limit = 10000000000000000
    ans = solve_limit(limit)
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;

public class Euler268 {

    static List<Integer> primesBelow100() {
        int max_n = 100;
        boolean[] composite = new boolean[max_n + 1];
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i < max_n; i++) {
            if (!composite[i])
                primes.add(i);
            for (int p : primes) {
                int v = i * p;
                if (v >= max_n)
                    break;
                composite[v] = true;
                if (i % p == 0)
                    break;
            }
        }
        return primes;
    }

    static long binom(int n, int k) {
        if (k < 0 || k > n)
            return 0;
        if (k == 0 || k == n)
            return 1;
        long numer = 1;
        long denom = 1;
        for (int i = 1; i <= k; i++) {
            numer *= (n - k + i);
            denom *= i;
        }
        return numer / denom;
    }

    static long[] subsetSums;
    static long maxValue;
    static List<Integer> primes;

    static void dfs(int index, int taken, long product) {
        if (taken >= 4) {
            subsetSums[taken] += maxValue / product;
        }

        if (index >= primes.size())
            return;

        for (int i = index; i < primes.size(); i++) {
            long p = primes.get(i);
            if (product > maxValue / p)
                continue;
            dfs(i + 1, taken + 1, product * p);
        }
    }

    static long solveLimit(long limit) {
        if (limit <= 1)
            return 0;

        primes = primesBelow100();
        subsetSums = new long[26];
        maxValue = limit - 1;

        dfs(0, 0, 1L);

        long answer = 0;
        for (int k = 4; k <= primes.size(); k++) {
            long coeff = (((k - 4) % 2 == 0) ? 1 : -1) * binom(k - 1, 3);
            answer += coeff * subsetSums[k];
        }
        return answer;
    }

    public static void main(String[] args) {
        long ans = solveLimit(10000000000000000L);
        System.out.println(ans);
    }
}