Project Euler Problem 944

Solution

Answer: 1228599511

Let $E \subseteq {1,2,\dots,n}$. An element $x\in E$ contributes to $\operatorname{sev}(E)$ iff another element of $E$ is divisible by $x$.

For fixed $x$, let

$$m=\left\lfloor \frac{n}{x}\right\rfloor$$

so there are $m-1$ other multiples of $x$ in ${1,\dots,n}$.

The number of subsets $E$ where $x$ is an elevisor is:

include $x$,
include at least one of the $m-1$ other multiples,
all remaining elements arbitrary.

Hence:

$$N_x = 2^{,n-m}\left(2^{m-1}-1\right) = 2^{n-1}-2^{n-m}$$

$$S(n) = \sum_{x=1}^{n} x\left(2^{n-1}-2^{n-\lfloor n/x\rfloor}\right)$$

Rearranging gives a much faster summation:

$$S(n) = \sum_{d=2}^{n} 2^{n-d}\cdot \frac{\lfloor n/d\rfloor(\lfloor n/d\rfloor+1)}{2}$$

This can be evaluated in $O(\sqrt n)$ using the standard floor-division interval trick:

values of $\lfloor n/d\rfloor$ are constant on intervals.

Clean Python implementation:

import math

MOD = 1234567891
INV2 = (MOD + 1) // 2

def solve(n):
    D = math.isqrt(n)
    res = 0

    # Small d: direct
    p = pow(2, (n - 2) % (MOD - 1), MOD)

    for d in range(2, D + 1):
        q = n // d
        tri = (q * (q + 1) // 2) % MOD
        res += p * tri
        p = (p * INV2) % MOD

    res %= MOD

    # Large d grouped by equal floor(n/d)
    cache = {1: INV2}

    cur = D + 1
    pcur = pow(2, (n - cur) % (MOD - 1), MOD)

    for q in range(D - 1, 0, -1):
        R = n // q
        if R < cur:
            continue

        length = R - cur + 1

        if length not in cache:
            cache[length] = pow(INV2, length, MOD)

        ip = cache[length]

        tri = (q * (q + 1) // 2) % MOD
        interval_sum = (2 * pcur * (1 - ip)) % MOD

        res += interval_sum * tri
        pcur = (pcur * ip) % MOD
        cur = R + 1

    return res % MOD

print(solve(10**14))

Checking the sample:

$$S(10)=4927$$

which matches the problem statement.

For $n=10^{14}$:

Answer: 1228599511