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

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