euler/969.py

#!/usr/bin/env python3

# Generated by GPT-5 on 2025-11-05. Tested on Ubuntu 25.04, pypy3 3.11.11. Runtime on Lenovo ThinkPad P14s: 0m0.077s
# Output: 412543690

# Computes sum_{n=1..N} S(n) for Project Euler 969 formulation, modulo 1e9+7.
# N is set to 10**18 per problem statement.
# This script is self-contained and uses Python's big integers.

MOD = 10**9 + 7

def primes_upto(n):
    sieve = bytearray(b"\x01") * (n+1)
    primes = []
    for i in range(2, n+1):
        if sieve[i]:
            primes.append(i)
            step = i
            start = i*i
            if start > n:
                continue
            sieve[start:n+1:step] = b"\x00" * (((n - start)//step) + 1)
    return primes

def sum_pows_lagrange(T, k, MOD):
    # returns sum_{t=1}^T t^k mod MOD using Lagrange interpolation on points 0..d (d=k+1)
    d = k + 1
    if T <= d:
        s = 0
        for i in range(1, T+1):
            s = (s + pow(i, k, MOD)) % MOD
        return s

    # y_i = sum_{t=1}^i t^k for i=0..d
    ys = [0] * (d + 1)
    for i in range(1, d + 1):
        ys[i] = (ys[i-1] + pow(i, k, MOD)) % MOD

    fact = [1] * (d + 1)
    for i in range(1, d + 1):
        fact[i] = fact[i-1] * i % MOD
    invfact = [1] * (d + 1)
    invfact[d] = pow(fact[d], MOD-2, MOD)
    for i in range(d, 0, -1):
        invfact[i-1] = invfact[i] * i % MOD

    Tmod = T % MOD

    # pre[i] = prod_{j=0..i-1} (T - j)
    pre = [1] * (d + 1)
    for i in range(1, d + 1):
        pre[i] = pre[i-1] * (Tmod - (i-1)) % MOD

    # suf[i] = prod_{j=i+1..d} (T - j)  => we build suf with one extra slot
    suf = [1] * (d + 2)
    for i in range(d, -1, -1):
        suf[i] = suf[i+1] * (Tmod - i) % MOD

    res = 0
    # Lagrange basis for x_i = i, i=0..d
    for i in range(0, d + 1):
        # numerator = pre[i] * suf[i+1] = prod_{j != i} (T - j)
        num = pre[i] * suf[i+1] % MOD
        # denom = ∏_{j != i} (i - j) = (-1)^{d-i} * i! * (d-i)!
        denom = fact[i] * fact[d - i] % MOD
        if ((d - i) & 1) == 1:
            denom = (-denom) % MOD
        invden = pow(denom, MOD-2, MOD)
        li = num * invden % MOD
        res = (res + ys[i] * li) % MOD

    return res

def compute_Mk(k, primes):
    # compute M_k = ∏_{p <= k} p^{ceil(v_p(k!)/k)}
    if k == 0:
        return 1
    M = 1
    for p in primes:
        if p > k:
            break
        v = 0
        pp = p
        # v_p(k!)
        while pp <= k:
            v += k // pp
            pp *= p
        e = (v + k - 1) // k  # ceil(v/k)
        if e > 0:
            M *= pow(p, e)
    return M

def compute_sum_S_upto_N(N):
    primes = primes_upto(1000)  # plenty for our k range
    MOD = 10**9 + 7

    # precompute factorials (for modular inverses) up to a safe bound
    MAXF = 200
    fact = [1] * (MAXF + 1)
    for i in range(1, MAXF + 1):
        fact[i] = fact[i-1] * i % MOD

    total = 0
    k = 0
    while True:
        M = compute_Mk(k, primes) if k > 0 else 1
        if M > N:
            break
        if k == 0:
            T = N  # m runs over 1..N (since M_0 = 1)
        else:
            if N < k:
                break
            T = (N - k) // M
            if T <= 0:
                k += 1
                continue

        # C_k = (-1)^k * M^k / k!  (mod MOD)
        if k == 0:
            C = 1
        else:
            C = pow(M % MOD, k, MOD) * pow(fact[k], MOD-2, MOD) % MOD
            if k & 1:
                C = (-C) % MOD

        sum_tk = sum_pows_lagrange(T, k, MOD)
        contrib = C * sum_tk % MOD
        total = (total + contrib) % MOD

        k += 1

    return total

if __name__ == "__main__":
    N = 10**18
    print(compute_sum_S_upto_N(N) % MOD)