Problem 712: Exponent Difference
View on Project EulerProject Euler Problem 712 Solution
EulerSolve provides an optimized solution for Project Euler Problem 712, Exponent Difference, with C++, Python, Java, and a step-by-step mathematical explanation.
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....
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
- Problem page: https://projecteuler.net/problem=712
- \(p\)-adic valuation: Wikipedia - \(p\)-adic valuation
- Prime-counting function: Wikipedia - Prime-counting function
- Meissel-Lehmer method: Wikipedia - Meissel-Lehmer algorithm
- Segmented sieve overview: cp-algorithms - Sieve of Eratosthenes
Problem 712 source code
C++
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <vector>
namespace {
using i64 = std::int64_t;
using u64 = std::uint64_t;
constexpr i64 kMod = 1'000'000'007LL;
class PrimeCounter {
public:
PrimeCounter() { sieve(); }
u64 count_primes(const u64 x) { return lehmer(x); }
private:
static constexpr int kSieveLimit = 5'000'000;
std::vector<int> primes_;
std::vector<int> pi_;
std::unordered_map<u64, u64> lehmer_cache_;
void sieve() {
std::vector<bool> is_prime(static_cast<std::size_t>(kSieveLimit + 1), true);
is_prime[0] = false;
is_prime[1] = false;
for (int i = 2; static_cast<i64>(i) * i <= kSieveLimit; ++i) {
if (!is_prime[static_cast<std::size_t>(i)]) {
continue;
}
for (int j = i * i; j <= kSieveLimit; j += i) {
is_prime[static_cast<std::size_t>(j)] = false;
}
}
pi_.assign(static_cast<std::size_t>(kSieveLimit + 1), 0);
for (int i = 2; i <= kSieveLimit; ++i) {
if (is_prime[static_cast<std::size_t>(i)]) {
primes_.push_back(i);
}
pi_[static_cast<std::size_t>(i)] =
pi_[static_cast<std::size_t>(i - 1)] + (is_prime[static_cast<std::size_t>(i)] ? 1 : 0);
}
}
static u64 isqrt(u64 x) {
u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
while ((r + 1) <= x / (r + 1)) {
++r;
}
while (r > x / r) {
--r;
}
return r;
}
static u64 icbrt(u64 x) {
u64 r = static_cast<u64>(std::cbrt(static_cast<long double>(x)));
while ((r + 1) <= x / (r + 1) / (r + 1)) {
++r;
}
while (r > 0 && r > x / r / r) {
--r;
}
return r;
}
static u64 iroot4(u64 x) {
u64 r = static_cast<u64>(std::sqrt(std::sqrt(static_cast<long double>(x))));
while ((r + 1) <= x / (r + 1) / (r + 1) / (r + 1)) {
++r;
}
while (r > 0 && r > x / r / r / r) {
--r;
}
return r;
}
u64 phi(const u64 x, const int s) {
if (s == 0) {
return x;
}
if (s == 1) {
return x - x / 2;
}
if (s == 2) {
return x - x / 2 - x / 3 + x / 6;
}
if (s == 3) {
return x - x / 2 - x / 3 - x / 5 + x / 6 + x / 10 + x / 15 - x / 30;
}
if (x <= kSieveLimit && static_cast<u64>(primes_[static_cast<std::size_t>(s - 1)]) *
static_cast<u64>(primes_[static_cast<std::size_t>(s - 1)]) >
x) {
return static_cast<u64>(pi_[static_cast<std::size_t>(x)] - s + 1);
}
return phi(x, s - 1) - phi(x / static_cast<u64>(primes_[static_cast<std::size_t>(s - 1)]), s - 1);
}
u64 lehmer(const u64 x) {
if (x <= kSieveLimit) {
return static_cast<u64>(pi_[static_cast<std::size_t>(x)]);
}
const auto cached = lehmer_cache_.find(x);
if (cached != lehmer_cache_.end()) {
return cached->second;
}
const u64 a = lehmer(iroot4(x));
const u64 b = lehmer(isqrt(x));
const u64 c = lehmer(icbrt(x));
u64 sum = phi(x, static_cast<int>(a));
sum += (b + a - 2) * (b - a + 1) / 2;
for (u64 i = a + 1; i <= b; ++i) {
const u64 p = static_cast<u64>(primes_[static_cast<std::size_t>(i - 1)]);
const u64 w = x / p;
sum -= lehmer(w);
if (i <= c) {
const u64 bi = lehmer(isqrt(w));
for (u64 j = i; j <= bi; ++j) {
const u64 pj = static_cast<u64>(primes_[static_cast<std::size_t>(j - 1)]);
sum -= lehmer(w / pj) - (j - 1);
}
}
}
lehmer_cache_.emplace(x, sum);
return sum;
}
};
i64 solve(const u64 n) {
PrimeCounter prime_counter;
const u64 y = std::min<u64>(n, 100'000'000ULL);
const i64 n_mod = static_cast<i64>(n % static_cast<u64>(kMod));
i64 ans = 0;
const auto add_q = [&](const u64 q, const u64 count, i64& total) {
const i64 q_mod = static_cast<i64>(q % static_cast<u64>(kMod));
i64 term = static_cast<i64>(((__int128)n_mod * q_mod - (__int128)q_mod * q_mod) % kMod);
if (term < 0) {
term += kMod;
}
const i64 count_mod = static_cast<i64>(count % static_cast<u64>(kMod));
total = static_cast<i64>((total + (__int128)2 * term % kMod * count_mod) % kMod);
};
{
const u64 segment_size = 1'000'000ULL;
const int root = static_cast<int>(std::sqrt(static_cast<long double>(y)));
std::vector<bool> base_mark(static_cast<std::size_t>(root + 1), true);
std::vector<int> base_primes;
for (int i = 2; i <= root; ++i) {
if (!base_mark[static_cast<std::size_t>(i)]) {
continue;
}
base_primes.push_back(i);
if (static_cast<i64>(i) * i <= root) {
for (int j = i * i; j <= root; j += i) {
base_mark[static_cast<std::size_t>(j)] = false;
}
}
}
for (u64 low = 2; low <= y; low += segment_size) {
const u64 high = std::min(y, low + segment_size - 1);
std::vector<bool> is_prime(static_cast<std::size_t>(high - low + 1), true);
for (const int p : base_primes) {
const u64 prime = static_cast<u64>(p);
u64 start = (low + prime - 1) / prime * prime;
const u64 p2 = prime * prime;
if (start < p2) {
start = p2;
}
if (start > high) {
continue;
}
for (u64 v = start; v <= high; v += prime) {
is_prime[static_cast<std::size_t>(v - low)] = false;
}
}
for (u64 p = low; p <= high; ++p) {
if (!is_prime[static_cast<std::size_t>(p - low)] || p < 2) {
continue;
}
u64 power = p;
while (power <= n) {
add_q(n / power, 1, ans);
if (power > n / p) {
break;
}
power *= p;
}
}
}
}
if (y < n) {
const u64 q_max = n / (y + 1);
for (u64 q = 1; q <= q_max; ++q) {
u64 l = n / (q + 1) + 1;
u64 r = n / q;
if (r <= y) {
continue;
}
if (l <= y) {
l = y + 1;
}
if (l > r) {
continue;
}
const u64 count = prime_counter.count_primes(r) - prime_counter.count_primes(l - 1);
if (count != 0) {
add_q(q, count, ans);
}
}
}
return ans;
}
} // namespace
int main() {
assert(solve(10) == 210);
assert(solve(100) == 37'018);
std::cout << solve(1'000'000'000'000ULL) << '\n';
return 0;
}
Python
import math
def solve():
MOD = 1000000007
n = 1000000000000
# Lehmer prime counting
SLIMIT = 5000000
is_p = bytearray(b'\x01') * (SLIMIT + 1); is_p[0] = is_p[1] = 0
for i in range(2, int(SLIMIT**0.5)+1):
if is_p[i]:
for j in range(i*i, SLIMIT+1, i): is_p[j] = 0
primes_list = [i for i in range(2, SLIMIT+1) if is_p[i]]
pi_arr = [0]*(SLIMIT+1)
for i in range(2, SLIMIT+1): pi_arr[i] = pi_arr[i-1] + is_p[i]
cache = {}
def phi_f(x, s):
if s == 0: return x
if s == 1: return x - x//2
if s == 2: return x - x//2 - x//3 + x//6
if s == 3: return x - x//2 - x//3 - x//5 + x//6 + x//10 + x//15 - x//30
if x <= SLIMIT and primes_list[s-1]**2 > x:
return pi_arr[x] - s + 1
key = (x, s)
if key in cache: return cache[key]
r = phi_f(x, s-1) - phi_f(x // primes_list[s-1], s-1)
cache[key] = r; return r
lc = {}
def isqrt(x):
r = int(math.isqrt(x))
while (r+1)*(r+1) <= x: r += 1
while r*r > x: r -= 1
return r
def icbrt(x):
r = int(round(x**(1/3)))
while (r+1)**3 <= x: r += 1
while r > 0 and r**3 > x: r -= 1
return r
def iroot4(x):
r = int(round(x**0.25))
while (r+1)**4 <= x: r += 1
while r > 0 and r**4 > x: r -= 1
return r
def lehmer(x):
if x <= SLIMIT: return pi_arr[x]
if x in lc: return lc[x]
a = lehmer(iroot4(x)); b = lehmer(isqrt(x)); c = lehmer(icbrt(x))
s = phi_f(x, a) + (b+a-2)*(b-a+1)//2
for i in range(a+1, b+1):
p = primes_list[i-1]; w = x // p; s -= lehmer(w)
if i <= c:
bi = lehmer(isqrt(w))
for j in range(i, bi+1):
pj = primes_list[j-1]; s -= lehmer(w//pj) - (j-1)
lc[x] = s; return s
n_mod = n % MOD; ans = 0
y = min(n, 100000000)
# Segmented sieve for primes up to y
root = int(y**0.5); base_primes = [p for p in primes_list if p <= root]
SEG = 1000000
for low in range(2, y+1, SEG):
high = min(y, low+SEG-1); mark = bytearray(b'\x01') * (high-low+1)
for p in base_primes:
start = max(p*p, ((low+p-1)//p)*p)
if start > high: continue
for v in range(start, high+1, p): mark[v-low] = 0
for p in range(low, high+1):
if p < 2 or not mark[p-low]: continue
pw = p
while pw <= n:
q = n // pw; q_mod = q % MOD
term = (n_mod * q_mod - q_mod * q_mod) % MOD
ans = (ans + 2 * term) % MOD
if pw > n // p: break
pw *= p
if y < n:
qmax = n // (y+1)
for q in range(1, qmax+1):
l = n//(q+1)+1; r = n//q
if r <= y: continue
if l <= y: l = y+1
if l > r: continue
cnt = lehmer(r) - lehmer(l-1)
if cnt:
q_mod = q % MOD
term = (n_mod * q_mod - q_mod * q_mod) % MOD
ans = (ans + 2 * term % MOD * (cnt % MOD)) % MOD
return str(ans % MOD)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
public class Euler712 {
static final long kMod = 1000000007L;
static class PrimeCounter {
static final int kSieveLimit = 5000000;
List<Integer> primes = new ArrayList<>();
int[] pi = new int[kSieveLimit + 1];
Map<Long, Long> lehmerCache = new HashMap<>();
PrimeCounter() {
sieve();
}
void sieve() {
byte[] isPrime = new byte[kSieveLimit + 1];
java.util.Arrays.fill(isPrime, (byte) 1);
isPrime[0] = isPrime[1] = 0;
for (int i = 2; (long) i * i <= kSieveLimit; ++i) {
if (isPrime[i] != 0) {
for (int j = i * i; j <= kSieveLimit; j += i) {
isPrime[j] = 0;
}
}
}
int count = 0;
for (int i = 2; i <= kSieveLimit; ++i) {
if (isPrime[i] != 0) {
primes.add(i);
count++;
}
pi[i] = count;
}
}
long isqrt(long x) {
long r = (long) Math.sqrt(x);
while ((r + 1) <= x / (r + 1))
++r;
while (r > 0 && r > x / r)
--r;
return r;
}
long icbrt(long x) {
long r = (long) Math.cbrt(x);
while ((r + 1) <= x / (r + 1) / (r + 1))
++r;
while (r > 0 && r > x / r / r)
--r;
return r;
}
long iroot4(long x) {
long r = (long) Math.sqrt(Math.sqrt(x));
while ((r + 1) <= x / (r + 1) / (r + 1) / (r + 1))
++r;
while (r > 0 && r > x / r / r / r)
--r;
return r;
}
long phi(long x, int s) {
if (s == 0)
return x;
if (s == 1)
return x - x / 2;
if (s == 2)
return x - x / 2 - x / 3 + x / 6;
if (s == 3)
return x - x / 2 - x / 3 - x / 5 + x / 6 + x / 10 + x / 15 - x / 30;
if (x <= kSieveLimit && (long) primes.get(s - 1) * primes.get(s - 1) > x) {
return pi[(int) x] - s + 1;
}
return phi(x, s - 1) - phi(x / primes.get(s - 1), s - 1);
}
long lehmer(long x) {
if (x <= kSieveLimit) {
return pi[(int) x];
}
if (lehmerCache.containsKey(x)) {
return lehmerCache.get(x);
}
long a = lehmer(iroot4(x));
long b = lehmer(isqrt(x));
long c = lehmer(icbrt(x));
long sum = phi(x, (int) a);
sum += (b + a - 2) * (b - a + 1) / 2;
for (long i = a + 1; i <= b; ++i) {
long p = primes.get((int) (i - 1));
long w = x / p;
sum -= lehmer(w);
if (i <= c) {
long bi = lehmer(isqrt(w));
for (long j = i; j <= bi; ++j) {
long pj = primes.get((int) (j - 1));
sum -= lehmer(w / pj) - (j - 1);
}
}
}
lehmerCache.put(x, sum);
return sum;
}
}
public static String solve() {
long n = 1000000000000L;
PrimeCounter primeCounter = new PrimeCounter();
long y = Math.min(n, 100000000L);
long nMod = n % kMod;
long[] ans = { 0 };
java.util.function.BiConsumer<Long, Long> addQ = (q, count) -> {
long qMod = q % kMod;
long term = (nMod * qMod % kMod - qMod * qMod % kMod + kMod) % kMod;
long countMod = count % kMod;
ans[0] = (ans[0] + 2 * term % kMod * countMod) % kMod;
};
long segmentSize = 1000000L;
int root = (int) Math.sqrt(y);
byte[] baseMark = new byte[root + 1];
java.util.Arrays.fill(baseMark, (byte) 1);
List<Integer> basePrimes = new ArrayList<>();
for (int i = 2; i <= root; ++i) {
if (baseMark[i] != 0) {
basePrimes.add(i);
if ((long) i * i <= root) {
for (int j = i * i; j <= root; j += i) {
baseMark[j] = 0;
}
}
}
}
for (long low = 2; low <= y; low += segmentSize) {
long high = Math.min(y, low + segmentSize - 1);
byte[] isPrime = new byte[(int) (high - low + 1)];
java.util.Arrays.fill(isPrime, (byte) 1);
for (int p : basePrimes) {
long prime = p;
long start = (low + prime - 1) / prime * prime;
long p2 = prime * prime;
if (start < p2)
start = p2;
if (start > high)
continue;
for (long v = start; v <= high; v += prime) {
isPrime[(int) (v - low)] = 0;
}
}
for (long p = low; p <= high; ++p) {
if (isPrime[(int) (p - low)] != 0) {
long power = p;
while (power <= n) {
addQ.accept(n / power, 1L);
if (power > n / p)
break;
power *= p;
}
}
}
}
if (y < n) {
long qMax = n / (y + 1);
for (long q = 1; q <= qMax; ++q) {
long l = n / (q + 1) + 1;
long r = n / q;
if (r <= y)
continue;
if (l <= y)
l = y + 1;
if (l > r)
continue;
long count = primeCounter.lehmer(r) - primeCounter.lehmer(l - 1);
if (count != 0) {
addQ.accept(q, count);
}
}
}
return Long.toString(ans[0]);
}
public static void main(String[] args) {
System.out.println(solve());
}
}