Project Euler Problem 956

Solution

Answer: 882086212

Let

$$F(n,m)=D(n,m)=\sum_{\substack{d\mid n\ m\mid \Omega(d)}} d.$$

For

$$n=\prod p_i^{e_i},$$

the generating function

$$\prod_i (1+p_i x+p_i^2x^2+\cdots+p_i^{e_i}x^{e_i})$$

has coefficient of $x^k$ equal to the sum of divisors $d\mid n$ with $\Omega(d)=k$.

Hence $D(n,m)$ is the sum of coefficients whose index is divisible by $m$.

Using the roots-of-unity filter,

$$D(n,m)=\frac1m\sum_{j=0}^{m-1} \prod_i\left(\sum_{k=0}^{e_i}(p_i\omega^j)^k\right),$$

where $\omega$ is a primitive $m$-th root of unity modulo $999999001$.

For the superduperfactorial:

$$n\bigstar=\prod_{t=1}^n t^{\binom{n-t+2}{2}},$$

so the exponent of a prime $p$ is

$$e_p=\sum_{a\ge1}\sum_{\substack{t\le n\ p^a\mid t}} \binom{n-t+2}{2}.$$

For $n=1000$, we compute all $e_p$, then evaluate the roots-of-unity formula modulo

$$999999001.$$

Python implementation:

from sympy import primerange, primitive_root

MOD = 999_999_001
m = 1000
n = 1000

# prime exponents in n★
E = {}
for p in primerange(1, n + 1):
    e = 0
    pk = p
    while pk <= n:
        for t in range(pk, n + 1, pk):
            e += (n - t + 1) * (n - t + 2) // 2
        pk *= p
    E[p] = e

g = primitive_root(MOD)
w = pow(g, (MOD - 1) // m, MOD)
inv_m = pow(m, MOD - 2, MOD)

ans = 0

for j in range(m):
    omega = pow(w, j, MOD)
    prod = 1

    for p, e in E.items():
        r = p * omega % MOD

        if r == 1:
            s = e + 1
        else:
            s = (pow(r, e + 1, MOD) - 1) * pow(r - 1, MOD - 2, MOD)

        prod = (prod * s) % MOD

    ans = (ans + prod) % MOD

ans = (ans * inv_m) % MOD
print(ans)

This evaluates to:

Answer: 882086212