Problem 590: Sets with a Given Least Common Multiple

Project Euler Problem 590 Solution

EulerSolve provides an optimized solution for Project Euler Problem 590, Sets with a Given Least Common Multiple, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a positive integer $n$, let $H(n)$ be the number of sets of positive integers whose least common multiple is exactly $n$. The problem then defines $$L(n)=\operatorname{lcm}(1,2,\dots,n),\qquad HL(n)=H(L(n)),$$ and asks for $HL(50000)\bmod 10^9$. If a set has lcm $N$, then every element of that set must divide $N$. So the problem is really about counting subsets of the divisor set of $N=L(50000)$ whose lcm is exactly $N$. The challenge is that $L(50000)$ has many prime factors and an enormous number of divisors, so direct enumeration is completely infeasible. Mathematical Approach Write a general target number as $$N=\prod_{i=1}^{k} p_i^{e_i},\qquad d_i=e_i+1.$$ Each divisor of $N$ is determined by choosing, for every prime $p_i$, an exponent from $0$ to $e_i$. Hence $d_i$ is the number of choices contributed by the $i$-th prime. Step 1: Turn the Problem into a Subset Count Let $\mathcal{D}(N)$ be the set of all positive divisors of $N$. Any valid set is a subset of $\mathcal{D}(N)$, and because we are counting sets rather than sequences, each divisor is either present or absent exactly once. For the lcm of a chosen subset to equal $N$, every prime power $p_i^{e_i}$ must appear at full strength in at least one selected divisor....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, let $H(n)$ be the number of sets of positive integers whose least common multiple is exactly $n$. The problem then defines

$$L(n)=\operatorname{lcm}(1,2,\dots,n),\qquad HL(n)=H(L(n)),$$

and asks for $HL(50000)\bmod 10^9$.

If a set has lcm $N$, then every element of that set must divide $N$. So the problem is really about counting subsets of the divisor set of $N=L(50000)$ whose lcm is exactly $N$. The challenge is that $L(50000)$ has many prime factors and an enormous number of divisors, so direct enumeration is completely infeasible.

Mathematical Approach

Write a general target number as

$$N=\prod_{i=1}^{k} p_i^{e_i},\qquad d_i=e_i+1.$$

Each divisor of $N$ is determined by choosing, for every prime $p_i$, an exponent from $0$ to $e_i$. Hence $d_i$ is the number of choices contributed by the $i$-th prime.

Step 1: Turn the Problem into a Subset Count

Let $\mathcal{D}(N)$ be the set of all positive divisors of $N$. Any valid set is a subset of $\mathcal{D}(N)$, and because we are counting sets rather than sequences, each divisor is either present or absent exactly once.

For the lcm of a chosen subset to equal $N$, every prime power $p_i^{e_i}$ must appear at full strength in at least one selected divisor. In other words, for each prime $p_i$, at least one chosen divisor must use exponent $e_i$ at that prime.

Step 2: Apply Inclusion-Exclusion to Force Every Maximal Exponent

For each $i$, let $A_i$ be the bad event that no selected divisor reaches exponent $e_i$ at prime $p_i$. If a set avoids all bad events, then its lcm is exactly $N$.

Now choose a subset $S\subseteq\{1,\dots,k\}$ of primes that are forced to miss their top exponent. For $i\in S$, a divisor may use only $0,1,\dots,e_i-1$, so it has $d_i-1$ choices. For $i\notin S$, it still has all $d_i$ choices. Therefore the number of divisors compatible with these restrictions is

$$M(S)=\prod_{i\in S}(d_i-1)\prod_{i\notin S}d_i.$$

Any subset of those $M(S)$ divisors satisfies all restrictions in $S$, so there are $2^{M(S)}$ such subsets. Inclusion-exclusion gives

$$\boxed{H(N)=\sum_{S\subseteq\{1,\dots,k\}}(-1)^{|S|}2^{M(S)}.}$$

The empty subset is harmless here: it lies in every bad event, so inclusion-exclusion cancels it automatically.

Step 3: Specialize to $L(n)$ and Group Equal Local Factors

For $N=L(n)$, the exponent of a prime $p\le n$ is

$$e_p=\max\{a\ge 0 : p^a\le n\},$$

$$d_p=e_p+1.$$

Many primes have the same $d_p$. That matters because the formula above depends on a prime only through $d_p$. Suppose one particular value $d$ occurs for exactly $c$ primes. If exactly $i$ of those $c$ primes belong to $S$, then:

$$(-1)^i\binom{c}{i}$$

is the signed combinatorial coefficient, and

$$(d-1)^i d^{\,c-i}$$

is the factor contributed to $M(S)$.

If the distinct $d$-values are $d_1,\dots,d_t$ with multiplicities $c_1,\dots,c_t$, then the subset sum collapses to

$$\boxed{HL(n)=\sum_{0\le i_j\le c_j}\left(\prod_{j=1}^{t}(-1)^{i_j}\binom{c_j}{i_j}\right)\,2^{\prod_{j=1}^{t}(d_j-1)^{i_j}d_j^{\,c_j-i_j}}.}$$

This is the key reduction: the naive formula has one binary decision per prime, while the grouped formula has only one integer choice per distinct $d$-value. Since $d_p=e_p+1$ and $e_p\le \lfloor \log_2 n\rfloor$, the number of groups is only $O(\log n)$.

Step 4: Reduce Huge Powers of Two Modulo $10^9$

The exponent

$$E=\prod_{j=1}^{t}(d_j-1)^{i_j}d_j^{\,c_j-i_j}$$

is enormous, so the implementation never stores it exactly. Instead it only stores the information needed to evaluate $2^E \bmod 10^9$.

Because

$$10^9=2^9\cdot 5^9,$$

we distinguish two cases. If $E<9$, then $2^E$ can be computed directly. If $E\ge 9$, then the residue modulo $2^9$ is already fixed, and modulo $5^9$ the powers of $2$ repeat with multiplicative order

$$\operatorname{ord}_{5^9}(2)=4\cdot 5^8=1{,}562{,}500.$$

Therefore, for $E\ge 9$,

$$2^E \bmod 10^9 = 2^{\,9+((E-9)\bmod 1{,}562{,}500)} \bmod 10^9.$$

So each grouped contribution only needs two pieces of data: the exponent modulo $1{,}562{,}500$, and a capped indicator telling whether the true exponent is still below $9$.

Step 5: Worked Example with $HL(4)$

We have

$$L(4)=\operatorname{lcm}(1,2,3,4)=12=2^2\cdot 3.$$

Hence the local choice counts are $d=3$ for prime $2$ and $d=2$ for prime $3$. There are four inclusion-exclusion cases:

$$\begin{aligned} S=\varnothing &:&& M=3\cdot 2=6,\\ S=\{2\} &:&& M=2\cdot 2=4,\\ S=\{3\} &:&& M=3\cdot 1=3,\\ S=\{2,3\} &:&& M=2\cdot 1=2. \end{aligned}$$

$$HL(4)=H(12)=2^6-2^4-2^3+2^2=64-16-8+4=44.$$

This is exactly the small checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the grouped inclusion-exclusion formula above. They first sieve the primes up to $n$, then compute $e_p$ and $d_p=e_p+1$ for each prime using repeated integer division, which avoids any floating-point issues.

Next, they count how many times each $d$-value occurs and build one choice table per group. For every possible value $i=0,\dots,c$ in a group of size $c$, the table stores the signed binomial coefficient, the corresponding exponent factor modulo $1{,}562{,}500$, and the capped information needed to detect whether the true exponent is still below $9$.

After that, the implementation combines all groups except the largest into aggregate states. Each aggregate state represents many original subsets, but it keeps only the data required for the final modular evaluation. The last accumulation step pairs those aggregate states with the choices from the largest group and converts the stored exponent information into the correct value of $2^E \bmod 10^9$.

The C++ implementation parallelizes that outer accumulation over several threads. The Python implementation performs the same arithmetic serially. The Java implementation acts as a thin wrapper around the same compiled computation, so all three language versions are numerically consistent.

Complexity Analysis

Let the distinct grouped values be $d_1,\dots,d_t$ with multiplicities $c_1,\dots,c_t$. A naive inclusion-exclusion over primes would require $2^{\pi(n)}$ cases, which is hopeless for $n=50000$. The grouped method reduces that to roughly

$$\prod_{j=1}^{t}(c_j+1)$$

combined choices, plus sieve and preprocessing work.

The prime sieve costs $O(n\log\log n)$ time and $O(n)$ memory. Building the per-group binomial data costs $O\!\left(\sum_j c_j^2\right)$ time. The final state-combination phase is proportional to the grouped state count above, with extra practical speedup in the multithreaded C++ version. Memory usage is $O(n)$ for the sieve plus the stored grouped states.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=590
Least common multiple: Wikipedia - Least common multiple
Divisor function: Wikipedia - Divisor function
Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
Multiplicative order: Wikipedia - Multiplicative order

Problem 590 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <string>
#include <thread>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using u8 = std::uint8_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

constexpr u32 MOD = 1'000'000'000U;
constexpr u32 POW2_THRESHOLD = 9U;
constexpr u32 POW2_CYCLE = 1'562'500U;  // ord_{5^9}(2) = 4 * 5^8

struct Options {
    int n = 50'000;
    bool run_checkpoints = true;
    unsigned requested_threads = 0U;
};

struct ChoiceData {
    u32 coef = 0U;      // (-1)^i * C(count, i) mod MOD
    u32 exp_mod = 0U;   // factor exponent mod POW2_CYCLE
    u8 exp_cap = 0U;    // min(real exponent factor, 9)
};

struct GroupData {
    int d = 0;          // d = e_p + 1
    int count = 0;      // how many primes have this d
    std::vector<ChoiceData> choices;
};

struct RestState {
    u32 coef = 0U;
    u32 exp_mod = 0U;
    u8 exp_cap = 1U;
};

struct Pow2Tables {
    std::array<u32, POW2_THRESHOLD> small{};
    std::vector<u32> cycle_shifted;  // cycle_shifted[r] = 2^(9+r) mod MOD
};

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

    long long parsed = 0LL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10LL + static_cast<long long>(c - '0');
        if (parsed > static_cast<long long>(std::numeric_limits<int>::max())) {
            return false;
        }
    }

    value = static_cast<int>(parsed);
    return true;
}

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

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
        if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
            return false;
        }
    }

    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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_int_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

unsigned pick_thread_count(const unsigned requested) {
    if (requested > 0U) {
        return requested;
    }
    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0U) {
        hw = 4U;
    }
    return hw;
}

u8 cap_mul_9(const u8 a, const u8 b) {
    const unsigned product = static_cast<unsigned>(a) * static_cast<unsigned>(b);
    if (product >= POW2_THRESHOLD) {
        return static_cast<u8>(POW2_THRESHOLD);
    }
    return static_cast<u8>(product);
}

u32 pow_mod_u32(u32 base, u64 exp, const u32 mod) {
    if (mod == 1U) {
        return 0U;
    }

    u64 result = 1ULL % mod;
    u64 cur = base % mod;
    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = (result * cur) % mod;
        }
        cur = (cur * cur) % mod;
        exp >>= 1ULL;
    }
    return static_cast<u32>(result);
}

const Pow2Tables& get_pow2_tables() {
    static const Pow2Tables tables = []() {
        Pow2Tables t;
        t.small[0] = 1U;
        for (u32 i = 1U; i < POW2_THRESHOLD; ++i) {
            t.small[i] = (2U * t.small[i - 1U]) % MOD;
        }

        t.cycle_shifted.assign(POW2_CYCLE, 0U);
        u32 value = static_cast<u32>((2ULL * t.small[POW2_THRESHOLD - 1U]) % MOD);  // 2^9 mod MOD
        for (u32 i = 0U; i < POW2_CYCLE; ++i) {
            t.cycle_shifted[i] = value;
            value = static_cast<u32>((2ULL * value) % MOD);
        }
        return t;
    }();
    return tables;
}

u32 pow2_mod_from_exp(const u32 exp_mod, const u8 exp_cap) {
    const Pow2Tables& tables = get_pow2_tables();
    if (exp_cap < POW2_THRESHOLD) {
        return tables.small[exp_cap];
    }

    // For E >= 9:
    // 2^E mod 10^9 = 2^(9 + ((E - 9) mod ord_{5^9}(2))) mod 10^9.
    u32 idx = exp_mod;
    if (idx >= POW2_THRESHOLD) {
        idx -= POW2_THRESHOLD;
    } else {
        idx += POW2_CYCLE - POW2_THRESHOLD;
    }
    return tables.cycle_shifted[idx];
}

std::vector<int> sieve_primes(const int limit) {
    if (limit < 2) {
        return {};
    }

    std::vector<bool> is_prime(static_cast<std::size_t>(limit) + 1ULL, true);
    is_prime[0] = false;
    is_prime[1] = false;

    for (int p = 2; p * p <= limit; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int q = p * p; q <= limit; q += p) {
            is_prime[static_cast<std::size_t>(q)] = false;
        }
    }

    std::vector<int> primes;
    for (int p = 2; p <= limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)]) {
            primes.push_back(p);
        }
    }
    return primes;
}

std::vector<int> exponent_plus_one_list_for_lcm_upto_n(const int n) {
    const std::vector<int> primes = sieve_primes(n);
    std::vector<int> result;
    result.reserve(primes.size());

    for (const int p : primes) {
        int d = 0;
        int x = n;
        while (x > 0) {
            x /= p;
            ++d;
        }
        result.push_back(d);
    }

    return result;
}

std::vector<u32> build_binomial_row_mod(const int n) {
    std::vector<u32> row(static_cast<std::size_t>(n) + 1ULL, 0U);
    row[0] = 1U;

    for (int i = 1; i <= n; ++i) {
        for (int k = i; k >= 1; --k) {
            u32 value = row[static_cast<std::size_t>(k)] +
                        row[static_cast<std::size_t>(k - 1)];
            if (value >= MOD) {
                value -= MOD;
            }
            row[static_cast<std::size_t>(k)] = value;
        }
    }

    return row;
}

GroupData build_group_data(const int d, const int count) {
    GroupData group;
    group.d = d;
    group.count = count;
    group.choices.assign(static_cast<std::size_t>(count) + 1ULL, ChoiceData{});

    const std::vector<u32> binom = build_binomial_row_mod(count);

    std::vector<u32> pow_d_mod(static_cast<std::size_t>(count) + 1ULL, 1U);
    std::vector<u32> pow_dm1_mod(static_cast<std::size_t>(count) + 1ULL, 1U);
    std::vector<u8> pow_d_cap(static_cast<std::size_t>(count) + 1ULL, 1U);
    std::vector<u8> pow_dm1_cap(static_cast<std::size_t>(count) + 1ULL, 1U);

    for (int k = 1; k <= count; ++k) {
        pow_d_mod[static_cast<std::size_t>(k)] = static_cast<u32>(
            (static_cast<u64>(pow_d_mod[static_cast<std::size_t>(k - 1)]) *
             static_cast<u64>(d)) %
            POW2_CYCLE);
        pow_dm1_mod[static_cast<std::size_t>(k)] = static_cast<u32>(
            (static_cast<u64>(pow_dm1_mod[static_cast<std::size_t>(k - 1)]) *
             static_cast<u64>(d - 1)) %
            POW2_CYCLE);

        pow_d_cap[static_cast<std::size_t>(k)] = cap_mul_9(
            pow_d_cap[static_cast<std::size_t>(k - 1)], static_cast<u8>(d));
        pow_dm1_cap[static_cast<std::size_t>(k)] = cap_mul_9(
            pow_dm1_cap[static_cast<std::size_t>(k - 1)], static_cast<u8>(d - 1));
    }

    for (int i = 0; i <= count; ++i) {
        ChoiceData choice;

        const u32 c = binom[static_cast<std::size_t>(i)];
        if ((i & 1) == 0) {
            choice.coef = c;
        } else {
            choice.coef = (c == 0U ? 0U : MOD - c);
        }

        const u32 left = pow_dm1_mod[static_cast<std::size_t>(i)];
        const u32 right = pow_d_mod[static_cast<std::size_t>(count - i)];
        choice.exp_mod = static_cast<u32>(
            (static_cast<u64>(left) * static_cast<u64>(right)) % POW2_CYCLE);

        choice.exp_cap = cap_mul_9(pow_dm1_cap[static_cast<std::size_t>(i)],
                                   pow_d_cap[static_cast<std::size_t>(count - i)]);

        group.choices[static_cast<std::size_t>(i)] = choice;
    }

    return group;
}

std::vector<GroupData> build_groups_for_hl(const int n) {
    const std::vector<int> d_list = exponent_plus_one_list_for_lcm_upto_n(n);

    std::unordered_map<int, int> count_by_d;
    count_by_d.reserve(d_list.size() * 2ULL + 1ULL);
    for (const int d : d_list) {
        ++count_by_d[d];
    }

    std::vector<std::pair<int, int>> sorted_pairs;
    sorted_pairs.reserve(count_by_d.size());
    for (const auto& kv : count_by_d) {
        sorted_pairs.push_back(kv);
    }

    // Sort by multiplicity descending so the largest bucket is parallelized.
    std::sort(sorted_pairs.begin(), sorted_pairs.end(),
              [](const std::pair<int, int>& a, const std::pair<int, int>& b) {
                  if (a.second != b.second) {
                      return a.second > b.second;
                  }
                  return a.first > b.first;
              });

    std::vector<GroupData> groups;
    groups.reserve(sorted_pairs.size());
    for (const auto& [d, count] : sorted_pairs) {
        groups.push_back(build_group_data(d, count));
    }

    return groups;
}

std::vector<RestState> build_rest_states(const std::vector<GroupData>& groups) {
    std::vector<RestState> states;
    states.push_back(RestState{1U, 1U, 1U});

    for (std::size_t g = 1; g < groups.size(); ++g) {
        const GroupData& group = groups[g];
        std::vector<RestState> next;
        next.reserve(states.size() * group.choices.size());

        for (const RestState& state : states) {
            for (const ChoiceData& choice : group.choices) {
                const u32 coef = static_cast<u32>(
                    (static_cast<u64>(state.coef) * static_cast<u64>(choice.coef)) %
                    MOD);
                if (coef == 0U) {
                    continue;
                }

                const u32 exp_mod = static_cast<u32>(
                    (static_cast<u64>(state.exp_mod) *
                     static_cast<u64>(choice.exp_mod)) %
                    POW2_CYCLE);
                const u8 exp_cap = cap_mul_9(state.exp_cap, choice.exp_cap);
                next.push_back(RestState{coef, exp_mod, exp_cap});
            }
        }

        states.swap(next);
    }

    return states;
}

u32 solve_hl_mod(const int n, const unsigned requested_threads) {
    const std::vector<GroupData> groups = build_groups_for_hl(n);
    if (groups.empty()) {
        return 0U;
    }

    const GroupData& lead = groups[0];
    const std::vector<RestState> rest_states = build_rest_states(groups);

    const std::size_t lead_choices = lead.choices.size();
    unsigned thread_count = requested_threads;
    if (thread_count == 0U) {
        thread_count = 1U;
    }
    thread_count = std::min<unsigned>(thread_count, static_cast<unsigned>(lead_choices));
    if (thread_count == 0U) {
        thread_count = 1U;
    }

    std::vector<u64> partial(static_cast<std::size_t>(thread_count), 0ULL);
    std::vector<std::thread> workers;
    workers.reserve(thread_count);

    const std::size_t block = (lead_choices + static_cast<std::size_t>(thread_count) - 1ULL) /
                              static_cast<std::size_t>(thread_count);

    auto worker = [&](const unsigned worker_id) {
        const std::size_t begin = static_cast<std::size_t>(worker_id) * block;
        const std::size_t end = std::min<std::size_t>(lead_choices, begin + block);
        u64 local = 0ULL;

        for (std::size_t i = begin; i < end; ++i) {
            const ChoiceData& prefix = lead.choices[i];
            for (const RestState& state : rest_states) {
                const u32 coef = static_cast<u32>(
                    (static_cast<u64>(prefix.coef) * static_cast<u64>(state.coef)) %
                    MOD);
                if (coef == 0U) {
                    continue;
                }

                const u8 exp_cap = cap_mul_9(prefix.exp_cap, state.exp_cap);
                const u32 exp_mod = static_cast<u32>(
                    (static_cast<u64>(prefix.exp_mod) * static_cast<u64>(state.exp_mod)) %
                    POW2_CYCLE);
                const u32 value_2exp = pow2_mod_from_exp(exp_mod, exp_cap);

                local += (static_cast<u64>(coef) * static_cast<u64>(value_2exp)) % MOD;
            }
        }

        partial[static_cast<std::size_t>(worker_id)] = local % MOD;
    };

    for (unsigned t = 0U; t < thread_count; ++t) {
        workers.emplace_back(worker, t);
    }
    for (std::thread& th : workers) {
        th.join();
    }

    u64 total = 0ULL;
    for (const u64 value : partial) {
        total += value;
    }
    return static_cast<u32>(total % MOD);
}

u64 brute_h(const u64 n) {
    std::vector<u64> divisors;
    for (u64 d = 1ULL; d <= n; ++d) {
        if (n % d == 0ULL) {
            divisors.push_back(d);
        }
    }

    const std::size_t m = divisors.size();
    if (m >= 63U) {
        throw std::runtime_error("brute_h requested for too many divisors");
    }

    const u64 subset_count = (1ULL << m);
    u64 answer = 0ULL;
    for (u64 mask = 1ULL; mask < subset_count; ++mask) {
        u64 current_lcm = 1ULL;
        for (std::size_t i = 0; i < m; ++i) {
            if (((mask >> i) & 1ULL) == 0ULL) {
                continue;
            }
            const u64 x = divisors[i];
            const u64 g = std::gcd(current_lcm, x);
            current_lcm = (current_lcm / g) * x;
        }
        if (current_lcm == n) {
            ++answer;
        }
    }

    return answer;
}

u32 solve_hl_direct_ie_small(const int n) {
    const std::vector<int> d_list = exponent_plus_one_list_for_lcm_upto_n(n);
    const int k = static_cast<int>(d_list.size());
    if (k >= 60) {
        throw std::runtime_error("Direct IE checkpoint is only for small k");
    }

    const u64 subset_count = (1ULL << static_cast<u64>(k));
    u64 answer = 0ULL;

    for (u64 mask = 0ULL; mask < subset_count; ++mask) {
        u64 exponent = 1ULL;
        for (int i = 0; i < k; ++i) {
            const bool in_subset = ((mask >> static_cast<u64>(i)) & 1ULL) != 0ULL;
            const int factor = in_subset ? (d_list[static_cast<std::size_t>(i)] - 1)
                                         : d_list[static_cast<std::size_t>(i)];
            exponent *= static_cast<u64>(factor);
        }

        const u32 term = pow_mod_u32(2U, exponent, MOD);
        if ((__builtin_popcountll(mask) & 1) != 0) {
            answer += (MOD - term);
        } else {
            answer += term;
        }
        answer %= MOD;
    }

    return static_cast<u32>(answer);
}

bool run_checkpoints(const unsigned thread_count) {
    if (brute_h(6ULL) != 10ULL) {
        std::cerr << "Checkpoint failed: H(6) must be 10\n";
        return false;
    }

    const u32 sample = solve_hl_mod(4, 1U);
    if (sample != 44U) {
        std::cerr << "Checkpoint failed: HL(4) must be 44, got " << sample << '\n';
        return false;
    }

    for (int n = 2; n <= 20; ++n) {
        const u32 fast = solve_hl_mod(n, 1U);
        const u32 direct = solve_hl_direct_ie_small(n);
        if (fast != direct) {
            std::cerr << "Checkpoint failed at n=" << n << ": fast=" << fast
                      << " direct=" << direct << '\n';
            return false;
        }
    }

    if (thread_count > 1U) {
        const unsigned check_threads = std::min<unsigned>(thread_count, 4U);
        const u32 single = solve_hl_mod(200, 1U);
        const u32 multi = solve_hl_mod(200, check_threads);
        if (single != multi) {
            std::cerr << "Checkpoint failed: thread consistency mismatch at n=200\n";
            return false;
        }
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.n < 2) {
        std::cerr << "n must be at least 2\n";
        return 1;
    }

    const unsigned thread_count = pick_thread_count(options.requested_threads);

    if (options.run_checkpoints) {
        if (!run_checkpoints(thread_count)) {
            return 1;
        }
    }

    const u32 answer = solve_hl_mod(options.n, thread_count);
    std::cout << answer << '\n';
    return 0;
}

Python

def solve():
    MOD = 10**9; N = 50000
    P2T = 9; P2C = 1562500  # ord_{5^9}(2)

    def sieve(n):
        ip=[True]*(n+1); ip[0]=ip[1]=False
        for i in range(2,int(n**0.5)+1):
            if ip[i]:
                for j in range(i*i,n+1,i): ip[j]=False
        return [i for i in range(2,n+1) if ip[i]]

    def pm(b,e,m):
        r=1; b%=m
        while e>0:
            if e&1: r=r*b%m
            b=b*b%m; e>>=1
        return r

    def cap9(a,b):
        p=a*b; return min(p,P2T)

    # Pre-build pow2 tables
    small=[1]*P2T
    for i in range(1,P2T): small[i]=2*small[i-1]%MOD
    cyc=[0]*P2C; v=2*small[P2T-1]%MOD
    for i in range(P2C): cyc[i]=v; v=2*v%MOD

    def pow2_mod(em, ec):
        if ec<P2T: return small[ec]
        idx=em
        if idx>=P2T: idx-=P2T
        else: idx+=P2C-P2T
        return cyc[idx]

    # d_list: for each prime p <= N, d = floor(log_p(N)) + 1
    primes = sieve(N)
    d_list = []
    for p in primes:
        d=0; x=N
        while x>0: x//=p; d+=1
        d_list.append(d)

    # Group by d value
    from collections import Counter
    cnt = Counter(d_list)
    groups = sorted(cnt.items(), key=lambda x: (-x[1], -x[0]))

    def build_binom(n2):
        row=[0]*(n2+1); row[0]=1
        for i in range(1,n2+1):
            for k in range(i,0,-1): row[k]=(row[k]+row[k-1])%MOD
        return row

    def build_group(d,count):
        binom=build_binom(count)
        pdm=[1]*(count+1); pd1m=[1]*(count+1)
        pdc=[1]*(count+1); pd1c=[1]*(count+1)
        for k in range(1,count+1):
            pdm[k]=pdm[k-1]*d%P2C; pd1m[k]=pd1m[k-1]*(d-1)%P2C
            pdc[k]=cap9(pdc[k-1],d); pd1c[k]=cap9(pd1c[k-1],d-1)
        choices=[]
        for i in range(count+1):
            c=binom[i]; coef=(MOD-c)%MOD if i&1 else c
            em=pd1m[i]*pdm[count-i]%P2C
            ec=cap9(pd1c[i],pdc[count-i])
            choices.append((coef,em,ec))
        return choices

    # Build rest states from groups[1:]
    group_choices = [build_group(d,c) for d,c in groups]
    states = [(1,1,1)]  # (coef, exp_mod, exp_cap)
    for gi in range(1,len(groups)):
        nxt=[]
        for sc,se,sec in states:
            for cc,ce,cec in group_choices[gi]:
                nc=sc*cc%MOD
                if nc==0: continue
                ne=se*ce%P2C; nec=cap9(sec,cec)
                nxt.append((nc,ne,nec))
        states=nxt

    # Combine with lead group
    total=0
    for pc,pe,pec in group_choices[0]:
        for sc,se,sec in states:
            c=pc*sc%MOD
            if c==0: continue
            ec=cap9(pec,sec); em=pe*se%P2C
            v=pow2_mod(em,ec)
            total=(total+c*v)%MOD
    return str(total%MOD)

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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler590 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler590.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler590_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler590 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler590 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler590 C++ bridge produced empty output.");
        }
        return parsed;
    }

    public static void main(String[] args) throws Exception {
        System.out.println(solveViaCppBridge());
    }
}