Detailed mathematical approach

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.

Step 2: Möbius Inversion Enforces \(\gcd(a,b,c)=1\)

Use the standard coprimality indicator

$\mathbf{1}_{\gcd(a,b,c)=1}=\sum_{d\mid \gcd(a,b,c)} \mu(d),$

where \(\mu\) is the Möbius function. Summing this identity over all triples in the box yields

$D(N)=\sum_{0\le a,b,c\le N}\mathbf{1}_{(a,b,c)\neq(0,0,0)}\sum_{d\mid \gcd(a,b,c)}\mu(d).$

Now exchange the order of summation. For a fixed divisor \(d\), the condition \(d\mid a,b,c\) means that each coordinate can be any multiple of \(d\) up to \(N\). That gives \(\left\lfloor N/d\right\rfloor+1\) choices per coordinate, including 0. The all-zero triple must then be removed exactly once.

Hence

$\boxed{D(N)=\sum_{d=1}^{N}\mu(d)\left(\left(\left\lfloor\frac{N}{d}\right\rfloor+1\right)^3-1\right).}$

This is the core formula implemented in the C++, Python, and Java files.

Step 3: Quotient Blocks Compress the Outer Sum

A naive loop over every \(d\le N\) is too slow when \(N=10^{10}\). The crucial observation is that the quotient

$q=\left\lfloor\frac{N}{d}\right\rfloor$

stays constant on whole intervals. If a block starts at \(l\), then every \(d\) in

$l\le d\le r=\left\lfloor\frac{N}{q}\right\rfloor$

has the same quotient \(q\). Over that interval the cubic factor is constant, so only the Möbius prefix sum changes. Introduce the Mertens function

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

Then for one block,

$\sum_{d=l}^{r}\mu(d)=M(r)-M(l-1),$

and the full sum becomes a block sum of the form

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

where the summation runs over all quotient blocks. The number of such blocks is only about \(2\sqrt N\), which is why the outer loop is practical.

Step 4: Fast Evaluation of the Mertens Prefix

The source files use the same two-level plan. First, they build \(\mu(d)\) with a linear sieve up to

$B=\max\left(10^6,\left\lfloor N^{2/3}\right\rfloor+1000\right).$

From this sieve they build the prefix array \(M(1),M(2),\dots,M(B)\).

For larger arguments, they use the classical identity

$\sum_{d=1}^{x}\mu(d)\left\lfloor\frac{x}{d}\right\rfloor=1,$

which can be rewritten as

$M(x)=1-\sum_{k=2}^{x} M\left(\left\lfloor\frac{x}{k}\right\rfloor\right).$

Again, equal floor quotients are grouped, so one recursive subtraction handles a whole interval at once. Large values are cached after the first computation, exactly as the `cache` structure does in the implementations.

Worked Checkpoints

For \(N=1\), every nonzero triple in \(\{0,1\}^3\) is primitive, so

$D(1)=7.$

For \(N=2\), the Möbius formula already shows the mechanism:

$D(2)=\mu(1)\left((2+1)^3-1\right)+\mu(2)\left((1+1)^3-1\right)=26-7=19.$

The C++ code also checks the optimized result against brute force at \(N=40\), where both methods give

$D(40)=56335.$

A larger hard-coded checkpoint is

$D(10^6)=831909254469114121,$

and only after these tests does the program attempt the target \(N=10^{10}\).

How the Code Works

The C++ file is the authoritative implementation. The function mu_sieve constructs the linear Möbius sieve, the class MertensPrefix serves \(M(x)\) either from the prefix array or from the memoized recurrence, and distinct_lines performs the quotient-block accumulation. The main program stores the running total in __int128, exposes --n= and --skip-checkpoints, and formats the final output through first9_last9_token.

The Python file is a compact fixed-\(N\) version of the same mathematics and relies on Python's arbitrary-precision integers. The Java file mirrors the same algorithm, but uses BigInteger for the cubic block value and the final accumulated answer because Java has no built-in signed 128-bit integer type.

Complexity Analysis

Let \(B=\max\left(10^6,\left\lfloor N^{2/3}\right\rfloor+1000\right)\). Building the linear sieve and the prefix Mertens table costs \(O(B)\) time and \(O(B)\) memory. The outer distinct-lines sum uses \(O(\sqrt N)\) quotient blocks. Large Mertens arguments are memoized and are themselves processed in grouped floor-quotient intervals, so the whole method is far below \(O(N)\) work and is practical for the implemented target \(N=10^{10}\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=388
2. Möbius inversion formula: Wikipedia — Möbius inversion formula
3. Möbius function: Wikipedia — Möbius function
4. Mertens function: Wikipedia — Mertens function
5. Quotient grouping and floor-sum methods: Wikipedia — Dirichlet hyperbola method