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 is a triangle with integer sides $a\le b\le c$ in geometric progression, that is $b^2=ac$. ^[2] The statement's example has $a=144$, $b=156$, $c=169$, and indeed $156^2=24336=144\cdot169$. Write $G(L)$ for the number of geometric triangles of perimeter $a+b+c\le L$. The statement gives $G(10^6)=861805$ and asks for $G(L)$ at $L=2.5\cdot10^{13}$. ^[1]

The whole solution in one sentence: every geometric triangle is a whole multiple of a primitive one of the shape $(p^2,pq,q^2)$ with $\gcd(p,q)=1$ and $q/p$ below the golden ratio, so $G(L)$ is a sum of floor terms $\lfloor L/(p^2+pq+q^2)\rfloor$ over those primitive pairs, and the sum is evaluated by sweeping $q$ in blocks on which that floor never changes.

Mathematical Approach

1. The smallest geometric triangles

Any equilateral triangle $(n,n,n)$ qualifies, with common ratio $1$, and up to perimeter $18$ there are no others. The first other one is $(4,6,9)$, of ratio $3/2$ and perimeter $19$: here $6^2=36=4\cdot9$, and $9\lt4+6$, so it really closes into a triangle. Up to perimeter $50$ a direct search finds exactly $19$ triangles: the sixteen equilateral ones from $(1,1,1)$ to $(16,16,16)$, plus $(4,6,9)$, its double $(8,12,18)$, and $(9,12,16)$ of ratio $4/3$. So $G(10)=3$, $G(20)=7$ and $G(50)=19$.

2. Every geometric triangle has the shape $(mp^2,mpq,mq^2)$

Let $(a,b,c)$ be a triple of positive integers with $a\le b\le c$ and $b^2=ac$. Put $m=\gcd(a,c)$ and write $a=mu$, $c=mv$ with $\gcd(u,v)=1$. Then $b^2=m^2uv$, so $m^2\mid b^2$. Writing $v_\pi(n)$ for the exponent of a prime $\pi$ in $n$, so that $v_2(12)=2$, this gives $2v_\pi(m)\le2v_\pi(b)$, hence $v_\pi(m)\le v_\pi(b)$ and $m\mid b$. Therefore $uv=(b/m)^2$ is a perfect square. Since $u$ and $v$ are coprime, no prime occurs in both, so in $uv$ every prime exponent already comes from one factor alone and must be even there: $u$ and $v$ are themselves squares. Write $u=p^2$, $v=q^2$ with $p,q\ge1$; then $\gcd(p,q)=1$, and $b/m=pq$ because $b/m$ is positive. We have proved

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

with $\gcd(p,q)=1$ and $p\le q$, the last inequality because $a\le c$. Conversely every such triple satisfies $b^2=m^2p^2q^2=ac$ and is correctly ordered, and the representation is unique: $\gcd(mp^2,mq^2)=m\gcd(p,q)^2=m$ recovers $m$, then $p=\sqrt{a/m}$, $q=\sqrt{c/m}$. Such triples are therefore in one-to-one correspondence with the parameters $(m,p,q)$, $m\ge1$, $1\le p\le q$, $\gcd(p,q)=1$ — the statement's example being $m=1$, $p=12$, $q=13$, and $(8,12,18)$ being $m=2$, $p=2$, $q=3$. Not all of them are triangles: $m=p=1$, $q=2$ gives $(1,2,4)$ with $4\gt1+2$, so the triangle inequality is still to be imposed.

3. The triangle inequality is a golden-ratio bound

Three positive lengths form a triangle of positive area exactly when each is smaller than the sum of the other two. ^[3] With $a\le b\le c$ only $c\lt a+b$ has to be imposed, the other two being automatic. Substituting the parameterisation and dividing by $m$ gives

$$q^2\lt p^2+pq.$$

Divide by $p^2$ and set $x=q/p\ge1$: the condition becomes $x^2-x-1\lt0$. The positive root of $x^2-x-1$ is the golden ratio ^[4]

$$\varphi=\frac{1+\sqrt5}{2}=1.6180339\ldots,$$

so the admissible ratios are exactly $1\le q/p\lt\varphi$. The scale $m$ has dropped out, as it must: scaling a triangle does not change its shape. Equality $q^2=p^2+pq$ is impossible in integers, since it would force the irrational ratio $q/p=\varphi$, so a geometric triangle is never degenerate and it makes no difference whether the triangle inequality is read strictly.

4. The counting formula

The perimeter of $(mp^2,mpq,mq^2)$ is $m\,N(p,q)$ with

$$N(p,q)=p^2+pq+q^2.$$

For a fixed admissible pair the perimeter condition $m\,N(p,q)\le L$ is satisfied by exactly $\lfloor L/N(p,q)\rfloor$ values of $m$. Let $\mathcal P(L)$ be the set of pairs $(p,q)$ with

$$1\le p\le q,\quad\gcd(p,q)=1,\quad q^2\lt p^2+pq,$$

and $N(p,q)\le L$. Since the representation of Section 2 is unique, summing over $\mathcal P(L)$ counts every triangle once:

$$G(L)=\sum_{(p,q)\in\mathcal P(L)}\left\lfloor\frac{L}{N(p,q)}\right\rfloor.$$

The quantity $N(p,q)$ is the quadratic form of the hexagonal lattice, and the integers it represents are the Loeschian numbers $0,1,3,4,7,9,12,13,16,19,\ldots$ ^[5] Every primitive perimeter in this problem is one of them.

5. Worked example: $G(100)=42$

Take $L=100$. The pair $(1,1)$ has $N=3$ and gives $\lfloor100/3\rfloor=33$ equilateral triangles. For $p=2$ the ratio bound allows $q\le3$; $q=2$ is not coprime to $p$, and $q=3$ gives $N=19$ and $5$ triangles, $(4,6,9)$ scaled by $m=1,2,3,4,5$. For $p=3$ it allows $q\le4$, and only $q=4$ is coprime: $N=37$, two triangles $(9,12,16)$ and $(18,24,32)$. For $p=4$ it allows $q\le6$, but $q=4$ and $q=6$ share a factor with $p$, so $q=5$ is left with $N=61$ and $(16,20,25)$. For $p=5$ it allows $q\le8$, yet $N(5,7)=109$ and $N(5,8)=129$ exceed $100$, so $q=6$ is left with $N=91$ and $(25,30,36)$. Beyond $p=5$ nothing survives, since $3p^2\gt100$ for $p\ge6$. Altogether

$$G(100)=33+5+2+1+1=42.$$

6. The two bounds that make the sum finite

For $q\ge p$ the norm is at least $N(p,p)=3p^2$, a lower bound that holds even where $q=p$ itself is barred by coprimality, so a pair of $\mathcal P(L)$ can only occur while $3p^2\le L$:

$$1\le p\le p_{\max}=\bigl\lfloor\sqrt{L/3}\bigr\rfloor.$$

For a fixed $p$ the form $N(p,q)$ grows strictly with $q$: raising $q$ by $1$ adds $p+2q+1\gt0$ to it. Solving $q^2+pq+p^2\le L$ for $q$ gives the perimeter bound

$$q\le Q(p,L)=\left\lfloor\frac{\sqrt{4L-3p^2}-p}{2}\right\rfloor,$$

while Section 3 gives the shape bound $q\le\lfloor\varphi p\rfloor$. The candidates for this $p$ are the integers in $[p,q_{\mathrm{end}}]$ with $q_{\mathrm{end}}=\min(\lfloor\varphi p\rfloor,Q(p,L))$, and those coprime to $p$ are its pairs in $\mathcal P(L)$. The interval is never empty, because $3p^2\le L$ makes $q=p$ pass both bounds — though $q=p$ is coprime to $p$ only for $p=1$. Which bound bites is easiest to see before the rounding: with $R(p,L)=\tfrac12(\sqrt{4L-3p^2}-p)$, so that $Q=\lfloor R\rfloor$, putting $q=\varphi p$ into the norm gives $(3+\sqrt5)p^2$, hence $R(p,L)\ge\varphi p$ exactly when $L\ge(3+\sqrt5)p^2$, that is for

$$p\le\sqrt{L/(3+\sqrt5)}\approx0.437\sqrt L.$$

Below this crossing the shape bound rules, above it the perimeter bound, which is therefore in charge on about $24\%$ of the range of $p$. The two integer bounds may still tie right beside the crossing: at $L=10^6$, where the crossing is $437.02$, both equal $707$ at $p=437$, and $p=438$ is the first $p$ with a strictly smaller perimeter bound, $706$ against $708$. At $L=10^{12}$, crossing $437016.02$, the first such $p$ is $437017$.

7. Why the pairs cannot be listed one by one

How large is $\mathcal P(L)$? Dropping the coprimality condition and scaling $(p,q)=\sqrt L\,(x,y)$ turns the conditions into $x\le y\lt\varphi x$ and $x^2+xy+y^2\le1$, a fixed region. Put $t=y/x$, so $1\le t\lt\varphi$ and $0\lt x\le(1+t+t^2)^{-1/2}$; a strip of width $dt$ at fixed $x$ is a strip of width $x\,dt$ in $y$, so the area element is $x\,dx\,dt$ and the inner integral over $x$ leaves

$$\begin{aligned}A&=\frac12\int_1^{\varphi}\frac{dt}{1+t+t^2}\\&=\frac{1}{\sqrt3}\left(\arctan\frac{2+\sqrt5}{\sqrt3}-\frac\pi3\right),\end{aligned}$$

numerically $A=0.0782074\ldots$, so there are about $A\,L$ pairs with $N\le L$. The Möbius function has $\mu(1)=1$, $\mu(n)=(-1)^k$ when $n$ is a product of $k$ distinct primes, and $\mu(n)=0$ when a prime square divides $n$. ^[6] A pair with $\gcd(p,q)=d$ is $d$ times a primitive pair whose norm is $N/d^2$, so sorting all pairs by their common factor and inverting with $\mu$ multiplies the density by $\sum_d\mu(d)/d^2=1/\zeta(2)$:

$$\#\mathcal P(L)\approx\kappa L,\qquad \kappa=\frac{6}{\pi^2}A=0.0475444\ldots$$

A count confirms both constants: at $L=10^8$ there are $7823606$ pairs before the coprimality filter and $4754410$ after it, densities $0.07823606$ and $0.04754410$, which miss $A$ and $\kappa$ by $2.9\cdot10^{-5}$ and $3.2\cdot10^{-7}$ — the primitive one to six decimal places, the raw one, with the larger boundary error, to four. Adding up the candidate intervals at the limit of the problem gives exactly $1955187452248$ pairs, roughly $2\cdot10^{12}$: that is how many steps a loop over every candidate $q$ would take. Far too many, so the pairs must be handled in groups.

8. Counting the coprime $q$ in an interval

The first group operation removes the coprimality test. The divisor sum $\sum_{d\mid n}\mu(d)$ is $1$ for $n=1$ and $0$ otherwise, ^[6] so applying it to $n=\gcd(p,q)$ and swapping the two sums, the number of $q\in[\ell,r]$ with $\gcd(p,q)=1$ is

$$\begin{aligned}C_p(\ell,r)&=\sum_{d\mid p}\mu(d)\,\delta_d,\\ \delta_d&=\Bigl\lfloor\frac rd\Bigr\rfloor-\Bigl\lfloor\frac{\ell-1}d\Bigr\rfloor,\end{aligned}$$

because $\delta_d$ counts the multiples of $d$ in $[\ell,r]$. Only squarefree $d$ contribute, so the sum has $2^{\omega(p)}$ terms, $\omega(p)$ being the number of distinct primes of $p$: this is inclusion–exclusion over those primes. For $p=30$ and the interval $[30,48]$ the eight terms are $+19$, $-10$, $-7$, $-4$, $+4$, $+2$, $+2$, $-1$ for $d=1,2,3,5,6,10,15,30$, and they add up to $5$ — the values $31,37,41,43,47$.

9. Blocks on which the floor is constant

The second group operation is the decisive one. Since $N(p,q)$ increases with $q$, the floor $t=\lfloor L/N(p,q)\rfloor$ never increases, and it stays constant over long runs. Let a run start at $q=\ell$ and put

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

For any $q$ with $N(p,q)\le M$ we have $N(p,q)\le L/t$, hence $L/N(p,q)\ge t$ and $\lfloor L/N(p,q)\rfloor\ge t$ because $t$ is an integer. ^[7] If in addition $q\ge\ell$ then $N(p,q)\ge N(p,\ell)$ gives $\lfloor L/N(p,q)\rfloor\le t$. So the floor equals $t$ on exactly the integers $q\in[\ell,r]$ with $r=\min(q_{\mathrm{end}},Q(p,M))$, the right end being computable with the same formula as in Section 6. The run contributes $t\cdot C_p(\ell,r)$, and the next run starts at $r+1$; since $M\ge N(p,\ell)$ we always have $r\ge\ell$, so the sweep advances.

An example at $L=10^6$ and $p=500$: the shape bound allows $q\le809$, the perimeter bound only $q\le651$, so $q_{\mathrm{end}}=651$. At $\ell=500$ we get $N=750000$, $t=1$ and $M=10^6$, whose own bound is $Q(500,10^6)=651$: the whole range is a single block of $152$ values. With $500=2^2\cdot5^3$ the four terms of $C_{500}(500,651)$ are $+152$, $-76$, $-31$, $+16$, so $61$ of them are coprime to $500$ and this $p$ contributes $61$ triangles — four Möbius terms instead of $152$ separate tests.

A second example shows the quotient really changing. At $L=1000$ and $p=10$ we have $q_{\mathrm{end}}=16$, and the sweep uses three blocks: $[10,11]$ with $t=3$ and $M=333$, then $[12,15]$ with $t=2$ and $M=500$, then $[16,16]$ with $t=1$ and $M=1000$. Their coprime values are $q=11$, $q=13$ and none at all, so the three contributions are $3$, $2$ and $0$, and this $p$ adds $5$ triangles.

10. How many blocks there are

Write $B(p)$ for the number of blocks used for a fixed $p$. It cannot exceed the number of candidates in $[p,q_{\mathrm{end}}]$, coprime or not, which is at most $(\varphi-1)p+1$, nor the number of values the floor takes; since $N\ge3p^2$, those lie in $[1,L/(3p^2)]$. Hence

$$B(p)\le\min\Bigl((\varphi-1)p+1,\ \frac{L}{3p^2}+1\Bigr).$$

Put $\alpha=\varphi-1$. The two bounds cross at $p_0=(L/(3\alpha))^{1/3}$; let $K=\min(p_{\max},\lfloor p_0\rfloor)$. Summing the first bound up to $K$ and the second beyond it, and bounding that tail by $\int_K^{\infty}x^{-2}dx=1/K$,

$$\sum_{p\le p_{\max}}B(p)\le\frac{\alpha K(K+1)}2+\frac L{3K}+p_{\max},$$

which is $c\,L^{2/3}+O(\sqrt L)$ with $c=(3\alpha)^{1/3}/2=0.6143\ldots$ The true count is smaller: at $L=2.5\cdot10^{13}$ the sweep uses $297351015$ blocks, that is $0.348\,L^{2/3}$, against the $1955187452248$ pairs of Section 7 — a saving by a factor of about $6575$.

How the Code Works

The C++ program ^[10] is the reference implementation. parse_arguments fills an Options record from --limit= and --skip-checkpoints, insisting on a limit of at least $3$, the smallest possible perimeter. isqrt_u128 is the integer square root used everywhere: a floating-point estimate corrected up and down until it is exact. to_string_u128 prints the result digit by digit, since the standard streams cannot format the 128-bit type.

Its constructor computes the largest $p$ the run will ever need and calls build_smallest_prime_factors, a sieve of Eratosthenes storing for every index its smallest prime factor instead of a Boolean flag. ^[8] count recomputes $p_{\max}=$isqrt_u64(limit/3) and adds up count_for_single_p over $1\le p\le p_{\max}$, the single loop of count_range; asked for a limit needing a larger sieve, it returns $0$ rather than read past its table, which main prevents by sizing the counter for the largest limit of the run.

q_max_triangle and q_max_for_norm implement the two bounds of Section 6. Both start from a long double estimate and walk to the exact integer with the predicates triangle_ok and norm_ok, which compare $q^2$ and $N(p,q)$ in 128-bit arithmetic although $N$ stays below $2^{45}$ here. extract_distinct_primes peels the distinct primes of $p$ off the sieve, and build_squarefree_divisors_mu doubles the list of squarefree divisors once per prime, flipping the sign of $\mu$, so divisors and Möbius values are produced together.

count_for_single_p is the sweep of Section 9: it computes $q_{\mathrm{end}}$, then repeatedly takes norm_left, the quotient, the derived limit norm_limit and the block end, adds coprime_count * quotient to a 128-bit accumulator ^[9] and moves left past the block. count_coprime_in_interval picks between the two ways of counting: a span of at most eight values is settled by $\gcd$ tests, cheaper than a divisor sum of up to $128$ terms, anything longer by $C_p$ from Section 8. Its data sit in three fixed-capacity arrays, $8$ slots in primes and $128$ each in divisors and mus, which fit because every $p\le2886751$ has at most $7$ distinct primes, the smallest number with $8$ being $2\cdot3\cdots19=9699690$: exactly right here, but not for an arbitrarily large --limit.

brute_force_count is the comparison implementation: same two bounds, but every single $q$ is tested with std::gcd and $\lfloor L/N\rfloor$ added directly, with no blocks and no Möbius sum. run_checkpoints compares the fast count against the known values and against that brute force, and check_equal_u128 reports a mismatch on standard error and stops the run.

The Java port follows the same plan with long arithmetic: qMaxTri and qMaxNorm are the two bounds, distinctPrimes reads the spf table, and the squarefree divisors with their signs come from a bit mask over the distinct primes instead of a doubled list. Java has no 128-bit integer, so total is a long, and a crude bound shows that this suffices without knowing the result: each $p$ offers at most $\alpha p+1$ candidates, each contributing at most $L/(3p^2)$, so with $P=p_{\max}$ and $H_P=\sum_{p\le P}1/p\le1+\ln P$,

$$G(L)\le\frac L3\Bigl(\alpha H_P+\frac{\pi^2}6\Bigr).$$

At the limit of the problem this stays below $10^{14}$, far below $2^{63}-1$, and no partial sum is larger, since no contribution is negative; the C++ accumulator is wider than any such bound requires. The Java port has no self-checks and no options; its limit is the constant LIMIT. The Python file is a bridge: it looks for solutionsCpp/Euler370.cpp beside its own directory, compiles it with clang++ or else g++ and -std=c++17 -O2 into a cached binary when that is missing or older than the source, and returns what parse_output reads from the run. should_skip_cpp_checkpoints looks for --skip-checkpoints in the C++ source and passes the flag if it is there, so the run prints nothing but the final number — which matters, because parse_output prefers a line matching ANSWER_RE, then the last line matching EQUAL_RE, and only then the last line of the output.

Complexity and Verification

The sieve costs $O(\sqrt L\log\log\sqrt L)$ time and $O(\sqrt L)$ memory: here its $p_{\max}+1=2886752$ 32-bit entries take about $11.5$ MB, which dominates the memory used, as the sweep needs only the three small arrays above. By Section 10 it handles $O(L^{2/3})$ blocks, $297351015$ of them at the limit of the problem, each costing two quotients and a corrected integer square root, then either at most eight greatest common divisors or $2^{\omega(p)}\le128$ Möbius terms of two divisions each. In the measured run $199479070$ blocks took the short greatest-common-divisor path and $97871945$ the divisor sum, needing $821310840$ terms in all, $8.4$ per divisor-sum block. One pass over all $p$ was measured at about $11$ seconds for the C++ program and $16$ for the Java port on one desktop; the operation counts are the machine-independent part of that.

The program checks itself in two stages. It reproduces $G(10^6)=861805$ from the statement together with $G(100)=42$, $G(1000)=532$, $G(10^4)=6427$, $G(10^5)=75243$, and it compares the fast count with brute_force_count at $L=3,17,111,999,12345,54321$, where the counts are $1,5,48,532,8063,39297$. All pass. Since the brute force shares the bound functions with the fast path, it tests the blocks and the Möbius sums but not the parameterisation.

So the parameterisation was checked for this article by two programs that know nothing about $(m,p,q)$. The first loops over all pairs $a\le c$ up to $L$, sets $b=\lfloor\sqrt{ac}\rfloor$ and keeps the triple when $b^2=ac$, $a+b+c\le L$ and $c\lt a+b$; it confirms $G(10)=3$, $G(20)=7$, $G(50)=19$, $G(100)=42$, $G(500)=249$, $G(1000)=532$. The second runs over the middle side $b$, enumerating the divisors $a\le b$ of $b^2$ with $c=b^2/a$; it reaches further: $G(10^4)=6427$, $G(10^5)=75243$, $G(10^6)=861805$, $G(10^7)=9712598$. These comparisons all use the C++ program, the only one of the three with an adjustable limit; the Java port and the Python bridge each compute the single count their fixed limit asks for, and both return the same value as it does.

One asymptotic prediction closes the verification. A lattice count of this shape is its area up to its perimeter, so $\#\mathcal P(X)=\kappa X+O(\sqrt X\log X)$ — the densities of Section 7 confirm the constant — and partial summation turns the counting formula into $G(L)=\kappa L\ln L+O(L)$. Since there are about $\kappa$ primitive pairs per unit of norm, each contributing $L/n$, this predicts that $G(L)/L$ grows by $\kappa\ln10=0.1094751$ per factor of $10$. The measured $0.42$, $0.532$, $0.6427$, $0.75243$, $0.861805$, $0.9712598$ at $L=10^2,\ldots,10^7$ grow by $0.112$, $0.1107$, $0.10973$, $0.109375$, $0.1094548$, matching that to four digits. The $O(L)$ remainder alone does not force those differences to converge, so this is a consistency check, not a proof.

Footnotes and References

The references below supply the background facts used above; the parameterisation, the golden-ratio bound, the counting formula, the block argument and the operation counts are derived here.

Project Euler 370 — Geometric Triangles. The official statement defines a geometric triangle as an integer sided triangle with sides a ≤ b ≤ c forming a geometric progression, b² = a·c, gives the example with sides 144, 156 and 169, states that there are 861805 geometric triangles with perimeter at most 10⁶, and asks how many exist with perimeter at most 2.5·10¹³.
Wikipedia — Geometric progression. A geometric progression is a sequence of non-zero numbers in which each term is the previous one times a fixed common ratio, and the article gives the recurrence aₙ = aₙ₋₁²/aₙ₋₂, which rearranged is the condition b² = a·c used throughout this article.
Wikipedia — Triangle inequality. For a non-degenerate triangle each side is strictly shorter than the sum of the other two, and conversely three lengths with that property are the sides of a triangle of positive area; with equality the triangle is degenerate, of zero area. Section 3 imposes this condition.
Wikipedia — Golden ratio. The golden ratio φ = (1+√5)/2 = 1.618033988749… is the positive root of φ² − φ − 1 = 0, the quadratic that Section 3 obtains from the triangle inequality.
OEIS A003136 — Loeschian numbers. The numbers of the form x² + xy + y², the norms of vectors in the hexagonal A₂ lattice and of the Eisenstein integers, beginning 0, 1, 3, 4, 7, 9, 12, 13, 16, 19, 21, 25, 27, 28, 31, 36, 37, … The primitive perimeters 3, 19, 37, 61, 91 of Section 5 are terms of this sequence.
Wikipedia — Möbius function. The Möbius function is multiplicative, equals (−1)ᵏ on a product of k distinct primes — so μ(2) = −1, μ(6) = 1 — and 0 on any number divisible by the square of a prime, for instance μ(4) = 0, and its divisor sum gives 1 for n = 1 and 0 for every n > 1. Sections 7 and 8 use that identity and the density factor 1/ζ(2) = 6/π².
Wikipedia — Floor and ceiling functions. The article lists the equivalence n ≤ x if and only if n ≤ ⌊x⌋ for an integer n, the step that turns N ≤ ⌊L/t⌋ into ⌊L/N⌋ ≥ t in Section 9.
Wikipedia — Sieve of Eratosthenes. The sieve marks the multiples of each prime in turn and runs in O(n log log n) operations, the cost quoted in the complexity section. The variant used by the code writes the smallest prime factor into each slot instead of a Boolean flag, so one lookup chain factorises any p.
GCC manual — 128-bit integers. The page documents the __int128 and unsigned __int128 extensions, available on targets with a wide enough integer mode; the C++ program uses the unsigned form for its accumulator and for the comparisons of the two bounds.
C++ — Euler370.cpp. The linked immutable revision contains q_max_triangle, q_max_for_norm, build_squarefree_divisors_mu, count_coprime_in_interval, count_for_single_p and brute_force_count with the packaged checkpoints. It is the implementation record for this article; the independent counts above are not part of it.

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