Problem 379: Least Common Multiple Count

View on Project Euler

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

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());
    }
}