Detailed mathematical approach

Problem Summary

For each \(n\), define

$B_n=\prod_{k=0}^{n}\binom{n}{k},\qquad D_n=\sigma(B_n),\qquad S(N)=\sum_{n=1}^{N} D_n \pmod{10^9+7}.$

The task is to compute \(S(20000)\). The number \(B_n\) grows far too quickly to build directly, so the useful object is not \(B_n\) itself but its prime factorization. Once the exponent of every prime in \(B_n\) is known, the divisor-sum function \(\sigma\) becomes a clean multiplicative product.

Mathematical Approach

The C++, Python, and Java implementations all follow the same strategy: rewrite \(B_n\) in factorial form, describe each prime exponent exactly, derive an update rule from \(n-1\) to \(n\), and then convert the resulting factorization into \(D_n\).

Step 1: Rewrite the product of binomial coefficients

Each factor satisfies

$\binom{n}{k}=\frac{n!}{k!(n-k)!}.$

Multiplying over all \(k=0,1,\dots,n\) gives

$B_n=\prod_{k=0}^{n}\binom{n}{k}=\frac{(n!)^{n+1}}{\left(\prod_{j=0}^{n} j!\right)^2}.$

The numerator contributes \(n+1\) copies of \(n!\). The denominator contains one copy of every \(k!\) and one copy of every \((n-k)!\), which is the same list reversed, so the whole denominator becomes the square of \(\prod_{j=0}^{n} j!\).

Step 2: Express \(B_n\) through prime valuations

Fix a prime \(p\). Let

$A_{n,p}=v_p(n!),\qquad T_{n,p}=\sum_{j=0}^{n} A_{j,p},$

where \(v_p(m)\) is the exponent of \(p\) in \(m\). Equivalently, \(A_{n,p}\) is the quantity described by Legendre's formula.

If \(e_{n,p}=v_p(B_n)\), then the factorial identity above yields

$e_{n,p}=(n+1)A_{n,p}-2T_{n,p}.$

So the entire problem reduces to keeping track of these exponents for all primes \(p\le n\).

Step 3: Derive the increment from \(n-1\) to \(n\)

The implementations do not recompute \(e_{n,p}\) from scratch. Instead they update it incrementally. Since

$A_{n,p}=A_{n-1,p}+v_p(n),\qquad T_{n,p}=T_{n-1,p}+A_{n,p},$

we get

$\begin{aligned} e_{n,p}-e_{n-1,p} &= (n+1)A_{n,p}-2T_{n,p}-\left(nA_{n-1,p}-2T_{n-1,p}\right) \\ &= (n-1)v_p(n)-A_{n-1,p}. \end{aligned}$

This formula is the heart of the solution. Every prime receives the subtractive term \(-A_{n-1,p}\), and only primes dividing \(n\) receive the additive correction \((n-1)v_p(n)\).

Step 4: Turn prime exponents into the divisor sum

Once

$B_n=\prod_{p\le n} p^{e_{n,p}}$

is known, multiplicativity of the divisor-sum function gives

$D_n=\sigma(B_n)=\prod_{p\le n}\sigma\!\left(p^{e_{n,p}}\right)=\prod_{p\le n}\frac{p^{e_{n,p}+1}-1}{p-1} \pmod{10^9+7}.$

That is why the implementations maintain \(p^{e_{n,p}+1}\) modulo \(10^9+7\) rather than the bare exponent itself. Because the modulus is prime and much larger than every prime used here, both \(p\) and \(p-1\) have modular inverses.

Step 5: Worked example for \(n=5\)

From the product definition,

$B_5=\binom{5}{0}\binom{5}{1}\binom{5}{2}\binom{5}{3}\binom{5}{4}\binom{5}{5}=1\cdot 5\cdot 10\cdot 10\cdot 5\cdot 1=2500.$

The valuation formula recovers the same factorization. For the relevant primes,

$\begin{aligned} A_{5,2}&=v_2(5!)=3, & T_{5,2}&=0+0+1+1+3+3=8, & e_{5,2}&=6\cdot 3-2\cdot 8=2, \\ A_{5,3}&=v_3(5!)=1, & T_{5,3}&=0+0+0+1+1+1=3, & e_{5,3}&=6\cdot 1-2\cdot 3=0, \\ A_{5,5}&=v_5(5!)=1, & T_{5,5}&=0+0+0+0+0+1=1, & e_{5,5}&=6\cdot 1-2\cdot 1=4. \end{aligned}$

So

$B_5=2^2\cdot 5^4,$

and therefore

$D_5=\frac{2^3-1}{2-1}\cdot\frac{5^5-1}{5-1}=7\cdot 781=5467.$

The earlier values are \(D_1=1\), \(D_2=3\), \(D_3=13\), and \(D_4=252\), so

$S(5)=1+3+13+252+5467=5736,$

which matches the checkpoint used by the implementations.

How the Code Works

The implementations begin with a smallest-prime-factor sieve up to \(N\). This provides the full prime list and allows each new integer \(n\) to be factorized quickly.

For every prime, the implementation stores two modular states. One state represents the inverse of the prime power already accumulated in the current factorial valuation, and the other represents the current value of \(p^{e_{n,p}+1}\). At the start of the transition from \(n-1\) to \(n\), every stored prime-power term is multiplied by the inverse factorial contribution corresponding to \(-A_{n-1,p}\).

Then the current integer \(n\) is factorized. If \(p^v\) divides \(n\), the stored value for that prime is multiplied by \(p^{(n-1)v}\), and the stored inverse factorial contribution is extended by another factor of \(p^{-v}\). After these updates, the maintained prime-power term is exactly \(p^{e_{n,p}+1}\).

Finally, the implementation evaluates

$D_n=\prod_{p}\frac{p^{e_{n,p}+1}-1}{p-1}\pmod{10^9+7}$

across the prime list and adds the result to the running total. Primes larger than the current \(n\) still have exponent \(0\), so their factor is \(1\) and they do not affect the answer.

Complexity Analysis

The smallest-prime-factor sieve is built in \(O(N)\) time and uses \(O(N)\) memory. Factoring each \(n\) through that sieve costs \(O(\Omega(n))\) in the straightforward bound, with total average-order work \(O(N\log\log N)\). In these implementations the dominant term is the full prime scan performed for each \(n\): once to apply the common recurrence factor and once to assemble \(D_n\). Therefore the overall running time is \(O(N\pi(N)+N\log\log N)\), usually summarized as \(O(N\pi(N))\), and the extra memory beyond the sieve is \(O(\pi(N))\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=650
2. Binomial coefficient: Wikipedia - Binomial coefficient
3. Legendre's formula: Wikipedia - Legendre's formula
4. Divisor function: Wikipedia - Divisor function
5. Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes