Problem 858: LCM

Project Euler Problem 858 Solution

EulerSolve provides an optimized solution for Project Euler Problem 858, LCM, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every subset $S\subseteq\{1,\dots,N\}$, including the empty set, we take its least common multiple and sum all those values: $$G(N)=\sum_{S\subseteq\{1,\dots,N\}}\operatorname{lcm}(S),\qquad \operatorname{lcm}(\varnothing)=1.$$ The goal is to evaluate $G(800)\bmod(10^9+7)$. Since there are $2^{800}$ subsets, the solution cannot enumerate subsets directly. Instead it groups subsets by their exact LCM and rewrites the problem as a weighted sum over divisors. Mathematical Approach Let $$L_N=\operatorname{lcm}(1,2,\dots,N)=\prod_{p\le N}p^{M_p},\qquad M_p=\left\lfloor\log_p N\right\rfloor.$$ Every subset LCM is a divisor of $L_N$, so the whole problem lives on the divisor lattice of $L_N$. Step 1: Count subsets whose LCM divides a fixed divisor For a divisor $d\mid L_N$, define $$A(d)=\#\{m\in\{1,\dots,N\}:m\mid d\}.$$ These are exactly the numbers that may appear in a subset whose LCM divides $d$. Therefore, if $$F(d)=\#\{S\subseteq\{1,\dots,N\}:\operatorname{lcm}(S)\mid d\},$$ then every subset of those $A(d)$ admissible numbers is valid, so $$F(d)=2^{A(d)}.$$ This converts subset enumeration into divisor counting....

Detailed mathematical approach

Problem Summary

For every subset $S\subseteq\{1,\dots,N\}$, including the empty set, we take its least common multiple and sum all those values:

$$G(N)=\sum_{S\subseteq\{1,\dots,N\}}\operatorname{lcm}(S),\qquad \operatorname{lcm}(\varnothing)=1.$$

The goal is to evaluate $G(800)\bmod(10^9+7)$. Since there are $2^{800}$ subsets, the solution cannot enumerate subsets directly. Instead it groups subsets by their exact LCM and rewrites the problem as a weighted sum over divisors.

Mathematical Approach

Let

$$L_N=\operatorname{lcm}(1,2,\dots,N)=\prod_{p\le N}p^{M_p},\qquad M_p=\left\lfloor\log_p N\right\rfloor.$$

Every subset LCM is a divisor of $L_N$, so the whole problem lives on the divisor lattice of $L_N$.

Step 1: Count subsets whose LCM divides a fixed divisor

For a divisor $d\mid L_N$, define

$$A(d)=\#\{m\in\{1,\dots,N\}:m\mid d\}.$$

These are exactly the numbers that may appear in a subset whose LCM divides $d$. Therefore, if

$$F(d)=\#\{S\subseteq\{1,\dots,N\}:\operatorname{lcm}(S)\mid d\},$$

then every subset of those $A(d)$ admissible numbers is valid, so

$$F(d)=2^{A(d)}.$$

This converts subset enumeration into divisor counting.

Step 2: Recover the exact LCM by Möbius inversion

Let

$$T(d)=\#\{S\subseteq\{1,\dots,N\}:\operatorname{lcm}(S)=d\}.$$

Since “the LCM divides $d$” means the exact LCM is some divisor of $d$, we have

$$F(d)=\sum_{c\mid d}T(c).$$

Möbius inversion on divisors gives

$$T(d)=\sum_{r\mid d}\mu(r)\,F\left(\frac{d}{r}\right).$$

Now substitute this into

$$G(N)=\sum_{d\mid L_N} d\,T(d).$$

After reindexing the sum, we obtain

$$G(N)=\sum_{d\mid L_N} W(d)\,2^{A(d)},$$

where the weight is

$$W(d)=d\sum_{r\mid L_N/d} r\,\mu(r).$$

Step 3: Factor the weight prime by prime

The Möbius function is zero on non-squarefree numbers, so only the choices $r=1$ and $r=p$ matter for each prime that is still missing from $d$ at full exponent. Hence

$$W(d)=d\prod_{p\mid L_N/d}(1-p).$$

If $v_p(d)=e$, the local factor is therefore

$$w_p(e)=\begin{cases} p^e(1-p), & e<M_p,\\ p^{M_p}, & e=M_p. \end{cases}$$

So the global weight factors multiplicatively as

$$W(d)=\prod_{p\le N} w_p\bigl(v_p(d)\bigr).$$

This is the mathematical reason the implementation multiplies by the prime power $p^e$ and appends the factor $(1-p)$ whenever the exponent is not maximal.

Step 4: Split the divisor into small and large primes

Write any divisor $d\mid L_N$ as

$$d=a\,b,$$

where $a$ uses only primes $q\le\sqrt N$ and $b$ uses only primes $p>\sqrt N$.

For a large prime $p>\sqrt N$, we have $p^2>N$, so its exponent in $L_N$ is only $1$. Thus the large-prime part $b$ is squarefree.

Define

$$B_a(M)=\#\{n\le M:n\mid a,\ \text{and every prime factor of }n\text{ is }\le\sqrt N\}.$$

Any number $m\le N$ dividing $d$ contains either no large prime or exactly one large prime. It cannot contain two distinct large primes, because their product would exceed $N$. Therefore

$$A(a b)=B_a(N)+\sum_{p\mid b} B_a\left(\left\lfloor\frac{N}{p}\right\rfloor\right).$$

This identity is the key simplification used by all three implementations.

Step 5: Sum all choices of large primes at once

Fix the small-prime part $a$. Then

$$2^{A(a b)}=2^{B_a(N)}\prod_{p\mid b}2^{B_a(\lfloor N/p\rfloor)}.$$

Combining this with the local weight factors shows that the sum over every possible large-prime subset factorizes independently:

$$\sum_{b} W(a b)\,2^{A(a b)} = W_{\mathrm{small}}(a)\,2^{B_a(N)} \prod_{p>\sqrt N}\left((1-p)+p\,2^{B_a(\lfloor N/p\rfloor)}\right),$$

where $W_{\mathrm{small}}(a)$ is the product of the local weights coming from the small primes only.

Many large primes share the same quotient

$$M=\left\lfloor\frac{N}{p}\right\rfloor.$$

Grouping primes by this common value lets the implementation precompute one product table per group.

Step 6: Encode the small-prime part as mixed-radix states

If the small primes are $q_1,\dots,q_k$, then $a$ is determined by its exponent vector

$$\bigl(e_1,\dots,e_k\bigr),\qquad 0\le e_i\le M_{q_i}.$$

The number of states is

$$\prod_{i=1}^k (M_{q_i}+1).$$

For each needed limit $M$, the implementation enumerates all small-only integers $n\le M$, records their exact exponent vectors, and applies multidimensional prefix sums. After this transform, each state stores $B_a(M)$: the number of small-only divisors of the corresponding $a$ that do not exceed $M$.

Once these tables are built, every state can be evaluated independently and all contributions are accumulated modulo $10^9+7$.

Worked Example: $N=5$

Here

$$L_5=\operatorname{lcm}(1,2,3,4,5)=60=2^2\cdot 3\cdot 5.$$

The only small prime is $2$, so the small-prime states are $a\in\{1,2,4\}$. The large primes are $3$ and $5$, and both satisfy

$$\left\lfloor\frac{5}{3}\right\rfloor=\left\lfloor\frac{5}{5}\right\rfloor=1.$$

The small-only divisor counts are

$$B_1(5)=1,\qquad B_2(5)=2,\qquad B_4(5)=3,$$

because the admissible small-only divisors up to $5$ are respectively $\{1\}$, $\{1,2\}$, and $\{1,2,4\}$. Also

$$B_1(1)=B_2(1)=B_4(1)=1.$$

The small-prime weights are

$$W_{\mathrm{small}}(1)=1\cdot(1-2)=-1,\qquad W_{\mathrm{small}}(2)=2\cdot(1-2)=-2,\qquad W_{\mathrm{small}}(4)=4.$$

The common large-prime factor is

$$\bigl((1-3)+3\cdot 2^1\bigr)\bigl((1-5)+5\cdot 2^1\bigr)=4\cdot 6=24.$$

So the three state contributions are

$$\begin{aligned} a=1&:&&-1\cdot 2^{1}\cdot 24=-48,\\ a=2&:&&-2\cdot 2^{2}\cdot 24=-192,\\ a=4&:&&4\cdot 2^{3}\cdot 24=768. \end{aligned}$$

Adding them gives

$$-48-192+768=528,$$

which matches the known checkpoint for $G(5)$.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They first generate all primes up to $N$ and split them at $\lfloor\sqrt N\rfloor$. For each small prime they record the maximum exponent occurring in $L_N$, which determines the mixed-radix state space.

Next, they enumerate every small-only integer up to $N$, attach it to its exponent state, and sort those values. For each needed limit $M\in\{N\}\cup\{\lfloor N/p\rfloor:p>\sqrt N\}$, they mark the eligible small-only integers, run a multidimensional prefix accumulation, and obtain the table $B_a(M)$ for every small state $a$.

They also precompute the powers $2^t\bmod(10^9+7)$ for $0\le t\le N$, together with one table for each group of large primes sharing the same quotient $\lfloor N/p\rfloor$. Finally, they iterate over all small-prime states, reconstruct the corresponding small weight, read the precomputed divisor counts, multiply the grouped large-prime factors, and accumulate the answer modulo $10^9+7$.

Complexity Analysis

Let

$$S=\prod_{q\le\sqrt N}(M_q+1)$$

be the number of small-prime states, and let $K$ be the number of distinct values of $\lfloor N/p\rfloor$ among primes $p>\sqrt N$. Prime generation costs $O(N\log\log N)$. After that, the main preprocessing and evaluation steps are near $O(S\cdot K)$, with only a small multiplicative factor from the number of small primes and the prefix transforms. Memory usage is $O(S\cdot K)$ for the stored count tables and grouped large-prime tables. For $N=800$, these sizes remain comfortably manageable.

Footnotes and References

Problem page: https://projecteuler.net/problem=858
Möbius inversion formula: Wikipedia — Möbius inversion formula
Möbius function: Wikipedia — Möbius function
Least common multiple: Wikipedia — Least common multiple
Prime factorization: Wikipedia — Fundamental theorem of arithmetic

Problem 858 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <utility>
#include <vector>

using i64 = std::int64_t;
using u64 = std::uint64_t;
using u16 = std::uint16_t;
using u128 = unsigned __int128;

static constexpr i64 MOD = 1'000'000'007LL;

static i64 mod_norm(i64 x) {
    x %= MOD;
    if (x < 0) x += MOD;
    return x;
}

static i64 mod_mul(i64 a, i64 b) {
    return static_cast<i64>((static_cast<u128>(mod_norm(a)) * static_cast<u128>(mod_norm(b))) % MOD);
}

static std::vector<int> primes_up_to(int n) {
    std::vector<bool> is_prime(n + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int p = 2; p * p <= n; ++p) {
        if (!is_prime[p]) continue;
        for (int q = p * p; q <= n; q += p) is_prime[q] = false;
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) if (is_prime[i]) primes.push_back(i);
    return primes;
}

struct Engine {
    int N;
    std::vector<int> small_primes;
    std::vector<int> large_primes;
    std::vector<int> max_exp;
    std::vector<int> radix;
    std::vector<int> stride;
    int states = 1;
    std::vector<std::vector<i64>> p_pow;
    std::vector<std::pair<int, int>> small_numbers;

    explicit Engine(int n) : N(n) {
        auto primes = primes_up_to(N);
        int root = 1;
        while ((root + 1) * (root + 1) <= N) ++root;

        for (int p : primes) {
            if (p <= root) small_primes.push_back(p);
            else large_primes.push_back(p);
        }

        int k = static_cast<int>(small_primes.size());
        max_exp.assign(k, 0);
        radix.assign(k, 0);

        for (int i = 0; i < k; ++i) {
            int p = small_primes[i];
            int e = 0;
            i64 v = 1;
            while (v * p <= N) {
                v *= p;
                ++e;
            }
            max_exp[i] = e;
            radix[i] = e + 1;
            states *= radix[i];
        }

        stride.assign(k, 1);
        for (int i = k - 2; i >= 0; --i) stride[i] = stride[i + 1] * radix[i + 1];

        p_pow.resize(k);
        for (int i = 0; i < k; ++i) {
            int emax = max_exp[i];
            p_pow[i].assign(emax + 1, 1);
            for (int e = 1; e <= emax; ++e) p_pow[i][e] = mod_mul(p_pow[i][e - 1], small_primes[i]);
        }

        std::vector<int> exps(k, 0);
        gen_small_numbers(0, 1, 0, exps);
        std::sort(small_numbers.begin(), small_numbers.end());
    }

    void gen_small_numbers(int i, i64 value, int idx, std::vector<int>& exps) {
        if (i == static_cast<int>(small_primes.size())) {
            small_numbers.emplace_back(static_cast<int>(value), idx);
            return;
        }

        i64 v = value;
        int p = small_primes[i];
        int base_idx = idx * radix[i];
        for (int e = 0; e <= max_exp[i]; ++e) {
            exps[i] = e;
            gen_small_numbers(i + 1, v, base_idx + e, exps);
            if (e == max_exp[i]) break;
            if (v > N / p) break;
            v *= p;
        }
    }

    void prefix_transform(std::vector<int>& arr) const {
        int k = static_cast<int>(small_primes.size());
        for (int d = 0; d < k; ++d) {
            int st = stride[d];
            int block = radix[d] * st;
            for (int base = 0; base < states; base += block) {
                for (int off = 0; off < st; ++off) {
                    int idx = base + off;
                    for (int e = 1; e < radix[d]; ++e) {
                        idx += st;
                        arr[idx] += arr[idx - st];
                    }
                }
            }
        }
    }

    std::vector<u16> counts_for_limit(int limit) const {
        std::vector<int> arr(states, 0);
        for (const auto& [val, idx] : small_numbers) {
            if (val > limit) break;
            arr[idx] = 1;
        }
        prefix_transform(arr);

        std::vector<u16> out(states);
        for (int i = 0; i < states; ++i) out[i] = static_cast<u16>(arr[i]);
        return out;
    }

    static int floor_div(int a, int b) {
        return a / b;
    }

    i64 solve_mod() {
        std::vector<int> needed_limits;
        needed_limits.push_back(N);

        std::unordered_map<int, std::vector<int>> groups;
        groups.reserve(64);
        for (int p : large_primes) {
            int L = floor_div(N, p);
            groups[L].push_back(p);
        }
        for (auto& kv : groups) needed_limits.push_back(kv.first);

        std::sort(needed_limits.begin(), needed_limits.end());
        needed_limits.erase(std::unique(needed_limits.begin(), needed_limits.end()), needed_limits.end());

        std::unordered_map<int, std::vector<u16>> counts;
        counts.reserve(needed_limits.size() * 2);
        for (int L : needed_limits) counts[L] = counts_for_limit(L);

        std::vector<i64> pow2(N + 1, 1);
        for (int i = 1; i <= N; ++i) pow2[i] = (pow2[i - 1] * 2) % MOD;

        std::unordered_map<int, std::vector<i64>> group_value;
        group_value.reserve(groups.size() * 2);
        for (const auto& kv : groups) {
            int L = kv.first;
            const auto& primes = kv.second;
            std::vector<i64> table(N + 1, 1);
            for (int t = 0; t <= N; ++t) {
                i64 v = 1;
                for (int p : primes) {
                    i64 term = mod_norm(1 - p + mod_mul(p, pow2[t]));
                    v = mod_mul(v, term);
                }
                table[t] = v;
            }
            group_value[L] = std::move(table);
        }

        i64 ans = 0;
        int k = static_cast<int>(small_primes.size());

        for (int idx = 0; idx < states; ++idx) {
            int x = idx;
            i64 d = 1;
            i64 cut = 1;

            for (int i = k - 1; i >= 0; --i) {
                int e = x % radix[i];
                x /= radix[i];

                d = mod_mul(d, p_pow[i][e]);
                if (e < max_exp[i]) cut = mod_mul(cut, mod_norm(1 - small_primes[i]));
            }

            i64 term = mod_mul(d, cut);
            int t0 = counts[N][idx];
            term = mod_mul(term, pow2[t0]);

            for (const auto& kv : group_value) {
                int L = kv.first;
                int t = counts[L][idx];
                term = mod_mul(term, kv.second[t]);
            }

            ans += term;
            if (ans >= MOD) ans -= MOD;
        }

        return ans;
    }
};

static i64 G_mod(int N) {
    Engine e(N);
    return e.solve_mod();
}

int main() {
    assert(G_mod(5) == 528);
    assert(G_mod(20) == 108589719);
    std::cout << G_mod(800) << '\n';
    return 0;
}

Python

import math

kMod = 1000000007

def mod_norm(x):
    x %= kMod
    if x < 0:
        x += kMod
    return x

def mod_mul(a, b):
    return (mod_norm(a) * mod_norm(b)) % kMod

def primes_up_to(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False
    for p in range(2, int(math.isqrt(n)) + 1):
        if is_prime[p]:
            for q in range(p * p, n + 1, p):
                is_prime[q] = False
    return [i for i, prime in enumerate(is_prime) if prime]

class Engine:
    def __init__(self, n):
        self.N = n
        primes = primes_up_to(n)
        root = 1
        while (root + 1) ** 2 <= n:
            root += 1

        self.small_primes = [p for p in primes if p <= root]
        self.large_primes = [p for p in primes if p > root]

        k = len(self.small_primes)
        self.max_exp = [0] * k
        self.radix = [0] * k
        self.states = 1

        for i in range(k):
            p = self.small_primes[i]
            e = 0
            v = 1
            while v * p <= self.N:
                v *= p
                e += 1
            self.max_exp[i] = e
            self.radix[i] = e + 1
            self.states *= self.radix[i]

        self.stride = [1] * k
        for i in range(k - 2, -1, -1):
            self.stride[i] = self.stride[i + 1] * self.radix[i + 1]

        self.p_pow = []
        for i in range(k):
            emax = self.max_exp[i]
            pw = [1] * (emax + 1)
            for e in range(1, emax + 1):
                pw[e] = mod_mul(pw[e - 1], self.small_primes[i])
            self.p_pow.append(pw)

        self.small_numbers = []
        exps = [0] * k
        self.gen_small_numbers(0, 1, 0, exps)
        self.small_numbers.sort()

    def gen_small_numbers(self, i, value, idx, exps):
        if i == len(self.small_primes):
            self.small_numbers.append((value, idx))
            return

        v = value
        p = self.small_primes[i]
        base_idx = idx * self.radix[i]
        for e in range(self.max_exp[i] + 1):
            exps[i] = e
            self.gen_small_numbers(i + 1, v, base_idx + e, exps)
            if e == self.max_exp[i]:
                break
            if v > self.N // p:
                break
            v *= p

    def prefix_transform(self, arr):
        k = len(self.small_primes)
        for d in range(k):
            st = self.stride[d]
            block = self.radix[d] * st
            for base in range(0, self.states, block):
                for off in range(st):
                    idx = base + off
                    for e in range(1, self.radix[d]):
                        idx += st
                        arr[idx] += arr[idx - st]

    def counts_for_limit(self, limit):
        arr = [0] * self.states
        for val, idx in self.small_numbers:
            if val > limit:
                break
            arr[idx] = 1
        
        self.prefix_transform(arr)
        return arr

    def solve_mod(self):
        needed_limits = [self.N]
        groups = {}
        for p in self.large_primes:
            L = self.N // p
            if L not in groups:
                groups[L] = []
            groups[L].append(p)
            
        for L in groups.keys():
            needed_limits.append(L)

        needed_limits = sorted(list(set(needed_limits)))

        counts = {}
        for L in needed_limits:
            counts[L] = self.counts_for_limit(L)

        pow2 = [1] * (self.N + 1)
        for i in range(1, self.N + 1):
            pow2[i] = (pow2[i - 1] * 2) % kMod

        group_value = {}
        for L, primes in groups.items():
            table = [1] * (self.N + 1)
            for t in range(self.N + 1):
                v = 1
                for p in primes:
                    term = mod_norm(1 - p + mod_mul(p, pow2[t]))
                    v = mod_mul(v, term)
                table[t] = v
            group_value[L] = table

        ans = 0
        k = len(self.small_primes)

        for idx in range(self.states):
            x = idx
            d = 1
            cut = 1

            for i in range(k - 1, -1, -1):
                e = x % self.radix[i]
                x //= self.radix[i]

                d = mod_mul(d, self.p_pow[i][e])
                if e < self.max_exp[i]:
                    cut = mod_mul(cut, mod_norm(1 - self.small_primes[i]))

            term = mod_mul(d, cut)
            t0 = counts[self.N][idx]
            term = mod_mul(term, pow2[t0])

            for L, table in group_value.items():
                t = counts[L][idx]
                term = mod_mul(term, table[t])

            ans += term
            if ans >= kMod:
                ans -= kMod

        return ans

def solve():
    engine = Engine(800)
    ans = engine.solve_mod()
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;

public class Euler858 {
    static final long MOD = 1000000007L;

    static long modNorm(long x) {
        x %= MOD;
        if (x < 0)
            x += MOD;
        return x;
    }

    static long modMul(long a, long b) {
        return (modNorm(a) * modNorm(b)) % MOD;
    }

    static ArrayList<Integer> primesUpTo(int n) {
        boolean[] isPrime = new boolean[n + 1];
        java.util.Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;
        for (int p = 2; p * p <= n; ++p) {
            if (!isPrime[p])
                continue;
            for (int q = p * p; q <= n; q += p) {
                isPrime[q] = false;
            }
        }
        ArrayList<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; ++i) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    static class Pair implements Comparable<Pair> {
        int val;
        int idx;

        Pair(int val, int idx) {
            this.val = val;
            this.idx = idx;
        }

        @Override
        public int compareTo(Pair o) {
            return Integer.compare(this.val, o.val);
        }
    }

    static class Engine {
        int N;
        ArrayList<Integer> smallPrimes = new ArrayList<>();
        ArrayList<Integer> largePrimes = new ArrayList<>();
        int[] maxExp;
        int[] radix;
        int[] stride;
        int states = 1;
        long[][] pPow;
        ArrayList<Pair> smallNumbers = new ArrayList<>();

        Engine(int n) {
            this.N = n;
            ArrayList<Integer> primes = primesUpTo(n);
            int root = 1;
            while ((root + 1) * (root + 1) <= n)
                root++;

            for (int p : primes) {
                if (p <= root)
                    smallPrimes.add(p);
                else
                    largePrimes.add(p);
            }

            int k = smallPrimes.size();
            maxExp = new int[k];
            radix = new int[k];

            for (int i = 0; i < k; ++i) {
                int p = smallPrimes.get(i);
                int e = 0;
                long v = 1;
                while (v * p <= N) {
                    v *= p;
                    e++;
                }
                maxExp[i] = e;
                radix[i] = e + 1;
                states *= radix[i];
            }

            stride = new int[k];
            if (k > 0)
                stride[k - 1] = 1;
            for (int i = k - 2; i >= 0; --i) {
                stride[i] = stride[i + 1] * radix[i + 1];
            }

            pPow = new long[k][];
            for (int i = 0; i < k; ++i) {
                int emax = maxExp[i];
                pPow[i] = new long[emax + 1];
                pPow[i][0] = 1;
                for (int e = 1; e <= emax; ++e) {
                    pPow[i][e] = modMul(pPow[i][e - 1], smallPrimes.get(i));
                }
            }

            genSmallNumbers(0, 1, 0);
            Collections.sort(smallNumbers);
        }

        void genSmallNumbers(int i, long value, int idx) {
            if (i == smallPrimes.size()) {
                smallNumbers.add(new Pair((int) value, idx));
                return;
            }

            long v = value;
            int p = smallPrimes.get(i);
            int baseIdx = idx * radix[i];
            for (int e = 0; e <= maxExp[i]; ++e) {
                genSmallNumbers(i + 1, v, baseIdx + e);
                if (e == maxExp[i])
                    break;
                if (v > N / p)
                    break;
                v *= p;
            }
        }

        void prefixTransform(int[] arr) {
            int k = smallPrimes.size();
            for (int d = 0; d < k; ++d) {
                int st = stride[d];
                int block = radix[d] * st;
                for (int base = 0; base < states; base += block) {
                    for (int off = 0; off < st; ++off) {
                        int idx = base + off;
                        for (int e = 1; e < radix[d]; ++e) {
                            idx += st;
                            arr[idx] += arr[idx - st];
                        }
                    }
                }
            }
        }

        short[] countsForLimit(int limit) {
            int[] arr = new int[states];
            for (Pair p : smallNumbers) {
                if (p.val > limit)
                    break;
                arr[p.idx] = 1;
            }
            prefixTransform(arr);

            short[] out = new short[states];
            for (int i = 0; i < states; ++i)
                out[i] = (short) arr[i];
            return out;
        }

        long solveMod() {
            ArrayList<Integer> neededLimitsList = new ArrayList<>();
            neededLimitsList.add(N);

            HashMap<Integer, ArrayList<Integer>> groups = new HashMap<>();
            for (int p : largePrimes) {
                int l = N / p;
                groups.computeIfAbsent(l, k -> new ArrayList<>()).add(p);
            }
            neededLimitsList.addAll(groups.keySet());

            Collections.sort(neededLimitsList);
            ArrayList<Integer> neededLimits = new ArrayList<>();
            if (!neededLimitsList.isEmpty()) {
                neededLimits.add(neededLimitsList.get(0));
                for (int i = 1; i < neededLimitsList.size(); ++i) {
                    if (!neededLimitsList.get(i).equals(neededLimitsList.get(i - 1))) {
                        neededLimits.add(neededLimitsList.get(i));
                    }
                }
            }

            HashMap<Integer, short[]> counts = new HashMap<>();
            for (int l : neededLimits) {
                counts.put(l, countsForLimit(l));
            }

            long[] pow2 = new long[N + 1];
            pow2[0] = 1;
            for (int i = 1; i <= N; ++i)
                pow2[i] = (pow2[i - 1] * 2) % MOD;

            HashMap<Integer, long[]> groupValue = new HashMap<>();
            for (java.util.Map.Entry<Integer, ArrayList<Integer>> entry : groups.entrySet()) {
                int l = entry.getKey();
                ArrayList<Integer> primes = entry.getValue();
                long[] table = new long[N + 1];
                for (int t = 0; t <= N; ++t) {
                    long v = 1;
                    for (int p : primes) {
                        long term = modNorm(1 - p + modMul(p, pow2[t]));
                        v = modMul(v, term);
                    }
                    table[t] = v;
                }
                groupValue.put(l, table);
            }

            long ans = 0;
            int k = smallPrimes.size();

            for (int idx = 0; idx < states; ++idx) {
                int x = idx;
                long d = 1;
                long cut = 1;

                for (int i = k - 1; i >= 0; --i) {
                    int e = x % radix[i];
                    x /= radix[i];

                    d = modMul(d, pPow[i][e]);
                    if (e < maxExp[i]) {
                        cut = modMul(cut, modNorm(1 - smallPrimes.get(i)));
                    }
                }

                long term = modMul(d, cut);
                int t0 = counts.get(N)[idx];
                term = modMul(term, pow2[t0]);

                for (java.util.Map.Entry<Integer, long[]> entry : groupValue.entrySet()) {
                    int l = entry.getKey();
                    long[] table = entry.getValue();
                    int t = counts.get(l)[idx];
                    term = modMul(term, table[t]);
                }

                ans += term;
                if (ans >= MOD)
                    ans -= MOD;
            }

            return ans;
        }
    }

    public static String solve() {
        Engine e = new Engine(800);
        return Long.toString(e.solveMod());
    }

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