Problem 956: Super Duper Sum

Project Euler Problem 956 Solution

EulerSolve provides an optimized solution for Project Euler Problem 956, Super Duper Sum, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The quantity behind Problem 956 can be collapsed to a single integer $$N_n=\prod_{t=1}^{n} t^{c_t},\qquad c_t=\binom{n-t+2}{2}=\frac{(n-t+1)(n-t+2)}{2}.$$ The required value is the sum of those divisors $d\mid N_n$ for which the total exponent count $$\Omega(d)=\sum_p \alpha_p$$ is divisible by $m$, where $d=\prod_p p^{\alpha_p}$. The implementations evaluate the case $n=m=1000$ and reduce the answer modulo $999999001$. A brute-force divisor enumeration is impossible. Even after prime factorization, the exponent bounds $E_p$ of $N_{1000}$ are enormous, so the number of divisors is astronomically large. The workable approach is to separate the contribution of each prime and track only the residue of the exponent sum modulo $m$. Mathematical Approach The entire solution is a prime-by-prime generating-function computation. The key is that the condition $\Omega(d)\equiv 0 \pmod m$ depends only on exponent residues, while the divisor value $d$ factors multiplicatively over primes. Collapsing the weighted product The code never expands the nested product directly. Instead it records how many times each integer $t$ appears after everything is flattened. That multiplicity is the triangular number $c_t$, so the whole object is exactly $N_n=\prod t^{c_t}$....

Detailed mathematical approach

Problem Summary

The quantity behind Problem 956 can be collapsed to a single integer

$$N_n=\prod_{t=1}^{n} t^{c_t},\qquad c_t=\binom{n-t+2}{2}=\frac{(n-t+1)(n-t+2)}{2}.$$

The required value is the sum of those divisors $d\mid N_n$ for which the total exponent count

$$\Omega(d)=\sum_p \alpha_p$$

is divisible by $m$, where $d=\prod_p p^{\alpha_p}$. The implementations evaluate the case $n=m=1000$ and reduce the answer modulo $999999001$.

A brute-force divisor enumeration is impossible. Even after prime factorization, the exponent bounds $E_p$ of $N_{1000}$ are enormous, so the number of divisors is astronomically large. The workable approach is to separate the contribution of each prime and track only the residue of the exponent sum modulo $m$.

Mathematical Approach

The entire solution is a prime-by-prime generating-function computation. The key is that the condition $\Omega(d)\equiv 0 \pmod m$ depends only on exponent residues, while the divisor value $d$ factors multiplicatively over primes.

Collapsing the weighted product

The code never expands the nested product directly. Instead it records how many times each integer $t$ appears after everything is flattened. That multiplicity is the triangular number $c_t$, so the whole object is exactly $N_n=\prod t^{c_t}$.

For any prime $p$, the exponent of $p$ in $N_n$ is therefore

$$E_p=v_p(N_n)=\sum_{t=1}^{n} c_t\,v_p(t).$$

Once all $E_p$ are known, the original problem has been reduced to a divisor-sum problem on the prime factorization

$$N_n=\prod_p p^{E_p}.$$

Divisors as exponent vectors

Every divisor of $N_n$ has the form

$$d=\prod_p p^{\alpha_p},\qquad 0\le \alpha_p\le E_p.$$

The quantity being summed is the value of the divisor itself, not just the number of admissible divisors. So a choice of exponent vector $(\alpha_p)$ contributes

$$\prod_p p^{\alpha_p}.$$

The only global constraint is

$$\Omega(d)=\sum_p \alpha_p\equiv 0 \pmod m.$$

This is exactly the kind of condition that can be handled by a residue-class convolution.

A generating function that matches the divisor sum

For a fixed prime $p$, introduce

$$F_p(x)=\sum_{a=0}^{E_p} p^a x^a.$$

Multiplying these factors over all primes gives

$$F(x)=\prod_p F_p(x).$$

The coefficient of $x^k$ in $F(x)$ is precisely the sum of all divisor values $d\mid N_n$ with $\Omega(d)=k$, because choosing the term $p^{\alpha_p}x^{\alpha_p}$ from each prime contributes the divisor value and records the total exponent count in the power of $x$.

We do not need the exact exponent total $k$; only its residue modulo $m$ matters. So we can work in the quotient ring where $x^m=1$. In that setting, each prime contributes only $m$ residue classes:

$$F_p(x)\equiv \sum_{r=0}^{m-1} C_{p,r}x^r \pmod{x^m-1},$$

with

$$C_{p,r}=\sum_{\substack{0\le a\le E_p\\ a\equiv r \pmod m}} p^a.$$

The desired answer is then the coefficient of $x^0$ in the product of these reduced polynomials.

Closed form for one residue class

If $E_p<r$, then $C_{p,r}=0$. Otherwise the admissible exponents are

$$a=r,\ r+m,\ r+2m,\ \dots,\ r+(N_{p,r}-1)m,$$

where

$$N_{p,r}=\left\lfloor\frac{E_p-r}{m}\right\rfloor+1.$$

Let $q=p^m$. Then

$$C_{p,r}=p^r\sum_{j=0}^{N_{p,r}-1} q^j.$$

Modulo $M=999999001$, this becomes a geometric series. If $q\not\equiv 1 \pmod M$, then

$$C_{p,r}\equiv p^r\,(q^{N_{p,r}}-1)\,(q-1)^{-1}\pmod M.$$

If $q\equiv 1 \pmod M$, the denominator would vanish, and the sum collapses to

$$C_{p,r}\equiv p^r\,N_{p,r}\pmod M.$$

This is the only subtle modular case in the entire computation.

Worked example: $n=m=3$

Here

$$c_1=6,\qquad c_2=3,\qquad c_3=1,$$

$$N_3=1^6\cdot 2^3\cdot 3=24=2^3\cdot 3.$$

We want divisors whose total exponent count is a multiple of $3$. For the prime $2$ with exponent $3$, the residue-class sums are

$$C_{2,0}=1+2^3=9,\qquad C_{2,1}=2,\qquad C_{2,2}=4.$$

For the prime $3$ with exponent $1$, they are

$$C_{3,0}=1,\qquad C_{3,1}=3,\qquad C_{3,2}=0.$$

To obtain total residue $0\pmod 3$, we combine residue pairs $(0,0)$, $(2,1)$, and $(1,2)$. Therefore the required sum is

$$9\cdot 1+4\cdot 3+2\cdot 0=21.$$

Indeed the valid divisors are $1$, $8$, and $12$, and their sum is $21$. This small case is exactly the same convolution that the full solution performs for $n=m=1000$.

Why only $m$ states are needed

After reducing exponents modulo $m$, the coefficient vector for each prime has length $m$. Multiplying in another prime means a cyclic convolution on those $m$ residues. No larger state space is needed, because the problem never distinguishes between exponent totals that differ by a multiple of $m$.

How the Code Works

Prime-exponent preprocessing

The C++, Python, and Java implementations first build a smallest-prime-factor table up to $n$. Then they scan $t=1,2,\dots,n$, factor each $t$, and add $c_t\,v_p(t)$ to the stored exponent of every prime $p$ appearing in $t$. At the end of this pass they know every nonzero $E_p$ in the factorization of $N_n$.

Residue tables for each prime

For each prime $p$, the implementation computes the $m$ numbers $C_{p,0},C_{p,1},\dots,C_{p,m-1}$ modulo $M$. Instead of summing powers one by one, it maintains the current factor $p^r$ and evaluates the remaining geometric series in closed form. The branch $p^m\equiv 1\pmod M$ is handled separately so that no invalid modular inverse is taken.

The dynamic-programming invariant

After processing any subset of the primes, the state vector satisfies the invariant

The state $\mathrm{dp}[r]$ is the sum of divisor values formed from the processed primes whose total exponent residue is $r$.

When a new prime is incorporated, the update is the cyclic convolution

$$\mathrm{next}[(r+t)\bmod m]\;{+}{=}\;\mathrm{dp}[r]\cdot C_{p,t}\pmod M.$$

After the final prime has been processed, $\mathrm{dp}[0]$ is exactly the required answer. The C++ and Java versions also compare the modular routine against a tiny exact computation on a small sample, which is a good sanity check for the formulae.

Complexity Analysis

The sieve used to build smallest prime factors is linear in $n$. Factoring every $t\le n$ by repeated division through that table is close to linear for these input sizes and is much cheaper than the final convolution phase.

The dominant cost is the residue DP. If $\pi(n)$ denotes the number of primes up to $n$, then the transition over all primes uses $O(\pi(n)m^2)$ modular operations, because each prime contributes an $m\times m$ cyclic convolution. Memory usage is $O(n+m)$: the sieve and exponent arrays take $O(n)$, and the algorithm needs only a constant number of DP buffers of length $m$.

Footnotes and References

Problem page: https://projecteuler.net/problem=956
Prime factorization: Wikipedia - Prime factorization
Divisor function: Wikipedia - Divisor function
Prime omega function: Wikipedia - Prime omega function
Generating function: Wikipedia - Generating function
Geometric series: Wikipedia - Geometric series
Dynamic programming: Wikipedia - Dynamic programming

Problem 956 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>
#include <string>
#include <utility>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using boost::multiprecision::cpp_int;

constexpr u64 kMod = 999'999'001ULL;

u64 mod_mul(u64 a, u64 b, u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}

u64 mod_pow(u64 base, u64 exp, u64 mod) {
    u64 result = 1ULL % mod;
    u64 cur = base % mod;
    while (exp > 0) {
        if ((exp & 1ULL) != 0ULL) {
            result = mod_mul(result, cur, mod);
        }
        cur = mod_mul(cur, cur, mod);
        exp >>= 1ULL;
    }
    return result;
}

u64 mod_inv(u64 x, u64 mod) {
    return mod_pow(x, mod - 2ULL, mod);
}

std::vector<int> build_spf(int n) {
    std::vector<int> spf(n + 1, 0);
    std::vector<int> primes;
    primes.reserve(n / 5);

    for (int i = 2; i <= n; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            const u64 v = static_cast<u64>(i) * static_cast<u64>(p);
            if (v > static_cast<u64>(n) || p > spf[i]) {
                break;
            }
            spf[static_cast<std::size_t>(v)] = p;
        }
    }
    return spf;
}

std::vector<std::pair<int, u64>> prime_exponents_superduper(int n) {
    const std::vector<int> spf = build_spf(n);
    std::vector<u64> exp(static_cast<std::size_t>(n + 1), 0ULL);

    for (int t = 1; t <= n; ++t) {
        const u64 coef = static_cast<u64>(n - t + 1) * static_cast<u64>(n - t + 2) / 2ULL;
        int x = t;
        while (x > 1) {
            const int p = spf[static_cast<std::size_t>(x)];
            int cnt = 0;
            while (x % p == 0) {
                x /= p;
                ++cnt;
            }
            exp[static_cast<std::size_t>(p)] += coef * static_cast<u64>(cnt);
        }
    }

    std::vector<std::pair<int, u64>> factors;
    for (int p = 2; p <= n; ++p) {
        if (exp[static_cast<std::size_t>(p)] != 0ULL) {
            factors.push_back({p, exp[static_cast<std::size_t>(p)]});
        }
    }
    return factors;
}

u64 D_mod(int n, int m, u64 mod) {
    const auto factors = prime_exponents_superduper(n);
    std::vector<u64> dp(static_cast<std::size_t>(m), 0ULL);
    std::vector<u64> next(static_cast<std::size_t>(m), 0ULL);
    std::vector<u64> coeff(static_cast<std::size_t>(m), 0ULL);
    dp[0] = 1ULL;

    for (const auto& [p_int, e] : factors) {
        const u64 p = static_cast<u64>(p_int);
        const u64 q = mod_pow(p, static_cast<u64>(m), mod);

        u64 p_pow_t = 1ULL;
        for (int t = 0; t < m; ++t) {
            if (e >= static_cast<u64>(t)) {
                const u64 cnt = (e - static_cast<u64>(t)) / static_cast<u64>(m) + 1ULL;
                u64 geom = 0ULL;
                if (q == 1ULL) {
                    geom = cnt % mod;
                } else {
                    const u64 num = (mod_pow(q, cnt, mod) + mod - 1ULL) % mod;
                    const u64 den_inv = mod_inv((q + mod - 1ULL) % mod, mod);
                    geom = mod_mul(num, den_inv, mod);
                }
                coeff[static_cast<std::size_t>(t)] = mod_mul(p_pow_t, geom, mod);
            } else {
                coeff[static_cast<std::size_t>(t)] = 0ULL;
            }
            p_pow_t = mod_mul(p_pow_t, p, mod);
        }

        std::fill(next.begin(), next.end(), 0ULL);
        for (int r = 0; r < m; ++r) {
            if (dp[static_cast<std::size_t>(r)] == 0ULL) {
                continue;
            }
            const u64 base = dp[static_cast<std::size_t>(r)];
            for (int t = 0; t < m; ++t) {
                const int idx = (r + t >= m) ? (r + t - m) : (r + t);
                next[static_cast<std::size_t>(idx)] =
                    (next[static_cast<std::size_t>(idx)] + mod_mul(base, coeff[static_cast<std::size_t>(t)], mod)) % mod;
            }
        }
        dp.swap(next);
    }

    return dp[0];
}

cpp_int D_exact_small(int n, int m) {
    const auto factors = prime_exponents_superduper(n);
    std::vector<cpp_int> dp(static_cast<std::size_t>(m), cpp_int(0));
    std::vector<cpp_int> next(static_cast<std::size_t>(m), cpp_int(0));
    std::vector<cpp_int> coeff(static_cast<std::size_t>(m), cpp_int(0));
    dp[0] = 1;

    for (const auto& [p_int, e] : factors) {
        const cpp_int p = p_int;

        for (int t = 0; t < m; ++t) {
            if (e < static_cast<u64>(t)) {
                coeff[static_cast<std::size_t>(t)] = 0;
                continue;
            }
            cpp_int term = 1;
            for (int i = 0; i < t; ++i) {
                term *= p;
            }
            cpp_int sum = 0;
            for (u64 a = static_cast<u64>(t); a <= e; a += static_cast<u64>(m)) {
                sum += term;
                for (int i = 0; i < m; ++i) {
                    term *= p;
                }
                if (e - a < static_cast<u64>(m)) {
                    break;
                }
            }
            coeff[static_cast<std::size_t>(t)] = sum;
        }

        std::fill(next.begin(), next.end(), cpp_int(0));
        for (int r = 0; r < m; ++r) {
            if (dp[static_cast<std::size_t>(r)] == 0) {
                continue;
            }
            for (int t = 0; t < m; ++t) {
                const int idx = (r + t >= m) ? (r + t - m) : (r + t);
                next[static_cast<std::size_t>(idx)] += dp[static_cast<std::size_t>(r)] * coeff[static_cast<std::size_t>(t)];
            }
        }
        dp.swap(next);
    }

    return dp[0];
}

void run_validations() {
    const cpp_int sample_exact = D_exact_small(6, 6);
    assert(sample_exact == cpp_int("6368195719791280"));
    const u64 sample_mod = D_mod(6, 6, kMod);
    assert(sample_mod == static_cast<u64>(sample_exact % kMod));
}

}  // namespace

int main() {
    run_validations();
    std::cout << D_mod(1'000, 1'000, kMod) << '\n';
    return 0;
}

Python

def solve():
    MOD = 999999001
    n_val = 1000; m_val = 1000

    def mod_pow(base, exp, mod=MOD):
        r = 1; base %= mod
        while exp > 0:
            if exp & 1: r = r * base % mod
            base = base * base % mod; exp >>= 1
        return r

    def mod_inv(x): return mod_pow(x, MOD - 2)

    # SPF sieve
    spf = list(range(n_val + 1))
    primes = []
    for i in range(2, n_val + 1):
        if spf[i] == i: primes.append(i)
        for p in primes:
            if i * p > n_val or p > spf[i]: break
            spf[i * p] = p

    # Prime exponents of superproduct
    exp_arr = [0] * (n_val + 1)
    for t in range(1, n_val + 1):
        coef = (n_val - t + 1) * (n_val - t + 2) // 2
        x = t
        while x > 1:
            p = spf[x]; cnt = 0
            while x % p == 0: x //= p; cnt += 1
            exp_arr[p] += coef * cnt

    factors = [(p, exp_arr[p]) for p in range(2, n_val + 1) if exp_arr[p] > 0]

    dp = [0] * m_val; dp[0] = 1
    for p_int, e in factors:
        p = p_int
        q = mod_pow(p, m_val)
        coeff = [0] * m_val
        pp = 1
        for t in range(m_val):
            if e >= t:
                cnt = (e - t) // m_val + 1
                if q == 1: geom = cnt % MOD
                else:
                    num = (mod_pow(q, cnt) - 1) % MOD
                    geom = num * mod_inv((q - 1) % MOD) % MOD
                coeff[t] = pp * geom % MOD
            pp = pp * p % MOD

        nxt = [0] * m_val
        for r in range(m_val):
            if dp[r] == 0: continue
            for t in range(m_val):
                idx = (r + t) % m_val
                nxt[idx] = (nxt[idx] + dp[r] * coeff[t]) % MOD
        dp = nxt

    return str(dp[0])

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public class Euler956 {

    static final long K_MOD = 999999001L;

    static long modPow(long base, long exp, long mod) {
        long result = 1L % mod;
        long cur = base % mod;
        while (exp > 0) {
            if ((exp & 1L) != 0L) {
                result = (result * cur) % mod;
            }
            cur = (cur * cur) % mod;
            exp >>= 1L;
        }
        return result;
    }

    static long modInv(long x, long mod) {
        return modPow(x, mod - 2L, mod);
    }

    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>(n / 5);

        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long v = (long) i * p;
                if (v > n || p > spf[i]) {
                    break;
                }
                spf[(int) v] = p;
            }
        }
        return spf;
    }

    static class Factor {
        int p;
        long e;

        Factor(int p, long e) {
            this.p = p;
            this.e = e;
        }
    }

    static List<Factor> primeExponentsSuperduper(int n) {
        int[] spf = buildSpf(n);
        long[] exp = new long[n + 1];

        for (int t = 1; t <= n; ++t) {
            long coef = (long) (n - t + 1) * (n - t + 2) / 2L;
            int x = t;
            while (x > 1) {
                int p = spf[x];
                int cnt = 0;
                while (x % p == 0) {
                    x /= p;
                    ++cnt;
                }
                exp[p] += coef * cnt;
            }
        }

        List<Factor> factors = new ArrayList<>();
        for (int p = 2; p <= n; ++p) {
            if (exp[p] != 0L) {
                factors.add(new Factor(p, exp[p]));
            }
        }
        return factors;
    }

    static long dMod(int n, int m, long mod) {
        List<Factor> factors = primeExponentsSuperduper(n);
        long[] dp = new long[m];
        long[] next = new long[m];
        long[] coeff = new long[m];
        dp[0] = 1L;

        for (Factor f : factors) {
            long p = f.p;
            long e = f.e;
            long q = modPow(p, m, mod);

            long pPowT = 1L;
            for (int t = 0; t < m; ++t) {
                if (e >= t) {
                    long cnt = (e - t) / m + 1L;
                    long geom;
                    if (q == 1L) {
                        geom = cnt % mod;
                    } else {
                        long num = (modPow(q, cnt, mod) + mod - 1L) % mod;
                        long denInv = modInv((q + mod - 1L) % mod, mod);
                        geom = (num * denInv) % mod;
                    }
                    coeff[t] = (pPowT * geom) % mod;
                } else {
                    coeff[t] = 0L;
                }
                pPowT = (pPowT * p) % mod;
            }

            Arrays.fill(next, 0L);
            for (int r = 0; r < m; ++r) {
                if (dp[r] == 0L)
                    continue;
                long base = dp[r];
                for (int t = 0; t < m; ++t) {
                    int idx = (r + t >= m) ? (r + t - m) : (r + t);
                    next[idx] = (next[idx] + (base * coeff[t]) % mod) % mod;
                }
            }
            System.arraycopy(next, 0, dp, 0, m);
        }

        return dp[0];
    }

    public static String solve() {
        return Long.toString(dMod(1000, 1000, K_MOD));
    }

    public static void main(String[] args) {
        BigInteger exact = new BigInteger("6368195719791280");
        long sampleMod = dMod(6, 6, K_MOD);
        if (sampleMod != exact.mod(BigInteger.valueOf(K_MOD)).longValue()) {
            System.out.println("Validation failed");
            return;
        }
        System.out.println(solve());
    }
}