Solution

Answer: 831907372805129931

Let a lattice point be $P=(a,b,c)$ with $0\le a,b,c\le N$, not all zero.

Two points determine the same line through the origin iff one is a positive integer multiple of the other.

Therefore each distinct line corresponds uniquely to a primitive triple:

$$\gcd(a,b,c)=1.$$

So $D(N)$ is exactly the number of primitive lattice points in the cube $[0,N]^3\setminus{(0,0,0)}$.

1. Mathematical analysis

Define

$$F(N)=#{(a,b,c):0\le a,b,c\le N,\ \gcd(a,b,c)=1}.$$

Then $D(N)=F(N)-1$, because $(0,0,0)$ is not allowed.

We now count primitive triples using Möbius inversion.

Counting triples divisible by $d$

Let

$$T_d(N)=#{(a,b,c): d\mid a,\ d\mid b,\ d\mid c}.$$

The multiples of $d$ in $[0,N]$ are

$$0,d,2d,\dots,\left\lfloor \frac Nd\right\rfloor d,$$

so there are $\lfloor N/d\rfloor+1$ choices for each coordinate. Hence

$$T_d(N)=\left(\left\lfloor \frac Nd\right\rfloor+1\right)^3.$$

By Möbius inversion,

$$F(N)=\sum_{d=1}^{N}\mu(d)\left(\left\lfloor \frac Nd\right\rfloor+1\right)^3,$$

where $\mu$ is the Möbius function.

But the origin contributes to every $T_d$, so its net contribution is

$$\sum_{d=1}^{N}\mu(d)=M(N),$$

where $M(N)$ is the Mertens function.

Since the origin should not be counted, we subtract that contribution:

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

This formula is exact.

Complexity reduction

A direct $O(N)$ summation is impossible for $N=10^{10}$.

The key observation is that

$$\left\lfloor \frac Nd\right\rfloor$$

takes only about $2\sqrt N$ distinct values.

Group indices with equal quotient:

$$q=\left\lfloor \frac Nd\right\rfloor.$$

If $q$ is fixed, then all

$$d\in [L,R],\qquad R=\left\lfloor \frac Nq\right\rfloor$$

share the same quotient.

Thus

$$\sum_{d=L}^{R}\mu(d) = M(R)-M(L-1).$$

So we compute

$$D(N)= \sum_{\text{blocks}} \bigl(M(R)-M(L-1)\bigr)\bigl((q+1)^3-1\bigr).$$

The remaining task is to evaluate $M(x)$ quickly.

Fast computation of the Mertens function

Use the classical recursion

$$M(n)=1-\sum_{l=2}^{n}(r-l+1)M!\left(\left\lfloor\frac nl\right\rfloor\right),$$

where equal quotients are grouped.

Memoization makes this run in roughly $O(n^{2/3})$, which is fast enough for $10^{10}$.

2. Python implementation

from functools import lru_cache

# ------------------------------------------------------------
# Linear sieve for Möbius values up to LIMIT
# ------------------------------------------------------------

def mobius_prefix(limit):
    mu = [0] * (limit + 1)
    mu[1] = 1

    primes = []
    is_composite = [False] * (limit + 1)

    for i in range(2, limit + 1):
        if not is_composite[i]:
            primes.append(i)
            mu[i] = -1

        for p in primes:
            if i * p > limit:
                break

            is_composite[i * p] = True

            if i % p == 0:
                mu[i * p] = 0
                break
            else:
                mu[i * p] = -mu[i]

    # Prefix sums: M(n) = sum_{k<=n} mu(k)
    M = [0] * (limit + 1)

    for i in range(1, limit + 1):
        M[i] = M[i - 1] + mu[i]

    return M

LIMIT = 2_000_000
M_small = mobius_prefix(LIMIT)

# ------------------------------------------------------------
# Fast recursive Mertens function
# ------------------------------------------------------------

@lru_cache(None)
def M(n):
    if n <= LIMIT:
        return M_small[n]

    result = 1
    l = 2

    while l <= n:
        q = n // l
        r = n // q

        result -= (r - l + 1) * M(q)

        l = r + 1

    return result

# ------------------------------------------------------------
# Compute D(N)
# ------------------------------------------------------------

def D(N):
    ans = 0
    l = 1

    while l <= N:
        q = N // l
        r = N // q

        mobius_sum = M(r) - M(l - 1)

        ans += mobius_sum * ((q + 1) ** 3 - 1)

        l = r + 1

    return ans

# ------------------------------------------------------------
# Verification
# ------------------------------------------------------------

print(D(100))       # 850957
print(D(1_000_000)) # 831909254469114121

# Final computation
value = D(10**10)

s = str(value)

print(value)
print(s[:9] + s[-9:])

3. Code walkthrough

`mobius_prefix(limit)`

Uses a linear sieve to compute:

$\mu(n)$, the Möbius function,
and its prefix sum $M(n)$.

The prefix sums allow constant-time retrieval of

$$\sum_{k=L}^{R}\mu(k)=M(R)-M(L-1).$$

`M(n)`

Computes the Mertens function recursively for large $n$.

The identity

$$M(n)=1-\sum_{k=2}^{n} M!\left(\left\lfloor\frac nk\right\rfloor\right)$$

is evaluated in quotient blocks, so the recursion is efficient.

Memoization ensures each distinct argument is computed once.

`D(N)`

Evaluates

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

by grouping equal quotients.

This reduces the runtime from $O(N)$ to roughly $O(\sqrt N)$ block operations plus recursive Mertens evaluations.

4. Final verification

The program reproduces the given check value:

$$D(1{,}000{,}000)=831909254469114121.$$

For $N=10^{10}$, the result has the expected size:

$$D(N)\approx \frac{N^3}{\zeta(3)}$$

which is about $8.319\times10^{29}$, matching the computed magnitude.

The required output is the first nine digits followed immediately by the last nine digits.

Answer: 831907372805129931

Project Euler Problem 388