Problem 299: Three Similar Triangles

Project Euler Problem 299 Solution

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

Problem Summary The problem starts with $$A(a,0),\qquad B(b,0),\qquad C(0,a),\qquad D(0,d),$$ where $0<a<b$ and $0<a<d$. We seek integer triplets $(a,b,d)$ for which there exists an integer point $P$ on $AC$ such that the three triangles $ABP$, $CDP$, and $BDP$ are all similar. The quantity to count is the number of such triplets with $$b+d<N.$$ Mathematical Approach 1. Why $a=c$ and why two cases appear Project Euler already states that similarity is possible only when $a=c$, so the segment $AC$ lies on the line $$x+y=a.$$ Because $OA$ and $OC$ are symmetric, the angles at $A$ and $C$ are both $135^\circ$. Matching the three similar triangles leaves two possible angle correspondences: Incenter case. The lines $PB$ and $PD$ are angle bisectors of the right triangle $OBD$, so $P$ is the incenter of $\triangle OBD$. Parallel case. The line $AC$ is parallel to $BD$. These two families are disjoint, and the code counts them separately. 2. Incenter case: reduction to a Pythagorean triple If $P$ is the incenter of the right triangle with legs $b$ and $d$, then its coordinates must be $$P=(i,i),$$ where $i$ is the inradius. Since $P\in AC$, we have $$a=2i.$$ For a right triangle, the inradius is $$i=\frac{b+d-\sqrt{b^2+d^2}}{2}.$$ Therefore $a$ is integral exactly when $\sqrt{b^2+d^2}$ is integral....

Detailed mathematical approach

Problem Summary

Four points with integer coordinates lie on the axes: $A(a,0)$, $B(b,0)$, $C(0,c)$ and $D(0,d)$, with $0\lt a\lt b$ and $0\lt c\lt d$. A point $P$, also with integer coordinates, is chosen on the line $AC$ so that the three triangles $ABP$, $CDP$ and $BDP$ are all similar, that is, have the same three angles. The statement notes that this is possible only if $a=c$, so a configuration is described by the triplet $(a,b,d)$, and triplets that differ only in the order of $b$ and $d$, such as $(2,3,4)$ and $(2,4,3)$, count as different. For a positive integer $N$, let $T(N)$ be the number of triplets with $b+d\lt N$ for which at least one such point $P$ exists. For $(2,3,4)$ the point $P(1,1)$ works, the statement gives $T(100)=92$ and $T(100\,000)=320471$, and we must find $T(100\,000\,000)$. ^[1]

The whole solution in one sentence: the three triangles are similar exactly when $P$ is the centre of the circle inscribed in the right triangle formed by the origin, $B$ and $D$, or when $b=d$ and $P$ lies on one particular circle through $B$ and $D$; in both cases the triplets correspond to right triangles with integer sides, which we count by running through the primitive ones and counting the multiples of each with one division.

Mathematical Approach

Start with the example of the statement, $(a,b,d)=(2,3,4)$ with $P(1,1)$. The squared side lengths of $ABP$, $CDP$ and $BDP$, in increasing order, are $1,2,5$, then $2,4,10$, then $5,10,25$. The second and third lists are $2$ and $5$ times the first, so the side lengths are proportional and the three triangles are similar. ^[2] To find all such configurations, we first locate $P$ and the angles of the triangles; only two shapes survive, and both lead to Pythagorean triangles, which one identity describes completely.

1. Where the point $P$ can lie

Let $O$ be the origin. The statement tells us that similarity forces $a=c$. For a point $P$ between $A$ and $C$ the reason is short: the angle of $ABP$ at $A$ is $180^\circ-\angle OAC$ and the angle of $CDP$ at $C$ is $180^\circ-\angle OCA$. Both are obtuse, a triangle has at most one obtuse angle, and similar triangles have equal angles, so these two angles are equal. Then the right triangle $OAC$ has equal angles at $A$ and $C$, and $a=c$.

From now on $a=c$. The points $A$ and $C$ lie on the line $x+y=a$, and $\angle OAC=\angle OCA=45^\circ$, so $\angle CAB=\angle ACD=135^\circ$. The point $P$ cannot lie outside the segment $AC$. If it lies beyond $A$, the angle of $ABP$ at $A$ is $180^\circ-135^\circ=45^\circ$, so its other two angles add up to $135^\circ$ and none of its angles is $135^\circ$; but $P$ is then on the ray from $C$ through $A$, and $CDP$ has the angle $\angle ACD=135^\circ$ at $C$. A point beyond $C$ fails in the same way, and $P=A$ or $P=C$ makes a triangle degenerate. Hence $P=(x,y)$ with integers $x,y\ge1$ and $x+y=a$.

2. The angles of the three triangles

For such a point the triangle $ABP$ has the angle $\angle CAB=135^\circ$ at $A$, and $CDP$ has the angle $135^\circ$ at $C$. So $BDP$ must have an angle of $135^\circ$ too. It is not at $B$ or $D$: the point $P$ lies inside the right triangle $OBD$, because $x/b+y/d\lt x/a+y/a=1$, so $\angle PBD\lt\angle OBD\lt90^\circ$ and likewise $\angle PDB\lt90^\circ$. Hence $\angle BPD=135^\circ$.

The quadrilateral $ABDC$, the triangle $OBD$ with the corner $OAC$ cut off, is convex: its angles are $135^\circ$, $\angle OBD$, $\angle ODB$ and $135^\circ$. The segments $PB$ and $PD$ cut it into the triangles $ABP$, $BDP$ and $DCP$, whose angles at $P$ fill the straight angle $APC$. With $\varphi=\angle APB$ and $\psi=\angle DPC$ this gives $\varphi+135^\circ+\psi=180^\circ$, that is

$$\varphi+\psi=45^\circ.$$

The angle sum of $ABP$ now gives $\angle ABP=45^\circ-\varphi=\psi$, and that of $CDP$ gives $\angle CDP=\varphi$. So $ABP$ and $CDP$ both have the angles $135^\circ$, $\varphi$, $\psi$, and only $BDP$ remains to be matched.

3. Only two shapes are possible

The triangle $BDP$ has the angle $135^\circ$ at $P$, and its angles at $B$ and $D$ add up to $45^\circ$. Similarity requires its angles to be $135^\circ$, $\varphi$, $\psi$, ^[2] and this happens exactly when one of two cases holds.

Case 1: $\angle PBD=\psi$ and $\angle PDB=\varphi$. Since $A$ lies on $OB$ and $C$ on $OD$, the angle $\angle ABP=\psi$ lies between $BO$ and $BP$, and $\angle CDP=\varphi$ between $DO$ and $DP$. So $BP$ and $DP$ bisect the angles of the triangle $OBD$ at $B$ and $D$, and $P$ is its incentre, the point where the angle bisectors meet. ^[3]
Case 2: $\angle PBD=\varphi$ and $\angle PDB=\psi$. The equal angles $\angle APB$ and $\angle PBD$ are alternate angles at the line $PB$, which crosses the lines $AC$ and $BD$ and separates $A$ from $D$. Hence $BD$ is parallel to $AC$. ^[4] The line $AC$ has slope $-1$ and $BD$ has slope $-d/b$, so $b=d$.

Conversely, both cases work. If $P$ is the incentre of $OBD$, then $\angle PBD+\angle PDB$ is half of $\angle OBD+\angle ODB=90^\circ$, so $\angle BPD=135^\circ$, Section 2 applies, and the bisectors give $\angle PBD=\angle ABP=\psi$ and $\angle PDB=\angle CDP=\varphi$. If $b=d$, the parallel lines give $\angle PBD=\angle APB=\varphi$ and $\angle PDB=\angle CPD=\psi$ for every $P$ between $A$ and $C$; both are angles of triangles with an angle of $135^\circ$, hence below $45^\circ$, so $BDP$ can have its angle of $135^\circ$ only at $P$. Thus for $b=d$ the triangles are similar exactly when $\angle BPD=135^\circ$, and then Section 2 gives all three the angles $135^\circ$, $\varphi$, $\psi$.

4. The incentre family

The incentre of $OBD$ has the same distance $r$ from both axes, so it is the point $(r,r)$, and the incircle touches the axes at $(r,0)$ and $(0,r)$. The two tangent segments from a vertex to the incircle are equally long, ^[3] so the tangent lengths from $B$ and $D$ are $b-r$ and $d-r$, and the hypotenuse $z=BD$ is their sum. Hence $r=(b+d-z)/2$, and $P=(r,r)$ lies on the line $x+y=a$ exactly when

$$a=b+d-z,\qquad z=\sqrt{b^2+d^2}.$$

Since $a$, $b$ and $d$ are integers, so is $z$: the numbers $b,d,z$ form a Pythagorean triple, the sides of a right triangle with integer sides. ^[5] Conversely, every Pythagorean triple gives a valid triplet. The number $a=b+d-z$ is positive by the triangle inequality, it is smaller than $b$ and $d$ because each leg is shorter than the hypotenuse, and it is even: every integer has the same parity as its square, so $b+d-z$ has the parity of $b^2+d^2-z^2=0$. So $P=(a/2,a/2)$ has integer coordinates. The incentre family therefore consists of one triplet $(b+d-z,\,b,\,d)$ for each ordered pair of legs $(b,d)$ of a Pythagorean triangle, and each of these triplets has exactly one point $P$, since Section 3 allows no other point when $b\ne d$. The opening example comes from the triangle $3,4,5$: it has $a=3+4-5=2$, inradius $1$ and incentre $P(1,1)$, and the ordered pairs $(3,4)$ and $(4,3)$ give the triplets $(2,3,4)$ and $(2,4,3)$.

5. The parallel family

Now let $b=d$ and $\angle BPD=135^\circ$. Then the triangles $APB$ and $PBD$ have the same angles, $135^\circ$, $\varphi$, $\psi$, at the corresponding vertices $A,P,B$ and $P,B,D$. Their sides are therefore proportional, $AP/PB=PB/BD$, which means $PB^2=AP\cdot BD$. Conversely, if $PB^2=AP\cdot BD$, the two triangles have proportional sides around the equal angles $\angle APB=\angle PBD$, so they are similar by the side–angle–side criterion, ^[2] and $\angle BPD=\angle PAB=135^\circ$. For $P=(x,y)$ with $x+y=a$ we have $AP=\sqrt2\,y$, $BD=\sqrt2\,b$ and $PB^2=(b-x)^2+y^2$, so the condition is $(b-x)^2+y^2=2by$. Adding $b^2-2by$ to both sides turns it into

$$(b-x)^2+(b-y)^2=b^2.$$

So $P$ lies on the circle with centre $(b,b)$ through $B$ and $D$. This matches the picture: the centre sees $B$ and $D$ at a right angle, so by the inscribed angle theorem every point of the arc facing $O$ sees $BD$ at half of $270^\circ$, which is $135^\circ$. ^[6]

With $k=b-x$ and $l=b-y$, both positive because $x,y\lt a\lt b$, the condition reads $k^2+l^2=b^2$: the numbers $k,l,b$ are the sides of a Pythagorean triangle with hypotenuse $b$, and $a=x+y=2b-k-l$. Conversely, every such triangle gives a valid point: $x=b-k$ and $y=b-l$ are positive because each leg is shorter than the hypotenuse, and $a\lt b$ because $k+l\gt b$. Exchanging $k$ and $l$ gives the mirror point with the same $a$, so each triangle gives one triplet $(a,b,b)$ with two points $P$. Different triangles with the same hypotenuse give different values of $a$, because $a$ fixes $k+l=2b-a$ and $kl=\frac12\bigl((k+l)^2-b^2\bigr)$, hence the pair $\{k,l\}$. For example, the triangle $3,4,5$ gives $P(2,1)$ and $P(1,2)$ and the triplet $(3,5,5)$; with $P(2,1)$ the squared sides of $ABP$, $CDP$ and $BDP$ are $2,4,10$, then $4,8,20$, then $10,20,50$. The hypotenuse $25$ belongs to the triangles $7,24,25$ and $15,20,25$, which give $(19,25,25)$ and $(15,25,25)$.

6. The two families do not overlap

A triplet with $b=d$ cannot belong to the incentre family, since $b^2+b^2=z^2$ would make $\sqrt2=z/b$ rational. Hence

$$T(N)=F_1(N)+F_2(N),$$

where $F_1(N)$ counts the ordered leg pairs $(b,d)$ of Pythagorean triangles with $b+d\lt N$, and $F_2(N)$ counts the Pythagorean triangles whose hypotenuse $b$ satisfies $2b\lt N$.

7. One identity for all Pythagorean triangles

To count these triangles we describe them as the programs do. For legs $b,d$ and hypotenuse $z$ put $a=b+d-z$, $u=b-a$ and $v=d-a$. Substituting $z=b+d-a$ into $b^2+d^2=z^2$ gives $a^2-2ab-2ad+2bd=0$, that is

$$a^2=2uv,$$

for example $2^2=2\cdot1\cdot2$ for the triangle $3,4,5$. Let $g=\gcd(u,v)$ be the greatest common divisor of $u$ and $v$, and write $u=gu'$ and $v=gv'$. Since $g^2$ divides $a^2$, $g$ divides $a$, and with $a=ga'$ we get $a'^2=2u'v'$ with coprime $u'$ and $v'$. So $a'$ is even, $a'=2t$, $u'v'=2t^2$, and exactly one of $u'$, $v'$ is even. If $v'$ is even, the coprime numbers $u'$ and $v'/2$ have the square product $t^2$, so both are squares: in the prime factorisation of a square every exponent is even, and coprime numbers share no prime. ^[7] This gives the first of the following two orientations, and an even $u'$ gives the second:

orientation 1: $u=gm^2$, $v=2gn^2$, $a=2gmn$ with $m$ odd, so $b=g(m^2+2mn)$, $d=g(2mn+2n^2)$ and $b+d=g\,s_1(m,n)$;
orientation 2: $u=2gm^2$, $v=gn^2$, $a=2gmn$ with $n$ odd, so $b=g(2m^2+2mn)$, $d=g(2mn+n^2)$ and $b+d=g\,s_2(m,n)$;

always with $m,n\ge1$, $\gcd(m,n)=1$ and $g\ge1$, where

$$\begin{aligned}s_1(m,n)&=m^2+4mn+2n^2,\\ s_2(m,n)&=2m^2+4mn+n^2.\end{aligned}$$

Every choice gives $a^2=2uv$ and hence a genuine triangle, and the choice is unique: in both orientations $\gcd(u,v)=g$, and $u/g$ is odd in the first and even in the second, so $u$ and $v$ determine $g$, the orientation, $m$ and $n$. In orientation 1 the hypotenuse is $z=b+d-a=g\,h(m,n)$ with $h(m,n)=m^2+2mn+2n^2$. With $\mu=m+n$ and $\nu=n$, the triangle of orientation 1 with $g=1$ becomes $(\mu^2-\nu^2,\,2\mu\nu,\,\mu^2+\nu^2)$ with coprime $\mu\gt\nu$ of opposite parity. This is Euclid's formula for the primitive triples, those whose sides have no common factor, so $g$ is the factor by which a primitive triangle is scaled. ^[5]

8. Counting with one division per primitive triangle

For fixed $m$, $n$ and orientation, the multiples with $b+d\lt N$ are those with $g\,s\le N-1$, and there are $\lfloor(N-1)/s\rfloor$ of them. Orientation 2 is orientation 1 with $b$ and $d$ exchanged and the letters $m$ and $n$ swapped, since $s_2(m,n)=s_1(n,m)$, so both orientations contribute the same sum. In the parallel family each triangle is counted once, whatever the order of its legs. Exactly one leg of a primitive triangle is odd, ^[5] so every triangle has exactly one representation of orientation 1, the one that lists first the leg that is odd in the primitive triangle. Its hypotenuse $g\,h(m,n)$ is the coordinate $b$ of the triplet, and $a=2b-k-l=g(m^2+2n^2)$. The condition $2b\lt N$ leaves $\lfloor(N-1)/(2h)\rfloor$ multiples. Hence, with all sums over the pairs of positive integers $(m,n)$ with $m$ odd and $\gcd(m,n)=1$,

$$F_1(N)=2\sum_{(m,n)}\left\lfloor\frac{N-1}{s_1(m,n)}\right\rfloor,$$

$$F_2(N)=\sum_{(m,n)}\left\lfloor\frac{N-1}{2h(m,n)}\right\rfloor.$$

Only finitely many terms are nonzero, because $s_1$ and $2h$ exceed both $m^2$ and $2n^2$.

9. Worked example: $N=100$

For $N=100$ the pairs $(m,n)$ of Section 8 with $s_1\le99$ are $(1,1)$, $(1,2)$, $(3,1)$, $(1,3)$, $(3,2)$, $(5,1)$, $(1,4)$, $(1,5)$, $(5,2)$, $(7,1)$, $(3,4)$ and $(1,6)$. Their leg sums $s_1$ are $7$, $17$, $23$, $31$, $41$, $47$, $49$, $71$, $73$, $79$, $89$ and $97$, those of the twelve primitive triangles from $3,4,5$ to $13,84,85$, ^[8] and the numbers of multiples $\lfloor99/s_1\rfloor$ are $14$, $5$, $4$, $3$, $2$, $2$, $2$, $1$, $1$, $1$, $1$ and $1$. They add up to $37$, so $F_1(100)=74$. The pair $(3,3)$, with $s_1=63$, is skipped because $\gcd(3,3)=3$: its triangle $27,36,45$ is the ninth multiple of $3,4,5$, which $(1,1)$ already counts.

For the parallel family the condition $2h\le99$ leaves the pairs $(1,1)$, $(1,2)$, $(3,1)$, $(1,3)$, $(3,2)$, $(5,1)$ and $(1,4)$, with the primitive hypotenuses $h=5$, $13$, $17$, $25$, $29$, $37$, $41$ ^[8] and $\lfloor99/(2h)\rfloor=9$, $3$, $2$, $1$, $1$, $1$, $1$, so $F_2(100)=18$. Together $T(100)=74+18=92$, the value in the statement. The same sums give $F_1(100\,000)=245280$ and $F_2(100\,000)=75191$, which add up to the second value of the statement, $320471$.

How the Code Works

The C++ program ^[10] reads the options --limit= (default $10^8$) and --skip-checkpoints, runs run_checkpoints unless told not to, and prints the value of count_triplets, the sum of count_family_x_eq_y, the incentre family, whose point $P$ lies on the line $x=y$, and count_family_u_eq_v, the parallel family, where $u=v$ means $b=d$. Both functions set numerator to $N-1$, so numerator / s is the floor $\lfloor(N-1)/s\rfloor$ of Section 8.

count_family_x_eq_y lets $m$ run from $1$ to $\lfloor\sqrt N\rfloor+3$. For odd $m$ it runs $n=1,2,\dots$, computes $s_1(m,n)$ and adds numerator / s when std::gcd(m, n) == 1; for every $m$ a second loop does the same with $s_2(m,n)$ for odd $n$, which is orientation 2. Each inner loop stops at the first $n$ with $s\ge N$, which is correct because $s$ grows with $n$. The range of $m$ suffices because $s_1\gt m^2$ and $s_2\gt2m^2$; the extra $3$ only guards against rounding in the floating-point square root. The program evaluates both orientations instead of doubling one sum, which costs time but gives the same result.

count_family_u_eq_v runs over the odd values $p=1,3,5,\dots$ up to $\lfloor\sqrt{N/2}\rfloor+3$ and over $q=1,2,\dots$; here $(p,q)$ is the pair $(m,n)$ of Section 8, and s $=p^2+2pq+2q^2$ is the primitive hypotenuse $h$. The loop over $q$ stops when denom $=2s$ reaches $N$, and the function adds numerator / denom when $\gcd(p,q)=1$; the range of $p$ suffices because $2h\gt2p^2$. The product $2s$ is formed in the 128-bit type u128, so that the comparison with $N$ cannot overflow. The counts are kept in the 64-bit type u64, which is ample for the target, since $T(N)\lt N^2=10^{16}$: a triplet is determined by $(b,d)$, or in the parallel family by $b$ and the shorter leg.

Both functions return $0$ for very small limits. The first guard, $N\le7$, is exact, because $3+4=7$ is the smallest leg sum. The second guard, $N\le13$, is stricter than necessary: the smallest member of the parallel family, $(3,5,5)$, has $b+d=10$, so for $N=11$, $12$ and $13$ the programs report $2$ instead of $3$. This does not affect the checkpoints or the requested limit. run_checkpoints compares count_triplets with $92$ for $N=100$ and with $320471$ for $N=100\,000$, and throws an exception, which main reports, if either value differs.

The Python and Java versions contain the same loops, in the functions count_x_eq_y_worker and count_u_eq_v_worker of Python and countXEqYWorker and countUEqVWorker of Java. Python uses exact integers, math.isqrt for the bounds and math.gcd, and its solve(limit) returns the count as a string. Java uses long arithmetic, Math.sqrt for the bounds and a recursive Euclidean gcd, and its solve() always uses $N=10^8$. Neither of them runs the checkpoints.

Complexity and Verification

In each of the three loops one parameter stays below $\sqrt N$ and the other below $\sqrt{N/2}$, because $s_1$, $s_2$ and $2h$ exceed the corresponding squares, so each loop examines fewer than $N/\sqrt2$ pairs. For $N=10^8$ the loops examine $15\,578\,139$, $31\,152\,770$ and $9\,815\,714$ pairs inside their bounds, $56\,546\,623$ in all, plus one final test that ends each inner loop. A pair costs a few multiplications and a gcd of two numbers below $10^4$, in the second loop only for odd $n$, and a gcd takes $O(\log N)$ division steps. The running time is therefore $O(N\log N)$, less than a second for the C++ program on the computer used for these checks. The memory is $O(1)$: the loops keep only a few integers, and only the recursive gcd of the Java version adds $O(\log N)$ stack frames.

The C++ program checks both values of the statement before the main count, unless --skip-checkpoints is given; the Python and Java versions contain no checks. For this article, a separate brute force compared the squared side lengths exactly, in integers, for every triplet with $b+d\lt400$ and every point $P=(x,a-x)$ with $0\le x\le a$, and for $b+d\lt160$ also with $-400\le x\le a+400$. It found no solution outside the segment $AC$, one point $P$ for each triplet of the incentre family and two for each triplet of the parallel family, and it agrees with the sums of Section 8 and with a direct enumeration of Pythagorean leg pairs for every $N\le400$; so do the three programs, except for $N=11$, $12$, $13$. Finally, an independent count that generates the $15\,915\,492$ primitive triples with hypotenuse below $10^8$ with the Berggren tree ^[9] and adds up their multiples agrees with the C++ program at $N=100$, $100\,000$ and $10^8$. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above. The classification into the incentre and parallel families, their reduction to Pythagorean triangles and the counting sums are derived explicitly in this article.

Project Euler 299 — Three Similar Triangles. The official statement asks for the number of triplets (a, b, d) with b + d < 100 000 000 for which a point P with integer coordinates on the line AC makes the triangles ABP, CDP and BDP similar. It notes that this forces a = c, counts (2, 3, 4) and (2, 4, 3) as distinct, and gives the example P(1, 1) and the counts 92 for b + d < 100 and 320471 for b + d < 100 000.
Wikipedia — Similarity (geometry). Two triangles are similar when their corresponding angles are equal, equivalently when all corresponding sides are proportional, or when two pairs of sides are proportional and the angles between them are equal, the SAS criterion used in Section 5.
Wikipedia — Incircle and excircles. The incentre is the intersection of the internal angle bisectors and has the same distance from all three sides, and the distances from a vertex to the two nearest points of contact are equal; Sections 3 and 4 use this for the triangle OBD.
Wikipedia — Transversal (geometry). By Euclid's Elements, Book I, Propositions 27 and 29, two lines cut by a transversal are parallel exactly when a pair of alternate angles is equal; Section 3 uses both directions.
Wikipedia — Pythagorean triple. Every Pythagorean triple is a multiple of a primitive one, and Euclid's formula (m² − n², 2mn, m² + n²) with coprime m > n > 0, exactly one of them even, gives every primitive triple exactly once up to the order of the legs. Exactly one leg of a primitive triple is even.
Wikipedia — Inscribed angle. An inscribed angle is half of the central angle on the same arc, and on a fixed side of a segment the points from which it is seen at a fixed angle form an arc of a circle; Section 5 uses the arc facing O of the circle with centre (b, b).
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes in exactly one way, up to the order of the factors; hence coprime positive integers whose product is a square are squares, as used in Section 7.
OEIS A198441 lists the sums of the two legs of the primitive Pythagorean triangles with multiplicity, 7, 17, 23, 31, 41, 47, 49, 71, 73, 79, 89, 97, …, and OEIS A020882 lists their hypotenuses with multiplicity, 5, 13, 17, 25, 29, 37, 41, 53, …; Section 9 uses the terms below 100 and below 50.
Wikipedia — Tree of primitive Pythagorean triples. Three matrices, found by B. Berggren in 1934, generate every primitive triple exactly once from (3, 4, 5), and the inverse step leads to a smaller hypotenuse, so the hypotenuse grows from parent to child; the independent check above uses this tree.
C++ — Euler299.cpp. The linked immutable C++ revision contains count_family_x_eq_y, count_family_u_eq_v, count_triplets and run_checkpoints with the two values of the statement. It is the implementation record for this article.

Problem 299 source code

C++

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

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

struct Options {
    u64 limit = 100000000ULL;
    bool run_checkpoints = true;
    unsigned requested_threads = 0U;
};

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 = 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) {
            throw std::overflow_error("u64 parse overflow");
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const std::string& 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())) {
        throw std::overflow_error("unsigned parse overflow");
    }

    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 (parse_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
            continue;
        }

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

unsigned pick_thread_count(const unsigned requested_threads) {
    if (requested_threads > 0U) {
        return requested_threads;
    }

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

template <typename Worker>
u64 parallel_sum_range(const int begin,
                       const int end,
                       unsigned threads,
                       const Worker& worker) {
    if (begin >= end) {
        return 0ULL;
    }

    const int work_items = end - begin;
    if (threads == 0U) {
        threads = 1U;
    }
    if (threads == 1U || work_items < 128) {
        return worker(begin, end);
    }

    threads = std::min<unsigned>(threads, static_cast<unsigned>(work_items));

    std::vector<std::thread> pool;
    pool.reserve(threads);

    std::vector<u64> partial(threads, 0ULL);

    for (unsigned tid = 0U; tid < threads; ++tid) {
        const int lo = begin + (work_items * static_cast<int>(tid)) /
                                   static_cast<int>(threads);
        const int hi = begin +
                       (work_items * static_cast<int>(tid + 1U)) /
                           static_cast<int>(threads);

        pool.emplace_back([&, tid, lo, hi]() { partial[tid] = worker(lo, hi); });
    }

    for (std::thread& th : pool) {
        th.join();
    }

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

u64 count_family_x_eq_y(const u64 limit, unsigned threads) {
    if (limit <= 7ULL) {
        return 0ULL;
    }

    const u64 numerator = limit - 1ULL;
    const int m_max = static_cast<int>(std::sqrt(static_cast<long double>(limit))) + 3;

    auto worker = [limit, numerator](const int lo, const int hi) -> u64 {
        u64 local = 0ULL;

        for (int m = lo; m < hi; ++m) {
            const u64 m64 = static_cast<u64>(m);
            const u64 m2 = m64 * m64;

            // Orientation 1: u = g*m^2, v = 2*g*n^2 (m must be odd).
            if ((m & 1) != 0) {
                for (int n = 1;; ++n) {
                    const u64 n64 = static_cast<u64>(n);
                    const u64 s = m2 + 4ULL * m64 * n64 + 2ULL * n64 * n64;
                    if (s >= limit) {
                        break;
                    }
                    if (std::gcd(m, n) == 1) {
                        local += numerator / s;
                    }
                }
            }

            // Orientation 2: u = 2*g*m^2, v = g*n^2 (n must be odd).
            for (int n = 1;; ++n) {
                const u64 n64 = static_cast<u64>(n);
                const u64 s = 2ULL * m2 + 4ULL * m64 * n64 + n64 * n64;
                if (s >= limit) {
                    break;
                }
                if ((n & 1) != 0 && std::gcd(m, n) == 1) {
                    local += numerator / s;
                }
            }
        }

        return local;
    };

    return parallel_sum_range(1, m_max + 1, threads, worker);
}

u64 count_family_u_eq_v(const u64 limit, unsigned threads) {
    if (limit <= 13ULL) {
        return 0ULL;
    }

    const u64 numerator = limit - 1ULL;
    const int p_max = static_cast<int>(std::sqrt(static_cast<long double>(limit / 2ULL))) + 3;
    const int odd_count = (p_max + 1) / 2;

    auto worker = [limit, numerator](const int lo, const int hi) -> u64 {
        u64 local = 0ULL;

        for (int idx = lo; idx < hi; ++idx) {
            const int p = 2 * idx + 1;
            const u64 p64 = static_cast<u64>(p);
            const u64 p2 = p64 * p64;

            for (int q = 1;; ++q) {
                const u64 q64 = static_cast<u64>(q);
                const u64 s = p2 + 2ULL * p64 * q64 + 2ULL * q64 * q64;
                const u128 denom = static_cast<u128>(2ULL) * static_cast<u128>(s);
                if (denom >= static_cast<u128>(limit)) {
                    break;
                }
                if (std::gcd(p, q) == 1) {
                    local += numerator / static_cast<u64>(denom);
                }
            }
        }

        return local;
    };

    return parallel_sum_range(0, odd_count, threads, worker);
}

u64 count_triplets(const u64 limit, unsigned threads) {
    if (limit <= 0ULL) {
        return 0ULL;
    }

    if (limit < 1000000ULL) {
        threads = 1U;
    }

    const u64 family_x_eq_y = count_family_x_eq_y(limit, threads);
    const u64 family_u_eq_v = count_family_u_eq_v(limit, threads);
    return family_x_eq_y + family_u_eq_v;
}

void run_checkpoints(const unsigned threads) {
    struct Checkpoint {
        u64 limit;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {100ULL, 92ULL},
        {100000ULL, 320471ULL},
    };

    for (const Checkpoint& checkpoint : checkpoints) {
        const u64 got = count_triplets(checkpoint.limit, threads);
        if (got != checkpoint.expected) {
            throw std::runtime_error(
                "Checkpoint failed for limit=" + std::to_string(checkpoint.limit) +
                ": expected " + std::to_string(checkpoint.expected) +
                ", got " + std::to_string(got));
        }
    }
}

} // namespace

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

        const unsigned threads = pick_thread_count(options.requested_threads);

        if (options.run_checkpoints) {
            run_checkpoints(threads);
        }

        const u64 answer = count_triplets(options.limit, threads);
        std::cout << answer << '\n';
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

import math
from concurrent.futures import ProcessPoolExecutor
import multiprocessing

def count_x_eq_y_worker(lo, hi, limit, numerator):
    local = 0
    for m in range(lo, hi):
        m2 = m * m
        if m % 2 != 0:
            n = 1
            while True:
                s = m2 + 4 * m * n + 2 * n * n
                if s >= limit: break
                if math.gcd(m, n) == 1:
                    local += numerator // s
                n += 1
                
        n = 1
        while True:
            s = 2 * m2 + 4 * m * n + n * n
            if s >= limit: break
            if n % 2 != 0 and math.gcd(m, n) == 1:
                local += numerator // s
            n += 1
            
    return local

def count_family_x_eq_y(limit, threads):
    if limit <= 7: return 0
    numerator = limit - 1
    m_max = math.isqrt(limit) + 3
    
    if threads <= 1:
        return count_x_eq_y_worker(1, m_max + 1, limit, numerator)
        
    chunk = (m_max + threads) // threads
    ranges = [(max(1, i * chunk), min(m_max + 1, (i + 1) * chunk)) for i in range(threads)]
    
    total = 0
    with ProcessPoolExecutor(max_workers=threads) as executor:
        futures = [executor.submit(count_x_eq_y_worker, r[0], r[1], limit, numerator) for r in ranges]
        for f in futures:
            total += f.result()
            
    return total

def count_u_eq_v_worker(lo, hi, limit, numerator):
    local = 0
    for idx in range(lo, hi):
        p = 2 * idx + 1
        p2 = p * p
        q = 1
        while True:
            s = p2 + 2 * p * q + 2 * q * q
            denom = 2 * s
            if denom >= limit: break
            if math.gcd(p, q) == 1:
                local += numerator // denom
            q += 1
            
    return local

def count_family_u_eq_v(limit, threads):
    if limit <= 13: return 0
    numerator = limit - 1
    p_max = math.isqrt(limit // 2) + 3
    odd_count = (p_max + 1) // 2
    
    if threads <= 1:
        return count_u_eq_v_worker(0, odd_count, limit, numerator)
        
    chunk = (odd_count + threads - 1) // threads
    ranges = [(i * chunk, min(odd_count, (i + 1) * chunk)) for i in range(threads)]
    
    total = 0
    with ProcessPoolExecutor(max_workers=threads) as executor:
        futures = [executor.submit(count_u_eq_v_worker, r[0], r[1], limit, numerator) for r in ranges]
        for f in futures:
            total += f.result()
            
    return total

def count_triplets(limit):
    threads = multiprocessing.cpu_count()
    if threads == 0: threads = 1
    
    if limit < 1000000:
        threads = 1
        
    total_x = count_family_x_eq_y(limit, threads)
    total_u = count_family_u_eq_v(limit, threads)
    return total_x + total_u

def solve(limit=100000000):
    return str(count_triplets(limit))

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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;

public class Euler299 {
    static int gcd(int a, int b) {
        if (b == 0)
            return a;
        return gcd(b, a % b);
    }

    static long countXEqYWorker(int lo, int hi, long limit, long numerator) {
        long local = 0;
        for (int m = lo; m < hi; ++m) {
            long m2 = (long) m * m;

            if ((m & 1) != 0) {
                for (int n = 1;; ++n) {
                    long s = m2 + 4L * m * n + 2L * n * n;
                    if (s >= limit)
                        break;
                    if (gcd(m, n) == 1) {
                        local += numerator / s;
                    }
                }
            }

            for (int n = 1;; ++n) {
                long s = 2L * m2 + 4L * m * n + (long) n * n;
                if (s >= limit)
                    break;
                if ((n & 1) != 0 && gcd(m, n) == 1) {
                    local += numerator / s;
                }
            }
        }
        return local;
    }

    static long countFamilyXEqY(long limit, int threads) {
        if (limit <= 7)
            return 0;
        long numerator = limit - 1;
        int mMax = (int) Math.sqrt(limit) + 3;

        if (threads <= 1) {
            return countXEqYWorker(1, mMax + 1, limit, numerator);
        }

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();
        int chunk = (mMax + threads) / threads;

        for (int i = 0; i < threads; ++i) {
            final int lo = Math.max(1, i * chunk);
            final int hi = Math.min(mMax + 1, (i + 1) * chunk);
            if (lo >= hi)
                continue;
            futures.add(executor.submit(() -> countXEqYWorker(lo, hi, limit, numerator)));
        }

        long total = 0;
        for (Future<Long> f : futures) {
            try {
                total += f.get();
            } catch (Exception e) {
                e.printStackTrace();
            }
        }
        executor.shutdown();
        return total;
    }

    static long countUEqVWorker(int lo, int hi, long limit, long numerator) {
        long local = 0;
        for (int idx = lo; idx < hi; ++idx) {
            long p = 2L * idx + 1;
            long p2 = p * p;
            for (long q = 1;; ++q) {
                long s = p2 + 2L * p * q + 2L * q * q;
                long denom = 2L * s;
                if (denom >= limit)
                    break;
                if (gcd((int) p, (int) q) == 1) {
                    local += numerator / denom;
                }
            }
        }
        return local;
    }

    static long countFamilyUEqV(long limit, int threads) {
        if (limit <= 13)
            return 0;
        long numerator = limit - 1;
        int pMax = (int) Math.sqrt(limit / 2) + 3;
        int oddCount = (pMax + 1) / 2;

        if (threads <= 1) {
            return countUEqVWorker(0, oddCount, limit, numerator);
        }

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();
        int chunk = (oddCount + threads - 1) / threads;

        for (int i = 0; i < threads; ++i) {
            final int lo = i * chunk;
            final int hi = Math.min(oddCount, (i + 1) * chunk);
            if (lo >= hi)
                continue;
            futures.add(executor.submit(() -> countUEqVWorker(lo, hi, limit, numerator)));
        }

        long total = 0;
        for (Future<Long> f : futures) {
            try {
                total += f.get();
            } catch (Exception e) {
                e.printStackTrace();
            }
        }
        executor.shutdown();
        return total;
    }

    public static String solve() {
        long limit = 100000000L;
        int threads = Runtime.getRuntime().availableProcessors();
        if (threads == 0)
            threads = 1;

        long totalX = countFamilyXEqY(limit, threads);
        long totalU = countFamilyUEqV(limit, threads);

        return String.valueOf(totalX + totalU);
    }

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