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

A positive integer qualifies if at least four different primes below $100$ divide it; how often each of these primes divides it, and which other prime factors it has, does not matter. Let $F(N)$ be the number of qualifying positive integers less than $N$. The statement gives $F(1000)=23$ and asks for $F(10^{16})$. ^[1]

The whole solution in one sentence: for every set of at least four primes below $100$ we count the multiples of their product below $10^{16}$ and add these counts with the weights $1,-4,10,-20,\dots$, which make every qualifying number count exactly once; only the sets whose product is below $10^{16}$ contribute, and a search that multiplies primes in increasing order lists all of them, about $9.6$ million.

Mathematical Approach

We count by hand below $1000$, see why this simple count fails for larger limits, find and prove weights that repair it, and list the sets of primes that contribute.

1. Eligible primes and multiples

There are $25$ primes below $100$, from $2,3,5,7,11$ to $89,97$; we call them the eligible primes, and for a set $S$ of eligible primes we write $q_S$ for the product of its elements, with $q_S=1$ for the empty set. The smallest qualifying number is $210=2\cdot3\cdot5\cdot7$. The number $420=2^2\cdot3\cdot5\cdot7$ also qualifies, since the prime $2$ counts once however often it divides, and so does $21210=210\cdot101$, since the prime $101$ is simply ignored; $3030=2\cdot3\cdot5\cdot101$ has only three eligible primes and does not qualify.

Two facts turn the problem into counting. First, $x$ is divisible by every prime of $S$ exactly when $q_S$ divides $x$. If $q_S$ divides $x$, so does every prime of $S$. Conversely, let the distinct primes $p_1,\dots,p_k$ all divide $x$, and suppose that $x=p_1\cdots p_j\,y$ for some $j\lt k$. The prime $p_{j+1}$ divides this product but none of the primes $p_1,\dots,p_j$, so by Euclid's lemma, applied to one factor after another, it divides $y$; hence $p_1\cdots p_{j+1}$ divides $x$, and induction on $j$ finishes the proof. ^[2] Second, the multiples of $q$ among $1,2,\dots,M$ are $q,2q,\dots,jq$, where $jq\le M\lt(j+1)q$, so there are $\lfloor M/q\rfloor$ of them. As the problem asks for integers less than $N$, we always take $M=N-1$. For example, $\lfloor 999/210\rfloor=4$ counts $210$, $420$, $630$ and $840$.

2. The checkpoint by hand

Below $1000$, every qualifying number is a multiple of $q_S$ for a set $S$ of four eligible primes with $q_S\le999$, and there are exactly $16$ such sets. The products $210$, $330$, $390$ and $462$ have $4$, $3$, $2$ and $2$ multiples up to $999$, and the other $12$ products, $510$, $546$, $570$, $690$, $714$, $770$, $798$, $858$, $870$, $910$, $930$ and $966$, have one each. No number below $1000$ is a multiple of two of these products, since it would then be divisible by at least five eligible primes, and five distinct primes have a product of at least $2\cdot3\cdot5\cdot7\cdot11=2310$. So all these multiples are different, and

$$F(1000)=4+3+2+2+12=23,$$

as the statement says. The $23$ numbers $210,330,390,420,\dots,966,990$ are also the first $23$ numbers with exactly four distinct prime factors; ^[3] below $1000$ such a number cannot have a prime factor above $100$, because $2\cdot3\cdot5\cdot101=3030$.

3. Overcounting from $2310$ on

For larger limits the same sum counts some numbers more than once. The number $2310=2\cdot3\cdot5\cdot7\cdot11$ has five eligible primes, and leaving out one of them in turn gives five sets of four, so $2310$ is counted five times. In general, a number divisible by exactly $r$ eligible primes is counted once for each set of four of them. The number of ways to choose $k$ of $r$ objects is the binomial coefficient $\binom rk$, ^[4] so this happens $\binom r4$ times; for $30030=2310\cdot13$ it happens $\binom64=15$ times. To correct this we also count the multiples of products of more primes, which is the idea of the inclusion–exclusion principle. ^[5] For $k\ge1$ let

$$T_k=\sum_{|S|=k}\left\lfloor\frac{M}{q_S}\right\rfloor,$$

where $S$ runs over all sets of $k$ eligible primes. By Section 1, $T_k$ counts every $x\le M$ once for each set of $k$ eligible primes that divide $x$: a number with exactly $r$ eligible primes is counted $\binom rk$ times, which is $0$ when $k\gt r$.

4. Finding the weights

We look for weights $a_4,a_5,\dots,a_{25}$ such that $a_4T_4+a_5T_5+\dots+a_{25}T_{25}$ counts every qualifying number exactly once. A number with exactly $r$ eligible primes contributes

$$w(r)=\sum_{k=4}^{r}a_k\binom{r}{k},$$

an empty sum, equal to $0$, when $r\lt4$, as it should be. We need $w(r)=1$ for $4\le r\le25$. Since $\binom rr=1$, the condition for $r$ fixes $a_r$ once the earlier weights are known:

$r=4$: $a_4=1$;
$r=5$: $5a_4+a_5=1$, so $a_5=-4$;
$r=6$: $15a_4+6a_5+a_6=1$, so $a_6=10$;
$r=7$: $35a_4+21a_5+7a_6+a_7=1$, so $a_7=-20$.

The numbers $1,4,10,20$ are $\binom33,\binom43,\binom53,\binom63$, and the signs alternate, which suggests

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

Solving all $22$ conditions by computer confirms this up to $a_{25}=-\binom{24}{3}=-2024$. For example, $2310$ now has the total weight $5\cdot1-4\cdot1=1$, and $30030$ has $15\cdot1-4\cdot6+10\cdot1=1$.

5. Proof for every $r$

Pascal's rule says that $\binom rk=\binom{r-1}{k}+\binom{r-1}{k-1}$. ^[4] With the weights above, $w(4)=\binom33\binom44=1$, and for $r\ge5$ the rule splits $w(r)$ into two sums. The first is $w(r-1)$, because its last term contains $\binom{r-1}{r}=0$. In the second we put $n=r-1\ge4$ and $i=k-1$:

$$\begin{aligned}&w(r)-w(r-1)\\&=\sum_{i=3}^{n}(-1)^{i-3}\binom{i}{3}\binom{n}{i}.\end{aligned}$$

Choosing $i$ of $n$ primes and then $3$ of these $i$ is the same as choosing the $3$ first and then the other $i-3$ from the remaining $n-3$, so ^[4]

$$\binom{n}{i}\binom{i}{3}=\binom{n}{3}\binom{n-3}{i-3}.$$

With $j=i-3$ and the binomial theorem, the difference becomes

$$\begin{aligned}&w(r)-w(r-1)\\&=\binom{n}{3}\sum_{j=0}^{n-3}(-1)^{j}\binom{n-3}{j}\\&=\binom{n}{3}(1-1)^{n-3}=0,\end{aligned}$$

because $n-3\ge1$. So $w(r)=w(r-1)$ for every $r\ge5$, hence $w(r)=1$ for every $r\ge4$, and

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

6. Where the weights come from

The coefficient $\binom{k-1}3$ also has a direct explanation. Group the qualifying numbers by their four smallest eligible primes, a set $Q$ with largest element $q$. A number belongs to the group of $Q$ when it is a multiple of $q_Q$ and is divisible by none of the eligible primes below $q$ outside $Q$; call the set of these excluded primes $B$. Inclusion–exclusion over $B$ gives the size of the group: ^[5]

$$\sum_{U\subseteq B}(-1)^{|U|}\left\lfloor\frac{M}{q_Q\,q_U}\right\rfloor.$$

Adding these sums over all $Q$, a set $S=Q\cup U$ of $k$ primes appears once for each way of choosing $Q$ inside it: $Q$ must contain the largest prime of $S$, because every prime of $U$ is smaller than $q$, and its other three primes can be any three of the remaining $k-1$. So $\lfloor M/q_S\rfloor$ appears $\binom{k-1}3$ times, always with the sign $(-1)^{k-4}$, and the formula of Section 5 follows again. With the $m$ smallest primes in place of four, the same argument counts the numbers divisible by at least $m$ eligible primes with the weights $(-1)^{k-m}\binom{k-1}{m-1}$; for events in place of divisibility this is one of Waring's formulae. ^[6]

7. A model to check by hand

The weights can be watched at work in a smaller model: only the primes $2,3,5,7$, the numbers below $100$, so $M=99$, and at least two primes, so $m=2$ and the weights $(-1)^{k-2}\binom{k-1}{1}$ are $1$, $-2$ and $3$ for $k=2,3,4$. The six products of two primes, $6,10,14,15,21,35$, have $16,9,7,6,4,2$ multiples up to $99$, so $T_2=44$. The products of three primes, $30,42,70,105$, have $3,2,1,0$ multiples, so $T_3=6$, and $T_4=0$ because $210\gt99$. The formula gives

$$T_2-2T_3+3T_4=44-12=32.$$

Six numbers below $100$ are divisible by three of the primes: $30,42,60,70,84,90$. Each of them is counted $\binom32=3$ times in $T_2$ and once in $T_3$, so it gets the weight $3-2=1$; the other $44-18=26$ terms of $T_2$ are the $26$ numbers with exactly two of the primes, and a direct check of all numbers below $100$ also gives $32$. For the real problem the formula gives

$$F(10^4)=943-4\cdot33=811,$$

$$\begin{aligned}F(10^5)&=15626-4\cdot1649\\&\quad+10\cdot25=9280,\end{aligned}$$

from the values of $T_4$, $T_5$ and $T_6$ at these limits, and a direct count agrees.

8. Which sets contribute, and how to list them

A set with $q_S\gt M$ adds $\lfloor M/q_S\rfloor=0$, so only sets with $q_S\le M$ matter. For $N=10^{16}$, every set of at most eight eligible primes has $q_S\le M$, since even the eight largest give $61\cdot67\cdots97\approx1.2\cdot10^{15}$. No set of $14$ or more primes does, since already the $14$ smallest give $2\cdot3\cdots43\approx1.31\cdot10^{16}$. In all, $9\,597\,723$ of the $2^{25}=33\,554\,432$ sets have $q_S\le M$, among them $25\,521$ sets of $13$ primes, and $9\,595\,097$ of them have at least four primes.

These sets are listed by backtracking, a depth-first search through a tree of partial sets. ^[7] It starts from the empty set and adds only primes larger than every prime already in the set, so each set is built exactly once, from its smallest prime upwards. A prime $p$ may be added to $S$ when $q_S\,p\le M$, which for integers is equivalent to $q_S\le\lfloor M/p\rfloor$: the first gives $q_S\le M/p$, hence $q_S\le\lfloor M/p\rfloor$ because $q_S$ is an integer, and the second gives $q_S\,p\le\lfloor M/p\rfloor\,p\le M$. So the test never forms a product larger than $M$. When it fails, every set built from $S\cup\{p\}$ has a product above $M$ too, and every set with $q_S\le M$ is reached, because the sets formed by its smallest primes have even smaller products.

How the Code Works

The C++ program ^[8] starts in main, which reads the options with parse_arguments: --limit= replaces the default limit $10^{16}$, and --skip-checkpoints turns off the checks. An unknown argument or a limit of $0$ ends the program with exit code 1; the digit parser does not check for overflow, so a limit of $2^{64}$ or more wraps around.

primes_below_100 builds the $25$ primes with a linear sieve: for each $i$ from $2$ to $99$ it records $i$ as a prime if $i$ has not been crossed out, and then crosses out $i\,p$ for the primes $p$ found so far, in increasing order, while $i\,p\lt100$, stopping after the first $p$ that divides $i$. Every composite number below $100$ is thus crossed out exactly once, as its smallest prime factor times $i$, which is the property of Euler's sieve. ^[9] binom(n, k) computes a binomial coefficient as a product of $k$ factors divided by $k!$; it is only needed for $\binom{k-1}3\le\binom{24}3=2024$.

solve(limit) returns $0$ for a limit of at most $1$ and sets max_value to $M$, that is limit - 1. The array subset_sums of $26$ unsigned 64-bit counters holds $T_k$ in entry $k$ for $k\ge4$. The recursive lambda dfs(index, taken, product) is the search of Section 8: product is $q_S$, taken is the number of primes in $S$, and index is the position of the first prime that may still be added. If taken is at least $4$, it adds max_value / product to subset_sums[taken]; then, for every later prime $p$, it skips $p$ when product > max_value / p and otherwise calls itself with $q_S\,p$. After the call that starts from the empty set, a loop over $k=4,\dots,25$ forms the sum of Section 5 in a signed 64-bit integer, with the coefficient $\pm$binom(k - 1, 3).

Unless --skip-checkpoints is given, run_checkpoints requires solve(1000) to be $23$ and compares solve(100000) with brute(100000), which tests every $x\lt100000$ by trial division by the $25$ primes; a failure ends the program with exit code 2. Then main prints the count for the chosen limit. The Python and Java versions, with solve_limit and solveLimit, contain the same sieve, binomial and search, with the limit fixed at $10^{16}$ and without options or checks. Python's integers have no size limit, and its recursion is at most $14$ calls deep; Java uses long and keeps the search state in static fields.

The 64-bit types are always wide enough. The search never forms a product above $M$, and $\lfloor M/q_S\rfloor\le M/q_S$ gives $T_k\le M\,e_k$, where $e_k$ is the sum of $1/q_S$ over all sets of $k$ eligible primes, that is, the coefficient of $z^k$ in the product of the $25$ factors $1+z/p$. Multiplying out this product with exact fractions shows that

$$\sum_{k=4}^{25}\binom{k-1}{3}\,e_k=0.41449\ldots$$

is below $0.415$, so every term $a_kT_k$ and every partial sum of the final loop is smaller than $0.415\,M$ in absolute value: below $2^{52}$ for the default limit, and below $2^{63}$ for every limit up to $2^{64}-1$.

Complexity and Verification

The search visits each of the $9\,597\,723$ sets with $q_S\le M$ once and tests every later prime there: $15\,999\,002$ tests in all, each with one division, of which $6\,401\,280$ fail. Since the primes increase, the loop could stop at its first failure, leaving $12\,587\,935$ tests; the code goes on with continue, which costs time but never changes the result. Another $9\,595\,097$ divisions compute the terms $\lfloor M/q_S\rfloor$. The memory is fixed: $25$ primes, $26$ counters and a recursion at most $14$ calls deep. However large the limit, there are at most $2^{25}$ sets. In our runs the C++ program took about $0.03$ seconds including its checkpoints, the Java program about $0.1$ seconds and the Python program about $2$ seconds.

The C++ checkpoints confirm $F(1000)=23$ and compare $F(10^5)$ with a brute-force count. For this article we also verified, with separate programs:

the weights $(-1)^{k-m}\binom{k-1}{m-1}$ give $w(r)=1$ for $r\ge m$ and $w(r)=0$ for $r\lt m$, for all $m\le9$ and $r\le39$;
the formula agrees with a sieve that counts the eligible primes of every $x\lt10^6$, for every $N\le3000$, for $300$ random $N$ up to $10^6$ and for $N=10^4,10^5,10^6$, with $F(10^6)=77\,579$; the C++ program with --limit= prints the same values;
the counts of sets, tests and divisions given above.

Finally, a recursion without any weights confirms the full-size count. Let $g(M,i,s)$ be the number of $y\le M$ divisible by at least $s$ of the eligible primes from the $i$-th one on, and let $p$ be the $i$-th prime. A multiple $p\,y$ of $p$, with $y\le\lfloor M/p\rfloor$, is divisible by the same larger eligible primes as $y$, so counting the multiples of $p$ with the last term and the other numbers with the first two gives

$$\begin{aligned}&g(M,i,s)=g(M,i+1,s)\\&\quad-g(\lfloor M/p\rfloor,i+1,s)\\&\quad+g(\lfloor M/p\rfloor,i+1,s-1),\end{aligned}$$

with $g(M,i,0)=M$; for $s\ge1$, $g=0$ when fewer than $s$ primes remain or $M$ is smaller than the product of the $s$ smallest remaining ones. A memoised evaluation of $g(10^{16}-1,1,4)$, which caches about $1.2\cdot10^7$ values, gives the same count as all three published programs, which agree with each other. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above. The weights and their two proofs, the model example, the bounds on the integer sizes and the operation counts are derived explicitly in this article.

Project Euler 268 — At Least Four Distinct Prime Factors Less Than 100. The official statement says that 23 positive integers less than 1000 are divisible by at least four distinct primes less than 100, and asks how many positive integers less than 10¹⁶ have this property.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes in exactly one way, up to the order of the factors, and Euclid's lemma says that a prime dividing a product of two integers divides one of them; Section 1 uses the lemma to show that distinct primes dividing x have a product dividing x.
OEIS A033993 — Numbers that are divisible by exactly four different primes. The sequence begins 210, 330, 390, 420, 462, …, 966, 990, 1020; its first 23 terms are the 23 numbers below 1000 of Section 2.
Wikipedia — Binomial coefficient. The binomial coefficient C(n, k) is the number of ways to choose k elements from n. The article gives Pascal's rule, the identity C(n, k)·C(k, j) = C(n, j)·C(n − j, k − j) and the alternating row sum Σⱼ (−1)ʲ C(n, j) = 0 for n > 0, the three facts used in Section 5.
Wikipedia — Inclusion–exclusion principle. The principle counts the elements of a union of finite sets, or the elements outside all of them, through the sizes of their intersections with alternating signs; the article applies it to the integers divisible by none of several given primes, as Section 6 does.
Wikipedia — Schuette–Nesbitt formula. Its history section lists Waring's formulae, among them: the probability that at least n of m events occur is Σₖ (−1)^(k − n) C(k − 1, n − 1) Sₖ, where Sₖ is the sum of the probabilities of all intersections of k of the events. For a number chosen uniformly from 1, …, M and n = 4 this is the formula of Section 5 divided by M.
Wikipedia — Backtracking. Backtracking builds candidates step by step, traverses the tree of partial candidates in depth-first order and skips the whole subtree of a partial candidate that cannot be completed, as the search of Section 8 does with sets whose product exceeds M.
C++ — Euler268.cpp. The linked immutable C++ revision contains primes_below_100, solve with its search, brute and run_checkpoints. It is the implementation record for this article.
Wikipedia — Sieve of Eratosthenes. The sieve finds the primes by marking the multiples of each prime as composite; in its variant Euler's sieve each composite number is eliminated exactly once, which takes linear time.

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);
    }
}