Project Euler

Distinct Lines through Lattice Points

Consider all lattice points (a, b, c) with 0 <= a, b, c <= N. From the origin O(0,0,0) a straight line is drawn to every point. Let D(N) be the number of distinct such lines. You are given that D(1...

Source sync Apr 19, 2026

Problem #0388

Level Level 19

Solved By 695

Languages C++, Python

Answer 831907372805129931

Length 345 words

number_theorygeometrylinear_algebra

Consider all lattice points $(a,b,c)$ with $0 \le a,b,c \le N$.

From the origin $O(0,0,0)$ all lines are drawn to the other lattice points.

Let $D(N)$ be the number of distinct such lines.

You are given that $D(1\,000\,000) = 831909254469114121$.

Find $D(10^{10})$. Give as your answer the first nine digits followed by the last nine digits.

Problem 388: Distinct Lines through Lattice Points

Mathematical Analysis

Two lattice points $(a_1, b_1, c_1)$ and $(a_2, b_2, c_2)$ define the same line through the origin if and only if one is a positive scalar multiple of the other. That is, $(a_1, b_1, c_1) = k(a_2, b_2, c_2)$ for some positive rational $k$ .

Equivalently, each distinct line corresponds to a unique primitive direction: a point $(a, b, c)$ in $\{0, \ldots, N\}^3 \setminus \{(0,0,0)\}$ with $\gcd(a, b, c) = 1$ .

Therefore:

D(N) = \#\{(a,b,c) \in \{0,\ldots,N\}^3 \setminus \{(0,0,0)\} : \gcd(a,b,c) = 1\}

Derivation / Algorithm

Mobius Inversion Formula

Using the identity $[\gcd(a,b,c) = 1] = \sum_{d \mid \gcd(a,b,c)} \mu(d)$ , we get:

D(N) = \sum_{\substack{(a,b,c) \in \{0,\ldots,N\}^3 \\ (a,b,c) \neq (0,0,0)}} \sum_{d \mid \gcd(a,b,c)} \mu(d)

Swapping summation order and substituting $a = da'$ , $b = db'$ , $c = dc'$ :

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

The $+1$ arises because each coordinate in $\{0, d, 2d, \ldots\}$ contributes $\lfloor N/d \rfloor + 1$ values. The $-1$ subtracts the origin.

Sub-linear Computation via Mertens Function

For large $N$ (e.g., $N = 10^{10}$ ), we cannot sieve $\mu(d)$ for all $d \le N$ . Instead, we group terms by $q = \lfloor N/d \rfloor$ :

D(N) = \sum_{\text{distinct } q} g(q) \cdot \bigl(M(\lfloor N/q \rfloor) - M(\lfloor N/(q+1) \rfloor)\bigr)

where $g(q) = (q+1)^3 - 1$ and $M(x) = \sum_{d=1}^{x} \mu(d)$ is the Mertens function.

There are only $O(\sqrt{N})$ distinct values of $\lfloor N/d \rfloor$ . The Mertens function is computed using the recursive identity:

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

with memoization and a sieve for small values up to threshold $N^{2/3}$ , giving an overall time complexity of $O(N^{2/3})$ .

Proof of Correctness

Bijection: Each distinct line through the origin and a point in $\{0,\ldots,N\}^3$ corresponds to a unique primitive direction vector $(a,b,c)$ with $\gcd(a,b,c) = 1$ . Every non-primitive point $(a,b,c)$ with $g = \gcd(a,b,c) > 1$ lies on the same line as $(a/g, b/g, c/g)$ , which is primitive and also in the cube since $a/g \le a \le N$ .
Mobius inversion: The standard identity $\sum_{d \mid n} \mu(d) = [n = 1]$ correctly extracts the coprimality condition.
Verification: $D(1\,000\,000) = 831\,909\,254\,469\,114\,121$ matches the problem’s given value. Small cases ( $D(1) = 7$ , $D(2) = 19$ , etc.) are verified against brute force.

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Complexity Analysis

Sieve phase: $O(N^{2/3})$ time and space to compute $\mu$ and prefix sums up to threshold.
Mertens recursion: $O(N^{2/3})$ time with memoization over $O(\sqrt{N})$ distinct values.
Main summation: $O(\sqrt{N})$ iterations.
Overall: $O(N^{2/3})$ time, $O(N^{2/3})$ space.

Answer

\boxed{831907372805129931}

C++ project_euler/problem_388/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Problem 388: Distinct Lines through Lattice Points
 *
 * D(N) = sum_{d=1}^{N} mu(d) * ((floor(N/d)+1)^3 - 1)
 *
 * Sub-linear computation via Mertens function M(x) = sum_{d<=x} mu(d).
 * Time: O(N^{2/3}), Space: O(N^{2/3}).
 */

using i64 = long long;
using i128 = __int128;

const i64 N = 10000000000LL;  // 10^10
int THRESH;

vector<int> mu_arr;
vector<i64> M_small;
unordered_map<i64, i64> M_cache;

void sieve_mu() {
    mu_arr.assign(THRESH + 1, 0);
    mu_arr[1] = 1;
    vector<bool> is_prime(THRESH + 1, true);
    vector<int> primes;
    for (int i = 2; i <= THRESH; i++) {
        if (is_prime[i]) {
            primes.push_back(i);
            mu_arr[i] = -1;
        }
        for (int p : primes) {
            if ((i64)i * p > THRESH) break;
            is_prime[i * p] = false;
            if (i % p == 0) {
                mu_arr[i * p] = 0;
                break;
            }
            mu_arr[i * p] = -mu_arr[i];
        }
    }
    M_small.assign(THRESH + 1, 0);
    for (int i = 1; i <= THRESH; i++)
        M_small[i] = M_small[i - 1] + mu_arr[i];
}

i64 mertens(i64 n) {
    if (n <= THRESH) return M_small[n];
    auto it = M_cache.find(n);
    if (it != M_cache.end()) return it->second;

    i64 result = 1;
    i64 k = 2;
    while (k <= n) {
        i64 q = n / k;
        i64 k_next = n / q + 1;
        result -= (k_next - k) * mertens(q);
        k = k_next;
    }
    M_cache[n] = result;
    return result;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    // Threshold ~ N^{2/3}
    THRESH = (int)cbrt((double)N * N) + 1;
    // More precise: N^{2/3}
    while ((i128)THRESH * THRESH * THRESH < (i128)N * N) THRESH++;

    sieve_mu();

    // D(N) = sum over blocks grouped by q = floor(N/d)
    i128 total = 0;
    i64 d = 1;
    while (d <= N) {
        i64 q = N / d;
        i64 d_max = N / q;
        i64 mu_sum = mertens(d_max) - mertens(d - 1);
        i128 g = (i128)(q + 1) * (q + 1) * (q + 1) - 1;
        total += (i128)mu_sum * g;
        d = d_max + 1;
    }

    // Format: first 9 digits + last 9 digits
    // Convert i128 to string
    string s;
    i128 val = total;
    if (val == 0) {
        s = "0";
    } else {
        bool neg = val < 0;
        if (neg) val = -val;
        while (val > 0) {
            s += '0' + (int)(val % 10);
            val /= 10;
        }
        if (neg) s += '-';
        reverse(s.begin(), s.end());
    }

    string answer = s.substr(0, 9) + s.substr(s.size() - 9);
    cout << answer << endl;
    return 0;
}

Python project_euler/problem_388/solution.py

"""
Problem 388: Distinct Lines through Lattice Points

Consider all lattice points (a, b, c) with 0 <= a, b, c <= N.
From the origin O(0,0,0) a straight line is drawn to every point.
Let D(N) be the number of distinct such lines.

Given: D(1,000,000) = 831909254469114121
Find: D(10^10), answer as first 9 digits followed by last 9 digits.

Key identity (Mobius inversion):
    D(N) = sum_{d=1}^{N} mu(d) * ((floor(N/d) + 1)^3 - 1)

For large N, use sub-linear Mertens function computation.
"""

from math import gcd

def solve_brute_force(n: int):
    """Brute-force D(N) by enumerating primitive directions. O(N^3)."""
    count = 0
    for a in range(n + 1):
        for b in range(n + 1):
            for c in range(n + 1):
                if a == 0 and b == 0 and c == 0:
                    continue
                if gcd(gcd(a, b), c) == 1:
                    count += 1
    return count

def solve_sieve(n: int):
    """Compute D(N) via Mobius sieve. O(N) time, works for N up to ~10^7."""
    mu = [0] * (n + 1)
    mu[1] = 1
    is_prime = bytearray(n + 1)
    primes = []
    for i in range(2, n + 1):
        if not is_prime[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            if i * p > n:
                break
            is_prime[i * p] = 1
            if i % p == 0:
                mu[i * p] = 0
                break
            mu[i * p] = -mu[i]

    result = 0
    for d in range(1, n + 1):
        if mu[d] == 0:
            continue
        q = n // d + 1
        result += mu[d] * (q * q * q - 1)
    return result

def solve_sublinear(n: int):
    """
    Compute D(N) in O(N^{2/3}) using sub-linear Mertens function.

    D(N) = sum_{d=1}^{N} mu(d) * ((floor(N/d)+1)^3 - 1)

    We group by q = floor(N/d) and use the Mertens function M(x) = sum_{d<=x} mu(d).
    """
    # Sieve threshold: N^{2/3}
    threshold = int(n ** (2 / 3)) + 1

    # Sieve mu and prefix sums up to threshold
    mu = [0] * (threshold + 1)
    mu[1] = 1
    is_prime = bytearray(threshold + 1)
    primes = []
    for i in range(2, threshold + 1):
        if not is_prime[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            if i * p > threshold:
                break
            is_prime[i * p] = 1
            if i % p == 0:
                mu[i * p] = 0
                break
            mu[i * p] = -mu[i]

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

    # Memoized Mertens for large values
    m_cache = {}

    def mertens(x):
        if x <= threshold:
            return m_small[x]
        if x in m_cache:
            return m_cache[x]
        result = 1
        k = 2
        while k <= x:
            q = x // k
            k_next = x // q + 1
            result -= (k_next - k) * mertens(q)
            k = k_next
        m_cache[x] = result
        return result

    # Main summation: group by q = floor(N/d)
    total = 0
    d = 1
    while d <= n:
        q = n // d
        d_max = n // q
        mu_sum = mertens(d_max) - mertens(d - 1)
        g = (q + 1) ** 3 - 1
        total += mu_sum * g
        d = d_max + 1

    return total

def format_answer(value: int) -> str:
    """Format as first 9 digits + last 9 digits."""
    s = str(value)
    return s[:9] + s[-9:]

# ---- Compute and verify ----

# Verify small cases against brute force
for test_n in [1, 2, 3, 5]:
    bf = solve_brute_force(test_n)
    sv = solve_sieve(test_n)
    assert bf == sv, f"Mismatch at N={test_n}: brute={bf}, sieve={sv}"

# Verify D(10^6) matches the given value
d_1m = solve_sieve(1_000_000)
assert d_1m == 831909254469114121, f"D(10^6) = {d_1m}, expected 831909254469114121"

# Compute D(10^10)
d_10b = solve_sublinear(10**10)
answer = format_answer(d_10b)
assert answer == "831907372805129931", f"Got {answer}, expected 831907372805129931"

print(answer)

# ---- Visualization ----