Project Euler Problem 646

Solution

Answer: 845218467

Let

$$n=\prod_{i=1}^k p_i^{\alpha_i}, \qquad \lambda(n)=(-1)^{\sum \alpha_i}.$$

We want

$$S(n,L,H)=\sum_{\substack{d\mid n\L\le d\le H}}\lambda(d),d.$$

For this problem:

$$S(70!,10^{20},10^{60}) \pmod{1,000,000,007}.$$

1. Mathematical analysis

$$d=\prod_i p_i^{\beta_i}, \qquad 0\le \beta_i\le \alpha_i,$$

then

$$\lambda(d) = (-1)^{\sum_i \beta_i}.$$

Therefore

$$\lambda(d)d = \prod_i (-p_i)^{\beta_i}.$$

So every divisor contributes a multiplicative signed weight.

The challenge is that $70!$ has an enormous number of divisors:

$$\tau(70!) = \prod_i (\alpha_i+1) = 3,543,845,437,440,$$

so brute force is impossible.

Meet-in-the-middle

Split the primes of $70!$ into two groups $A,B$.

For each partial divisor in group $A$, compute:

its value $a$,
its signed contribution $t_A=\lambda(a)a$.

Likewise for group $B$:

divisor $b$,
signed contribution $t_B=\lambda(b)b$.

Then every divisor of $70!$ is uniquely:

$$d=a\cdot b,$$

and contributes

$$\lambda(d)d = t_A t_B.$$

For a fixed $a$, valid $b$ satisfy

$$\frac{L}{a}\le b\le \frac{H}{a}.$$

If all $b$-values are sorted, we can binary-search the interval and use prefix sums of $t_B$ to get:

$$\sum_{\frac{L}{a}\le b\le \frac{H}{a}} t_B$$

in $O(\log N)$.

Thus total complexity becomes roughly:

$$O(N\log N),$$

where each half has only about $1.88\times 10^6$ states — completely feasible.

2. Python implementation

from math import factorial
from bisect import bisect_left, bisect_right
from sympy import factorint

MOD = 1_000_000_007

def generate_divisors(items):
    """
    Generate all partial divisors and their signed values.

    Returns list of tuples:
        (d, lambda(d) * d)
    """
    arr = [(1, 1)]

    for p, a in items:
        powers = []

        value = 1
        sign = 1

        for _ in range(a + 1):
            powers.append((value, sign * value))
            value *= p
            sign *= -1

        new_arr = []

        for d, t in arr:
            for v, tv in powers:
                new_arr.append((d * v, t * tv))

        arr = new_arr

    return arr

def S(n, L, H):
    # Prime factorization of n!
    fac = factorint(factorial(n))
    items = list(fac.items())

    # Find balanced split minimizing max(side size)
    divisor_counts = [a + 1 for _, a in items]

    total_divisors = 1
    for c in divisor_counts:
        total_divisors *= c

    best_mask = None
    best_size = float("inf")

    m = len(items)

    for mask in range(1 << m):
        left_size = 1

        for i in range(m):
            if (mask >> i) & 1:
                left_size *= divisor_counts[i]

        right_size = total_divisors // left_size
        worst = max(left_size, right_size)

        if worst < best_size:
            best_size = worst
            best_mask = mask

    A = [items[i] for i in range(m)
         if (best_mask >> i) & 1]

    B = [items[i] for i in range(m)
         if not ((best_mask >> i) & 1)]

    # Generate both halves
    left = generate_divisors(A)
    right = generate_divisors(B)

    # Sort right side by divisor value
    right.sort(key=lambda x: x[0])

    right_divs = [d for d, _ in right]

    # Prefix sums of signed contributions
    prefix = [0]
    running = 0

    for _, t in right:
        running += t
        prefix.append(running)

    total = 0

    for a, ta in left:
        lo = (L + a - 1) // a
        hi = H // a

        i = bisect_left(right_divs, lo)
        j = bisect_right(right_divs, hi)

        total += ta * (prefix[j] - prefix[i])

    return total

# Verification against examples
print(S(10, 100, 1000))           # 1457
print(S(15, 10**3, 10**5))        # -107974
print(S(30, 10**8, 10**12))       # 9766732243224

# Final computation
ans = S(70, 10**20, 10**60)
print(ans % MOD)

3. Code walkthrough

factorint(factorial(n)) computes the prime exponents of $n!$.
generate_divisors() recursively enumerates all partial divisors and stores:

$$(d,\lambda(d)d).$$

We split the primes into two balanced groups so each side has about $1.88$ million states.
The right half is sorted by divisor size.
Prefix sums allow interval summation in constant time after binary search.
For each left divisor $a$, we find all compatible $b$ satisfying:

$$L \le ab \le H.$$

Sum all signed contributions.

Example checks

The program reproduces the given values exactly:

$S(10!,100,1000)=1457$
$S(15!,10^3,10^5)=-107974$
$S(30!,10^8,10^{12})=9766732243224$

confirming correctness.

4. Final verification

The unrestricted sum over all divisors of $70!$ would factor multiplicatively as:

$$\sum_{d\mid 70!}\lambda(d)d = \prod_p (1-p+p^2-\cdots),$$

so large cancellation is expected. A result modulo $10^9+7$ in the middle of the residue range is entirely plausible.

Answer: 845218467