Detailed mathematical approach

Problem Summary

Write the factorial as

$n!=\prod_{p\le n}p^{a_p},\qquad a_p=v_p(n!)=\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor.$

The quantity to compute is the signed divisor sum

$S(n;L,H)=\sum_{\substack{d\mid n!\\L\le d\le H}}\lambda(d)\,d,\qquad \lambda(d)=(-1)^{\Omega(d)},$

reduced modulo \(10^9+7\). The interval \([L,H]\) is enormous in the target case, so the solution cannot enumerate divisors naively; it must exploit the multiplicative structure of \(n!\) and the fact that only divisors up to \(H\) can matter.

Mathematical Approach

The main idea is to represent each divisor by its prime exponents, split those primes into two independent groups, and turn the interval condition into two prefix queries.

Step 1: Prime Exponents of the Factorial

For every prime \(p\le n\), Legendre's formula gives the exponent \(a_p\) of \(p\) in \(n!\). Therefore every divisor of \(n!\) has a unique form

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

So the problem is equivalent to iterating over all allowable exponent vectors \((e_p)\), but doing so intelligently enough that the upper bound \(H\) removes most impossible branches before they are fully built.

Step 2: The Signed Weight is Completely Multiplicative on Prime Powers

If \(d=\prod p^{e_p}\), then

$\lambda(d)\,d=(-1)^{\Omega(d)}\prod p^{e_p}=\prod (-p)^{e_p}.$

Without the interval restriction, the total signed sum over all divisors would factor as

$\sum_{d\mid n!}\lambda(d)\,d=\prod_{p\le n}\left(\sum_{e=0}^{a_p}(-p)^e\right).$

The interval \(L\le d\le H\) destroys that simple product formula, because different primes now interact through the size of the full divisor. That is exactly why the implementation switches to a meet-in-the-middle strategy.

Step 3: Split the Prime Factors into Two Independent Groups

Partition the prime powers of \(n!\) into two disjoint groups. For the first group define partial states

$u=\prod_{p\in A}p^{e_p},\qquad \alpha=\prod_{p\in A}(-p)^{e_p},$

and for the second group define

$v=\prod_{p\in B}p^{e_p},\qquad \beta=\prod_{p\in B}(-p)^{e_p}.$

Because the groups use disjoint primes, every divisor corresponds to exactly one pair of partial states, with

$d=u\,v,\qquad \lambda(d)\,d=\alpha\beta.$

This reduces one very large search over all prime exponents to a merge of two precomputed lists.

Step 4: Enumerate Only Partial Products that Stay at Most \(H\)

During recursion, if the current partial product already exceeds \(H\), or if multiplying by one more copy of the current prime would exceed \(H\), then no continuation can ever contribute to a divisor in \([L,H]\) or to a prefix sum up to \(H\). Every extra prime power only makes the divisor larger.

So each side generates only feasible pairs \((u,\alpha)\) or \((v,\beta)\) with numeric value at most \(H\). The sign is carried through the factor \(-p\), while the divisor magnitude itself is stored exactly for later ordering and comparison.

Step 5: Turn the Interval into Prefix Queries

Define the bounded prefix function

$P(X)=\sum_{\substack{d\mid n!\\d\le X}}\lambda(d)\,d.$

After splitting, this becomes

$P(X)=\sum_{(u,\alpha)}\alpha\left(\sum_{\substack{(v,\beta)\\v\le X/u}}\beta\right).$

If the second list is sorted by \(v\) and its signed weights are accumulated into prefix sums, then the inner quantity is available immediately once we know how far the threshold \(X/u\) reaches. Sorting the first list as well allows a monotone sweep: as \(u\) grows, the allowed bound \(X/u\) shrinks, so one pointer moves only in one direction.

Step 6: Recover the Required Interval Sum

The target interval sum is just the difference of two prefixes:

$S(n;L,H)=P(H)-P(L-1)\pmod{10^9+7}.$

That identity is the final reduction. The algorithm therefore computes two threshold queries, one at \(H\) and one at \(L-1\), and subtracts them modulo \(10^9+7\).

Worked Example: \(n=5\), \(L=6\), \(H=30\)

Here

$5!=2^3\cdot 3\cdot 5.$

Split the prime factors into \(A=\{2^3,3\}\) and \(B=\{5\}\). The first side produces the sorted partial states

$\{(1,1),(2,-2),(3,-3),(4,4),(6,6),(8,-8),(12,-12),(24,24)\}.$

The second side produces

$\{(1,1),(5,-5)\},$

whose prefix sums of weights are \(1\) and \(-4\).

For \(P(30)\), the bound \(30/u\) is at least \(5\) when \(u=1,2,3,4,6\), so the inner sum is \(-4\). For \(u=8,12,24\), only the state \(v=1\) is allowed, so the inner sum is \(1\). Thus

$P(30)=-4+8+12-16-24-8-12+24=-20.$

For \(P(5)\), only \(u=1,2,3,4\) contribute, giving

$P(5)=-4-2-3+4=-5.$

Therefore

$S(5;6,30)=P(30)-P(5)=-15\equiv 10^9+7-15 \pmod{10^9+7}.$

This small example is exactly the same merge logic used in the full computation.

How the Code Works

The C++, Python, and Java implementations first generate all primes up to \(n\) and compute their exponents in \(n!\) by repeated division. They then split the resulting prime-factor list after the first five primes, which keeps both recursive enumerations manageable for the target instance.

Each side recursively enumerates every partial divisor whose exact value does not exceed \(H\). Alongside the exact divisor value, the implementation carries the signed multiplicative weight modulo \(10^9+7\). The divisor values themselves use exact wide integers: 256-bit arithmetic in C++ and arbitrary-precision integers in Python and Java, so comparisons against \(10^{60}\) remain exact.

After sorting both partial-state lists by divisor value, the implementation builds prefix sums of the signed weights on the second list. A linear sweep over the first list then evaluates \(P(H)\) and \(P(L-1)\). If \(L\le 1\), the lower prefix is zero; otherwise the algorithm evaluates the second threshold normally and subtracts the two results modulo \(10^9+7\).

Complexity Analysis

Generating primes up to \(n\) costs \(O(n\log\log n)\) time and \(O(n)\) memory. Computing the exponents \(a_p\) by repeated division contributes \(O(\pi(n)\log n)\), which is small here. If the two recursive enumerations produce \(A\) and \(B\) feasible partial states after pruning, then enumeration costs \(O(A+B)\), sorting costs \(O(A\log A+B\log B)\), and each threshold sweep costs \(O(A+B)\) because the pointer over the second list moves monotonically. The memory usage is \(O(A+B)\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=646
2. Legendre's formula: Wikipedia — Legendre's formula
3. Liouville function: Wikipedia — Liouville function
4. Prefix sum: Wikipedia — Prefix sum
5. Divisor: Wikipedia — Divisor