Detailed mathematical approach

Problem Summary

We must evaluate

$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$

for \(m=200!\) and \(n=10^{12}\), modulo \(10^9+7\). Here \(\sigma_0=\tau\) is the divisor-counting function. A direct double loop is hopeless, so the solution rewrites the divisor sum into an inclusion-exclusion over the primes of \(200!\), plus fast queries to the summatory divisor function.

Mathematical Approach

Write \(m=f!\) with \(f=200\), and rename the upper limit as \(N=n\). Let \(\mathbb{P}\) denote the set of prime numbers and define

$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$

For each \(p\in\mathcal{P}\), the exponent of \(p\) in \(f!\) is

$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor.$

The implementations exploit this prime-power structure very directly.

Step 1: Factor the divisor sum prime by prime

Every divisor of \(f!\) has the form

$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p.$

For a fixed integer \(k\), write \(b_p=v_p(k)\) for primes \(p\in\mathcal{P}\). Since \(\tau\) is multiplicative on prime powers, the sum over all divisors of \(f!\) separates into local contributions:

$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$

So the entire problem reduces to understanding the one-prime expression \(\sum_{a=0}^{e}(a+b+1)\).

Step 2: Rewrite the local sum with triangular numbers

Define the triangular-number function

$T(t)=\frac{t(t+1)}{2}.$

Then

$\sum_{a=0}^{e}(a+b+1)=(e+1)(b+1)+\frac{e(e+1)}{2}=T(e+1)(b+1)-T(e)b.$

This identity is the key algebraic step. The term \(T(e+1)(b+1)\) keeps the current exponent of \(p\) inside \(k\), while the correction term \(T(e)b\) corresponds to removing one copy of \(p\) when \(p\mid k\).

Step 3: Expand the product by inclusion-exclusion

Now define

$K=\prod_{p\in\mathcal{P}}T(e_p+1),\qquad \beta_p=\frac{T(e_p)}{T(e_p+1)}.$

Choosing the correction term at a prime \(p\) introduces a factor \(-\beta_p\) and can only happen when \(p\) already divides \(k\). If \(Q\subseteq\mathcal{P}\) is the set of primes where we choose that correction, and

$P_Q=\prod_{p\in Q}p,$

then one copy of every prime in \(Q\) is removed from \(k\). Therefore, for each fixed \(k\),

$\sum_{d \mid f!}\tau(kd)=K\sum_{Q\subseteq\mathcal{P},\ P_Q\mid k}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\tau\left(\frac{k}{P_Q}\right).$

This is exactly the inclusion-exclusion that the implementations realize through a recursive traversal of prime subsets.

Step 4: Swap the order of summation

Now sum over \(1\le k\le N\). For a fixed subset \(Q\), the condition \(P_Q\mid k\) lets us write \(k=P_Qr\), so

$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\sum_{r\le N/P_Q}\tau(r).$

Introduce the summatory divisor function

$S_\tau(x)=\sum_{r=1}^{x}\tau(r).$

Then the whole problem becomes

$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right).$

If \(P_Q>N\), the corresponding term is zero, which explains the branch pruning in the recursive search.

Step 5: Evaluate the summatory divisor function quickly

The identity

$S_\tau(x)=\sum_{d=1}^{x}\left\lfloor\frac{x}{d}\right\rfloor$

counts factor pairs \((a,b)\) with \(ab\le x\). Splitting those lattice points across the hyperbola \(ab=x\) gives the standard Dirichlet-hyperbola formula

$S_\tau(x)=2\sum_{d=1}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\frac{x}{d}\right\rfloor-\lfloor\sqrt{x}\rfloor^2.$

This is the large-\(x\) formula used by the implementation. Small values are handled by a precomputed prefix table.

Worked Example: \(m=3!=6\) and \(N=5\)

Here \(\mathcal{P}=\{2,3\}\) and \(e_2=e_3=1\). Hence

$T(2)=3,\qquad T(1)=1,\qquad K=3\cdot 3=9,\qquad \beta_2=\beta_3=\frac{1}{3}.$

The subset formula becomes

$D(6,5)=9S_\tau(5)-3S_\tau(2)-3S_\tau(1)+S_\tau(0).$

Now

$S_\tau(5)=1+2+2+3+2=10,\qquad S_\tau(2)=3,\qquad S_\tau(1)=1,\qquad S_\tau(0)=0,$

so

$D(6,5)=9\cdot 10-3\cdot 3-3\cdot 1=78.$

A direct check agrees:

$\sum_{k=1}^{5}\tau(k)=10,\qquad \sum_{k=1}^{5}\tau(2k)=17,\qquad \sum_{k=1}^{5}\tau(3k)=19,\qquad \sum_{k=1}^{5}\tau(6k)=32,$

and therefore

$10+17+19+32=78.$

This small case shows exactly how the inclusion-exclusion reproduces the brute-force sum.

How the Code Works

The C++, Python, and Java implementations first enumerate the primes up to \(200\) and compute each factorial exponent \(e_p\) with Legendre's formula. From those exponents they build the global factor \(K\) and the prime-specific ratios \(\beta_p\), always working modulo \(10^9+7\).

Next, they precompute \(S_\tau(x)\) for all \(x<10^6\) with a divisor sieve followed by prefix sums. For larger arguments they use the hyperbola formula above and store the results in a cache, so repeated requests for the same value are answered immediately.

The main recursive traversal walks through the primes in increasing order, maintains the current squarefree product \(P_Q\), the alternating sign, and the multiplicative weight \(K\prod_{p\in Q}\beta_p\). At each visited subset it adds the corresponding term

$K\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right)$

with the appropriate sign. Because multiplying by another prime only increases \(P_Q\), the recursion stops as soon as the next product would exceed \(N\).

Complexity Analysis

Let \(B=10^6\). Building the small prefix table costs \(O(B\log B)\) time and \(O(B)\) memory, because every divisor updates all of its multiples. The prime list and factorial exponents up to \(200\) are tiny by comparison.

Let \(R\) be the number of squarefree products of primes \(\le 200\) that do not exceed \(N\). The subset traversal uses \(O(R)\) arithmetic steps. If \(X\) is the set of distinct large arguments passed to \(S_\tau\), then the uncached large-\(x\) cost is

$O\left(\sum_{x\in X}\sqrt{x}\right),$

because each such value is computed once with the hyperbola identity and then memoized. In practice, the pruning \(P_Q\le N\) and the cache are what make the computation feasible.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=608
2. Divisor function: Wikipedia - Divisor function
3. Divisor summatory function: Wikipedia - Divisor summatory function
4. Legendre's formula: Wikipedia - Legendre's formula
5. Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
6. Dirichlet hyperbola method: Wikipedia - Dirichlet hyperbola method