Problem 370: Geometric Triangles

Project Euler Problem 370 Solution

EulerSolve provides an optimized solution for Project Euler Problem 370, Geometric Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A geometric triangle here means an integer-sided triangle whose side lengths are in geometric progression. If the sides are ordered as $a \le b \le c$, then the condition is $b^2 = ac$. For a perimeter limit $L$, we must count all such triangles with $a+b+c \le L$. Mathematical Approach The local solution files reduce the problem to a sum over primitive parameter pairs $(p,q)$. The full argument has three layers: an exact parameterization of all integer geometric triangles, a triangle-inequality bound on the ratio $q/p$, and a scaling count for all multiples of each primitive triangle. Step 1: Parameterize Every Integer Geometric Triangle Let $(a,b,c)$ be an integer geometric triangle with $a \le b \le c$. Since the sides are in geometric progression, $$b^2 = ac.$$ Set $m=\gcd(a,c)$, and write $$a=mu,\qquad c=mv,\qquad \gcd(u,v)=1.$$ Then $b^2 = m^2uv$, so $uv$ is a square. Because $u$ and $v$ are coprime, each must itself be a square. Hence there exist coprime integers $p,q$ such that $$u=p^2,\qquad v=q^2,\qquad \gcd(p,q)=1.$$ Therefore every integer geometric triangle has the form $$a=mp^2,\qquad b=mpq,\qquad c=mq^2,$$ with $m \ge 1$, $1 \le p \le q$, and $\gcd(p,q)=1$. Conversely, every triple of this form satisfies $b^2=ac$, so the parameterization is exact and unique once $p \le q$ and $\gcd(p,q)=1$ are imposed....

Detailed mathematical approach

Problem Summary

A geometric triangle here means an integer-sided triangle whose side lengths are in geometric progression. If the sides are ordered as $a \le b \le c$, then the condition is $b^2 = ac$. For a perimeter limit $L$, we must count all such triangles with $a+b+c \le L$.

Mathematical Approach

The local solution files reduce the problem to a sum over primitive parameter pairs $(p,q)$. The full argument has three layers: an exact parameterization of all integer geometric triangles, a triangle-inequality bound on the ratio $q/p$, and a scaling count for all multiples of each primitive triangle.

Step 1: Parameterize Every Integer Geometric Triangle

Let $(a,b,c)$ be an integer geometric triangle with $a \le b \le c$. Since the sides are in geometric progression,

$$b^2 = ac.$$

Set $m=\gcd(a,c)$, and write

$$a=mu,\qquad c=mv,\qquad \gcd(u,v)=1.$$

Then $b^2 = m^2uv$, so $uv$ is a square. Because $u$ and $v$ are coprime, each must itself be a square. Hence there exist coprime integers $p,q$ such that

$$u=p^2,\qquad v=q^2,\qquad \gcd(p,q)=1.$$

Therefore every integer geometric triangle has the form

$$a=mp^2,\qquad b=mpq,\qquad c=mq^2,$$

with $m \ge 1$, $1 \le p \le q$, and $\gcd(p,q)=1$. Conversely, every triple of this form satisfies $b^2=ac$, so the parameterization is exact and unique once $p \le q$ and $\gcd(p,q)=1$ are imposed.

Step 2: Triangle Inequality Gives the Golden-Ratio Bound

The only nontrivial triangle inequality is $c \lt a+b$. Substituting the parameterization gives

$$mq^2 \lt mp^2 + mpq \iff q^2 \lt p^2 + pq.$$

Divide by $p^2$ and set $x=q/p$. Then

$$x^2 - x - 1 \lt 0.$$

The positive root is the golden ratio

$$\varphi = \frac{1+\sqrt{5}}{2},$$

so admissible ratios satisfy

$$1 \le \frac{q}{p} \lt \varphi,\qquad q \lt \varphi p.$$

This is why the code computes an integer upper bound `q_max_triangle(p)` by starting from $\varphi p$ and then correcting the boundary with exact arithmetic.

Step 3: Perimeter Bound and the Primitive Contribution

The perimeter of the triangle $(mp^2,mpq,mq^2)$ is

$$P = m(p^2 + pq + q^2).$$

For a fixed primitive pair $(p,q)$, the allowed scale factors $m$ are exactly

$$1 \le m \le \left\lfloor\frac{L}{p^2+pq+q^2}\right\rfloor.$$

So the primitive pair $(p,q)$ contributes

$$\left\lfloor\frac{L}{p^2+pq+q^2}\right\rfloor$$

different triangles. Summing over all primitive ratios yields the counting formula used by the solver:

$$\boxed{G(L)=\sum_{\substack{p\ge 1\\ q\ge p\\ \gcd(p,q)=1\\ q^2 \lt p^2+pq}} \left\lfloor\frac{L}{p^2+pq+q^2}\right\rfloor.}$$

Step 4: Finite Search Bounds

Because $q \ge p$, the smallest primitive perimeter for a fixed $p$ occurs at $q=p$, namely $p^2+pq+q^2=3p^2$. Therefore no contribution is possible once

$$3p^2 \gt L,$$

so the outer loop only needs

$$p \le p_{\max} = \left\lfloor\sqrt{\frac{L}{3}}\right\rfloor.$$

For fixed $p$, the perimeter inequality

$$p^2+pq+q^2 \le L$$

is quadratic in $q$. Solving it gives

$$q \le \left\lfloor\frac{\sqrt{4L-3p^2}-p}{2}\right\rfloor.$$

The implementation names this second bound `q_max_for_norm(p, L)`, and finally uses

$$q_{\max}(p)=\min\!\left(\left\lfloor \varphi p \right\rfloor,\left\lfloor\frac{\sqrt{4L-3p^2}-p}{2}\right\rfloor\right).$$

Worked Example: $L=100$

The primitive pairs $(p,q)$ that survive both bounds are

$$ (1,1),\ (2,3),\ (3,4),\ (4,5),\ (5,6). $$

Their primitive perimeters are

$$3,\ 19,\ 37,\ 61,\ 91.$$

Hence

$$G(100)=\left\lfloor\frac{100}{3}\right\rfloor+\left\lfloor\frac{100}{19}\right\rfloor+\left\lfloor\frac{100}{37}\right\rfloor+\left\lfloor\frac{100}{61}\right\rfloor+\left\lfloor\frac{100}{91}\right\rfloor=33+5+2+1+1=42.$$

This matches the statement checkpoint $G(100)=42$ that the C++ file verifies automatically.

Step 5: Counting Coprime $q$ by Möbius Inversion

A naive implementation would inspect every $q$ individually. The optimized solver instead counts coprime values in an interval $[A,B]$ by inclusion-exclusion over the prime divisors of $p$:

$$C_p(A,B)=\left|\{q\in[A,B]:\gcd(p,q)=1\}\right|=\sum_{d\mid p}\mu(d)\left(\left\lfloor\frac{B}{d}\right\rfloor-\left\lfloor\frac{A-1}{d}\right\rfloor\right).$$

To evaluate this quickly, the code factors $p$ with a smallest-prime-factor sieve, extracts the distinct primes of $p$, and enumerates the squarefree divisors together with their Möbius signs $\mu(d)\in\{-1,1\}$. For very short intervals, both the C++ and Java implementations switch to direct $\gcd$ tests because that is cheaper than running the divisor sum.

Step 6: Quotient Blocking

For fixed $p$, define the primitive perimeter form

$$N(p,q)=p^2+pq+q^2,\qquad t(q)=\left\lfloor\frac{L}{N(p,q)}\right\rfloor.$$

Since $N(p,q)$ is increasing in $q$, the quotient $t(q)$ stays constant on contiguous ranges of $q$. If a block starts at $q=\ell$, the code sets

$$t=\left\lfloor\frac{L}{N(p,\ell)}\right\rfloor,\qquad M=\left\lfloor\frac{L}{t}\right\rfloor.$$

Then the same quotient holds precisely while $N(p,q)\le M$. So the block ends at the largest $r$ satisfying that inequality, which is found by another call to `q_max_for_norm(p, M)`. The whole block contributes

$$t\cdot C_p(\ell,r).$$

This quotient-blocking step is the main reason the optimized solution avoids a full scan over all $q$.

How the Code Works

The C++ file is the main optimized implementation. It precomputes smallest prime factors up to $p_{\max}=\lfloor\sqrt{L/3}\rfloor$, processes each $p$, and optionally parallelizes independent $p$-values across threads. It also validates the mathematics by checking the statement values $G(100)=42$, $G(1000)=532$, $G(10^4)=6427$, $G(10^5)=75243$, $G(10^6)=861805$, and by comparing fast and brute-force counts on several small limits. The Java file implements the same mathematics in a simpler single-threaded form with `long`. The Python file is intentionally a bridge: it compiles and runs the C++ solver so that Python returns exactly the same validated result.

Complexity Analysis

The sieve up to $p_{\max}$ costs $O(p_{\max}\log\log p_{\max})$ time and $O(p_{\max})$ memory. A naive double loop over all admissible $(p,q)$ pairs is essentially linear in $L$, because $p$ ranges up to $O(\sqrt L)$ and each $p$ has $O(p)$ possible $q$-values. The implemented solver replaces that raw scan by quotient blocks and Möbius-weighted interval counts. The exact number of blocks depends on the floor-division structure, so the sharp asymptotic is subtle, but in practice it is dramatically smaller than the naive scan, and each block touches only the squarefree divisors of one $p$.

Footnotes and References

Problem page: https://projecteuler.net/problem=370
Geometric progression: Wikipedia — Geometric progression
Golden ratio: Wikipedia — Golden ratio
Möbius inversion formula: Wikipedia — Möbius inversion formula
Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
Smallest-prime-factor sieve: Wikipedia — Sieve of Eratosthenes

Problem 370 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <thread>
#include <vector>

namespace {

using i64 = std::int64_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

constexpr u64 kDefaultLimit = 25'000'000'000'000ULL;
constexpr u64 kThreadConsistencyLimit = 100'000'000'000ULL;

struct Checkpoint {
    u64 limit;
    u64 expected;
    const char* label;
};

constexpr std::array<Checkpoint, 5> kStatementCheckpoints = {{
    {100ULL, 42ULL, "G(100)"},
    {1'000ULL, 532ULL, "G(1000)"},
    {10'000ULL, 6'427ULL, "G(10^4)"},
    {100'000ULL, 75'243ULL, "G(10^5)"},
    {1'000'000ULL, 861'805ULL, "G(10^6)"},
}};

struct Options {
    u64 limit = kDefaultLimit;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

u64 isqrt_u128(const u128 x) {
    long double root = std::sqrt(static_cast<long double>(x));
    u64 value = static_cast<u64>(root);

    while (static_cast<u128>(value + 1ULL) * static_cast<u128>(value + 1ULL) <= x) {
        ++value;
    }
    while (static_cast<u128>(value) * static_cast<u128>(value) > x) {
        --value;
    }
    return value;
}

u64 isqrt_u64(const u64 x) {
    return isqrt_u128(static_cast<u128>(x));
}

std::string to_string_u128(u128 value) {
    if (value == 0U) {
        return "0";
    }

    std::string digits;
    while (value > 0U) {
        const int digit = static_cast<int>(value % 10U);
        digits.push_back(static_cast<char>('0' + digit));
        value /= 10U;
    }

    std::reverse(digits.begin(), digits.end());
    return digits;
}

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(ch - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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 (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--limit=", parsed_u64)) {
            options.limit = parsed_u64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.limit < 3ULL) {
        std::cerr << "--limit must be at least 3.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const u64 work_items) {
    if (!allow_multithreading || work_items <= 1ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    threads = std::min<unsigned>(threads, static_cast<unsigned>(work_items));
    return std::max(1U, threads);
}

class GeometricTriangleCounter {
public:
    explicit GeometricTriangleCounter(const u64 max_limit)
        : p_max_capacity_(isqrt_u64(max_limit / 3ULL)) {
        build_smallest_prime_factors(static_cast<int>(p_max_capacity_));
    }

    u128 count(const u64 limit, const unsigned requested_threads) const {
        const u64 p_max = isqrt_u64(limit / 3ULL);
        if (p_max == 0ULL) {
            return 0U;
        }
        if (p_max > p_max_capacity_) {
            return 0U;
        }

        const unsigned threads =
            choose_thread_count(true, requested_threads, p_max);

        if (threads == 1U) {
            return count_range(limit, 1ULL, p_max);
        }

        std::atomic<u64> next_p{1ULL};
        std::vector<u128> partial(threads, 0U);
        std::vector<std::thread> workers;
        workers.reserve(threads);

        for (unsigned tid = 0U; tid < threads; ++tid) {
            workers.emplace_back([&, tid]() {
                u128 local = 0U;
                while (true) {
                    const u64 p = next_p.fetch_add(1ULL, std::memory_order_relaxed);
                    if (p > p_max) {
                        break;
                    }
                    local += count_for_single_p(limit, p);
                }
                partial[tid] = local;
            });
        }

        for (std::thread& worker : workers) {
            worker.join();
        }

        u128 total = 0U;
        for (const u128 part : partial) {
            total += part;
        }
        return total;
    }

    static u64 max_q_triangle(const u64 p) {
        return q_max_triangle(p);
    }

    static u64 max_q_for_norm(const u64 p, const u64 limit) {
        return q_max_for_norm(p, limit);
    }

private:
    std::vector<int> smallest_prime_factor_;
    u64 p_max_capacity_ = 0ULL;

    static bool triangle_ok(const u64 p, const u64 q) {
        return static_cast<u128>(q) * static_cast<u128>(q) <
               static_cast<u128>(p) * static_cast<u128>(p) +
                   static_cast<u128>(p) * static_cast<u128>(q);
    }

    static bool norm_ok(const u64 p, const u64 q, const u64 limit) {
        return static_cast<u128>(q) * static_cast<u128>(q) +
                   static_cast<u128>(p) * static_cast<u128>(q) +
                   static_cast<u128>(p) * static_cast<u128>(p) <=
               static_cast<u128>(limit);
    }

    static long double phi_value() {
        static const long double kPhi =
            (1.0L + std::sqrt(5.0L)) / 2.0L;
        return kPhi;
    }

    static u64 q_max_triangle(const u64 p) {
        u64 q = static_cast<u64>(phi_value() * static_cast<long double>(p));
        while (q > 0ULL && !triangle_ok(p, q)) {
            --q;
        }
        while (triangle_ok(p, q + 1ULL)) {
            ++q;
        }
        return q;
    }

    static u64 q_max_for_norm(const u64 p, const u64 limit) {
        if (static_cast<u128>(3U) * static_cast<u128>(p) * static_cast<u128>(p) >
            static_cast<u128>(limit)) {
            return 0ULL;
        }

        const u128 discriminant = static_cast<u128>(4U) * static_cast<u128>(limit) -
                                  static_cast<u128>(3U) * static_cast<u128>(p) *
                                      static_cast<u128>(p);
        const u64 root = isqrt_u128(discriminant);

        i64 q = (static_cast<i64>(root) - static_cast<i64>(p)) / 2LL;
        if (q < 0LL) {
            q = 0LL;
        }

        while (q > 0LL && !norm_ok(p, static_cast<u64>(q), limit)) {
            --q;
        }
        while (norm_ok(p, static_cast<u64>(q) + 1ULL, limit)) {
            ++q;
        }
        return static_cast<u64>(q);
    }

    void build_smallest_prime_factors(const int limit) {
        smallest_prime_factor_.assign(static_cast<std::size_t>(limit) + 1ULL, 0);
        if (limit >= 1) {
            smallest_prime_factor_[1] = 1;
        }

        for (int i = 2; i <= limit; ++i) {
            if (smallest_prime_factor_[static_cast<std::size_t>(i)] == 0) {
                smallest_prime_factor_[static_cast<std::size_t>(i)] = i;
                if (static_cast<i64>(i) * static_cast<i64>(i) <= limit) {
                    for (int j = i * i; j <= limit; j += i) {
                        if (smallest_prime_factor_[static_cast<std::size_t>(j)] == 0) {
                            smallest_prime_factor_[static_cast<std::size_t>(j)] = i;
                        }
                    }
                }
            }
        }
    }

    int extract_distinct_primes(u64 value, std::array<u32, 8>& primes) const {
        int count = 0;
        while (value > 1ULL) {
            const u32 prime = static_cast<u32>(
                smallest_prime_factor_[static_cast<std::size_t>(value)]);
            primes[static_cast<std::size_t>(count)] = prime;
            ++count;
            while (value % static_cast<u64>(prime) == 0ULL) {
                value /= static_cast<u64>(prime);
            }
        }
        return count;
    }

    static int build_squarefree_divisors_mu(const std::array<u32, 8>& primes,
                                            const int prime_count,
                                            std::array<u32, 128>& divisors,
                                            std::array<int8_t, 128>& mus) {
        divisors[0] = 1U;
        mus[0] = 1;
        int count = 1;

        for (int i = 0; i < prime_count; ++i) {
            const int current = count;
            const u32 p = primes[static_cast<std::size_t>(i)];
            for (int j = 0; j < current; ++j) {
                divisors[static_cast<std::size_t>(count)] =
                    divisors[static_cast<std::size_t>(j)] * p;
                mus[static_cast<std::size_t>(count)] =
                    static_cast<int8_t>(-mus[static_cast<std::size_t>(j)]);
                ++count;
            }
        }

        return count;
    }

    static u64 count_coprime_in_interval(
        const u64 p,
        const u64 left,
        const u64 right,
        const std::array<u32, 128>& divisors,
        const std::array<int8_t, 128>& mus,
        const int divisor_count) {
        if (left > right) {
            return 0ULL;
        }

        const u64 span = right - left + 1ULL;
        if (span <= 8ULL) {
            u64 count = 0ULL;
            for (u64 q = left; q <= right; ++q) {
                if (std::gcd(p, q) == 1ULL) {
                    ++count;
                }
            }
            return count;
        }

        i64 total = 0LL;
        for (int i = 0; i < divisor_count; ++i) {
            const u64 d = static_cast<u64>(divisors[static_cast<std::size_t>(i)]);
            total += static_cast<i64>(mus[static_cast<std::size_t>(i)]) *
                     static_cast<i64>(right / d - (left - 1ULL) / d);
        }
        return static_cast<u64>(total);
    }

    u128 count_for_single_p(const u64 limit, const u64 p) const {
        const u64 q_start = p;
        const u64 q_end =
            std::min(q_max_triangle(p), q_max_for_norm(p, limit));
        if (q_end < q_start) {
            return 0U;
        }

        std::array<u32, 8> primes{};
        std::array<u32, 128> divisors{};
        std::array<int8_t, 128> mus{};

        const int prime_count = extract_distinct_primes(p, primes);
        const int divisor_count =
            build_squarefree_divisors_mu(primes, prime_count, divisors, mus);

        u128 total = 0U;
        u64 left = q_start;

        while (left <= q_end) {
            const u64 norm_left = static_cast<u64>(
                static_cast<u128>(left) * static_cast<u128>(left) +
                static_cast<u128>(p) * static_cast<u128>(left) +
                static_cast<u128>(p) * static_cast<u128>(p));
            const u64 quotient = limit / norm_left;
            if (quotient == 0ULL) {
                break;
            }

            const u64 norm_limit = limit / quotient;
            const u64 right =
                std::min(q_end, q_max_for_norm(p, norm_limit));

            const u64 coprime_count = count_coprime_in_interval(
                p, left, right, divisors, mus, divisor_count);
            total += static_cast<u128>(coprime_count) *
                     static_cast<u128>(quotient);

            left = right + 1ULL;
        }

        return total;
    }

    u128 count_range(const u64 limit, const u64 start_p, const u64 end_p) const {
        u128 total = 0U;
        for (u64 p = start_p; p <= end_p; ++p) {
            total += count_for_single_p(limit, p);
        }
        return total;
    }
};

u64 brute_force_count(const u64 limit) {
    const u64 p_max = isqrt_u64(limit / 3ULL);
    u64 total = 0ULL;

    for (u64 p = 1ULL; p <= p_max; ++p) {
        const u64 q_end =
            std::min(GeometricTriangleCounter::max_q_triangle(p),
                     GeometricTriangleCounter::max_q_for_norm(p, limit));
        if (q_end < p) {
            continue;
        }

        for (u64 q = p; q <= q_end; ++q) {
            if (std::gcd(p, q) != 1ULL) {
                continue;
            }

            const u64 norm = static_cast<u64>(
                static_cast<u128>(p) * static_cast<u128>(p) +
                static_cast<u128>(p) * static_cast<u128>(q) +
                static_cast<u128>(q) * static_cast<u128>(q));
            total += limit / norm;
        }
    }

    return total;
}

bool check_equal_u128(const char* label, const u128 got, const u128 expected) {
    if (got != expected) {
        std::cerr << "Checkpoint failed for " << label << ": expected "
                  << to_string_u128(expected) << ", got " << to_string_u128(got)
                  << '\n';
        return false;
    }
    std::cout << "Checkpoint OK: " << label << " = " << to_string_u128(got)
              << '\n';
    return true;
}

bool run_checkpoints(const GeometricTriangleCounter& counter,
                     const unsigned threads) {
    for (const Checkpoint checkpoint : kStatementCheckpoints) {
        const u128 got = counter.count(checkpoint.limit, 1U);
        if (!check_equal_u128(checkpoint.label, got,
                              static_cast<u128>(checkpoint.expected))) {
            return false;
        }
    }

    constexpr std::array<u64, 6> kBruteLimits = {
        3ULL, 17ULL, 111ULL, 999ULL, 12'345ULL, 54'321ULL};
    for (const u64 limit : kBruteLimits) {
        const u128 fast = counter.count(limit, 1U);
        const u128 brute = static_cast<u128>(brute_force_count(limit));
        const std::string label = "fast-vs-brute G(" + std::to_string(limit) + ")";
        if (!check_equal_u128(label.c_str(), fast, brute)) {
            return false;
        }
    }

    if (threads > 1U) {
        const u128 single = counter.count(kThreadConsistencyLimit, 1U);
        const u128 multi = counter.count(kThreadConsistencyLimit, threads);
        if (!check_equal_u128("thread consistency",
                              multi,
                              single)) {
            return false;
        }
    }

    std::cout << "Validation checkpoints passed.\n";
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    u64 max_limit_for_solver = options.limit;
    if (options.run_checkpoints) {
        for (const Checkpoint checkpoint : kStatementCheckpoints) {
            max_limit_for_solver = std::max(max_limit_for_solver, checkpoint.limit);
        }
        max_limit_for_solver =
            std::max(max_limit_for_solver, kThreadConsistencyLimit);
    }

    const GeometricTriangleCounter counter(max_limit_for_solver);
    const u64 p_max = isqrt_u64(options.limit / 3ULL);
    const unsigned threads =
        choose_thread_count(options.allow_multithreading,
                            options.requested_threads,
                            p_max);

    if (options.run_checkpoints) {
        if (!run_checkpoints(counter, threads)) {
            return 1;
        }
    }

    const u128 answer = counter.count(options.limit, threads);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.*;

public class Euler370 {
    public static String solve() {
        long LIMIT = 25_000_000_000_000L;
        int pMax = (int) Math.sqrt(LIMIT / 3.0);
        double phi = (1.0 + Math.sqrt(5.0)) / 2.0;

        int[] spf = new int[pMax + 1];
        for (int i = 0; i <= pMax; i++)
            spf[i] = i;
        for (int i = 2; (long) i * i <= pMax; i++)
            if (spf[i] == i)
                for (int j = i * i; j <= pMax; j += i)
                    if (spf[j] == j)
                        spf[j] = i;

        long total = 0;
        for (int p = 1; p <= pMax; p++) {
            int qTriMax = qMaxTri(p, phi);
            int qNormMax = qMaxNorm(p, LIMIT);
            int qEnd = Math.min(qTriMax, qNormMax);
            if (qEnd < p)
                continue;

            int[] dprimes = distinctPrimes(p, spf);
            int nd = dprimes.length;
            int divCount = 1 << nd;
            int[] divs = new int[divCount];
            int[] mus = new int[divCount];
            for (int mask = 0; mask < divCount; mask++) {
                int d = 1, mu = 1;
                for (int i = 0; i < nd; i++)
                    if ((mask & (1 << i)) != 0) {
                        d *= dprimes[i];
                        mu = -mu;
                    }
                divs[mask] = d;
                mus[mask] = mu;
            }

            int left = p;
            while (left <= qEnd) {
                long normLeft = (long) left * left + (long) p * left + (long) p * p;
                long quotient = LIMIT / normLeft;
                if (quotient == 0)
                    break;
                long normLim = LIMIT / quotient;
                int right = Math.min(qEnd, qMaxNorm(p, normLim));

                int cop;
                int span = right - left + 1;
                if (span <= 8) {
                    cop = 0;
                    for (int q = left; q <= right; q++)
                        if (gcd(p, q) == 1)
                            cop++;
                } else {
                    cop = 0;
                    for (int i = 0; i < divCount; i++)
                        cop += mus[i] * (right / divs[i] - (left - 1) / divs[i]);
                }
                total += cop * quotient;
                left = right + 1;
            }
        }
        return String.valueOf(total);
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static int qMaxTri(int p, double phi) {
        int q = (int) (phi * p);
        while (q > 0 && (long) q * q >= (long) p * p + (long) p * q)
            q--;
        while ((long) (q + 1) * (q + 1) < (long) p * p + (long) p * (q + 1))
            q++;
        return q;
    }

    static int qMaxNorm(int p, long limit) {
        if (3L * p * p > limit)
            return 0;
        long disc = 4 * limit - 3L * p * p;
        long root = (long) Math.sqrt((double) disc);
        while ((root + 1) * (root + 1) <= disc)
            root++;
        while (root * root > disc)
            root--;
        int q = (int) ((root - p) / 2);
        if (q < 0)
            q = 0;
        while (q > 0 && (long) q * q + (long) p * q + (long) p * p > limit)
            q--;
        while ((long) (q + 1) * (q + 1) + (long) p * (q + 1) + (long) p * p <= limit)
            q++;
        return q;
    }

    static int[] distinctPrimes(int val, int[] spf) {
        List<Integer> ps = new ArrayList<>();
        while (val > 1) {
            int p = spf[val];
            ps.add(p);
            while (val % p == 0)
                val /= p;
        }
        return ps.stream().mapToInt(Integer::intValue).toArray();
    }

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