Problem 195: Inscribed Circles of Triangles with One Angle of 60 Degrees

Project Euler Problem 195 Solution

EulerSolve provides an optimized solution for Project Euler Problem 195, Inscribed Circles of Triangles with One Angle of 60 Degrees, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $T(N)$ denote the number of integer-sided triangles that contain a $60^\circ$ angle and whose incircle radius does not exceed $N$. The specific task is to evaluate $T(1{,}053{,}779)$. A brute-force search over side triples would be far too large, because the geometric condition mixes a Diophantine constraint from the $60^\circ$ angle with a second constraint on the inradius. The implementations therefore replace direct triangle enumeration by a parametrization of primitive objects and then count all admissible integer scalings. Mathematical Approach The key observation is that a $60^\circ$ triangle can be converted into a $120^\circ$ triangle with integer sides, and the latter has a classical coprime parametrization. Once that primitive parametrization is known, the inradius becomes a simple linear expression, so counting triangles reduces to summing the number of allowed scale factors. From a $60^\circ$ triangle to a $120^\circ$ triple Take a triangle whose $60^\circ$ angle lies between the integer sides $a$ and $b$, and let $c$ be the side opposite that angle. Without loss of generality order the adjacent sides so that $a \le b$. The law of cosines gives $$c^2=a^2+b^2-ab.$$ Now write $x=a$ and $y=b-a$, so $b=x+y$....

Detailed mathematical approach

Problem Summary

A triangle with whole-number sides and exactly one angle of $60^\circ$ is called a 60-degree triangle; since the angles add up to $180^\circ$, this excludes only the equilateral triangles. With $r$ the radius of the inscribed circle, $T(n)$ is the number of 60-degree triangles with $r\le n$, a triangle being determined by its three side lengths. The statement gives $T(100)=1234$, $T(1000)=22767$ and $T(10000)=359912$ and asks for $T(1053779)$. ^[1] We write $N$ for the bound on the radius, because the letter $n$ is needed for a parameter.

The whole solution in one sentence: generate the smallest whole-number triangle of each shape with a $120^\circ$ angle from a classical formula in two whole numbers $m$ and $n$, turn each into two 60-degree triangles, whose inradius is $\sqrt3/6$ or $\sqrt3/2$ times a simple product of $m$ and $n$, and count the whole-number enlargements of each that keep the inradius at most $N$.

Mathematical Approach

1. The 60-degree condition

Let the angle of $60^\circ$ lie between the sides $a$ and $b$, and let $c$ be the side opposite to it. The law of cosines with $\cos60^\circ=\tfrac12$ gives ^[2]

$$c^2=a^2+b^2-ab.$$

If $a=b$, this gives $c=a$, the equilateral triangle, so we may assume $a\lt b$. Then $c^2-a^2=b(b-a)\gt0$ and $b^2-c^2=a(b-a)\gt0$, so $a\lt c\lt b$. Conversely, whole numbers $0\lt a\lt b$ and $c$ with $c^2=a^2+b^2-ab$ always form a 60-degree triangle: $c^2$ lies strictly between $(b-a)^2$ and $(a+b)^2$, so the triangle inequality holds, and the law of cosines gives the angle $60^\circ$ between $a$ and $b$. We write such a triangle as $(a,b,c)$: the two sides at the $60^\circ$ angle, the shorter one first, then the opposite side. The smallest examples are $(3,8,7)$, because $9+64-24=49=7^2$, and $(5,8,7)$, because $25+64-40=49$.

2. The radius in whole numbers

Call the vertices $A$, $B$, $C$, with the $60^\circ$ angle at $C$, $a=BC$, $b=CA$ and $c=AB$. The inscribed circle touches each side once, and the two tangent segments from a vertex to the circle are equally long; ^[3] call these lengths $u$, $v$, $w$ at $A$, $B$, $C$. Each side consists of two tangent segments, so $a=v+w$, $b=u+w$, $c=u+v$, and therefore $w=(a+b-c)/2$. The centre $I$ of the circle lies on the bisector of the angle at $C$, and the radius to a touch point is perpendicular to the side, so $C$, $I$ and the touch point on $CA$ form a right triangle with the angle $30^\circ$ at $C$ and the legs $w$ and $r$. Hence $r=w\tan30^\circ=w/\sqrt3$, that is,

$$r=\frac{g}{2\sqrt3},\qquad g=a+b-c.$$

For $(3,8,7)$ we get $g=4$ and $r=2/\sqrt3\approx1.155$; as a check, the area $\tfrac12\cdot3\cdot8\cdot\sin60^\circ=6\sqrt3$ divided by half the perimeter, $9$, gives the same value. ^[3] For $(5,8,7)$, $g=6$ and $r=\sqrt3\approx1.732$.

Two consequences drive the whole count: $r\le N$ means $g\le2\sqrt3\,N$, and enlarging a triangle by a whole factor $k$ multiplies its sides, and with them $g$ and $r$, by $k$.

3. Cutting off an equilateral triangle

In the triangle $(a,b,c)$ mark the point $D$ on $CA$ with $CD=a$. The triangle $CDB$ has two sides $a$ around an angle of $60^\circ$, so it is equilateral and $DB=a$, and what remains is the triangle $DAB$ with the sides $x=a$, $y=b-a$, $z=c$ and the angle $120^\circ$ at $D$. Indeed, substituting $b=x+y$ into $c^2=a^2+b^2-ab$ gives

$$z^2=x^2+xy+y^2,$$

the law of cosines for $120^\circ$. Conversely, by Section 1 a solution in positive whole numbers gives the 60-degree triangles $(x,x+y,z)$ and $(y,x+y,z)$, because $x^2+(x+y)^2-x(x+y)$${}=x^2+xy+y^2$, and likewise with $x$ and $y$ exchanged; here $x\ne y$, since $3x^2=z^2$ is impossible, $\sqrt3$ being irrational. So every $120^\circ$ triangle gives exactly two 60-degree triangles, one for each order of its two shorter sides, and every 60-degree triangle $(a,b,c)$ arises exactly once, from the ordered pair $(x,y)=(a,b-a)$. For example, $(3,5,7)$, with $9+15+25=49$, gives $(3,8,7)$ and $(5,8,7)$.

A triangle is primitive if its sides have no common factor. Every triangle is a whole multiple of exactly one primitive triangle, and since $\gcd(a,b,c)=\gcd(a,b-a,c)$, a 60-degree triangle is primitive exactly when its $120^\circ$ triangle is.

4. Rational points on a curve

Dividing $z^2=x^2+xy+y^2$ by $z^2$ turns a solution into the point $(X,Y)=(x/z,y/z)$ with rational coordinates on the curve $X^2+XY+Y^2=1$. Primitive solutions and rational points with $X,Y\gt0$ correspond one to one: from a point, write both coordinates over their smallest common denominator. The curve passes through $(-1,0)$. The line $Y=s(X+1)$ through that point with slope $s$ meets the curve in a second point; substituting and dividing by $X+1$ gives

$$X=\frac{1-s^2}{1+s+s^2},\quad Y=\frac{s(2+s)}{1+s+s^2}.$$

A rational slope gives a rational point, a rational point gives the rational slope $s=Y/(X+1)$, and $X,Y\gt0$ holds exactly when $0\lt s\lt1$. This is the classical method for Pythagorean triples on the unit circle. ^[5] With $s=n/m$ in lowest terms, $m\gt n\ge1$, the point is $P(m,n)$ divided by its last entry, where

$$\begin{aligned}P(m,n)=\bigl(&m^2-n^2,\;2mn+n^2,\\ &m^2+mn+n^2\bigr).\end{aligned}$$

So every primitive solution is $P(m,n)$ divided by the greatest common divisor of its entries, for exactly one coprime pair $m\gt n\ge1$; the slope $s=1/2$ gives $P(2,1)=(3,5,7)$. These are the classical formulas for integer triangles with a $120^\circ$ angle. ^[4]

5. When the formula is primitive

Call a coprime pair $(m,n)$ with $m\gt n\ge1$ a pair of type 1 if $3$ does not divide $m-n$, and of type 3 if it does. The names are the greatest common divisors: the entries of $P(m,n)$ have no common factor for a pair of type 1, and exactly the common factor $3$ for a pair of type 3.

Proof. Let a prime $p$ divide both $x=(m-n)(m+n)$ and $z=m^2+mn+n^2$. If $p$ divides $m+n$, then $m\equiv-n\pmod p$ and $z\equiv m^2\pmod p$, so $p$ divides $m$ and then $n$, which is impossible. If $p$ divides $m-n$, then $z\equiv3m^2\pmod p$, so $p=3$, because $p$ cannot divide $m$. Hence $3$ is the only possible common prime, and only when $3\mid m-n$. In that case put $m=n+3j$: then $x=3j(m+n)$, $y=3n(n+2j)$ and $z=9j^2+3mn$ are multiples of $3$, while $z/3=3j^2+mn$ is not, because $3$ divides neither $m$ nor $n$. For example $P(4,1)=(15,9,21)=3\cdot(5,3,7)$.

6. Each primitive triangle exactly once

Exchanging the two shorter sides corresponds to a simple change of the pair. If $P(m,n)=(x,y,z)$, then

$$P(m+2n,\;m-n)=3\,(y,\;x,\;z).$$

Indeed, put $M=m+2n$ and $K=m-n$, so that $M-K=3n$, $M+K=2m+n$, $2M+K=3(m+n)$ and $MK=m^2+mn-2n^2$. The entries of $P(M,K)$ are

$$\begin{aligned}M^2-K^2&=3y,\\ 2MK+K^2&=3x,\\ M^2+MK+K^2&=3z,\end{aligned}$$

because $M^2-K^2=(M-K)(M+K)$, $2MK+K^2=K(2M+K)$ and $M^2+MK+K^2=(M-K)^2+3MK$.

If $(m,n)$ is of type 1, its partner $(m+2n,m-n)$ is coprime, since a common divisor divides $3m$ and $3n$, hence $3$, and $3\nmid m-n$; it is of type 3, since its entries differ by $3n$. Conversely, a pair $(m',n')$ of type 3 is the partner of the pair with the entries $(m'+2n')/3$ and $(m'-n')/3$, which is of type 1: these entries are whole because $m'\equiv n'\pmod3$, coprime because a common divisor would divide $m'$ and $n'$, and their difference $n'$ is not a multiple of $3$, as $m'$ and $n'$ are not both multiples of $3$. By Section 4, each of the two orders of the shorter sides of a primitive $120^\circ$ triangle comes from exactly one coprime pair, and by the identity these two pairs are partners, one of type 1 and one of type 3. Hence every primitive $120^\circ$ triangle arises from exactly one pair of type 1. ^[4] For example, $(2,1)$ and $(4,1)$ are partners, and $P(4,1)$ is three times $(3,5,7)$ with its shorter sides exchanged.

7. The two families

Let $(m,n)$ be of type 1 and $P(m,n)=(x,y,z)$. By Sections 3 and 6, the primitive 60-degree triangles are exactly the triangles $(x,x+y,z)$, family A, and $(y,x+y,z)$, family B, of these pairs, each obtained once. Their values of $g=a+b-c$ are $g_A=2x+y-z$ and $g_B=x+2y-z$, and inserting $P(m,n)$ turns them into products:

$$\begin{aligned}g_A&=(m-n)(m+2n),\\ g_B&=3mn.\end{aligned}$$

For instance, $2x+y-z=m^2+mn-2n^2$, which factors as shown. By Section 2 the radii are $r_A=\sqrt3\,(m-n)(m+2n)/6$ and $r_B=\sqrt3\,mn/2$. The pair $(2,1)$ gives back $(3,8,7)$ with $g_A=1\cdot4$ and $(5,8,7)$ with $g_B=3\cdot2$; the pair $(3,1)$, with $P(3,1)=(8,7,13)$, gives $(8,15,13)$ with $g_A=10$ and $(7,15,13)$ with $g_B=9$.

8. Counting the enlargements

The 60-degree triangles with $r\le N$ are the enlargements $k\cdot(a,b,c)$ of the primitive ones with $kg\le2\sqrt3\,N$. Writing $\lfloor u\rfloor$ for the largest whole number not above $u$, a primitive triangle with the value $g$ therefore contributes $F(g)$ triangles, and

$$\begin{aligned}T(N)&=\sum_{(m,n)}\bigl(F(g_A)+F(g_B)\bigr),\\ F(g)&=\Bigl\lfloor\frac{2\sqrt3\,N}{g}\Bigr\rfloor,\end{aligned}$$

with the sum over the pairs of type 1. The programs write $m=n+t$ with $t=m-n\ge1$. Then $\gcd(m,n)=\gcd(n,t)$, type 1 means $3\nmid t$, and

$$g_A=t\,(3n+t),\quad g_B=3n\,(n+t).$$

Each value $F(g)$ can be computed without rounding errors: for whole numbers $k,g\ge1$, both sides of $kg\le2\sqrt3\,N$ are positive, so this condition is the same as $k^2g^2\le12N^2$. The programs state this test as $3q^2k^2\le C^2$, with $q=g_A$ and $C=6N$ in family A, and with $q=n(n+t)=g_B/3$ and $C=2N$ in family B; multiplying out turns both into $k^2g^2\le12N^2$. The number of enlargements is the largest $k$ that passes, $K_C(q)=\lfloor C/(q\sqrt3)\rfloor$, which equals $F(g)$. Since $\lfloor u/g\rfloor=\lfloor\lfloor u\rfloor/g\rfloor$ for whole $g\ge1$, ^[7] also $F(g)=\lfloor L/g\rfloor$ with the single whole number $L=\lfloor2\sqrt3\,N\rfloor$.

9. Where the sums stop

A pair contributes only if $k=1$ passes, that is, if $3q^2\le C^2$, or $q\le Q$ with $Q=\lfloor C/\sqrt3\rfloor$:

$$Q_B=\Bigl\lfloor\frac{2N}{\sqrt3}\Bigr\rfloor,\quad Q_A=\Bigl\lfloor\frac{6N}{\sqrt3}\Bigr\rfloor=L.$$

In family B, $q=n(n+t)\le Q_B$ requires $n(n+1)\le Q_B$, and for such $n$ it allows $t=1,\dots,\lfloor Q_B/n\rfloor-n$. In family A, $q=t(3n+t)\le Q_A$ requires $t(t+3)\le Q_A$, and then allows $n=1,\dots,\lfloor(\lfloor Q_A/t\rfloor-t)/3\rfloor$. Rounding $Q_B/n$ and $Q_A/t$ down loses nothing, because a whole number is at most $u$ exactly when it is at most $\lfloor u\rfloor$. ^[7]

10. Worked example: $T(3)=7$

For $N=3$ we have $2\sqrt3\,N\approx10.39$, so $Q_A=10$ and $Q_B=\lfloor6/\sqrt3\rfloor=3$. Family B allows only $n=1$, because $2\cdot3\gt3$, and then $t=1,2$:

$(n,t)=(1,1)$, the pair $(m,n)=(2,1)$: $q=2$, $g_B=6$, the triangle $(5,8,7)$;
$(n,t)=(1,2)$, the pair $(3,1)$: $q=3$, $g_B=9$, the triangle $(7,15,13)$.

Family A allows $t=1,2$, because $3\cdot6\gt10$; the bound on $n$ is $(10-1)/3=3$ for $t=1$ and $(5-2)/3=1$ for $t=2$:

$(t,n)=(1,1)$, the pair $(2,1)$: $g_A=4$ and $\lfloor10/4\rfloor=2$, the triangles $(3,8,7)$ and $(6,16,14)$;
$(t,n)=(1,2)$, the pair $(3,2)$: $g_A=7$, the triangle $(5,21,19)$;
$(t,n)=(1,3)$, the pair $(4,3)$: $g_A=10$, the triangle $(7,40,37)$;
$(t,n)=(2,1)$, the pair $(3,1)$: $g_A=10$, the triangle $(8,15,13)$.

Together $2+5=7$ triangles. Pairs of type 3, such as $(4,1)$ and $(5,2)$, have $3\mid t$ and are always skipped; they would only repeat $(3,5,7)$ and $(8,7,13)$ with the shorter sides exchanged. A direct search over all pairs $a\lt b$ finds the same seven triangles, and also $T(1)=0$, $T(2)=2$, $T(4)=10$, $T(5)=14$, $T(10)=46$ and the given $T(100)=1234$.

How the Code Works

The C++ program ^[9] follows Sections 8 and 9 with two lookup tables. main reads the optional arguments --n=, the bound $N$ with default $1053779$, and --skip-checkpoints; unless the latter is given, run_validation_checkpoints first computes $T(100)$, $T(1000)$ and $T(10000)$ and stops with an error message if one differs from the statement. Then count_triangles computes $T(N)$, which is printed.

count_triangles returns $0$ for $N=0$ and refuses $N$ if $6N$ would not fit into 64 bits. It builds the table for $C=2N$ for count_type_b and the table for $C=6N$ for count_type_a, and adds the two counts. build_multiplier_table computes $C^2$ once in unsigned 128-bit arithmetic, finds $Q$ as the value for $q=1$ and stores $K_C(q)$ for $q=1,\dots,Q$ as 32-bit integers. Each value comes from floor_ratio_exact: it takes the long double estimate $C/(q\sqrt3)$, raises it while $k+1$ passes the test of Section 8 and lowers it while $k$ fails, where within_radius_bound evaluates $3(kq)^2\le C^2$ in 128-bit arithmetic. So the floating-point value only proposes and whole numbers decide, even where $C/(q\sqrt3)$ lies very close to a whole number.

count_type_b finds the largest $n$ with $n\le\lfloor Q_B/(n+1)\rfloor$, that is, $n(n+1)\le Q_B$. For each $n$ up to it, contribution_for_n runs over $t=1,\dots,\lfloor Q_B/n\rfloor-n$, skips the multiples of $3$ and the pairs for which std::gcd(n, t) is not $1$, and adds the table value at $q=n(n+t)$. count_type_a finds the largest $t$ with $t(t+3)\le Q_A$; for each $t$ not divisible by $3$, contribution_for_t runs over $n=1,\dots,\lfloor(\lfloor Q_A/t\rfloor-t)/3\rfloor$, skips the pairs that are not coprime and adds the value at $q=t(3n+t)$. At the target the table values stay below $2^{32}$, which would fail only for $N$ above about $1.24\cdot10^9$, and the totals are unsigned 64-bit integers.

The Python and Java programs run the same loops, family B first, with the same bounds and filters, but instead of tables they compute $K_C(q)$ on the spot for every coprime pair, in floor_c_div_q_sqrt3 or floorCDivQSqrt3. Python uses integers of unlimited size and starts from isqrt(C2 // (3 * q * q)) with C2 $=C^2$, which is already exact: for $k=\lfloor\sqrt u\rfloor$ we have $k^2\le\lfloor u\rfloor\le u\lt(k+1)^2$, so $\lfloor\sqrt{\lfloor u\rfloor}\rfloor=k$; its correction loops, the same as in C++, never change the value. Java corrects Math.sqrt of the integer quotient with the same loops in long arithmetic, which is safe because $C^2\lt4\cdot10^{13}$ and the products stay below $4C^2$; its gcd works on int values, enough for $n,t\lt1.3\cdot10^6$. Neither program checks the given values.

Complexity and Verification

Building the tables takes $O(N)$ time, a constant amount of work per entry. At the target $Q_B=1216799$ and $Q_A=3650397$, so the C++ tables hold $4\,867\,198$ values of four bytes, about $19.5$ MB, both kept until count_triangles returns; Python and Java need only constant extra memory. In family B the ranges of Section 9 contain $\sum_n(\lfloor Q_B/n\rfloor-n)$ pairs, with $n$ up to about $\sqrt{Q_B}$. Since $1+\tfrac12+\dots+\tfrac1x$ is about $\ln x$, ^[8] this is about $\tfrac12Q_B\ln Q_B$, and family A gives about $\tfrac16Q_A\ln Q_A$, so the work grows like $N\log N$. At the target the ranges contain $8\,618\,132$ pairs in family B and $9\,286\,092$ in family A; skipping $t$ divisible by $3$ leaves $12\,382\,096$ pairs that need a greatest common divisor, and $8\,948\,846$ of them are coprime and read a table value. Euclid's algorithm needs at most five steps per decimal digit of the smaller number, ^[6] so the time is $O(N\log^2N)$; on a current desktop computer the C++ and Java programs take a fraction of a second and the Python program a few seconds.

The C++ program checks the three given values unless --skip-checkpoints is given. For this article the following checks were run in addition; none of them is part of the published solutions.

A direct search over all pairs $a\lt b$ agrees with the C++ program for every $N\le300$. It stops at $a\lt2\sqrt3\,N$, because $g\gt a$, and at $2b-a\lt3a^2$, because $4c^2=(2b-a)^2+3a^2$ gives $(2c-2b+a)(2c+2b-a)=3a^2$ with two positive factors.
A second count runs over all $a\le L$ and all factorisations of $3a^2$ into these two factors, solves the two linear equations for $b$ and $c$ and keeps the whole-number solutions with $a\lt b$, without using $P(m,n)$. It reproduces the three given values and agrees with all three programs at $N=1053779$.
A search over all $x,y$ finds $828$ ordered primitive solutions of $z^2=x^2+xy+y^2$ with $z\le3000$, that is, $414$ triangles up to exchanging $x$ and $y$; each triangle comes from exactly one pair of type 1, in one of its two orders, as Section 6 claims.
At the target every table value equals $\lfloor L/g\rfloor$, and Python's start value never needs a correction; copies of the Python and Java programs with $N=100$, $1000$, $10000$ give the values of the statement.

Footnotes and References

The references below supply the background facts used above. The radius formula $r=g/(2\sqrt3)$, the partner identity of Section 6, the two families and the operation counts are derived explicitly in this article.

Project Euler 195 — 60-degree Triangle Inscribed Circles. The official statement calls an integer-sided triangle with exactly one angle of 60 degrees a 60-degree triangle, defines T(n) as the number of them whose inscribed circle has radius r ≤ n, gives T(100) = 1234, T(1000) = 22767 and T(10000) = 359912, and asks for T(1053779).
Wikipedia — Law of cosines. For a triangle with sides a, b, c opposite the angles α, β, γ, c² = a² + b² − 2ab cos γ; Sections 1 and 3 use it for angles of 60° and 120°.
Wikipedia — Incircle and excircles. The incentre lies on the internal angle bisectors, the two tangent segments from a vertex are equal (s − a from A, with s half the perimeter), and r = Δ/s, the area divided by s.
Wikipedia — Integer triangle. For integer triangles with a 120° angle it gives the sides m² + mn + n², 2mn + n² and m² − n² with coprime 0 < n < m and states that their greatest common divisor is 3 if m ≡ n (mod 3) and 1 otherwise, and that the extra condition m ≢ n (mod 3) generates every primitive triple exactly once; Sections 5 and 6 prove both statements.
Wikipedia — Pythagorean triple. The section on rational points of the unit circle writes them as x = (1 − t²)/(1 + t²), y = 2t/(1 + t²) with t = y/(1 + x), and states that, except for (−1, 0), a point is rational exactly when t is.
Wikipedia — Euclidean algorithm. By Lamé's theorem the algorithm never needs more steps than five times the number of decimal digits of the smaller number.
Wikipedia — Floor and ceiling functions. For an integer n and a real x, n ≤ x if and only if n ≤ ⌊x⌋, and for a positive integer n, ⌊x/n⌋ = ⌊⌊x⌋/n⌋.
Wikipedia — Harmonic number. The harmonic number 1 + 1/2 + … + 1/n is about as large as the natural logarithm of n; more precisely it behaves like ln n + γ + 1/(2n) − ….
C++ — Euler195.cpp. The linked immutable C++ revision contains floor_ratio_exact, count_type_b, count_type_a and run_validation_checkpoints. It is the implementation record for this article.

Problem 195 source code

C++

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

namespace {

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

constexpr u64 kTargetN = 1'053'779ULL;

struct Options {
    u64 n_limit = kTargetN;
    unsigned threads = 0;
    bool run_checkpoints = true;
};

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

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

    u64 parsed = 0;
    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_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;
        }

        u64 parsed = 0;
        if (parse_u64_after_prefix(arg, "--n=", parsed)) {
            options.n_limit = parsed;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--threads=", parsed)) {
            if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
                std::cerr << "Invalid thread count: " << parsed << "\n";
                return false;
            }
            options.threads = static_cast<unsigned>(parsed);
            continue;
        }

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

    return true;
}

unsigned resolve_thread_count(unsigned requested) {
    if (requested > 0) {
        return requested;
    }

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

inline bool within_radius_bound(const u64 multiplier,
                                const u64 q_value,
                                const u128 numerator_square) {
    const u128 scaled = static_cast<u128>(multiplier) * static_cast<u128>(q_value);
    return static_cast<u128>(3) * scaled * scaled <= numerator_square;
}

u64 floor_ratio_exact(const u64 numerator, const u64 q_value, const u128 numerator_square) {
    static const long double kSqrt3 = std::sqrt(3.0L);

    u64 estimate = static_cast<u64>(
        static_cast<long double>(numerator) /
        (static_cast<long double>(q_value) * kSqrt3));

    while (estimate < std::numeric_limits<u64>::max() &&
           within_radius_bound(estimate + 1, q_value, numerator_square)) {
        ++estimate;
    }
    while (!within_radius_bound(estimate, q_value, numerator_square)) {
        --estimate;
    }

    return estimate;
}

std::vector<u32> build_multiplier_table(const u64 numerator) {
    const u128 numerator_square = static_cast<u128>(numerator) * numerator;
    const u64 q_limit = floor_ratio_exact(numerator, 1, numerator_square);

    std::vector<u32> table(static_cast<std::size_t>(q_limit + 1), 0U);
    for (u64 q = 1; q <= q_limit; ++q) {
        table[static_cast<std::size_t>(q)] =
            static_cast<u32>(floor_ratio_exact(numerator, q, numerator_square));
    }
    return table;
}

template <typename Function>
u64 parallel_sum_over_range(const u64 first,
                            const u64 last,
                            unsigned threads,
                            const Function& contribution) {
    if (first > last) {
        return 0;
    }

    const u64 job_count = last - first + 1;
    if (threads == 0) {
        threads = 1;
    }
    threads = static_cast<unsigned>(std::min<u64>(threads, job_count));

    if (threads <= 1) {
        u64 total = 0;
        for (u64 i = first; i <= last; ++i) {
            total += contribution(i);
        }
        return total;
    }

    std::atomic<u64> next_index{first};
    std::vector<u64> partial(threads, 0);
    std::vector<std::thread> workers;
    workers.reserve(threads);

    for (unsigned tid = 0; tid < threads; ++tid) {
        workers.emplace_back([&, tid]() {
            u64 local_sum = 0;
            while (true) {
                const u64 index = next_index.fetch_add(1, std::memory_order_relaxed);
                if (index > last) {
                    break;
                }
                local_sum += contribution(index);
            }
            partial[tid] = local_sum;
        });
    }

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

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

u64 count_type_b(const std::vector<u32>& multiplier_table, const unsigned threads) {
    if (multiplier_table.size() <= 1) {
        return 0;
    }

    const u64 q_limit = static_cast<u64>(multiplier_table.size() - 1);

    u64 n_max = 0;
    for (u64 n = 1; n <= q_limit / (n + 1); ++n) {
        n_max = n;
    }

    auto contribution_for_n = [&](const u64 n) -> u64 {
        const u64 t_max = q_limit / n - n;
        u64 local = 0;

        for (u64 t = 1; t <= t_max; ++t) {
            if (t % 3 == 0) {
                continue;
            }
            if (std::gcd(n, t) != 1) {
                continue;
            }

            const u64 q = n * (n + t);  // q = m*n with m = n + t.
            local += multiplier_table[static_cast<std::size_t>(q)];
        }

        return local;
    };

    return parallel_sum_over_range(1, n_max, threads, contribution_for_n);
}

// Primitive 60-degree integer triangles come from primitive 120-degree triples:
// x = m^2 - n^2, y = 2mn + n^2, z = m^2 + mn + n^2 with gcd(m,n)=1 and (m-n)%3!=0.
// This yields two 60-degree families, parameterized with t = m-n:
// Type B: q = m*n = n*(n+t), r = k*sqrt(3)*q/2.
// Type A: q = (m-n)*(m+2n) = t*(3n+t), r = k*sqrt(3)*q/6.
// For fixed q, valid scales are k <= floor(C / (q*sqrt(3))) with C in {2N, 6N}.
u64 count_type_a(const std::vector<u32>& multiplier_table, const unsigned threads) {
    if (multiplier_table.size() <= 1) {
        return 0;
    }

    const u64 q_limit = static_cast<u64>(multiplier_table.size() - 1);

    u64 t_max = 0;
    for (u64 t = 1; t <= q_limit / (t + 3); ++t) {
        t_max = t;
    }

    auto contribution_for_t = [&](const u64 t) -> u64 {
        if (t % 3 == 0) {
            return 0;
        }

        const u64 n_max = (q_limit / t - t) / 3;
        u64 local = 0;

        for (u64 n = 1; n <= n_max; ++n) {
            if (std::gcd(n, t) != 1) {
                continue;
            }

            const u64 q = t * (3 * n + t);  // q = (m-n)*(m+2n).
            local += multiplier_table[static_cast<std::size_t>(q)];
        }

        return local;
    };

    return parallel_sum_over_range(1, t_max, threads, contribution_for_t);
}

u64 count_triangles(const u64 n_limit, const unsigned threads) {
    if (n_limit == 0) {
        return 0;
    }
    if (n_limit > std::numeric_limits<u64>::max() / 6ULL) {
        throw std::overflow_error("n is too large for safe arithmetic.");
    }

    const u64 numerator_b = 2ULL * n_limit;
    const std::vector<u32> table_b = build_multiplier_table(numerator_b);
    const u64 count_b = count_type_b(table_b, threads);

    const u64 numerator_a = 6ULL * n_limit;
    const std::vector<u32> table_a = build_multiplier_table(numerator_a);
    const u64 count_a = count_type_a(table_a, threads);

    return count_a + count_b;
}

bool run_validation_checkpoints() {
    struct Checkpoint {
        u64 n_limit;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {100, 1234},
        {1'000, 22'767},
        {10'000, 359'912},
    };

    for (const Checkpoint cp : checkpoints) {
        const u64 got = count_triangles(cp.n_limit, 1);
        if (got != cp.expected) {
            std::cerr << "Checkpoint failed for T(" << cp.n_limit << "): got " << got
                      << ", expected " << cp.expected << "\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_validation_checkpoints()) {
        return 1;
    }

    const unsigned threads = resolve_thread_count(options.threads);
    const u64 answer = count_triangles(options.n_limit, threads);
    std::cout << answer << '\n';
    return 0;
}

Python

# Problem 195: Inscribed circles of triangles with one angle of 60 degrees
# Count triangles with inscribed circle radius <= N = 1053779.

from math import gcd, isqrt

def solve():
    N = 1053779

    def floor_c_div_q_sqrt3(C, q):
        """floor(C / (q * sqrt(3))), i.e., max k such that 3*q*q*k*k <= C*C"""
        C2 = C * C
        est = isqrt(C2 // (3 * q * q))
        while (est + 1) * (est + 1) * 3 * q * q <= C2:
            est += 1
        while est > 0 and est * est * 3 * q * q > C2:
            est -= 1
        return est

    # Type B: q = n*(n+t) where t = m-n, gcd(n,t)=1, t%3 != 0
    # radius bound: k*sqrt(3)*q/2 <= N => k <= 2N/(q*sqrt(3))
    C_b = 2 * N
    q_limit_b = floor_c_div_q_sqrt3(C_b, 1)

    total = 0

    # Type B
    n = 1
    while n * (n + 1) <= q_limit_b:
        t_max = q_limit_b // n - n
        for t in range(1, t_max + 1):
            if t % 3 == 0: continue
            if gcd(n, t) != 1: continue
            q = n * (n + t)
            k = floor_c_div_q_sqrt3(C_b, q)
            total += k
        n += 1

    # Type A: q = t*(3n+t), gcd(n,t)=1, t%3 != 0
    # radius bound: k*sqrt(3)*q/6 <= N => k <= 6N/(q*sqrt(3))
    C_a = 6 * N
    q_limit_a = floor_c_div_q_sqrt3(C_a, 1)

    t = 1
    while t * (t + 3) <= q_limit_a:
        if t % 3 != 0:
            n_max = (q_limit_a // t - t) // 3
            for n in range(1, n_max + 1):
                if gcd(n, t) != 1: continue
                q = t * (3 * n + t)
                k = floor_c_div_q_sqrt3(C_a, q)
                total += k
        t += 1

    print(total)

solve()

Java

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

    static long floorCDivQSqrt3(long C, long q) {
        long C2 = C * C;
        long denom = 3 * q * q;
        long est = (long) Math.sqrt((double) (C2 / denom));
        while ((est + 1) * (est + 1) * 3 * q * q <= C2)
            est++;
        while (est > 0 && est * est * 3 * q * q > C2)
            est--;
        return est;
    }

    public static void main(String[] args) {
        long N = 1053779L;
        long total = 0;
        long CB = 2 * N, CA = 6 * N;
        long qLimB = floorCDivQSqrt3(CB, 1);
        // Type B: q = n*(n+t), t%3!=0, gcd(n,t)=1
        for (long n = 1; n * (n + 1) <= qLimB; n++) {
            long tMax = qLimB / n - n;
            for (long t = 1; t <= tMax; t++) {
                if (t % 3 == 0)
                    continue;
                if (gcd((int) n, (int) t) != 1)
                    continue;
                long q = n * (n + t);
                total += floorCDivQSqrt3(CB, q);
            }
        }
        long qLimA = floorCDivQSqrt3(CA, 1);
        // Type A: q = t*(3n+t), t%3!=0, gcd(n,t)=1
        for (long t = 1; t * (t + 3) <= qLimA; t++) {
            if (t % 3 == 0)
                continue;
            long nMax = (qLimA / t - t) / 3;
            for (long n2 = 1; n2 <= nMax; n2++) {
                if (gcd((int) n2, (int) t) != 1)
                    continue;
                long q = t * (3 * n2 + t);
                total += floorCDivQSqrt3(CA, q);
            }
        }
        System.out.println(total);
    }
}