Problem 388: Distinct Lines

View on Project Euler

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

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

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