Detailed mathematical approach

Problem Summary

For positive integers \(a=\prod p^{\alpha_p}\) and \(b=\prod p^{\beta_p}\), define their exponent difference by

$D(a,b)=\sum_{p}\left|\alpha_p-\beta_p\right|=\sum_{p}\left|\nu_p(a)-\nu_p(b)\right|.$

The task is to evaluate

$S(n)=\sum_{1\le a,b\le n} D(a,b)$

for the Project Euler instance \(n=10^{12}\), with the final result taken modulo \(10^9+7\). A direct scan of all ordered pairs would require \(10^{24}\) pair checks before even factoring the numbers, so the only viable route is to reorganize the sum prime by prime and prime-power by prime-power.

Mathematical Approach

The key simplification is that the contribution of each prime is independent, and the absolute difference of two exponents can be rewritten as a count of divisibility levels.

Step 1: Separate the contribution of each prime

Because only finitely many primes divide either \(a\) or \(b\), we may interchange the order of summation:

$S(n)=\sum_{p\le n}\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$

For a fixed prime \(p\), define

$T_p(n)=\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$

Then

$S(n)=\sum_{p\le n} T_p(n).$

So the entire problem reduces to understanding the contribution of one prime at a time.

Step 2: Rewrite the absolute difference as level indicators

For any nonnegative integers \(x\) and \(y\),

$|x-y|=\sum_{k\ge 1}\left|\mathbf{1}_{x\ge k}-\mathbf{1}_{y\ge k}\right|.$

Applied to \(x=\nu_p(a)\) and \(y=\nu_p(b)\), this says that level \(k\) contributes \(1\) exactly when one of \(a,b\) is divisible by \(p^k\) and the other is not. Therefore

$T_p(n)=\sum_{k\ge 1} N_{p,k},$

where \(N_{p,k}\) counts ordered pairs with different divisibility status by \(p^k\).

Step 3: Count one divisibility level

Let

$q_{p,k}=\#\{m\le n: p^k\mid m\}=\left\lfloor\frac{n}{p^k}\right\rfloor.$

There are \(q_{p,k}\) integers up to \(n\) divisible by \(p^k\), and \(n-q_{p,k}\) that are not. Since the pairs are ordered, the number of mismatches is

$N_{p,k}=q_{p,k}(n-q_{p,k})+(n-q_{p,k})q_{p,k}=2q_{p,k}(n-q_{p,k}).$

Hence

$T_p(n)=\sum_{k\ge 1} 2q_{p,k}(n-q_{p,k}),\qquad q_{p,k}=\left\lfloor\frac{n}{p^k}\right\rfloor.$

Only finitely many terms are nonzero, because \(q_{p,k}=0\) once \(p^k\gt n\).

Step 4: Collapse the problem into prime-power contributions

Substituting the previous identity yields

$\boxed{S(n)=2\sum_{p\le n}\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor\left(n-\left\lfloor\frac{n}{p^k}\right\rfloor\right).}$

So every prime power \(p^k\le n\) contributes one term that depends only on the quotient \(\left\lfloor n/p^k\right\rfloor\). The algorithm is therefore an efficient way to enumerate those prime powers, or to group them when many primes share the same quotient.

Step 5: Split small primes and large primes

The implementations choose the cutoff

$B=\min(n,10^8).$

For each prime \(p\le B\), all powers \(p,p^2,p^3,\dots\le n\) are visited explicitly. For the Project Euler instance \(n=10^{12}\), any prime \(p\gt B\) automatically satisfies \(p^2\gt n\), so such primes contribute only through the first power \(p\) itself.

That observation removes all higher-power work from the large-prime side.

Step 6: Group large primes by equal quotients

For a large prime \(p\gt B\), the only relevant quantity is

$q=\left\lfloor\frac{n}{p}\right\rfloor.$

All primes producing the same quotient \(q\) lie in the interval

$\frac{n}{q+1}\lt p\le \frac{n}{q}.$

Instead of iterating through those primes individually, we count how many primes lie in that interval and multiply by their common contribution

$2q(n-q).$

The number of primes in any interval \([L,R]\) is

$\pi(R)-\pi(L-1),$

so fast evaluation of the prime-counting function \(\pi(x)\) completes the large-prime part.

Worked Example: \(n=10\)

The primes up to \(10\) are \(2,3,5,7\).

For \(p=2\), the nonzero quotients are \(q_{2,1}=5\), \(q_{2,2}=2\), \(q_{2,3}=1\), giving

$2\cdot 5\cdot 5+2\cdot 2\cdot 8+2\cdot 1\cdot 9=50+32+18=100.$

For \(p=3\), we get \(q_{3,1}=3\) and \(q_{3,2}=1\), so

$2\cdot 3\cdot 7+2\cdot 1\cdot 9=42+18=60.$

For \(p=5\), only \(q_{5,1}=2\) survives, contributing

$2\cdot 2\cdot 8=32.$

For \(p=7\), only \(q_{7,1}=1\) survives, contributing

$2\cdot 1\cdot 9=18.$

Therefore

$S(10)=100+60+32+18=210,$

which matches the small checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. First, they keep all arithmetic modulo \(10^9+7\) and reduce each contribution immediately using

$2q(n-q)=2(nq-q^2).$

Second, they generate all primes up to the cutoff \(B\) with a segmented sieve, which avoids storing a full primality table up to \(10^8\). For every small prime, the implementation repeatedly multiplies by the same prime to enumerate \(p^1,p^2,p^3,\dots\) until the next power would exceed \(n\), and each power contributes one quotient term.

Third, the remaining primes \(p\gt B\) are processed in quotient blocks. The implementation loops over possible values of \(\left\lfloor n/p\right\rfloor\), derives the corresponding interval of primes, counts how many primes lie there, and adds that many copies of the same contribution.

Finally, those interval counts come from a Lehmer-style prime-counting routine backed by a small precomputed sieve and memoization. That makes \(\pi(x)\) queries fast enough that the large-prime phase is tiny compared with any direct enumeration of pairs.

Complexity Analysis

Let \(B=\min(n,10^8)\). The segmented sieve over small primes costs \(O(B\log\log B)\) time and uses memory proportional to the segment size, plus the small base sieve needed to mark composites. Enumerating the powers of each small prime adds one step for each prime power \(p^k\le n\) with \(p\le B\).

The large-prime phase iterates quotients up to \(n/(B+1)\). For the actual Project Euler input \(n=10^{12}\), that is only about \(10^4\) outer iterations. Each iteration performs prime-count queries at interval endpoints, and the memoized Lehmer method keeps those queries practical.

Overall, the method is many orders of magnitude faster than the naive \(O(n^2)\) scan over ordered pairs and is easily feasible for the required value \(10^{12}\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=712
2. \(p\)-adic valuation: Wikipedia - \(p\)-adic valuation
3. Prime-counting function: Wikipedia - Prime-counting function
4. Meissel-Lehmer method: Wikipedia - Meissel-Lehmer algorithm
5. Segmented sieve overview: cp-algorithms - Sieve of Eratosthenes