Detailed mathematical approach

Problem Summary

For each positive integer \(n\), define

$f(n)=\left|\left\{(x,y)\in\mathbb{Z}_{>0}^2 : x\le y,\ \operatorname{lcm}(x,y)=n\right\}\right|,$

and then define the summatory function

$g(N)=\sum_{n=1}^{N} f(n).$

The Project Euler task asks for \(g(10^{12})\). A direct search over pairs \((x,y)\) is hopeless, so the solution first derives a closed form for \(f(n)\), then transforms the summatory problem into a Möbius-weighted sum involving the ternary divisor function.

Mathematical Approach

Step 1: Count \(f(n)\) from the Prime Factorization of \(n\)

Write

$n=\prod_{p} p^{a_p}.$

If \(\operatorname{lcm}(x,y)=n\), then both \(x\) and \(y\) divide \(n\). For each prime \(p^{a_p}\parallel n\), write the \(p\)-adic exponents of \(x\) and \(y\) as \(\alpha_p\) and \(\beta_p\). They must satisfy

$0\le \alpha_p,\beta_p\le a_p,\qquad \max(\alpha_p,\beta_p)=a_p.$

For one fixed prime power \(p^{a}\), the number of ordered exponent pairs \((\alpha,\beta)\) is

$\underbrace{(a+1)^2}_{0\le \alpha,\beta\le a}-\underbrace{a^2}_{0\le \alpha,\beta\le a-1}=2a+1.$

Different primes act independently, so the total number of ordered pairs \((x,y)\) with \(\operatorname{lcm}(x,y)=n\) is

$\prod_{p\mid n}(2a_p+1)=\tau(n^2),$

because \(n^2=\prod p^{2a_p}\) and therefore \(\tau(n^2)=\prod (2a_p+1)\).

The problem does not count ordered pairs; it counts only pairs with \(x\le y\). Every non-diagonal ordered pair appears twice, once as \((x,y)\) and once as \((y,x)\). The only diagonal pair with least common multiple equal to \(n\) is \((n,n)\). Hence

$\boxed{f(n)=\frac{\tau(n^2)+1}{2}.}$

Step 2: Reduce \(g(N)\) to a Summatory Divisor Problem

Summing the closed form gives

$g(N)=\sum_{n\le N}\frac{\tau(n^2)+1}{2}=\frac{1}{2}\left(N+\sum_{n\le N}\tau(n^2)\right).$

So the core task is to evaluate

$S(N)=\sum_{n\le N}\tau(n^2).$

Once \(S(N)\) is known, the required value is simply

$g(N)=\frac{S(N)+N}{2}.$

Step 3: Möbius Inversion and the Ternary Divisor Function

Let \(d_3(m)\) denote the ternary divisor function, i.e. the number of ordered factorizations \(m=abc\). Its Dirichlet series is

$\sum_{m\ge 1}\frac{d_3(m)}{m^s}=\zeta(s)^3.$

For the square-divisor count we have the classical identity

$\sum_{n\ge 1}\frac{\tau(n^2)}{n^s}=\frac{\zeta(s)^3}{\zeta(2s)}.$

Using

$\frac{1}{\zeta(2s)}=\sum_{d\ge 1}\frac{\mu(d)}{d^{2s}},$

and comparing coefficients, we obtain

$\tau(n^2)=\sum_{d^2m=n}\mu(d)\,d_3(m)=\sum_{d^2\mid n}\mu(d)\,d_3\!\left(\frac{n}{d^2}\right).$

Now sum over \(n\le N\) and interchange the order:

$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$

where

$D_3(X)=\sum_{m\le X} d_3(m)=\left|\left\{(a,b,c)\in\mathbb{Z}_{>0}^3: abc\le X\right\}\right|.$

This identity is the main bridge from the arithmetic function \(\tau(n^2)\) to something the program can compute quickly.

Step 4: Fast Evaluation of \(D_3(X)\)

The helper function summatory_d3 computes \(D_3(X)\). Instead of counting all triples directly, it sorts them by their smallest factor. Let

$c=\left\lfloor X^{1/3}\right\rfloor,\qquad r_z=\left\lfloor\sqrt{\frac{X}{z}}\right\rfloor.$

After separating the cases “all three equal”, “exactly two equal”, and “all distinct”, the ordered triple count becomes

$D_3(X)=3\sum_{z=1}^{c}\left(2\sum_{x=z+1}^{r_z}\left\lfloor\frac{X}{zx}\right\rfloor-r_z^2+\left\lfloor\frac{X}{z^2}\right\rfloor\right)+c^3.$

The final term \(c^3\) is the compact form of the telescoping contribution coming from triples with repeated coordinates. The inner floor-sum

$\sum_{x=x_1}^{x_2}\left\lfloor\frac{n}{x}\right\rfloor$

is exactly what the C++ and Java helper sq_sum_floor evaluates. It uses quotient-difference updates rather than recomputing every division from scratch, which is why the routine stays fast even when called many times.

Step 5: Split the Möbius Sum into Small and Large Parts

A direct evaluation of

$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)$

still runs over too many values of \(d\). The implementation uses the two-scale choice

$I=\left\lfloor N^{1/3}\right\rfloor,\qquad D=\left\lfloor\sqrt{\frac{N}{I}}\right\rfloor.$

For \(d\le D\), the code evaluates the terms directly. For \(d>D\), the quantity \(\left\lfloor N/d^2\right\rfloor\) is smaller than \(I\), so many consecutive \(d\)-values share the same small argument. Introduce the Mertens function

$M(x)=\sum_{n\le x}\mu(n),$

and the first differences

$\Delta D_3(i)=D_3(i)-D_3(i-1)=d_3(i).$

Then the large-\(d\) tail can be regrouped as

$\sum_{d>D}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)=\sum_{i=1}^{I-1}\Delta D_3(i)\left(M\!\left(\left\lfloor\sqrt{\frac{N}{i}}\right\rfloor\right)-M(D)\right).$

This is exactly what the arrays small_diff, mertens_small, and mertens_big implement. The subtraction by M(D) * D3(I-1) in the code is the algebraic device that turns prefix values of \(D_3\) into the needed differences \(\Delta D_3(i)\).

Step 6: Recursive Computation of Large Mertens Values

The program does not sieve \(\mu(d)\) all the way to \(\sqrt N\). It only sieves up to \(D\) and computes larger Mertens values \(M(v)\) on demand for the arguments \(v=\left\lfloor\sqrt{N/i}\right\rfloor\). The recurrence comes from the identity

$\sum_{d\le v}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor=1.$

After separating the range at \(\lfloor\sqrt v\rfloor\), one obtains the form used in the code:

$M(v)=1-v+\lfloor\sqrt v\rfloor\,M(\lfloor\sqrt v\rfloor)-\sum_{d=2}^{\lfloor\sqrt v\rfloor} M\!\left(\left\lfloor\frac{v}{d}\right\rfloor\right)-\sum_{d=2}^{\lfloor\sqrt v\rfloor}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor.$

Because the needed arguments are highly structured, storing them in mertens_big[i] is enough; no full sieve up to \(\sqrt N\) is necessary.

Worked Checkpoints

Take \(n=12=2^2\cdot 3\). Then

$\tau(12^2)=\tau(2^4\cdot 3^2)=(4+1)(2+1)=15,$

so

$f(12)=\frac{15+1}{2}=8.$

The checkpoint values embedded in the C++ program are

$g(10)=29,\qquad g(100)=647,\qquad g(1000)=11751,\qquad g(10^4)=186991,\qquad g(10^6)=37429395,$

and the optimized method matches all of them before attempting \(g(10^{12})\).

How the Code Works

The C++ source is the authoritative implementation. It provides exact integer square-root and cube-root helpers, a fast linear Möbius sieve up to \(D\), the floor-sum helper sq_sum_floor, and the triple-count routine summatory_d3. The main function solve computes the direct small-\(d\) contribution, then adds the grouped tail through Mertens values, finally using

$g(N)=\frac{S(N)+N}{2}.$

The functions brute_f_from_lcm and formula_f_from_factorization are used only for checkpoint verification on small inputs. The Java file is a direct translation of the optimized algorithm. The Python file is intentionally different: it is a lightweight bridge that compiles and runs the C++ solver and then parses the printed answer.

Complexity Analysis

The memory footprint is dominated by the Möbius table up to \(D\) and the grouped Mertens arrays up to \(I\), so it is about \(O(N^{1/3})\). The total runtime is far below naive enumeration and, in the implementation-oriented estimate used for this approach, is roughly \(O(N^{2/3})\) arithmetic with fast floor-sum subroutines. That is what makes \(N=10^{12}\) feasible.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=379
2. Divisor function and generalized divisor functions: Wikipedia — Divisor function
3. Möbius function and inversion: Wikipedia — Möbius inversion formula
4. Dirichlet hyperbola method: Wikipedia — Dirichlet hyperbola method
5. Mertens function: Wikipedia — Mertens function