Problem 373: Circumscribed Circles

Project Euler Problem 373 Solution

EulerSolve provides an optimized solution for Project Euler Problem 373, Circumscribed Circles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer circumradius $r$, let $T(r)$ be the number of triangles with integer side lengths and circumradius exactly $r$. The program computes $$S(N)=\sum_{r=1}^{N} r\,T(r),\qquad N=10^7.$$ A brute-force search over all triples $(a,b,c)$ up to $2r$ would be far too slow, so the implementation first characterizes which side lengths can occur for a fixed radius and then checks only those candidates exactly. Mathematical Approach Step 1: Convert the circumradius condition into an exact integer identity For a triangle with sides $a,b,c$ and area $\Delta$, Heron's formula gives $$16\Delta^2=(a+b+c)(a+b-c)(a-b+c)(-a+b+c).$$ The circumradius also satisfies $$R=\frac{abc}{4\Delta}.$$ Squaring and eliminating $\Delta$ yields $$a^2b^2c^2=R^2(a+b+c)(a+b-c)(a-b+c)(-a+b+c).$$ This is exactly the predicate used in the solver. It avoids floating-point roundoff entirely: once $a$, $b$, $c$, and $r$ are integers, the question “is the circumradius equal to $r$?” becomes a pure integer comparison. Step 2: Describe every possible side length for a fixed radius If side $a$ is opposite angle $A$, then the extended law of sines gives $$a=2r\sin A.$$ Because $a$ and $r$ are integers, $\sin A=a/(2r)$ is rational....

Detailed mathematical approach

Problem Summary

For each integer circumradius $r$, let $T(r)$ be the number of triangles with integer side lengths and circumradius exactly $r$. The program computes

$$S(N)=\sum_{r=1}^{N} r\,T(r),\qquad N=10^7.$$

A brute-force search over all triples $(a,b,c)$ up to $2r$ would be far too slow, so the implementation first characterizes which side lengths can occur for a fixed radius and then checks only those candidates exactly.

Mathematical Approach

Step 1: Convert the circumradius condition into an exact integer identity

For a triangle with sides $a,b,c$ and area $\Delta$, Heron's formula gives

$$16\Delta^2=(a+b+c)(a+b-c)(a-b+c)(-a+b+c).$$

The circumradius also satisfies

$$R=\frac{abc}{4\Delta}.$$

Squaring and eliminating $\Delta$ yields

$$a^2b^2c^2=R^2(a+b+c)(a+b-c)(a-b+c)(-a+b+c).$$

This is exactly the predicate used in the solver. It avoids floating-point roundoff entirely: once $a$, $b$, $c$, and $r$ are integers, the question “is the circumradius equal to $r$?” becomes a pure integer comparison.

Step 2: Describe every possible side length for a fixed radius

If side $a$ is opposite angle $A$, then the extended law of sines gives

$$a=2r\sin A.$$

Because $a$ and $r$ are integers, $\sin A=a/(2r)$ is rational. Rational points on the unit circle are parameterized by primitive Pythagorean data: for coprime integers $m>n\ge 1$ with opposite parity,

$$\left(\frac{m^2-n^2}{m^2+n^2},\frac{2mn}{m^2+n^2}\right)$$

is a primitive rational point on $x^2+y^2=1$. Therefore every rational sine in $(0,1)$ can be written as one of

$$\sin A\in\left\{\frac{m^2-n^2}{m^2+n^2},\frac{2mn}{m^2+n^2}\right\}.$$

If we write $d=m^2+n^2$ and $u$ for one of the two numerators, then

$$a=2r\frac{u}{d}.$$

So $a$ is guaranteed to be integral whenever $d\mid r$. The special value $\sin A=1$ is not produced by finite $(m,n)$, so the code inserts the diameter case $a=2r$ explicitly.

Step 3: Precompute useful denominators once, then use only divisors of $r$

The C++ solver builds a table indexed by $d\le N$. For each primitive pair $(m,n)$ with $d=m^2+n^2\le N$, it stores the admissible numerators

$$u\in\{m^2-n^2,\ 2mn\}.$$

For a fixed radius $r$, only denominators dividing $r$ can contribute integer sides. If $r=kd$, then each stored numerator produces the candidate side

$$a=2ku.$$

This is why the implementation first factors $r$ using a smallest-prime-factor sieve, enumerates all divisors of $r$, and looks up only those divisor entries. The resulting candidate set is

$$\mathcal{S}_r=\{2r\}\cup \left\{2\frac{r}{d}u : d\mid r,\ u\in U(d)\right\},$$

where $U(d)$ is the precomputed list of numerators attached to denominator $d$. The list is sorted and deduplicated because distinct parameter pairs can lead to the same side length.

Step 4: Count valid triples from the candidate side set

After building $\mathcal{S}_r$, the solver enumerates

$$a\le b\le c,\qquad a,b,c\in \mathcal{S}_r.$$

The sorted order gives an immediate pruning rule: as soon as $a+b\le c$, the triangle inequality fails and the inner loop can stop for that pair $(a,b)$. For each remaining triple, the solver applies the exact identity from Step 1. Every triple that passes contributes $1$ to $T(r)$.

A simple example is $r=5$. The primitive pair $(m,n)=(2,1)$ gives $d=5$ and numerators $3$ and $4$. Because $5\mid r$, the candidate sides include

$$2\cdot \frac{5}{5}\cdot 3=6,\qquad 2\cdot \frac{5}{5}\cdot 4=8,\qquad 2r=10.$$

So $(6,8,10)$ appears naturally, and the exact check confirms that its circumradius is indeed $5$.

How the Code Works

The class Euler373Solver precomputes two structures: a smallest-prime-factor array for fast divisor generation, and a denominator-to-numerators table for rational-sine candidates. The method count_triangles_for_radius(r) builds $\mathcal{S}_r$, loops over ordered triples, and calls circumradius_equals(a,b,c,r) for the exact test.

The C++ implementation uses a custom 192-bit accumulator U192 because the products in the circumradius identity are too large for ordinary 64-bit arithmetic. The outer sum over radii is embarrassingly parallel, so the solver distributes radii across worker threads with an atomic counter and combines per-thread partial sums. The Python and Java files are thin bridges that compile and run this same C++ solver, so all three language solutions share identical mathematics and validation.

Before solving the full problem, the program checks the optimized method against brute force for $r\le 60$ and also verifies the checkpoints

$$S(100)=4950,\qquad S(1200)=1653605.$$

Complexity Analysis

Let $N$ be the radius limit. Building the smallest-prime-factor table costs $O(N\log\log N)$ time and $O(N)$ memory. The denominator table iterates over primitive pairs with $m^2+n^2\le N$, which is a one-time precomputation.

For one radius $r$, divisor enumeration is proportional to the number of divisors of $r$, and the dominant work is the triple scan over the deduplicated candidate set $\mathcal{S}_r$. In the worst case this stage is $O(|\mathcal{S}_r|^3)$, but in practice it is kept manageable because only divisors of $r$ contribute candidates, duplicates are removed, and the triangle inequality stops many inner-loop iterations early. The final summation over radii is parallelized.

References

Problem page: https://projecteuler.net/problem=373
Circumradius formulas and the extended law of sines: Wikipedia — Circumscribed circle
Heron's formula: Wikipedia — Heron's formula
Primitive Pythagorean triples and rational points on the unit circle: Wikipedia — Pythagorean triple

Problem 373 source code

C++

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

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr int kTargetRadius = 10'000'000;

struct U192 {
    u64 lo = 0;
    u64 mid = 0;
    u64 hi = 0;

    static U192 one() {
        U192 value;
        value.lo = 1;
        return value;
    }

    void multiply_small(u64 factor) {
        const u128 t0 = static_cast<u128>(lo) * factor;
        lo = static_cast<u64>(t0);
        u64 carry = static_cast<u64>(t0 >> 64);

        const u128 t1 = static_cast<u128>(mid) * factor + carry;
        mid = static_cast<u64>(t1);
        carry = static_cast<u64>(t1 >> 64);

        const u128 t2 = static_cast<u128>(hi) * factor + carry;
        hi = static_cast<u64>(t2);
    }

    bool operator==(const U192& other) const {
        return lo == other.lo && mid == other.mid && hi == other.hi;
    }
};

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

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

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

bool circumradius_equals(int a, int b, int c, int radius) {
    // Exact check of:
    // a^2 b^2 c^2 == R^2 (a+b+c)(a+b-c)(a-b+c)(-a+b+c)
    U192 lhs = U192::one();
    lhs.multiply_small(static_cast<u64>(radius));
    lhs.multiply_small(static_cast<u64>(radius));
    lhs.multiply_small(static_cast<u64>(a + b + c));
    lhs.multiply_small(static_cast<u64>(a + b - c));
    lhs.multiply_small(static_cast<u64>(a - b + c));
    lhs.multiply_small(static_cast<u64>(-a + b + c));

    U192 rhs = U192::one();
    rhs.multiply_small(static_cast<u64>(a));
    rhs.multiply_small(static_cast<u64>(a));
    rhs.multiply_small(static_cast<u64>(b));
    rhs.multiply_small(static_cast<u64>(b));
    rhs.multiply_small(static_cast<u64>(c));
    rhs.multiply_small(static_cast<u64>(c));

    return lhs == rhs;
}

class Euler373Solver {
public:
    explicit Euler373Solver(int radius_limit)
        : radius_limit_(radius_limit),
          smallest_prime_factor_(static_cast<std::size_t>(radius_limit) + 1),
          denominator_to_numerators_(static_cast<std::size_t>(radius_limit) + 1) {
        build_smallest_prime_factor();
        build_denominator_table();
    }

    u128 solve(bool allow_multithreading, unsigned requested_threads = 0) const {
        return sum_up_to(radius_limit_, allow_multithreading, requested_threads);
    }

    u128 sum_up_to(int max_radius,
                   bool allow_multithreading,
                   unsigned requested_threads = 0) const {
        if (max_radius > radius_limit_) {
            return 0;
        }

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

        if (!allow_multithreading || threads <= 1 || max_radius < 50'000) {
            threads = 1;
        } else {
            threads = std::min<unsigned>(threads, static_cast<unsigned>(max_radius));
        }

        if (threads == 1) {
            std::vector<int> divisors;
            std::vector<int> sides;
            u128 total = 0;
            for (int radius = 1; radius <= max_radius; ++radius) {
                total += static_cast<u128>(radius) *
                         count_triangles_for_radius(radius, divisors, sides);
            }
            return total;
        }

        std::atomic<int> next_radius{1};
        std::vector<u128> partial(threads, 0);
        std::vector<std::thread> workers;
        workers.reserve(threads);

        for (unsigned tid = 0; tid < threads; ++tid) {
            workers.emplace_back([&, tid]() {
                std::vector<int> divisors;
                std::vector<int> sides;
                u128 local = 0;

                while (true) {
                    const int radius = next_radius.fetch_add(1, std::memory_order_relaxed);
                    if (radius > max_radius) {
                        break;
                    }

                    local += static_cast<u128>(radius) *
                             count_triangles_for_radius(radius, divisors, sides);
                }

                partial[tid] = local;
            });
        }

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

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

    u64 count_triangles_for_radius(int radius) const {
        std::vector<int> divisors;
        std::vector<int> sides;
        return count_triangles_for_radius(radius, divisors, sides);
    }

private:
    int radius_limit_ = 0;
    std::vector<int> smallest_prime_factor_;
    std::vector<std::vector<int>> denominator_to_numerators_;

    void build_smallest_prime_factor() {
        for (int value = 0; value <= radius_limit_; ++value) {
            smallest_prime_factor_[static_cast<std::size_t>(value)] = value;
        }
        if (radius_limit_ >= 1) {
            smallest_prime_factor_[1] = 1;
        }

        for (int p = 2; static_cast<std::int64_t>(p) * p <= radius_limit_; ++p) {
            if (smallest_prime_factor_[static_cast<std::size_t>(p)] != p) {
                continue;
            }

            for (int multiple = p * p; multiple <= radius_limit_; multiple += p) {
                if (smallest_prime_factor_[static_cast<std::size_t>(multiple)] == multiple) {
                    smallest_prime_factor_[static_cast<std::size_t>(multiple)] = p;
                }
            }
        }
    }

    void build_denominator_table() {
        const int m_limit = static_cast<int>(std::sqrt(static_cast<long double>(radius_limit_)));

        for (int m = 2; m <= m_limit; ++m) {
            for (int n = 1; n < m; ++n) {
                if (((m + n) & 1) == 0) {
                    continue;
                }
                if (std::gcd(m, n) != 1) {
                    continue;
                }

                const int denom = m * m + n * n;
                if (denom > radius_limit_) {
                    break;
                }

                auto& numerators = denominator_to_numerators_[static_cast<std::size_t>(denom)];
                numerators.push_back(m * m - n * n);
                numerators.push_back(2 * m * n);
            }
        }

        for (int denom = 1; denom <= radius_limit_; ++denom) {
            auto& numerators = denominator_to_numerators_[static_cast<std::size_t>(denom)];
            if (numerators.empty()) {
                continue;
            }
            std::sort(numerators.begin(), numerators.end());
            numerators.erase(std::unique(numerators.begin(), numerators.end()),
                             numerators.end());
        }
    }

    void divisors_of(int value, std::vector<int>& out) const {
        out.clear();
        out.push_back(1);

        while (value > 1) {
            const int prime = smallest_prime_factor_[static_cast<std::size_t>(value)];
            int exponent = 0;
            while (value % prime == 0) {
                value /= prime;
                ++exponent;
            }

            const std::size_t base_size = out.size();
            int multiplier = 1;
            for (int i = 1; i <= exponent; ++i) {
                multiplier *= prime;
                for (std::size_t idx = 0; idx < base_size; ++idx) {
                    out.push_back(out[idx] * multiplier);
                }
            }
        }
    }

    u64 count_triangles_for_radius(int radius,
                                   std::vector<int>& divisors,
                                   std::vector<int>& sides) const {
        divisors_of(radius, divisors);

        sides.clear();
        sides.push_back(2 * radius); // Diameter (sin(theta)=1).

        for (const int denom : divisors) {
            const auto& numerators =
                denominator_to_numerators_[static_cast<std::size_t>(denom)];
            if (numerators.empty()) {
                continue;
            }

            const int scale = 2 * (radius / denom);
            for (const int numer : numerators) {
                sides.push_back(scale * numer);
            }
        }

        std::sort(sides.begin(), sides.end());
        sides.erase(std::unique(sides.begin(), sides.end()), sides.end());

        const int count = static_cast<int>(sides.size());
        u64 triangles = 0;

        for (int i = 0; i < count; ++i) {
            const int a = sides[static_cast<std::size_t>(i)];
            for (int j = i; j < count; ++j) {
                const int b = sides[static_cast<std::size_t>(j)];
                const int ab = a + b;
                for (int k = j; k < count; ++k) {
                    const int c = sides[static_cast<std::size_t>(k)];
                    if (ab <= c) {
                        break;
                    }

                    if (circumradius_equals(a, b, c, radius)) {
                        ++triangles;
                    }
                }
            }
        }

        return triangles;
    }
};

u64 brute_force_triangle_count_for_radius(int radius) {
    const int max_side = 2 * radius;
    u64 count = 0;

    for (int a = 1; a <= max_side; ++a) {
        for (int b = a; b <= max_side; ++b) {
            const int c_max = std::min(max_side, a + b - 1);
            for (int c = b; c <= c_max; ++c) {
                if (circumradius_equals(a, b, c, radius)) {
                    ++count;
                }
            }
        }
    }

    return count;
}

bool run_validation_checkpoints() {
    {
        // Structural validation: compare the optimized candidate enumeration
        // against brute force on small radii.
        Euler373Solver solver(60);
        for (int radius = 1; radius <= 60; ++radius) {
            const u64 fast = solver.count_triangles_for_radius(radius);
            const u64 brute = brute_force_triangle_count_for_radius(radius);
            if (fast != brute) {
                std::cerr << "Validation failed for radius " << radius << ": fast=" << fast
                          << ", brute=" << brute << '\n';
                return false;
            }
        }
    }

    struct Checkpoint {
        int radius_limit;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {100, 4'950},
        {1200, 1'653'605},
    };

    for (const Checkpoint cp : checkpoints) {
        const Euler373Solver solver(cp.radius_limit);
        const u64 got = static_cast<u64>(solver.solve(false, 1));
        if (got != cp.expected) {
            std::cerr << "Checkpoint failed for S(" << cp.radius_limit << "): got " << got
                      << ", expected " << cp.expected << '\n';
            return false;
        }
    }

    return true;
}

} // namespace

int main() {
    if (!run_validation_checkpoints()) {
        return 1;
    }

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 4;
    }

    const Euler373Solver solver(kTargetRadius);
    const u128 answer = solver.solve(true, 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.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler373 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler373.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler373_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler373 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler373 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler373 C++ bridge produced empty output.");
        }
        return parsed;
    }

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