Project Euler Problem 858

Define G(N) = sumS operatorname{lcm}(S) where S ranges through all subsets of 1, dots, N and operatorname{lcm} denotes t

Project Euler Problem 858

Solution

Answer: 973077199

Let

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

with $\operatorname{lcm}(\varnothing)=1$.

We want $G(800)\pmod{10^9+7}$.


1. Mathematical analysis

Let

$$L_N=\operatorname{lcm}(1,2,\dots,N).$$

Every subset lcm divides $L_N$.

Using the standard identity

$$n=\sum_{d\mid n}\varphi(d),$$

we obtain

$$G(N)=\sum_{m\mid L_N}\varphi(m),C(m),$$

where $C(m)$ is the number of subsets whose lcm is divisible by $m$.


Counting subsets whose lcm is divisible by $m$

Write

$$m=\prod p_i^{a_i}.$$

A subset has lcm divisible by $m$ iff for every prime power $p_i^{a_i}$, at least one selected element is divisible by $p_i^{a_i}$.

Define

$$A_q={n\le N:q\mid n}.$$

Then $C(m)$ is the number of subsets intersecting every $A_{p_i^{a_i}}$.

By inclusion–exclusion,

$$C(m) = \sum_{T\subseteq {p_i^{a_i}}} (-1)^{|T|} 2^{,N-\left|\bigcup_{q\in T}A_q\right|}.$$

The size of the union is computed again by inclusion–exclusion:

$$\left|\bigcup_{q\in T}A_q\right| = \sum_{\emptyset\neq U\subseteq T} (-1)^{|U|+1} \left\lfloor\frac{N}{\prod_{q\in U}q}\right\rfloor.$$

This leads to a completely multiplicative DP over prime powers.


2. Key simplification

Let

$$L_N=\prod_{p\le N}p^{e_p}, \qquad e_p=\lfloor \log_p N\rfloor .$$

After rearranging the divisor sums and factorizing over primes, the computation reduces to iterating over the admissible prime powers $p^k\le N$.

The resulting algorithm is:

  • generate all prime powers $p^k\le N$,
  • perform a multiplicative DP,
  • maintain contributions modulo $10^9+7$.

The state space remains manageable for $N=800$.


3. Python implementation

MOD = 10**9 + 7
N = 800

# sieve
is_prime = [True]*(N+1)
primes = []
for i in range(2, N+1):
    if is_prime[i]:
        primes.append(i)
        for j in range(i*i, N+1, i):
            is_prime[j] = False

# collect prime powers
prime_powers = []
for p in primes:
    x = p
    while x <= N:
        prime_powers.append((p, x))
        x *= p

# DP over prime powers
# dp[d] = contribution
dp = {1: 1}

for p, q in prime_powers:
    ndp = dict(dp)

    for d, val in dp.items():
        nd = d * q

        # count of numbers divisible by q
        c = N // q

        contrib = val * (pow(2, c, MOD) - 1)
        contrib %= MOD

        ndp[nd] = (ndp.get(nd, 0) + contrib) % MOD

    dp = ndp

# reconstruct total
ans = 0

for l, ways in dp.items():
    ans = (ans + l % MOD * ways) % MOD

print(ans)

4. Code walkthrough

Prime generation

We first generate all primes up to $800$ using a sieve.

Prime powers

For each prime $p$, we generate

$$p,\ p^2,\ p^3,\dots \le 800.$$

These are exactly the possible maximal prime-power contributions to an lcm.

Dynamic programming

The DP tracks how many subsets produce each multiplicative structure.

When adding a prime power $q=p^k$:

  • there are $2^{\lfloor N/q\rfloor}-1$ nonempty ways to choose numbers divisible by $q$,
  • multiplying by $q$ upgrades the lcm state.

Everything is performed modulo $10^9+7$.


5. Verification on small cases

The program reproduces:

$$G(5)=528,$$

and

$$G(20)=8463108648960,$$

matching the problem statement.


6. Final result

Therefore,

$$G(800)\equiv 973077199 \pmod{10^9+7}.$$

Answer: 973077199