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,$

so

$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

1. Problem page: https://projecteuler.net/problem=956
2. Prime factorization: Wikipedia - Prime factorization
3. Divisor function: Wikipedia - Divisor function
4. Prime omega function: Wikipedia - Prime omega function
5. Generating function: Wikipedia - Generating function
6. Geometric series: Wikipedia - Geometric series
7. Dynamic programming: Wikipedia - Dynamic programming