Problem 388: Distinct Lines

Project Euler Problem 388 Solution

EulerSolve provides an optimized solution for Project Euler Problem 388, Distinct Lines, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary In the implementation used in this repository, a distinct line is represented by a nonzero lattice direction $P=(a,b,c)$ with $0\le a,b,c\le N$. The line is the one passing through the origin and $P$. Two lattice points describe the same line exactly when one is a positive integer multiple of the other, so each line has a unique primitive representative whose coordinates have gcd equal to 1. Therefore the quantity being computed is $$D(N)=\#\left\{(a,b,c)\in \{0,\dots,N\}^3\setminus\{(0,0,0)\}:\gcd(a,b,c)=1\right\}.$$ The solutions then print the exact decimal value if it has at most 18 digits, or the concatenation of the first 9 and last 9 digits when the full answer is much longer. Mathematical Approach Step 1: Distinct Lines Become Primitive Triples If $g=\gcd(a,b,c)\gt 1$, then $(a,b,c)=g(a',b',c')$ with $\gcd(a',b',c')=1$. The points $(a,b,c)$ and $(a',b',c')$ lie on the same line through the origin, so the non-primitive vector is redundant. Conversely, every primitive triple inside the box gives one distinct line. Because the coordinates are restricted to $\{0,\dots,N\}$, there is no sign ambiguity: each line contributes exactly one primitive direction vector....

Detailed mathematical approach

Problem Summary

Draw the line from the origin $O(0,0,0)$ to every lattice point $(a,b,c)\ne O$ with $0\le a,b,c\le N$. Many of these lines coincide; $D(N)$ counts the distinct ones. The statement gives $D(1\,000\,000)=831909254469114121$ and asks for $D(10^{10})$, reported as its first nine digits followed by its last nine. ^[1]

The whole solution in one sentence: each line is named by the one point on it whose coordinates share no factor above $1$, the Möbius function counts those points by adding and subtracting the easy counts of points whose coordinates all share a factor $d$, and since $\lfloor N/d\rfloor$ takes at most $2\sqrt N$ values, that sum collapses into a few hundred thousand blocks weighted by prefix sums of $\mu$.

Mathematical Approach

1. Each line carries exactly one primitive point of the box

Call a lattice point $(a,b,c)\ne O$ primitive when $\gcd(a,b,c)=1$. Here $\gcd(0,0,c)=c$, so $(0,0,1)$ is primitive while $(0,0,2)$ is not.

Let $P\ne O$ lie in $B=\{0,1,\dots,N\}^3$ and put $g=\gcd(a,b,c)\ge1$. Then $P=gP'$ with $P'$ primitive, and no coordinate of $P'$ exceeds its counterpart in $P$, so $P'$ lies in $B$ too, on the line $OP$. It is the nearest nonzero lattice point of that line with non-negative coordinates: if $tP'$ is a lattice point and $t=u/v$ in lowest terms, $v$ divides every coordinate of $P'$, so $v=1$ and $t$ is an integer.

A line carries no second one. Let $P,Q\ne O$ lie in $B$ on one line through $O$, so $Q=tP$ with real $t\gt0$ (no coordinate is negative and neither point is $O$). With $P=gP'$ and $Q=hQ'$ primitive, $Q'=(tg/h)P'$; writing that positive rational as $u/v$ in lowest terms gives $vQ'=uP'$, whose coordinate gcd is $v$ on the left and $u$ on the right. So $u=v$ and $Q'=P'$, giving

$$D(N)=\#\{\,(a,b,c)\in B:\gcd(a,b,c)=1\,\}.$$

This already excludes $(0,0,0)$, whose gcd is $0$; the divisor identity of Section 2 needs a positive argument, so that point is removed there by hand.

2. The Möbius function as a coprimality test

The Möbius function is $\mu(1)=1$, $\mu(n)=(-1)^k$ for a product of $k$ distinct primes, and $\mu(n)=0$ when a square above $1$ divides $n$. ^[2] Its first twelve values are $1,-1,-1,0,-1,1,-1,0,0,1,-1,0$. We need the divisor sum

$$\sum_{d\mid n}\mu(d)=\begin{cases}1,&n=1\\0,&n\gt1\end{cases}$$

For instance $6$ has divisors $1,2,3,6$ with Möbius values $1,-1,-1,1$, summing to $0$. The reason is a pairing: only square-free divisors contribute, so fix a prime $p\mid n$ and match each square-free $e$ with $p\nmid e$ to $pe$; the two Möbius values differ only in sign, so all cancel. For $n=1$ only $\mu(1)=1$ is left.

So the divisor sum indicates $n=1$. Every point of $B$ but $O$ has a positive gcd, so applying the identity to it and summing gives

$$D(N)=\sum_{P\in B,\,P\ne O}\ \sum_{d\mid\gcd(a,b,c)}\mu(d).$$

3. Exchanging the two sums

Now count by $d$ instead of by $P$. The condition $d\mid\gcd(a,b,c)$ just says $d$ divides each coordinate, and the multiples of $d$ in $\{0,1,\dots,N\}$ number $\lfloor N/d\rfloor+1$, the value $0$ included. Removing the all-zero triple from each sub-box gives

$$D(N)=\sum_{d=1}^{N}\mu(d)\bigl((\lfloor N/d\rfloor+1)^3-1\bigr).$$

For $N=3$ the three terms are $\mu(1)(4^3-1)=63$, $\mu(2)(2^3-1)=-7$ and $\mu(3)(2^3-1)=-7$, so $D(3)=49$. For $N=2$ only two terms survive: $3^3-1=26$ nonzero points, minus the $2^3-1=7$ points whose coordinates are all even, so $D(2)=19$. Those $7$ discarded points are the doubles of the $7$ nonzero points of $\{0,1\}^3$, that is $D(1)=7$ again.

4. What the three parts of the summand count

Expanding $(q+1)^3-1=q^3+3q^2+3q$ with $q=\lfloor N/d\rfloor$ splits the formula into three Möbius sums, each with its own meaning. By the same exchange, $\sum_d\mu(d)\lfloor N/d\rfloor^{k}$ counts the $k$-tuples in $\{1,\dots,N\}^{k}$ of gcd $1$. So the pieces count the primitive points of $B$ with no zero coordinate, with one (the $3$ coordinate planes) and with two (the $3$ axes) — and the last is $1$ for every $N$, the only positive integer coprime to itself.

At $N=2$ the split is $19=7+9+3$: the axes give $(1,0,0)$, $(0,1,0)$, $(0,0,1)$; the planes give $(1,1,0)$, $(1,2,0)$, $(2,1,0)$ and their counterparts in the other two; the interior gives $(1,1,1)$, $(1,1,2)$, $(1,2,1)$, $(2,1,1)$, $(1,2,2)$, $(2,1,2)$, $(2,2,1)$. At $N=10$ it is $1033=841+189+3$, the middle part being $3$ times the $63$ coprime pairs in $\{1,\dots,10\}^2$. The split is an independent check, not part of the algorithm.

5. Only about $2\sqrt N$ different quotients

A loop over all $d\le N$ is hopeless at $N=10^{10}$, but $(q+1)^3-1$ depends on $d$ only through $q=\lfloor N/d\rfloor$, which is constant on long stretches: if $q=\lfloor N/l\rfloor$, every $d$ with $l\le d\le r=\lfloor N/q\rfloor$ has that quotient, and $r+1$ starts the next stretch. Stretches are few — at most $\sqrt N$ values $d\le\sqrt N$, and for $d\gt\sqrt N$ at most $\sqrt N$ values of $q\lt\sqrt N$ — so at most $2\sqrt N$ blocks, exactly $199\,999$ at $N=10^{10}$, against $10^{10}$ single terms.

Inside a block only the Möbius values change, and their sum is a difference of prefix sums of $\mu$ — the Mertens function ^[3]

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

so the formula becomes one term per block,

$$D(N)=\sum\bigl((q+1)^3-1\bigr)\bigl(M(r)-M(l-1)\bigr),$$

with $l$, $r$, $q$ running over the blocks just described.

At $N=10$ there are five blocks $l..r$ with quotient $q$: $1..1$, $q=10$; $2..2$, $q=5$; $3..3$, $q=3$; $4..5$, $q=2$; $6..10$, $q=1$. Their Möbius weights are $1,-1,-1,-1,1$, so

$$D(10)=1330-215-63-26+7=1033,$$

matching the brute-force count. What remains is to evaluate $M$ at the block ends: at most $2\sqrt N$ positive arguments, plus $M(0)=0$.

6. The Mertens function from its own divisor identity

Summing the indicator of Section 2 over $n\le x$ gives $\sum_{d\le x}\mu(d)\lfloor x/d\rfloor=1$, since only $n=1$ contributes. Reading that double sum by the size of the cofactor instead turns it into ^[3]

$$\sum_{q=1}^{x}M(\lfloor x/q\rfloor)=1,$$

and isolating $q=1$ leaves a recurrence in strictly smaller arguments:

$$M(x)=1-\sum_{q=2}^{x}M(\lfloor x/q\rfloor).$$

Take $x=10$: the values $M(\lfloor 10/q\rfloor)$ for $q=2,\dots,10$ are $-2,-1,0,0,1,1,1,1,1$, summing to $2$, so $M(10)=-1$, as $\mu(1)+\dots+\mu(10)$ confirms. The block trick applies here too, the arguments repeating: four grouped subtractions suffice, $1\cdot M(5)$, $1\cdot M(3)$, $2\cdot M(2)$ and $5\cdot M(1)$, giving $1-(-2-1+0+5)=-1$. Each $M(x)$ thus costs about $2\sqrt x$ grouped steps, and since $\lfloor\lfloor N/a\rfloor/b\rfloor=\lfloor N/(ab)\rfloor$, every argument appearing is again some $\lfloor N/k\rfloor$; memoising charges each one those steps once.

7. A linear sieve for the small values

The small arguments are very many, and each still costs its grouped steps. Producing them together is far cheaper: tabulate $\mu$ up to a cut-off $L$, accumulate its prefix sums, and answer every $M(x)$ with $x\le L$ by one array lookup. The table comes from a linear sieve, recording every $i\le L$ with its smallest prime factor $\mathrm{lp}(i)$ in $O(L)$ total time: each composite is written once, as $i\cdot p$ with $p\le\mathrm{lp}(i)$. ^[4] The same pass carries $\mu$: if $p\lt\mathrm{lp}(i)$ then $p\nmid i$, so $ip$ gains one new prime factor and $\mu(ip)=-\mu(i)$; if $p=\mathrm{lp}(i)$ then $p^2\mid ip$ and $\mu(ip)=0$, which is also where the inner loop stops. Primes get $\mu=-1$, and $\mu(1)=1$.

8. Where to put the cut-off

The sieve costs $O(L)$ steps. Above $L$ the values needed are $M(\lfloor N/k\rfloor)$ for $k\le N/L$, the $k$-th costing about $2\sqrt{N/k}$ grouped steps, so the recursion costs about

$$\sum_{k\le N/L}2\sqrt{N/k}\approx4\sqrt{N\cdot N/L}=\frac{4N}{\sqrt L}$$

steps: comparing the decreasing $t^{-1/2}$ with its integral gives $\sum_{k\le K}k^{-1/2}\le2\sqrt K$, so the recursion takes $O(N/\sqrt L)$ arithmetic operations at constant average cost per cache lookup. The total $L+4N/\sqrt L$ is smallest when $L^{3/2}=2N$, so any $L$ of order $N^{2/3}$ is within a constant factor of optimal and runs in $O(N^{2/3})$. The solutions take $L=\max(10^6,\lfloor N^{2/3}\rfloor+1000)$, the power evaluated in floating point, so a whole-number $N^{2/3}$ truncates one too low; neighbouring cut-offs serve as well. At $N=10^{10}$ the power is $4641588.83\ldots$ and $L=4\,642\,588$, with $325\,215$ primes below it, and the predicted $4N/\sqrt L\approx1.86\times10^7$ recursion steps match the $18\,264\,879$ performed. The floor of $10^6$ is why the checkpoint $N=10^6$ needs no recursion: there $L\ge N$, so all $1999$ blocks come from the sieve.

9. How large the numbers become

Write $\zeta(s)=\sum_{d\ge1}d^{-s}$ for $s\gt1$; the Möbius series gives $\sum_{d\ge1}\mu(d)/d^3=1/\zeta(3)$. ^[2] To reach our finite count from it, put $t=N/d\ge1$ and $q=\lfloor t\rfloor$. Then

$$t^3-1\le(q+1)^3-1\le t^3+7t^2,$$

because $q\le t\lt q+1$ and $(t+1)^3\le t^3+7t^2+1$ once $t\ge1$. Summed over $d\le N$ the error is at most $7\zeta(2)N^2$, and extending the cubic main term to all $d$ costs at most $N/2$, since $\sum_{d\gt N}d^{-3}\le1/(2N^2)$. Hence

$$D(N)=\frac{N^3}{\zeta(3)}+O(N^2).$$

The limiting density $1/\zeta(3)\approx0.8319$ is the probability that three integers drawn uniformly from $\{1,\dots,N\}$ share no prime factor, as $N$ grows. ^[5] ^[6] The given $D(10^6)$ lies within a relative $2.3\times10^{-6}$ of $10^{18}/\zeta(3)$, as predicted, and at $N=10^{10}$ the main term is near $0.83\times10^{30}$ while the error is only $O(N^2)$, so the value has thirty decimal digits. A signed 64-bit integer stops near $9.22\times10^{18}$ and cannot hold even the single term $d=1$, about $10^{30}$; the same inequality bounds the sum of absolute terms by $\zeta(3)N^3+7\zeta(2)N^2$, below $1.21\times10^{30}$ at the target and well inside the signed 128-bit range of about $1.70\times10^{38}$. Mertens arithmetic stays small, though: the largest grouped product $(r-l+1)M(q)$ at $N=10^{10}$ is $5\times10^{9}$ and the largest partial sum in the recurrence is $4\,999\,966\,278$, so 64 bits are ample there.

How the Code Works

The C++ program ^[7] keeps its settings in Options: the target $N$ and a checkpoint flag. parse_arguments accepts --skip-checkpoints and, through parse_i64_after_prefix, a --n= value, which is how the small cases above are reproduced.

mu_sieve is the linear sieve of Section 7, returning $\mu(1),\dots,\mu(L)$ in a vector whose unused entry $0$ stays zero; it keeps lp and a growing primes list, and breaks the inner loop when $i\cdot p$ leaves the range or p == lp[i]. MertensPrefix wraps the two regimes of Sections 6 and 7: the constructor fixes sieve_limit_ as $\max(10^6,\lfloor N^{2/3}\rfloor+1000)$, runs the sieve and accumulates prefix_, the running sums of $\mu$; M answers from prefix_ up to that limit, otherwise from cache_, and on a miss evaluates the recurrence block by block with ans -= (r - l + 1) * M(q) before storing it. The recursion stays shallow: the first block already halves the argument.

distinct_lines is Section 5 verbatim: it walks the blocks $l\dots r=\lfloor N/q\rfloor$, takes mertens.M(r) - mertens.M(l - 1) as the Möbius weight, and adds $((q+1)^3-1)$ times it to an __int128 accumulator, the cube itself built in 128-bit arithmetic since $(q+1)^3$ alone overflows 64 bits. ^[8] run_checkpoints tests the machinery on $D(1)=7$, on agreement with brute_distinct_lines at $N=40$ — a triple loop over all $41^3$ points that skips the origin and calls std::gcd twice on the rest — and on $D(10^6)$ against the statement. Finally to_string_i128 converts the accumulator digit by digit, the standard streams being unable to print a 128-bit integer, and first9_last9_token keeps the first and last nine characters of anything longer than eighteen.

Python packs the same algorithm into one function solve, with $N$ fixed, a nested M closure and a cache dictionary; its integers are unbounded, so the accumulator needs no care. ^[8] Java mirrors the C++ structure with static fields, a HashMap cache and an int[] prefix array — safe here, as $|M(x)|$ never exceeds $683$ for $x\le L$ — but has no 128-bit primitive, so the cube, the product and the accumulator are BigInteger while the recurrence stays in long. ^[8] Neither port runs the checkpoints.

Complexity and Verification

Time is $O(L+N/\sqrt L)=O(N^{2/3})$ for $L\asymp N^{2/3}$, as derived in Section 8: $199\,999$ outer blocks, and $18\,264\,879$ grouped subtractions for the $2153$ Mertens values above $L$. Memory is $O(L+N/L+\log N)$ — two $L$-sized arrays in the sieve, one for the prefix sums, $O(N/L)$ memoised values and recursion depth $O(\log N)$ — hence $O(L)$ at $L\asymp N^{2/3}$. The table holds those $2153$ entries although the constructor calls cache_.reserve for $2^{20}$; and prefix entries are 64-bit in C++, 32-bit in Java and object references in Python, so equal asymptotic space is not equal bytes. Measured on one machine — C++ built by GCC at -O2, timed with /usr/bin/time — the peak resident size was about $83$ MiB, $98$ MiB with the checkpoints first, and the times $0.11$ s for C++ with checkpoints, $0.17$ s for Java and $2.1$ s for Python, which the sieve dominates.

The packaged checks are the three in run_checkpoints: $D(1)=7$, the brute-force $D(40)=56335$, and $D(10^6)=831909254469114121$ from the statement. Verified independently for this article: a triple loop over the box agrees with the Möbius formula for every $N\le20$, giving $7,19,49,91,175,$${}253,415,571,805,1033$ at $N=1,\dots,10$, the known sequence for a cube of side $N$ ^[9]; the block-and-Mertens method equals the term-by-term sum for every $N\le200$ and for $N=500,997,1000,2000,4096,5000$; the split of Section 4 matches a brute force classifying each primitive point by its zero coordinates, at $N=1,2,3,5,10,20,40$, with $\sum_{d\le N}\mu(d)\lfloor N/d\rfloor=1$ for $N\le12$; and at $N=10^{10}$ the three programs give the same eighteen-character token, while C++ with --n=1000000 reproduces the statement value.

Two limitations deserve recording. Only the C++ program validates itself, so a reader running the Python or Java file gets no check; comparing them against C++ at a small $N$ is the remedy. And the --n= option has unchecked range limits: its decimal loop forms parsed * 10 + digit with no guard, so a long enough argument overflows long long and is accepted with a wrong value, and the estimated cut-off is then narrowed to an int, which leaves range near $N\approx10^{14}$ — where the arrays would need some $26$ GB anyway. Accepting digits is no promise that any positive input works; the target and the checkpoints sit far inside both limits.

Footnotes and References

The references below supply the background facts used above. The bijection with primitive points, the exchange of sums, the meaning of the summand's three parts, the block decomposition, the recurrence for $M$, the error estimate and the cost $L+4N/\sqrt L$ are derived here.

Project Euler 388 — Distinct Lines. The official statement takes all lattice points (a, b, c) with 0 ≤ a, b, c ≤ N, draws every line from the origin to the others, defines D(N) as the number of distinct such lines, gives D(1 000 000) = 831909254469114121, and asks for D(10^10) as its first nine digits then its last nine.
Wikipedia — Möbius function. Defines μ(1) = 1, μ(n) = (−1)^k for a product of k distinct primes and μ(n) = 0 when a square above 1 divides n; states that Σ μ(d) over the divisors of n is 1 for n = 1 and 0 otherwise; and gives Σ μ(n)/n^s = 1/ζ(s), used in Section 9.
Wikipedia — Mertens function. Defines M(x) as Σ μ(k) for k ≤ x, lists M(1) to M(12) as 1, 0, −1, −1, −2, −1, −2, −2, −2, −1, −2, −2, and states the identity of Section 6, that Σ M(n/d) over d ≤ n equals 1. It also gives O(x^(2/3)(log log x)^(1/3)) for isolated values by the methods of Lagarias–Miller–Odlyzko and Deléglise–Rivat — the exponent 2/3 reached here too.
cp-algorithms — Linear Sieve. Gives the sieve that records the minimum prime factor lp(i) of every i up to n in linear total time, and shows why the inner loop stops at the prime equal to lp(i): every composite x is written once, as i·p with p the smallest prime factor of x and p ≤ lp(i). Carrying μ along is not covered there; Section 7 justifies it.
Wikipedia — Coprime integers. States, citing Nymann, that k randomly chosen integers are setwise coprime with probability 1/ζ(k), with 6/π² = 1/ζ(2) at k = 2. The case k = 3 is the density of Section 9.
Wikipedia — Apéry's constant. Gives ζ(3) ≈ 1.2020569031595942854, notes Apéry's 1978 proof of its irrationality, and gives its reciprocal, about 0.8319, as the probability that three random positive integers are relatively prime.
C++ — Euler388.cpp. The linked immutable revision contains mu_sieve, MertensPrefix, distinct_lines, brute_distinct_lines and run_checkpoints. It is the implementation record here; the extra comparisons above are not packaged in it.
GCC manual — 128-bit integers documents __int128 as an extension on targets with an integer mode wide enough for 128 bits, and notes that where long long is narrower no 128-bit constant can be written; so the code widens 64-bit values before cubing and converts to decimal itself, stream insertion not accepting the extension type. Java — BigInteger documents immutable arbitrary-precision integers with the add, multiply and toString the port uses, and Python — numeric types states that integers have unlimited precision.
OEIS A090025. "Number of distinct lines through the origin in 3-dimensional cube of side length n", with terms 0, 7, 19, 49, 91, 175, 253, 415, 571, 805, 1033, … from n = 0, a comment identifying them as the lattice points of coordinate gcd 1, and the formula of Section 3.

Problem 388 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <unordered_map>
#include <vector>

namespace {

using i64 = long long;
using i128 = __int128_t;

struct Options {
    i64 n = 10000000000LL;
    bool run_checkpoints = true;
};

bool parse_i64_after_prefix(const std::string& arg, const std::string& prefix, i64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    i64 parsed = 0;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<i64>(ch - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(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;
        }
        if (parse_i64_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1;
}

std::string to_string_i128(i128 value) {
    if (value == 0) {
        return "0";
    }
    bool neg = false;
    if (value < 0) {
        neg = true;
        value = -value;
    }
    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    if (neg) {
        s.push_back('-');
    }
    std::reverse(s.begin(), s.end());
    return s;
}

std::string first9_last9_token(const std::string& s) {
    if (s.size() <= 18U) {
        return s;
    }
    return s.substr(0, 9) + s.substr(s.size() - 9);
}

std::vector<int> mu_sieve(const int limit) {
    std::vector<int> primes;
    primes.reserve(limit / 10);
    std::vector<int> lp(static_cast<std::size_t>(limit + 1), 0);
    std::vector<int> mu(static_cast<std::size_t>(limit + 1), 0);
    mu[1] = 1;

    for (int i = 2; i <= limit; ++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 i64 x = 1LL * i * p;
            if (x > limit) {
                break;
            }
            lp[static_cast<std::size_t>(x)] = p;
            if (p == lp[static_cast<std::size_t>(i)]) {
                mu[static_cast<std::size_t>(x)] = 0;
                break;
            }
            mu[static_cast<std::size_t>(x)] = -mu[static_cast<std::size_t>(i)];
        }
    }
    return mu;
}

class MertensPrefix {
   public:
    explicit MertensPrefix(const i64 n) {
        const i64 estimated = static_cast<i64>(std::pow(static_cast<long double>(n), 2.0L / 3.0L));
        sieve_limit_ = static_cast<int>(std::max<i64>(1000000LL, estimated + 1000));

        const std::vector<int> mu = mu_sieve(sieve_limit_);
        prefix_.assign(static_cast<std::size_t>(sieve_limit_ + 1), 0LL);
        for (int i = 1; i <= sieve_limit_; ++i) {
            prefix_[static_cast<std::size_t>(i)] =
                prefix_[static_cast<std::size_t>(i - 1)] + static_cast<i64>(mu[static_cast<std::size_t>(i)]);
        }
        cache_.reserve(1 << 20);
    }

    i64 M(const i64 n) {
        if (n <= sieve_limit_) {
            return prefix_[static_cast<std::size_t>(n)];
        }
        const auto it = cache_.find(n);
        if (it != cache_.end()) {
            return it->second;
        }

        i64 ans = 1;
        for (i64 l = 2; l <= n;) {
            const i64 q = n / l;
            const i64 r = n / q;
            ans -= (r - l + 1) * M(q);
            l = r + 1;
        }

        cache_[n] = ans;
        return ans;
    }

   private:
    int sieve_limit_{0};
    std::vector<i64> prefix_;
    std::unordered_map<i64, i64> cache_;
};

i128 distinct_lines(const i64 n) {
    MertensPrefix mertens(n);
    i128 answer = 0;

    for (i64 l = 1; l <= n;) {
        const i64 q = n / l;
        const i64 r = n / q;
        const i64 mu_block = mertens.M(r) - mertens.M(l - 1);

        const i128 qq = static_cast<i128>(q);
        const i128 f = (qq + 1) * (qq + 1) * (qq + 1) - 1;
        answer += f * static_cast<i128>(mu_block);

        l = r + 1;
    }

    return answer;
}

i128 brute_distinct_lines(const int n) {
    i128 count = 0;
    for (int a = 0; a <= n; ++a) {
        for (int b = 0; b <= n; ++b) {
            for (int c = 0; c <= n; ++c) {
                if (a == 0 && b == 0 && c == 0) {
                    continue;
                }
                const int g = std::gcd(a, std::gcd(b, c));
                if (g == 1) {
                    ++count;
                }
            }
        }
    }
    return count;
}

bool run_checkpoints() {
    if (distinct_lines(1) != 7) {
        std::cerr << "Checkpoint failed: D(1)\n";
        return false;
    }

    const int brute_n = 40;
    const i128 brute = brute_distinct_lines(brute_n);
    const i128 fast = distinct_lines(brute_n);
    if (fast != brute) {
        std::cerr << "Checkpoint failed against brute force at N=" << brute_n << '\n';
        return false;
    }

    const i128 known = distinct_lines(1000000);
    if (to_string_i128(known) != "831909254469114121") {
        std::cerr << "Checkpoint failed: D(10^6), got " << to_string_i128(known) << '\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 2;
    }

    const i128 answer = distinct_lines(options.n);
    const std::string full = to_string_i128(answer);
    std::cout << first9_last9_token(full) << '\n';
    return 0;
}

Python

def solve():
    import math

    n = 10000000000
    sieve_limit = max(1000000, int(math.pow(n, 2.0 / 3.0)) + 1000)

    primes = []
    lp = [0] * (sieve_limit + 1)
    mu = [0] * (sieve_limit + 1)
    mu[1] = 1

    for i in range(2, sieve_limit + 1):
        if lp[i] == 0:
            lp[i] = i
            primes.append(i)
            mu[i] = -1
        for p in primes:
            x = i * p
            if x > sieve_limit:
                break
            lp[x] = p
            if p == lp[i]:
                mu[x] = 0
                break
            mu[x] = -mu[i]

    prefix = [0] * (sieve_limit + 1)
    for i in range(1, sieve_limit + 1):
        prefix[i] = prefix[i - 1] + mu[i]

    cache = {}

    def M(x):
        if x <= sieve_limit:
            return prefix[x]
        if x in cache:
            return cache[x]

        ans = 1
        l = 2
        while l <= x:
            q = x // l
            r = x // q
            ans -= (r - l + 1) * M(q)
            l = r + 1

        cache[x] = ans
        return ans

    answer = 0
    l = 1
    while l <= n:
        q = n // l
        r = n // q
        mu_block = M(r) - M(l - 1)

        f = (q + 1) ** 3 - 1
        answer += f * mu_block

        l = r + 1

    s = str(answer)
    if len(s) <= 18:
        return s
    return s[:9] + s[-9:]

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

Java

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.math.BigInteger;

public class Euler388 {
    private static final long n = 10000000000L;
    private static int sieveLimit;
    private static int[] prefix;
    private static Map<Long, Long> cache = new HashMap<>();

    public static String solve() {
        sieveLimit = (int) Math.max(1000000L, (long) Math.pow(n, 2.0 / 3.0) + 1000);

        int[] lp = new int[sieveLimit + 1];
        int[] mu = new int[sieveLimit + 1];
        List<Integer> primes = new ArrayList<>();
        mu[1] = 1;

        for (int i = 2; i <= sieveLimit; ++i) {
            if (lp[i] == 0) {
                lp[i] = i;
                primes.add(i);
                mu[i] = -1;
            }
            for (int p : primes) {
                long x = (long) i * p;
                if (x > sieveLimit)
                    break;
                lp[(int) x] = p;
                if (p == lp[i]) {
                    mu[(int) x] = 0;
                    break;
                }
                mu[(int) x] = -mu[i];
            }
        }

        prefix = new int[sieveLimit + 1];
        long current = 0;
        for (int i = 1; i <= sieveLimit; ++i) {
            current += mu[i];
            prefix[i] = (int) current;
        }

        BigInteger answer = BigInteger.ZERO;
        long l = 1;
        while (l <= n) {
            long q = n / l;
            long r = n / q;
            long muBlock = M(r) - M(l - 1);

            BigInteger bq = BigInteger.valueOf(q + 1);
            BigInteger f = bq.multiply(bq).multiply(bq).subtract(BigInteger.ONE);
            answer = answer.add(f.multiply(BigInteger.valueOf(muBlock)));

            l = r + 1;
        }

        String s = answer.toString();
        if (s.length() <= 18)
            return s;
        return s.substring(0, 9) + s.substring(s.length() - 9);
    }

    private static long M(long x) {
        if (x <= sieveLimit)
            return prefix[(int) x];
        Long cached = cache.get(x);
        if (cached != null)
            return cached;

        long ans = 1;
        long l = 2;
        while (l <= x) {
            long q = x / l;
            long r = x / q;
            ans -= (r - l + 1) * M(q);
            l = r + 1;
        }

        cache.put(x, ans);
        return ans;
    }

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