Problem 379: Least Common Multiple Count

Project Euler Problem 379 Solution

EulerSolve provides an optimized solution for Project Euler Problem 379, Least Common Multiple Count, with C++, Python, Java, and a step-by-step mathematical explanation.

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$....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, let $f(n)$ be the number of couples $(x,y)$ of positive integers with $x\le y$ whose least common multiple is $n$, and let $g$ be the summatory function of $f$,

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

The statement gives $g(10^6)=37429395$ and asks for $g(10^{12})$. ^[1]

The whole solution in one sentence: the couples with least common multiple $n$ can be counted from the prime exponents of $n$ alone, which turns $g(N)$ into one sum of divisor counts, and that sum is then assembled from a few thousand counts of triples of positive integers whose product stays below a bound.

Mathematical Approach

1. The first few values

If the least common multiple of $x$ and $y$ is $n$, then $x$ and $y$ both divide $n$, so the search never leaves the divisors of $n$ and the small values can be written out by hand.

$n=1$: only $(1,1)$, so $f(1)=1$.
$n=2$: $(1,2)$ and $(2,2)$, so $f(2)=2$.
$n=4$: $(1,4)$, $(2,4)$, $(4,4)$, so $f(4)=3$.
$n=6$: $(1,6)$, $(2,3)$, $(2,6)$, $(3,6)$, $(6,6)$, so $f(6)=5$.
$n=12$: $(1,12)$, $(2,12)$, $(3,4)$, $(3,12)$, $(4,6)$, $(4,12)$, $(6,12)$, $(12,12)$, so $f(12)=8$.

Continuing in this way, $f(1),\dots,f(15)$ equals

$$1,2,2,3,2,5,2,4,3,5,2,8,2,5,5,$$

and the running totals $g(1),\dots,g(15)$ are

$$\begin{aligned}&1,3,5,8,10,15,17,21,\\&24,29,31,39,41,46,51.\end{aligned}$$

In particular $g(10)=29$, one of the checkpoints the program uses. The value $2$ appears at every prime and the large values at numbers with several prime factors, which suggests that $f(n)$ depends only on the exponents in the prime factorisation of $n$.

2. One prime at a time

Write $n=\prod_p p^{a_p}$. Since $x$ and $y$ divide $n$, we may write $x=\prod_p p^{\alpha_p}$ and $y=\prod_p p^{\beta_p}$ with $0\le\alpha_p,\beta_p\le a_p$. The least common multiple takes the larger exponent at each prime, so the condition becomes

$$\max(\alpha_p,\beta_p)=a_p$$

for every prime $p$, and the exponents at different primes may be chosen independently. For a single prime power $p^a$, the number of ordered pairs $(\alpha,\beta)$ with $0\le\alpha,\beta\le a$ and $\max(\alpha,\beta)=a$ is the full square minus the smaller square that avoids the exponent $a$,

$$(a+1)^2-a^2=2a+1.$$

Multiplying over the primes, the number of ordered pairs with least common multiple $n$ is $\prod_p(2a_p+1)$. The number of divisors of a number is the product of its exponents each increased by one, ^[3] and $n^2=\prod_p p^{2a_p}$, so that product is exactly $d(n^2)$:

$$\#\{(x,y):\mathrm{lcm}(x,y)=n\}=d(n^2).$$

For $n=12=2^2\cdot3$ this gives $(2\cdot2+1)(2\cdot1+1)=15$, and indeed $d(144)=15$. ^[2]

3. From ordered pairs to couples

Swapping the two entries permutes the ordered solutions and fixes exactly those with $x=y$. Since the least common multiple of $x$ with itself is $x$, the only fixed solution is $(n,n)$; all others fall into pairs with $x\lt y$ in one member. Hence

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

The halving is exact, because $d(n^2)=\prod_p(2a_p+1)$ is a product of odd numbers and therefore odd. For $n=12$ we get $f(12)=(15+1)/2=8$, the eight couples of Section 1.

The same formula counts something else. For positive integers $u\le v$ the equation $1/n=1/u+1/v$ forces $u,v\gt n$ and rearranges to $(u-n)(v-n)=n^2$, and adding $n$ to both factors of a factorisation of $n^2$ inverts that, so there are again $(d(n^2)+1)/2$ of them. ^[4] For $n=6$ the five decompositions

$$\begin{aligned}&\tfrac17+\tfrac1{42},\quad\tfrac18+\tfrac1{24},\\&\tfrac19+\tfrac1{18},\quad\tfrac1{10}+\tfrac1{15},\\&\tfrac1{12}+\tfrac1{12}\end{aligned}$$

match $f(6)=5$. Summing the closed form over $n\le N$ reduces the task to a single divisor sum,

$$g(N)=\frac{N+D(N)}{2},$$

where

$$D(N)=\sum_{n\le N}d(n^2).$$

4. Removing the square part

$D(N)$ still asks for a factorisation of every $n\le10^{12}$, so we need an identity that rewrites $d(n^2)$ through a function whose summatory values can be counted geometrically. Let $d_3(m)$ be the number of ordered triples $(a,b,c)$ of positive integers with $abc=m$; it is the three-factor member of the family of generalised divisor functions. ^[7] For real $s\gt1$ put $\zeta(s)=\sum_{a\ge1}a^{-s}$; all series below converge absolutely there, so they may be multiplied and regrouped by the index. Expanding the cube of $\zeta(s)$ and grouping each triple by its product shows

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

For $d(n^2)$, unique factorisation lets the exponent at each prime be chosen independently, so multiplying the series of the single prime powers gives the series of all integers — a product over primes called an Euler product. With $x=p^{-s}$, the geometric sum $\sum_{a\ge0}(2a+1)x^a=(1+x)/(1-x)^2$ and the factorisation $1+x=(1-x^2)/(1-x)$ give

$$\sum_{n\ge1}\frac{d(n^2)}{n^s}=\prod_p\frac{1-x^2}{(1-x)^3}=\frac{\zeta(s)^3}{\zeta(2s)}.$$

This is the Dirichlet series recorded for the sequence $d(n^2)$. ^[2] An integer is squarefree if no prime square divides it. The Möbius function is $\mu(1)=1$, equals $(-1)^r$ on a product of $r$ distinct primes and $0$ on all other integers, so $\mu(2)=-1$, $\mu(4)=0$ and $\mu(6)=1$; its Dirichlet series is $1/\zeta(s)$. ^[5] Hence $1/\zeta(2s)=\sum_{d\ge1}\mu(d)d^{-2s}$, a contribution to $n^{-s}$ has $n=d^2m$, and comparing coefficients yields

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

The identity can also be checked without series. Both sides are multiplicative, since $\mu$ and $d_3$ multiply over coprime factors, so it is enough to compare them at a prime power $p^a$, where the left side is $2a+1$. If a triple of exponents summing to $a$ starts with $j$, the remaining $a-j$ can be split in $a-j+1$ ways, so $d_3(p^a)=(a+1)(a+2)/2$. For $a\le1$ only $d=1$ contributes, giving $1$ and $3$ as required; for $a\ge2$ only $d=1$ and $d=p$ do, and

$$\begin{aligned}&d_3(p^a)-d_3(p^{a-2})\\&=\tfrac{(a+1)(a+2)}{2}-\tfrac{(a-1)a}{2}\\&=2a+1 .\end{aligned}$$

For $n=12$ the sum reads $d_3(12)-d_3(3)=18-3$, which is $d(144)=15$.

Now sum over $n\le N$ and group by $d$. For a fixed $d$ the admissible $n$ are $n=d^2m$ with $d^2m\le N$, that is $m\le\lfloor N/d^2\rfloor$, so

$$D(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)$ counts the triples with $abc\le x$. Section 8 evaluates this sum at $N=10$ and recovers the hand count $g(10)=29$.

5. Counting triples under a hyperbola

Only $D_3$ is left. For two factors the count is classical: the pairs with $ab\le m$ lie under a hyperbola, the region is symmetric, and every pair has $\min(a,b)\le\sqrt m$, so counting both halves and removing the doubly counted square gives the Dirichlet hyperbola method ^[6]

$$D_2(m)=2\sum_{a\le\sqrt m}\left\lfloor\frac{m}{a}\right\rfloor-\lfloor\sqrt m\rfloor^2 ,$$

which needs only $O(\sqrt m)$ divisions. For three factors we classify a triple by its smallest entry $v$. From $v^3\le abc\le x$ we get $v\le c=\lfloor x^{1/3}\rfloor$, and for fixed $v$ there are three shapes: all three entries equal $v$, one triple; exactly two equal to $v$, the third being one of the $\lfloor x/v^2\rfloor-v$ values above $v$ in any of three positions; and exactly one entry equal to $v$, the other two forming a pair counted by

$$\begin{aligned}P(v)=\#\{(a,b):{}&ab\le m_v,\\&\min(a,b)\gt v\},\end{aligned}$$

with $m_v=\lfloor x/v\rfloor$, again in three positions. Adding the three shapes over all $v\le c$,

$$\begin{aligned}D_3(x)&=3\sum_{v=1}^{c}\left(P(v)+\left\lfloor\frac{x}{v^2}\right\rfloor-v\right)\\&\quad+c .\end{aligned}$$

Here $\lfloor m_v/v\rfloor=\lfloor x/v^2\rfloor$ by the nested-division rule for floors, ^[9] and $\lfloor x/v^2\rfloor\ge v$ because $v^3\le x$. To evaluate $P(v)$ we remove from $D_2(m_v)$ the strips $a\le v$ and $b\le v$ and add back their intersection, the full square of side $v$ because $v^2\le m_v$:

$$\begin{aligned}P(v)&=D_2(m_v)-2\sum_{a\le v}\left\lfloor\frac{m_v}{a}\right\rfloor+v^2\\&=2\!\!\sum_{a=v+1}^{\lfloor\sqrt{m_v}\rfloor}\!\!\left\lfloor\frac{m_v}{a}\right\rfloor-\lfloor\sqrt{m_v}\rfloor^2+v^2 .\end{aligned}$$

Substituting this expression for $P(v)$ leaves the extra term $v^2-v$ in every summand. The published code omits that term from the bracket and replaces the final $c$ by $c^3$; the two changes cancel, because

$$\begin{aligned}&3\sum_{v=1}^{c}v(v-1)+c\\&=c(c+1)(c-1)+c=c^3 .\end{aligned}$$

A worked example: $x=10$, so $c=2$. For $v=1$ we have $m_1=10$ and $\lfloor\sqrt{10}\rfloor=3$, hence $P(1)=2(5+3)-9+1=8$ — the eight pairs $(2,2)$, $(2,3)$, $(2,4)$, $(2,5)$, $(3,2)$, $(3,3)$, $(4,2)$, $(5,2)$ — so $v=1$ contributes $3(8+10-1)+1=52$. For $v=2$ we have $m_2=5$, $\lfloor\sqrt5\rfloor=2$, an empty inner sum, $P(2)=0-4+4=0$ and $\lfloor10/4\rfloor=2$, so $v=2$ contributes $3(0+2-2)+1=1$. The total is $D_3(10)=53$, which matches the sum of the first ten values of $d_3$,

$$\begin{aligned}&1+3+3+6+3\\&\quad+9+3+10+6+9=53 .\end{aligned}$$

6. Two scales for the Möbius sum

Summing Section 4 directly would run $d$ up to $\lfloor\sqrt N\rfloor=10^6$ and need $\mu(d)$ on that whole range. But $\lfloor N/d^2\rfloor$ soon becomes small and then repeats, so the tail can be grouped. Fix the two cut-offs the code uses,

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

which for $N=10^{12}$ are $I=D=10^4$. For $d\le D$ we have $\lfloor N/d^2\rfloor\ge I$ and the terms are evaluated one by one; only the $6083$ squarefree values of $d$ contribute. For $d\gt D$ the definition of $D$ gives $(D+1)^2\gt N/I$ and hence $N/d^2\lt I$, so the argument is one of $1,\dots,I-1$. Writing $v_i=\lfloor\sqrt{N/i}\rfloor$, the condition $\lfloor N/d^2\rfloor\ge i$ is equivalent to $d\le v_i$, so the $d$ with $\lfloor N/d^2\rfloor=i$ are exactly those in $v_{i+1}\lt d\le v_i$ and their Möbius weights add up to $M(v_i)-M(v_{i+1})$, where

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

is the Mertens function. ^[8] These intervals partition $D\lt d\le\lfloor\sqrt N\rfloor$, since $v_1=\lfloor\sqrt N\rfloor$ and $v_I=D$. The tail is therefore the sum of $D_3(i)\bigl(M(v_i)-M(v_{i+1})\bigr)$ over $1\le i\lt I$. Expanding the brackets and collecting the coefficient of each $M(v_i)$ turns it into $D_3(i)-D_3(i-1)=d_3(i)$, where $D_3(0)=0$, and leaves one uncancelled term $-D_3(I-1)M(v_I)$. This reindexing, called summation by parts, gives

$$\sum_{i=1}^{I-1}d_3(i)\,M(v_i)\;-\;D_3(I-1)\,M(D).$$

Almost a million Möbius terms have become $9999$ Mertens values and the first differences $d_3(i)=D_3(i)-D_3(i-1)$, which are free once the prefix values $D_3(i)$ are known for $i\lt I$.

7. Mertens values on demand

The arguments $v_i$ reach $10^6$, so sieving $\mu$ that far is just what we wanted to avoid; each $M(v)$ is computed from smaller values instead. Since the divisor sums of $\mu$ vanish except at $1$, ^[5] counting the pairs $(k,e)$ with $ke\le v$ by each coordinate in turn gives two identities

$$\sum_{k\le v}\mu(k)\left\lfloor\frac{v}{k}\right\rfloor=1=\sum_{e\le v}M\!\left(\left\lfloor\frac{v}{e}\right\rfloor\right),$$

the second of which is already a recursion for the Mertens function. ^[8] A faster one follows from splitting the same region at $K=\lfloor\sqrt v\rfloor$ with the pair $(k,e)$ weighted by $\mu(k)$: every pair with $ke\le v$ has $k\le K$ or $e\le K$, the first strip carries $\sum_{k\le K}\mu(k)\lfloor v/k\rfloor$, the second $\sum_{e\le K}M(\lfloor v/e\rfloor)$, and their overlap, the square $k,e\le K$, carries $K\,M(K)$, as every $k\le K$ appears there with $K$ partners. Strips minus overlap is the total weight $1$, and the terms at index $1$ are $v$ and $M(v)$, so isolating $M(v)$ gives ^[6]

$$\begin{aligned}M(v)&=1-v+K\,M(K)\\&\quad-\sum_{d=2}^{K}\mu(d)\lfloor v/d\rfloor\\&\quad-\sum_{d=2}^{K}M(\lfloor v/d\rfloor).\end{aligned}$$

For $v=10$ we have $K=3$, $M(3)=-1$ and $M(5)=-2$, so the right-hand side is

$$\begin{aligned}&1-10-3-(-5-2)\\&\quad-(-3-1)=-1 ,\end{aligned}$$

which agrees with $\mu(1)+\dots+\mu(10)=-1$.

Each recursive argument $\lfloor v/d\rfloor$ is either at most $D$, where a sieved table of $M$ is available, or it is again one of the $v_j$. Indeed $v_i=\lfloor\sqrt{N/i}\rfloor$, and the nested-division rule ^[9] together with $\lfloor\sqrt{\lfloor y\rfloor}\rfloor=\lfloor\sqrt y\rfloor$ gives

$$\left\lfloor\frac{v_i}{d}\right\rfloor=\left\lfloor\sqrt{\frac{N}{id^2}}\,\right\rfloor=v_{id^2}.$$

The square-root rule holds because $k=\lfloor\sqrt y\rfloor$ satisfies $k^2\le\lfloor y\rfloor\le y\lt(k+1)^2$, the first inequality because $k^2$ is an integer. Moreover $\lfloor v_i/d\rfloor\gt D$ forces $(D+1)^2\le N/(id^2)$, hence $id^2\lt I$: the index is in range and larger than $i$, so a loop running $i$ downwards always finds the value it needs already stored, with no dictionary of visited arguments.

8. The whole method on a small case

Take $N=10$, where $I=D=2$. The direct part runs over $d\le2$ and contributes $D_3(10)-D_3(2)=53-4=49$. The tail has a single index, $i=1$: its first difference is $d_3(1)=1$ and $v_1=\lfloor\sqrt{10}\rfloor=3$, so it contributes $M(3)=-1$, while the correction term $D_3(I-1)\,M(D)$ vanishes, because $D_3(1)=1$ and $M(2)=0$. The tail is therefore $-1$, which is right: the only $d$ in $2\lt d\le3$ is $d=3$, with $\mu(3)D_3(\lfloor10/9\rfloor)=-1$. Altogether $D(10)=49-1=48$ and

$$g(10)=\frac{10+48}{2}=29 ,$$

the value obtained by listing couples in Section 1.

How the Code Works

The C++ program ^[10] reads its options first: parse_arguments accepts --skip-checkpoints and --n= followed by digits, which parse_u64_after_prefix rejects on a non-digit or on overflow. Unless that flag is given, run_checkpoints runs first and main stops with a message on standard error if a check fails. Then solve is called with the chosen bound, whose default is $10^{12}$. The method needs a positive bound and the parser does not enforce that: --n=0 is accepted and makes solve divide by zero when it forms n / I.

Five helpers carry the mathematics. Both root helpers start from a floating-point estimate and correct it with integer comparisons: isqrt_u64 compares $r\gt x/r$ instead of forming $r^2$, while icbrt_u64 compares cubes in unsigned __int128, the only use of a 128-bit type in the program. sq_sum_floor(n,x_1,x_2) evaluates $\sum_{x=x_1}^{x_2}\lfloor n/x\rfloor$ walking $x$ downwards. At the start of a step beta is $\lfloor n/(x+1)\rfloor$ with remainder epsilon, delta predicts how much the quotient grows when the denominator drops to $x$, and gamma holds $\beta-x\delta$, so that epsilon + gamma is exactly the remainder $n-x(\beta+\delta)$ of the predicted quotient. Correcting it into $0\le r\lt x$ therefore produces $\lfloor n/x\rfloor$, normally with additions and comparisons alone; when more than two corrections are needed the loop breaks into a second one that recovers the jump with a single division epsilon / x. A third loop would fall back to plain division, but the second has no exit of its own and always reaches the end of the range, so it is unreachable. In the evaluation of $D_3(10^{12})$ the first loop takes $147\,281\,575$ of the $148\,534\,687$ steps, the second the other $1\,253\,112$. summatory_d3 is Section 5 as written there, with the bracket accumulated over $z\le$ icbrt, then tripled, then $c^3$ added. mobius_sieve is a linear sieve that stores the least prime factor of every integer up to $D$ and fills $\mu$ from it.

solve follows Sections 4, 6 and 7 in order. It sets $I$ and $D$, sieves $\mu$ up to $D$, and in one pass adds $\mu(d)\,D_3(\lfloor n/d^2\rfloor)$ for every $d\le D$ while building mertens_small as the prefix sums of $\mu$. If $I\ge2$, it fills small_diff[i] with $D_3(i)$ for $i\lt I$, subtracts mertens_small[D] times small_diff[I-1], then differences small_diff in place from the top down, leaving $d_3(i)$. The final loop runs $i$ from $I-1$ down to $1$, evaluates the recursion of Section 7 for $v=$isqrt(n/i) — reading mertens_small[q] when $q=\lfloor v/d\rfloor\le D$ and the cached mertens_big[i*d*d] otherwise — stores it in mertens_big[i] and adds small_diff[i] times it. Finally $n$ is added and the total halved.

The self-checks are independent of all this machinery:

For every $n\le200$, the double loop brute_f_from_lcm over $x\le y\le n$ must agree with formula_f_from_factorization, which trial divides and multiplies the factors $2a+1$.
summatory_g_via_spf sums $f(n)$ over a smallest-prime-factor sieve and must give $g(10)=29$, $g(100)=647$, $g(1000)=11751$ and $g(10^4)=186991$.
solve itself must reproduce the published checkpoint $g(10^6)=37429395$.

The largest quantity handled is $D_3(10^{12})$, about $4.0\cdot10^{14}$, far inside the signed 64-bit range of roughly $9.2\cdot10^{18}$, so std::int64_t accumulators suffice for the sums. The Java program is the same algorithm in long arithmetic, with isqrt and icbrt correcting Math.sqrt and Math.cbrt; it hardcodes $N=10^{12}$ and has neither options nor checkpoints, so it reproduces only the main computation. The Python file is no independent implementation: it locates the C++ source relative to its own path, compiles it with clang++ or g++ at -std=c++17 -O2 into a cached binary, recompiles it only when it is newer, runs it and reads the value with parse_output. Since it passes --skip-checkpoints whenever that string occurs in the source, the Python route skips the self-checks above.

Complexity and Verification

One evaluation of $D_3(x)$ performs $\lfloor\sqrt{x/v}\rfloor-v$ steps for each $v\le c=\lfloor x^{1/3}\rfloor$. Since $\sum_{v\le c}\sqrt{x/v}\approx2\sqrt{xc}=2x^{2/3}$ and $\sum_{v\le c}v\approx c^2/2=x^{2/3}/2$, the cost is about $1.5\,x^{2/3}$ steps; for $x=10^{12}$ that estimate is $1.5\cdot10^8$ against the $148\,534\,687$ steps actually taken. Over the small-$d$ part the costs behave like

$$1.5\,N^{2/3}\sum_{d\ge1}\frac{\mu(d)^2}{d^{4/3}}=1.5\,N^{2/3}\frac{\zeta(4/3)}{\zeta(8/3)},$$

and the series converges to about $2.80$, so the large evaluations cost $O(N^{2/3})$ arithmetic operations. The other phases are cheaper: the separate values $D_3(i)$ for $i\lt I$ cost $O(I^{5/3})=O(N^{5/9})$; the inner Mertens loop at index $i$ runs at most $(N/i)^{1/4}$ times, which sums to $O(N^{1/4}I^{3/4})=O(N^{1/2})$; the sieve and the array passes are $O(N^{1/3})$. As the single term $d=1$ already needs $\Theta(N^{2/3})$ steps, the whole method is $\Theta(N^{2/3})$ arithmetic operations. Counting inner-loop iterations only, the computation at $N=10^{12}$ uses $398\,020\,178$ floor-sum steps in the $6083$ large evaluations of $D_3$ — $148\,534\,687$ of them in the single call with $d=1$ — then $3\,126\,877$ in the $9999$ small evaluations and $1\,317\,502$ in the inner Mertens loop, $402\,464\,557$ in all, and both the C++ and the Java program finish in well under a second. Memory is $O(N^{1/3})$: the arrays of length about $10^4$ — the sieved $\mu$ values and least prime factors, the two tables of Mertens values and the table of $D_3$ differences — together need well under a megabyte.

Beyond the three self-checks above, the following were verified independently for this article.

A brute-force loop over all couples agrees with $(d(n^2)+1)/2$ for every $n\le300$, as does the number of decompositions $1/n=1/u+1/v$ with $u\le v$.
The identity of Section 4 holds for every $n\lt600$, and $D(N)=\sum_{d\le\sqrt N}\mu(d)D_3(\lfloor N/d^2\rfloor)$ matches the direct sum of $d(n^2)$ at $N=10,12,100,1000,10^4$.
summatory_d3 agrees with a direct count of the triples $abc\le n$ for every $n\lt400$ and at $n=1000,5000,12345$; its value at $10^{12}$ also agrees with an independent routine that groups equal quotients instead of using incremental differences.
The incremental sq_sum_floor agrees with plain division sums on all $19068$ ranges that summatory_d3 uses for the bounds $n\le2000$ and $n=5000,12345,10^6,10^9+7$, and its state invariants for beta, epsilon and gamma hold at every step of those runs.
A Python port of the whole method reproduces the direct sum of $f(n)$ for the twelve bounds $N=2,3,5,10,12,50$ and $100,999,1000,5000,10^4,10^5$, and gives $37429395$ at $N=10^6$.
At $N=10^{12}$, the cache rule $\lfloor v_i/d\rfloor=v_{id^2}$ and the bound $id^2\le I-1$ hold for every index pair the recursion of Section 7 uses, and the C++ and Java programs produce identical output.

Footnotes and References

The references below supply the background facts used above. The closed form for $f(n)$, the triple classification for $D_3$, the two-scale split of the Möbius sum and the operation counts are derived explicitly in this article.

Project Euler 379 — Least Common Multiple Count. The official statement defines f(n) as the number of couples (x, y) of positive integers with x ≤ y and least common multiple n, defines g as the summatory function of f, gives g(10⁶) = 37429395 and asks for g(10¹²).
OEIS A048691 — d(n²). The entry states that the sequence counts the pairs (x, y) with lcm(x, y) = n, that it is multiplicative with value 2e + 1 at a prime power pᵉ, and records the Dirichlet generating function ζ(s)³/ζ(2s) used in Section 4.
Wikipedia — Divisor function. The article states that the divisor function is multiplicative, that the number of divisors of a number equals the product of its prime exponents each increased by one, and that its Dirichlet series is ζ(s)².
OEIS A018892 — two unit fractions. The entry counts the ways of writing 1/n as a sum of exactly two unit fractions, gives the formula (τ(n²) + 1)/2 and records the comment that this is equally the number of pairs (x, y) with lcm(x, y) = n; its first terms are the values f(n) of Section 1.
Wikipedia — Möbius function. The article defines μ, proves that it is multiplicative, states that the sum of μ(d) over the divisors of n is 1 for n = 1 and 0 otherwise, and gives the Dirichlet series 1/ζ(s) = Σ μ(n)n⁻ˢ, used at 2s in Section 4.
Wikipedia — Dirichlet hyperbola method. The article splits a sum over the lattice points under a hyperbola into two overlapping strips minus their common rectangle, derives D(n) = 2 Σ ⌊n/x⌋ − ⌊√n⌋² in O(√n) operations, and notes that Dirichlet introduced the technique in 1849.
Wikipedia — Divisor summatory function. The article introduces the generalised divisor function d_k(n), which counts the ways of writing n as a product of k numbers, and its summatory function D_k(x); the case k = 3 supplies the notation used in Sections 4 and 5.
Wikipedia — Mertens function. The article defines M(x) as the sum of μ(k) for k ≤ x, gives the identity that the sum of M(⌊n/d⌋) over d ≤ n equals 1, and states that sieving all values up to x costs O(x log log x), which is what the recursion of Section 7 avoids.
Wikipedia — Floor and ceiling functions. Its section on nested divisions states ⌊⌊x/m⌋/n⌋ = ⌊x/(mn)⌋ for a positive integer n and real m and x, the rule that lets the code replace repeated floor divisions by a single one in Sections 5 and 7.
C++ — Euler379.cpp. The linked immutable revision contains sq_sum_floor, summatory_d3, mobius_sieve, solve and the checks run by run_checkpoints. It is the implementation record for this article.

Problem 379 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;

constexpr u64 kDefaultN = 1'000'000'000'000ULL;
constexpr u64 kCheckpointN = 1'000'000ULL;
constexpr u64 kCheckpointExpected = 37'429'395ULL;

struct Options {
    u64 n = kDefaultN;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--n=", parsed_u64)) {
            options.n = parsed_u64;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return true;
}

u64 isqrt_u64(const u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1ULL) > 0ULL && (r + 1ULL) <= x / (r + 1ULL)) {
        ++r;
    }
    while (r > x / r) {
        --r;
    }
    return r;
}

u64 icbrt_u64(const u64 x) {
    u64 r = static_cast<u64>(std::cbrt(static_cast<long double>(x)));
    while ((r + 1ULL) > 0ULL &&
           static_cast<unsigned __int128>(r + 1ULL) * static_cast<unsigned __int128>(r + 1ULL) *
                   static_cast<unsigned __int128>(r + 1ULL) <=
               static_cast<unsigned __int128>(x)) {
        ++r;
    }
    while (static_cast<unsigned __int128>(r) * static_cast<unsigned __int128>(r) *
               static_cast<unsigned __int128>(r) >
           static_cast<unsigned __int128>(x)) {
        --r;
    }
    return r;
}

i64 sq_sum_floor(const i64 n, const i64 x1, const i64 x2) {
    if (x1 > x2) {
        return 0;
    }

    i64 x = x2;
    i64 s = 0;
    i64 beta = n / (x + 1);
    i64 epsilon = n % (x + 1);
    i64 delta = (n / x) - beta;
    i64 gamma = beta - x * delta;

    while (x >= x1) {
        epsilon += gamma;
        if (epsilon >= x) {
            delta += 1;
            gamma -= x;
            epsilon -= x;
            if (epsilon >= x) {
                delta += 1;
                gamma -= x;
                epsilon -= x;
                if (epsilon >= x) {
                    break;
                }
            }
        } else if (epsilon < 0) {
            delta -= 1;
            gamma += x;
            epsilon += x;
        }

        gamma += 2 * delta;
        beta += delta;
        s += beta;
        --x;
    }

    if (x < x1) {
        return s;
    }

    epsilon = n % (x + 1);
    delta = (n / x) - beta;
    gamma = beta - x * delta;

    while (x >= x1) {
        epsilon += gamma;
        const i64 delta2 = epsilon / x;
        delta += delta2;
        epsilon -= x * delta2;
        gamma += 2 * delta - x * delta2;
        beta += delta;
        s += beta;
        --x;
    }

    while (x >= x1) {
        s += n / x;
        --x;
    }

    return s;
}

i64 summatory_d3(const i64 n) {
    const i64 cbrt_n = static_cast<i64>(icbrt_u64(static_cast<u64>(n)));
    i64 acc = 0;

    for (i64 z = 1; z <= cbrt_n; ++z) {
        const i64 nz = n / z;
        const i64 sqrt_nz = static_cast<i64>(isqrt_u64(static_cast<u64>(nz)));
        acc += 2 * sq_sum_floor(nz, z + 1, sqrt_nz) - sqrt_nz * sqrt_nz + nz / z;
    }

    acc *= 3;
    acc += cbrt_n * cbrt_n * cbrt_n;
    return acc;
}

std::vector<int> mobius_sieve(const int n) {
    std::vector<int> mu(static_cast<std::size_t>(n) + 1ULL, 0);
    std::vector<int> lp(static_cast<std::size_t>(n) + 1ULL, 0);
    std::vector<int> primes;

    mu[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (lp[static_cast<std::size_t>(i)] == 0) {
            lp[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
            mu[static_cast<std::size_t>(i)] = -1;
        }
        for (const int p : primes) {
            const u64 v = static_cast<u64>(i) * static_cast<u64>(p);
            if (v > static_cast<u64>(n)) {
                break;
            }
            lp[static_cast<std::size_t>(v)] = p;
            if (i % p == 0) {
                mu[static_cast<std::size_t>(v)] = 0;
                break;
            }
            mu[static_cast<std::size_t>(v)] = -mu[static_cast<std::size_t>(i)];
        }
    }

    return mu;
}

u64 gcd_u64(u64 a, u64 b) {
    while (b != 0ULL) {
        const u64 t = a % b;
        a = b;
        b = t;
    }
    return a;
}

u64 brute_f_from_lcm(const u64 n) {
    u64 count = 0ULL;
    for (u64 x = 1ULL; x <= n; ++x) {
        for (u64 y = x; y <= n; ++y) {
            const u64 g = gcd_u64(x, y);
            const u64 l = (x / g) * y;
            if (l == n) {
                ++count;
            }
        }
    }
    return count;
}

u64 formula_f_from_factorization(u64 n) {
    u64 divisor_count_n2 = 1ULL;
    for (u64 p = 2ULL; p * p <= n; ++p) {
        if (n % p != 0ULL) {
            continue;
        }

        int exponent = 0;
        while (n % p == 0ULL) {
            n /= p;
            ++exponent;
        }
        divisor_count_n2 *= static_cast<u64>(2 * exponent + 1);
    }
    if (n > 1ULL) {
        divisor_count_n2 *= 3ULL;
    }
    return (divisor_count_n2 + 1ULL) / 2ULL;
}

u64 summatory_g_via_spf(const int n) {
    std::vector<int> spf(static_cast<std::size_t>(n) + 1ULL, 0);
    for (int i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] != 0) {
            continue;
        }
        for (u64 j = static_cast<u64>(i); j <= static_cast<u64>(n); j += static_cast<u64>(i)) {
            if (spf[static_cast<std::size_t>(j)] == 0) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }

    u64 total = 0ULL;
    for (int x = 1; x <= n; ++x) {
        int value = x;
        u64 divisor_count_n2 = 1ULL;
        while (value > 1) {
            const int p = spf[static_cast<std::size_t>(value)];
            int exponent = 0;
            while (value % p == 0) {
                value /= p;
                ++exponent;
            }
            divisor_count_n2 *= static_cast<u64>(2 * exponent + 1);
        }
        total += (divisor_count_n2 + 1ULL) / 2ULL;
    }
    return total;
}

u64 solve(const u64 n) {
    const i64 I = static_cast<i64>(icbrt_u64(n));
    const i64 D = static_cast<i64>(isqrt_u64(n / static_cast<u64>(I)));

    i64 res = 0;
    std::vector<int> mu = mobius_sieve(static_cast<int>(D));
    std::vector<i64> mertens_small(static_cast<std::size_t>(D) + 1ULL, 0);

    for (i64 d = 1; d <= D; ++d) {
        const int mud = mu[static_cast<std::size_t>(d)];
        if (mud != 0) {
            const u64 div = static_cast<u64>(d) * static_cast<u64>(d);
            res += static_cast<i64>(mud) * summatory_d3(static_cast<i64>(n / div));
        }
        mertens_small[static_cast<std::size_t>(d)] =
            mertens_small[static_cast<std::size_t>(d - 1)] + static_cast<i64>(mud);
    }

    if (I >= 2) {
        std::vector<i64> small_diff(static_cast<std::size_t>(I), 0);
        for (i64 i = 1; i < I; ++i) {
            small_diff[static_cast<std::size_t>(i)] = summatory_d3(i);
        }

        res -= mertens_small[static_cast<std::size_t>(D)] *
               small_diff[static_cast<std::size_t>(I - 1)];

        for (i64 i = I - 1; i >= 2; --i) {
            small_diff[static_cast<std::size_t>(i)] -= small_diff[static_cast<std::size_t>(i - 1)];
        }

        std::vector<i64> mertens_big(static_cast<std::size_t>(I), 0);
        for (i64 i = I - 1; i >= 1; --i) {
            const u64 v = isqrt_u64(n / static_cast<u64>(i));
            const u64 v_sqrt = isqrt_u64(v);

            i64 m = 1 - static_cast<i64>(v) +
                    static_cast<i64>(v_sqrt) * mertens_small[static_cast<std::size_t>(v_sqrt)];

            for (u64 d = 2; d <= v_sqrt; ++d) {
                const u64 q = v / d;
                if (q <= static_cast<u64>(D)) {
                    m -= mertens_small[static_cast<std::size_t>(q)];
                } else {
                    const u64 idx = static_cast<u64>(i) * d * d;
                    m -= mertens_big[static_cast<std::size_t>(idx)];
                }
                m -= (mertens_small[static_cast<std::size_t>(d)] -
                      mertens_small[static_cast<std::size_t>(d - 1)]) *
                     static_cast<i64>(q);
            }

            mertens_big[static_cast<std::size_t>(i)] = m;
            res += small_diff[static_cast<std::size_t>(i)] * m;
        }
    }

    res += static_cast<i64>(n);
    return static_cast<u64>(res / 2);
}

bool run_checkpoints() {
    for (u64 n = 1ULL; n <= 200ULL; ++n) {
        const u64 brute = brute_f_from_lcm(n);
        const u64 formula = formula_f_from_factorization(n);
        if (brute != formula) {
            std::cerr << "Identity check failed at n=" << n << ": brute f(n)=" << brute
                      << ", formula f(n)=" << formula << '\n';
            return false;
        }
    }

    struct Checkpoint {
        int n;
        u64 expected;
    };
    const std::vector<Checkpoint> small = {
        {10, 29ULL},
        {100, 647ULL},
        {1'000, 11'751ULL},
        {10'000, 186'991ULL},
    };

    for (const Checkpoint& cp : small) {
        const u64 got = summatory_g_via_spf(cp.n);
        if (got != cp.expected) {
            std::cerr << "Small checkpoint failed for g(" << cp.n << "): expected " << cp.expected
                      << ", got " << got << '\n';
            return false;
        }
    }

    const u64 got = solve(kCheckpointN);
    if (got != kCheckpointExpected) {
        std::cerr << "Main checkpoint failed for g(" << kCheckpointN << "): expected "
                  << kCheckpointExpected << ", got " << got << '\n';
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 1;
    }

    std::cout << solve(options.n) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


if __name__ == "__main__":
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;

public class Euler379 {

    static long isqrt(long x) {
        if (x < 0)
            return 0;
        long r = (long) Math.sqrt(x);
        while ((r + 1) * (r + 1) <= x)
            r++;
        while (r * r > x)
            r--;
        return r;
    }

    static long icbrt(long x) {
        if (x < 0)
            return 0;
        long r = (long) Math.cbrt(x);
        while ((r + 1) * (r + 1) * (r + 1) <= x)
            r++;
        while (r * r * r > x)
            r--;
        return r;
    }

    static long sqSumFloor(long n, long x1, long x2) {
        if (x1 > x2)
            return 0;

        long x = x2;
        long s = 0;
        long beta = n / (x + 1);
        long epsilon = n % (x + 1);
        long delta = (n / x) - beta;
        long gamma = beta - x * delta;

        while (x >= x1) {
            epsilon += gamma;
            if (epsilon >= x) {
                delta += 1;
                gamma -= x;
                epsilon -= x;
                if (epsilon >= x) {
                    delta += 1;
                    gamma -= x;
                    epsilon -= x;
                    if (epsilon >= x) {
                        break;
                    }
                }
            } else if (epsilon < 0) {
                delta -= 1;
                gamma += x;
                epsilon += x;
            }

            gamma += 2 * delta;
            beta += delta;
            s += beta;
            x -= 1;
        }

        if (x < x1)
            return s;

        epsilon = n % (x + 1);
        delta = (n / x) - beta;
        gamma = beta - x * delta;

        while (x >= x1) {
            epsilon += gamma;
            long delta2 = epsilon / x;
            delta += delta2;
            epsilon -= x * delta2;
            gamma += 2 * delta - x * delta2;
            beta += delta;
            s += beta;
            x -= 1;
        }

        while (x >= x1) {
            s += n / x;
            x -= 1;
        }

        return s;
    }

    static long summatoryD3(long n) {
        long cbrtN = icbrt(n);
        long acc = 0;

        for (long z = 1; z <= cbrtN; z++) {
            long nz = n / z;
            long sqrtNz = isqrt(nz);
            acc += 2 * sqSumFloor(nz, z + 1, sqrtNz) - sqrtNz * sqrtNz + nz / z;
        }

        acc *= 3;
        acc += cbrtN * cbrtN * cbrtN;
        return acc;
    }

    static int[] mobiusSieve(int n) {
        int[] mu = new int[n + 1];
        int[] lp = new int[n + 1];
        int[] primes = new int[n / 2 + 10];
        int primeCount = 0;

        if (n >= 1)
            mu[1] = 1;
        for (int i = 2; i <= n; i++) {
            if (lp[i] == 0) {
                lp[i] = i;
                primes[primeCount++] = i;
                mu[i] = -1;
            }
            for (int j = 0; j < primeCount; j++) {
                int p = primes[j];
                long v = (long) i * p;
                if (v > n)
                    break;
                lp[(int) v] = p;
                if (i % p == 0) {
                    mu[(int) v] = 0;
                    break;
                }
                mu[(int) v] = -mu[i];
            }
        }
        return mu;
    }

    static String solve() {
        long n = 1000000000000L;
        long I = icbrt(n);
        int D = (int) isqrt(n / I);

        long res = 0;
        int[] mu = mobiusSieve(D);
        long[] mertensSmall = new long[D + 1];

        for (int d = 1; d <= D; d++) {
            int mud = mu[d];
            if (mud != 0) {
                long div = (long) d * d;
                res += mud * summatoryD3(n / div);
            }
            mertensSmall[d] = mertensSmall[d - 1] + mud;
        }

        if (I >= 2) {
            long[] smallDiff = new long[(int) I];
            for (int i = 1; i < I; i++) {
                smallDiff[i] = summatoryD3(i);
            }

            res -= mertensSmall[D] * smallDiff[(int) I - 1];

            for (int i = (int) I - 1; i >= 2; i--) {
                smallDiff[i] -= smallDiff[i - 1];
            }

            long[] mertensBig = new long[(int) I];
            for (int i = (int) I - 1; i >= 1; i--) {
                long v = isqrt(n / i);
                long vSqrt = isqrt(v);

                long m = 1 - v + vSqrt * mertensSmall[(int) vSqrt];
                for (long d = 2; d <= vSqrt; d++) {
                    long q = v / d;
                    if (q <= D) {
                        m -= mertensSmall[(int) q];
                    } else {
                        int idx = (int) (i * d * d);
                        m -= mertensBig[idx];
                    }
                    m -= (mertensSmall[(int) d] - mertensSmall[(int) d - 1]) * q;
                }

                mertensBig[i] = m;
                res += smallDiff[i] * m;
            }
        }

        res += n;
        long ans = res / 2;
        return Long.toString(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}