Problem 279: Triangles with Integral Sides and an Integral Angle

View on Project Euler

Project Euler Problem 279 Solution

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

Problem Summary We must count all integer-sided triangles with perimeter at most \(L\) that contain an angle of exactly \(60^\circ\), \(90^\circ\), or \(120^\circ\). The final Project Euler total is intentionally omitted; this page explains the parameterizations and counting logic used by the C++ solution. Mathematical Approach 1) Three disjoint angle classes Let $$N(L)=N_{60}(L)+N_{90}(L)+N_{120}(L).$$ For integer-sided triangles these classes are disjoint: \(90^\circ\) and \(120^\circ\) cannot both occur because their sum already exceeds \(180^\circ\). \(60^\circ\) and \(120^\circ\) also cannot both occur because that would leave a third angle of \(0^\circ\). \(60^\circ\) and \(90^\circ\) would force a \(30\)-\(60\)-\(90\) triangle, whose side ratio is \(1:\sqrt3:2\). Since \(\sqrt3\) is irrational, no positive integer scaling can make all three sides integers. So every valid triangle belongs to exactly one of the three counters. Equilateral triangles belong only to the \(60^\circ\) family. 2) Law of cosines gives three Diophantine equations If side \(c\) is opposite the special angle \(\theta\), then $$c^2=a^2+b^2-2ab\cos\theta.$$ For the three angles in the problem: $$90^\circ:\ c^2=a^2+b^2,$$ $$60^\circ:\ c^2=a^2+b^2-ab,$$ $$120^\circ:\ c^2=a^2+b^2+ab.$$ So Project Euler 279 is really three separate integer-parameterization problems plus a scaling count....

Detailed mathematical approach

Problem Summary

We must count the triangles whose three sides have integer lengths, whose perimeter is at most \(10^8\), and which have at least one angle that is a whole number of degrees. A triangle is determined by its side lengths, so we count the triples \(a\le b\le c\) of positive integers with \(a+b\gt c\) and \(a+b+c\le L\), where \(L=10^8\). The statement gives no values for smaller limits; Section 11 works out the count for \(L=100\), which the published C++ program also checks. [1]

The whole solution in one sentence: every angle of a triangle with integer sides has a rational cosine, and the only whole-degree angles with a rational cosine are \(60^\circ\), \(90^\circ\) and \(120^\circ\), so we list the primitive right triangles with Euclid's formula and the primitive \(120^\circ\) triangles with a similar formula, obtain the \(60^\circ\) triangles from the \(120^\circ\) ones by attaching an equilateral triangle, and count every primitive triangle together with all its multiples by one division.

Mathematical Approach

We first find out which angles can occur at all, then sort the triangles into four classes, find a formula for the primitive triangles of each class, and finally count all their multiples at once.

1. Small perimeters

For perimeters up to \(15\) the triangles with a whole-degree angle can be listed by hand. The equilateral triangles \((1,1,1)\) to \((5,5,5)\) have three angles of \(60^\circ\). The triangle \((3,4,5)\) with perimeter \(12\) has a right angle, and \((3,5,7)\) with perimeter \(15\) has an angle of \(120^\circ\). To see this we use the law of cosines: for the angle \(\gamma\) opposite the side \(c\), [2]

$$\cos\gamma=\frac{a^2+b^2-c^2}{2ab},$$

and for \((3,5,7)\) the angle opposite \(7\) has \(\cos\gamma=(9+25-49)/30=-1/2\). So the counts for the limits \(10\), \(12\) and \(15\) are \(3\), \(5\) and \(7\). The first \(60^\circ\) angle outside an equilateral triangle appears in \((3,7,8)\), perimeter \(18\), where the angle opposite \(7\) has cosine \((9+64-49)/48=1/2\).

A floating-point scan of all \(6\,986\,111\) triangles with perimeter at most \(1000\) finds no angle within \(10^{-9}\) degree of a whole number other than \(60^\circ\), \(90^\circ\) and \(120^\circ\), except for a single one that Section 3 unmasks.

2. Only \(60^\circ\), \(90^\circ\) and \(120^\circ\) can occur

If the sides are integers, the formula above makes the cosine of every angle a rational number. [3] Now suppose that \(\gamma\) is a whole number \(d\) of degrees, \(0\lt d\lt 180\), and put \(x=2\cos\gamma\). Define polynomials by

$$\begin{aligned}P_0&=2,\qquad P_1=x,\\ P_{k+1}&=x\,P_k-P_{k-1}.\end{aligned}$$

Then \(P_k(2\cos\gamma)=2\cos(k\gamma)\) for every \(k\ge0\): this holds for \(k=0\) and \(k=1\), and the step from \(k\) to \(k+1\) is the product formula \(\cos(A+B)+\cos(A-B)\)\({}=2\cos A\cos B\) with \(A=k\gamma\) and \(B=\gamma\). For example \(P_2=x^2-2\) and \(P_3=x^3-3x\); these are rescaled Chebyshev polynomials, \(P_k(x)=2T_k(x/2)\). [4] By induction every \(P_k\) with \(k\ge1\) has integer coefficients and leading coefficient \(1\). Since \(360\gamma\) is \(d\) full turns, \(P_{360}(x)=2\cos(360\gamma)=2\), so \(x\) is a root of \(P_{360}-2\), a polynomial with integer coefficients and leading coefficient \(1\).

Such a polynomial has no rational roots except integers. Indeed, if \(p/q\) in lowest terms with \(q\ge1\) is a root of \(x^n+c_{n-1}x^{n-1}+\dots+c_0\), multiplying by \(q^n\) gives

$$p^n=-q\left(c_{n-1}p^{n-1}+\dots+c_0\,q^{n-1}\right),$$

so every prime factor of \(q\) would divide \(p\), and therefore \(q=1\). This is the rational root theorem. [5] Hence \(2\cos\gamma\) is an integer, and \(-2\lt 2\cos\gamma\lt 2\) leaves only \(1\), \(0\) and \(-1\): the angle is \(60^\circ\), \(90^\circ\) or \(120^\circ\). This is the cosine case of Niven's theorem. [6] With \(\cos\gamma=1/2\), \(0\), \(-1/2\) the law of cosines turns the three angles into three equations, in which \(a\) and \(b\) are the sides next to the angle and \(c\) is the side opposite it, not necessarily the longest:

$$\begin{aligned}60^\circ&:\quad c^2=a^2-ab+b^2,\\ 90^\circ&:\quad c^2=a^2+b^2,\\ 120^\circ&:\quad c^2=a^2+ab+b^2.\end{aligned}$$

3. Angles that are almost whole

The exception of Section 1 is the triangle \((209,362,418)\). Its smallest angle is about \(29.9999999993^\circ\), less than \(10^{-9}\) degree below \(30^\circ\), so a floating-point test with that tolerance accepts it, although Section 2 rules out \(30^\circ\). The triangle is almost a \(30^\circ\)–\(60^\circ\)–\(90^\circ\) triangle, whose sides have the ratio \(1:\sqrt3:2\), because \(362/209\) is a continued fraction convergent of \(\sqrt3\), one of its best approximations by fractions. [7] Later convergents come much closer: the triangle \((7865521,13623482,15731042)\), with perimeter \(37\,220\,045\), has an angle that differs from \(30^\circ\) by about \(4\cdot10^{-28}\) degree, and double-precision arithmetic returns \(30.000000000000004^\circ\). So angles must be decided by the exact integer equations above, as the programs do.

4. At most one special angle

The angles of a triangle add up to \(180^\circ\). Two angles from \(90^\circ\) and \(120^\circ\), or \(60^\circ\) together with \(120^\circ\), already use up \(180^\circ\) or more. A \(60^\circ\) angle and a \(90^\circ\) angle would leave \(30^\circ\) for the third angle, which Section 2 excludes. So the only triangles with two special angles have two angles of \(60^\circ\), and then the third is \(60^\circ\) as well: they are equilateral. Every triangle we count therefore belongs to exactly one of four classes: equilateral, exactly one angle of \(60^\circ\), exactly one angle of \(90^\circ\), exactly one angle of \(120^\circ\). With \(E(L)\), \(S(L)\), \(R(L)\) and \(O(L)\) the sizes of the classes,

$$\begin{aligned}N(L)&=E(L)+S(L)\\ &\quad+R(L)+O(L),\\ E(L)&=\left\lfloor\frac L3\right\rfloor,\end{aligned}$$

because the equilateral triangles are \((k,k,k)\) with \(3k\le L\). The programs keep \(E+S\) as their \(60^\circ\) count, \(R\) as the \(90^\circ\) count and \(O\) as the \(120^\circ\) count.

5. Primitive triangles and their multiples

A triangle is primitive if \(\gcd(a,b,c)=1\). Multiplying all sides by \(k\) keeps the angles and multiplies the perimeter by \(k\), and every triangle is \(k\) times exactly one primitive triangle, with \(k=\gcd(a,b,c)\). A primitive triangle with perimeter \(p\) therefore has \(\lfloor L/p\rfloor\) multiples with perimeter at most \(L\), and each class is counted as

$$R(L)=\sum_{p}\left\lfloor\frac{L}{p}\right\rfloor,$$

with one term for each primitive right triangle, \(p\) being its perimeter, and likewise for \(S\) and \(O\). For example \((3,4,5)\) has the \(\lfloor100/12\rfloor=8\) multiples \((3,4,5)\), \((6,8,10)\), \(\dots\), \((24,32,40)\) with perimeter at most \(100\). It remains to list the primitive triangles of each class exactly once.

6. Right angles: Euclid's formula

The primitive right triangles are the primitive Pythagorean triples. Euclid's formula says that for coprime integers \(m\gt n\ge1\) with \(m-n\) odd,

$$\begin{aligned}a&=m^2-n^2,\qquad b=2mn,\\ c&=m^2+n^2\end{aligned}$$

is a primitive triple, and that every primitive triple arises from exactly one such pair, with \(b\) the even leg. [8] Every primitive triple does have this form. Exactly one leg is even, because two odd squares add up to \(2\) modulo \(4\), which no square does, and two even legs would make \(c\) even. With \(b\) even, \(u=(c+a)/2\) and \(v=(c-a)/2\) are positive integers with \(uv=(b/2)^2\), and they are coprime, since a common prime would divide \(u+v=c\), \(u-v=a\) and then \(b\). Coprime numbers whose product is a square are squares themselves, so \(u=m^2\) and \(v=n^2\) with coprime \(m\gt n\), of opposite parity because \(a=m^2-n^2\) is odd; then \(c=u+v\), \(a=u-v\) and \(b=2mn\) are Euclid's formula, and the triple determines the pair. The perimeter is \(2m(m+n)\); for example \((m,n)=(2,1)\) gives \((3,4,5)\) with perimeter \(12\).

7. A formula for \(120^\circ\) triangles

Let \(a\), \(b\) and \(c\) be positive integers with \(c^2=a^2+ab+b^2\), which we write as \((c-a)(c+a)=b(a+b)\). Then \(c\gt a\), and \(c\lt a+b\) because \((a+b)^2=c^2+ab\). So the number

$$r=\frac{c-a}{a+b}=\frac{b}{c+a}$$

lies strictly between \(0\) and \(1\), and in lowest terms \(r=n/m\) with coprime \(m\gt n\ge1\). Consider the two equations \(m(c-a)=n(a+b)\) and \(mb=n(c+a)\). Multiplying the first by \(n\) and replacing \(nc\) by \(mb-na\) from the second gives \(b(m^2-n^2)=a(2mn+n^2)\), and then \(nc=mb-na\) gives \(c\). Hence

$$\begin{aligned}a:b:c&=(m^2-n^2):(2mn+n^2)\\ &\quad:(m^2+mn+n^2).\end{aligned}$$

Conversely, for coprime \(m\gt n\ge1\) these three numbers satisfy the equation, as expanding both sides shows, and their ratio \(r\) is \(n(m+2n)/(m(m+2n))=n/m\). So the ordered solutions \((a,b,c)\) without a common factor correspond one-to-one to the fractions \(n/m\) between \(0\) and \(1\): a solution determines \(r\), and \(r=n/m\) gives back the solution as the triple above divided by the greatest common divisor \(g\) of its three entries. The perimeter of the triple is \((m+n)(2m+n)\).

This is the \(120^\circ\) analogue of Euclid's formula. For the triple above, \((m-n\omega)^2=a-b\omega\), where \(\omega\) is a primitive cube root of unity, and the norm \(N(x+y\omega)=x^2-xy+y^2\) of the Eisenstein integers, the squared absolute value, is multiplicative; so \(N(a-b\omega)=N(m-n\omega)^2\), which is the identity \(a^2+ab+b^2=c^2\). [9]

8. Every primitive \(120^\circ\) triangle exactly once

A triangle is unordered, but Section 7 lists ordered solutions: \((a,b,c)\) and \((b,a,c)\) describe the same triangle, and they differ, because \(a=b\) would give \(c^2=3a^2\), which is impossible since \(\sqrt3\) is irrational. Two facts decide which orientation the formula produces without a common factor.

  • The common factor. A prime \(p\) dividing \(a=(m-n)(m+n)\), \(b=n(2m+n)\) and \(c\) also divides \(a+b=m(m+2n)\). If \(p\mid m\), then \(p\mid b\) forces \(p\mid n\); if \(p\mid n\), then \(p\mid a+b\) forces \(p\mid m\); both contradict \(\gcd(m,n)=1\). So \(p\mid 2m+n\) and \(p\mid m+2n\), hence \(p\mid 2(m+2n)-(2m+n)=3n\), and \(p=3\). As \(m+2n\equiv m-n\pmod 3\), this happens only if \(m\equiv n\pmod 3\); conversely, \(m\equiv n\pmod 3\) makes \(3\) divide \(m-n\), \(2m+n\) and \(c=(m-n)^2+3mn\). So \(g\gt1\) exactly when \(m\equiv n\pmod 3\). In that case \(3\) divides neither \(m\) nor \(n\), so \(mn\equiv n^2\equiv1\pmod 3\) and \(c=(m-n)^2+3mn\equiv3\pmod 9\); hence \(g=3\). [3]
  • The other orientation. For the triple of the pair \((m,n)\), the swapped solution has the ratio \((c-b)/(a+b)\)\({}=(m-n)/(m+2n)\). Since \(\gcd(m-n,m+2n)\)\({}=\gcd(m-n,3n)\)\({}=\gcd(m-n,3)\), the swapped solution comes from the pair \((m+2n,\,m-n)\) if \(m\not\equiv n\pmod 3\), and from \(((m+2n)/3,\,(m-n)/3)\) if \(m\equiv n\pmod 3\). In the first case the new pair has the difference \(3n\), so its entries are congruent modulo \(3\). In the second case its difference is \(n\), which is not divisible by \(3\), because \(3\mid n\) and \(3\mid m-n\) would make \(3\) divide \(m\).

So of the two orientations of a primitive \(120^\circ\) triangle exactly one comes from a pair with \(m\not\equiv n\pmod 3\), and for that pair \(g=1\): the formula gives the triangle itself. The primitive \(120^\circ\) triangles therefore correspond one-to-one to the coprime pairs \(m\gt n\ge1\) with \(m\not\equiv n\pmod 3\). For instance \((4,1)\) gives \((15,9,21)=3\cdot(5,3,7)\), while the other orientation \((3,5,7)\) has the ratio \(4/8=1/2\) and comes directly from \((2,1)\).

9. \(60^\circ\) triangles from \(120^\circ\) triangles

Take a triangle \(ABC\) with a \(120^\circ\) angle at \(C\), \(CB=a\), \(CA=b\) and \(AB=c\). Extend the side \(AC\) beyond \(C\) to the point \(D\) with \(CD=a\). The angle \(BCD\) is \(180^\circ-120^\circ=60^\circ\) and \(CB=CD\), so \(BCD\) is equilateral, and the triangle \(ABD\) has the sides \(BD=a\), \(AD=a+b\), \(AB=c\) and a \(60^\circ\) angle at \(D\). In formulas, \(a^2+ab+b^2=a^2-a(a+b)+(a+b)^2\). Extending the other side instead gives the \(60^\circ\) triangle \((b,a+b,c)\).

Conversely, in a \(60^\circ\) triangle that is not equilateral, the two sides next to the \(60^\circ\) angle are different, since equal sides \(x\) would give \(c^2=x^2\); call them \(x\lt y\). Cutting off the equilateral triangle with side \(x\) leaves the \(120^\circ\) triangle \((x,\,y-x,\,c)\), and attaching it again returns \((x,y,c)\). So every non-equilateral \(60^\circ\) triangle comes from exactly one \(120^\circ\) triangle together with the choice of the side that is kept, and both have the same common factor, because \(\gcd(a,a+b,c)=\gcd(a,b,c)\). Summing over the primitive \(120^\circ\) triangles,

$$\begin{aligned}S(L)&=\sum\left\lfloor\frac{L}{2a+b+c}\right\rfloor\\ &\quad+\sum\left\lfloor\frac{L}{a+2b+c}\right\rfloor.\end{aligned}$$

For example \((3,5,7)\) yields \((3,8,7)\) and \((5,8,7)\), the two smallest non-equilateral \(60^\circ\) triangles, with perimeters \(18\) and \(20\).

10. The \(60^\circ\) formulas of the programs

The programs do not reuse the \(120^\circ\) pairs for this class. They loop over coprime \(m\) and \(n\) with \(1\le n\le(m-1)/2\) and use

$$\begin{aligned}a&=m^2-n^2,\quad c=m^2-mn+n^2,\\ b_1&=2mn-n^2,\quad b_2=m^2-2mn.\end{aligned}$$

Both \((a,b_1,c)\) and \((a,b_2,c)\) satisfy \(c^2=a^2-ab+b^2\); as \(b_1+b_2=a\), they are the two roots \(b\) of the quadratic \(b^2-ab+a^2-c^2=0\). The substitution \(m=m'+n'\), \(n=n'\) explains them: it turns \(a\), \(b_1\), \(b_2\), \(c\) into \(a'+b'\), \(b'\), \(a'\), \(c'\), where \((a',b',c')\) is the \(120^\circ\) triple of the pair \((m',n')\). So the two branches are exactly the two extensions of Section 9. The condition \(n\le(m-1)/2\) becomes \(n'\lt m'\), \(\gcd(m,n)=\gcd(m',n')\), and both branches have the common factor of \((a',b',c')\). Keeping the branches with \(\gcd(a,b,c)=1\) therefore meets every non-equilateral primitive \(60^\circ\) triangle exactly once. Their perimeters are \((m+n)(2m-n)\) and \(3m(m-n)\).

11. Worked example: \(L=100\)

There are \(E(100)=33\) equilateral triangles. For the other classes we list the primitive triangles with perimeter at most \(100\) and add the quotients \(\lfloor100/p\rfloor\):

  • Right angles. The pairs \((2,1)\), \((3,2)\), \((4,1)\), \((4,3)\), \((5,2)\), \((6,1)\) and \((5,4)\) give \((3,4,5)\), \((5,12,13)\), \((15,8,17)\), \((7,24,25)\), \((21,20,29)\), \((35,12,37)\) and \((9,40,41)\), with perimeters \(12\), \(30\), \(40\), \(56\), \(70\), \(84\) and \(90\). So \(R(100)=8+3+2+1\)\({}+1+1+1=17\).
  • \(120^\circ\). The pairs \((2,1)\), \((3,1)\), \((3,2)\), \((5,1)\), \((4,3)\) and \((6,1)\) give \((3,5,7)\), \((8,7,13)\), \((5,16,19)\), \((24,11,31)\), \((7,33,37)\) and \((35,13,43)\), with perimeters \(15\), \(28\), \(40\), \(66\), \(77\) and \(91\). So \(O(100)=6+3+2\)\({}+1+1+1=14\). The pairs \((4,1)\) and \((5,2)\) are skipped, since they give \(3\cdot(5,3,7)\) and \(3\cdot(7,8,13)\).
  • \(60^\circ\), not equilateral. The two extensions of these six triangles have the perimeters \(18\) and \(20\), \(36\) and \(35\), \(45\) and \(56\), \(90\) and \(77\), \(84\) and \(110\), \(126\) and \(104\). So \(S(100)=5+5+2+2+2\)\({}+1+1+1+1=20\).

Altogether \(N(100)=33+20+17+14=84\), of which \(53=33+20\) have a \(60^\circ\) angle. The C++ program checks exactly these values.

12. How far \(m\) has to run

For fixed \(m\) the perimeters \(2m(m+n)\) and \((m+n)(2m+n)\) grow with \(n\), so they are at least their values \(2m(m+1)\) and \(2m^2+3m+1\) at \(n=1\). In the \(60^\circ\) loop, \((m+n)(2m-n)\) grows on \(1\le n\le(m-1)/2\), since one step in \(n\) adds \(m-2n-1\ge0\), so its smallest value is \(2m^2+m-1\); and \(3m(m-n)\) is smallest at the largest \(n\). The code bounds the latter from below by \(3m\lfloor(m+1)/2\rfloor\). This is the exact minimum for odd \(m\); for even \(m\) the minimum is \(3m(m/2+1)\), and the smaller bound is harmless. All these bounds grow with \(m\), so once they exceed \(L\) no larger \(m\) can contribute. For \(L=10^8\) the three loops end at \(m=7070\), \(7070\) and \(8164\).

How the Code Works

The C++ program [10] keeps the three counts of Section 4 in a struct Counts with the fields with_60, with_90 and with_120. main reads the options --limit= (default \(10^8\), at least \(3\)) and --skip-checkpoints, runs run_checkpoints unless told to skip them, and prints the total of count_by_parameterization(limit).

count_by_parameterization finds the ends of the three loops with max_m_90, max_m_120 and max_m_60, which increase \(m\) while the bounds of Section 12 are at most the limit, sets with_60 to limit / 3 and adds the results of three range functions. These take a first value, a last value and a step for \(m\): the driver cuts the range of \(m\) into interleaved slices \(m=2+t,\,2+t+s,\,\dots\) for \(t=0,\dots,s-1\) and adds their partial counts, the same sum as one sweep with step \(1\), since every \(m\) lies in exactly one slice.

  • count_90_range skips the pairs with \(m-n\) even or \(\gcd(m,n)\ne1\) and adds limit / (2*m*(m+n)) (Sections 5 and 6).
  • count_120_range skips the pairs with \(\gcd(m,n)\ne1\), forms \(a\), \(b\), \(c\) of Section 7, skips the triple when gcd3(a, b, c) is not \(1\), which by Section 8 is a disguised test of \(m\equiv n\pmod 3\), and adds limit divided by the perimeter.
  • count_60_non_eq_range runs over the \(n\) with \(1\le n\le\lfloor(m-1)/2\rfloor\) and \(\gcd(m,n)=1\), forms a, c, b_a \(=b_1\) and b_b \(=b_2\) of Section 10, and adds a term for each branch whose triple passes the same gcd3 test. Its extra test b_b > 0 always holds, because \(b_2=m(m-2n)\ge m\).

The inner loops always run over all \(n\); a pair whose perimeter exceeds the limit adds the quotient \(0\). Sides and perimeters are 64-bit integers i64 below \(3\cdot10^8\), the counts are u64, and max_m_60 evaluates its bounds in __int128, a margin this limit does not need.

count_by_bruteforce is the independent reference: for every triangle \(a\le b\le c\) with \(a+b\gt c\) within the limit it takes each side in turn as \(z\), compares \(z^2\) exactly with \(x^2+y^2\), \(x^2+y^2-xy\) and \(x^2+y^2+xy\) for the other two sides \(x\), \(y\), and counts the triangle once for each kind of angle found. run_checkpoints demands the values \(84\), \(53\), \(17\) and \(14\) of Section 11 at the limit \(100\) and compares the fast method with the brute force, count by count, at \(120\), \(200\) and \(300\); a failure ends the program with exit code \(2\).

The Python and Java versions have the same range functions, bounds and slices, with their own Euclidean gcd. Python's solve(limit) uses arbitrary-precision integers and returns a string; Java's solve() uses long with the limit fixed at \(10^8\). Only the C++ program contains the checkpoints and the brute force.

Complexity and Verification

For each \(m\) the right-angle and \(120^\circ\) loops run over \(m-1\) values of \(n\), and the \(60^\circ\) loop over \(\lfloor(m-1)/2\rfloor\). With the loop ends of Section 12 this gives \(24\,988\,915\), \(24\,988\,915\) and \(16\,658\,642\) inner iterations for \(L=10^8\), close to \(L/4+L/4+L/6=2L/3\), because the loop ends are close to \(\sqrt{L/2}\) and \(\sqrt{2L/3}\). Each iteration computes at most a few greatest common divisors of numbers of size \(O(L)\), each in \(O(\log L)\) steps (Euclid's algorithm in Python and Java, std::gcd in C++). With machine-word operations at unit cost, the time is \(O(L\log L)\), and only a constant number of integers is stored. Counting both \(60^\circ\) extensions inside the \(120^\circ\) loop, as in Section 9, would remove the separate \(60^\circ\) loop, whose second perimeter \(3m(m-n)\) falls as \(n\) grows, and let the two remaining loops stop at the first perimeter above \(L\), roughly halving the iterations; the published programs keep the full loops.

Besides the checkpoints of the C++ program, the following checks were run for this article; they are not part of the published solutions:

  • An exact search that solves each of the three equations of Section 2 for \(c\) agrees with the C++ program for every limit \(L\) from \(3\) to \(10^4\). Up to that perimeter the three loops produce each of the \(703\) primitive right triangles, the \(655\) primitive \(120^\circ\) triangles and the \(1053\) non-equilateral primitive \(60^\circ\) triangles exactly once.
  • For all \(683\,781\) coprime pairs with \(m\lt1500\), the common factor \(g\) of the \(120^\circ\) triple is \(1\) or \(3\), and it is \(1\) exactly when \(m\not\equiv n\pmod 3\).
  • A separate implementation with these shortcuts and the test \(m\not\equiv n\pmod 3\) agrees with the C++ program at \(L=10^3\), \(12\,345\), \(10^6\), \(10^7\) and \(10^8\), and the Python and Java versions agree with the C++ program at \(L=10^8\).

Footnotes and References

The references below supply the background facts used above; the proofs of Sections 2 and 6 to 10 and the iteration counts are derived explicitly in this article.

  1. Project Euler 279 — Triangles with Integral Sides and an Integral Angle. The official statement asks how many triangles there are with integral sides, at least one integral angle measured in degrees, and a perimeter that does not exceed 10⁸.
  2. Wikipedia — Law of cosines. The law c² = a² + b² − 2ab cos γ, with γ the angle opposite the side c; for a right angle it is the Pythagorean theorem.
  3. Wikipedia — Integer triangle. Every angle of an integer triangle has a rational cosine, and a rational cosine belongs to a rational angle only when it is 0, ±1/2 or ±1. For the 120° triangle with sides m² − n², 2mn + n² and m² + mn + n², with coprime m > n, the greatest common divisor of the sides is 3 if m ≡ n (mod 3) and 1 otherwise.
  4. Wikipedia — Chebyshev polynomials. Tₙ(cos θ) = cos(nθ), with T₀ = 1, T₁ = x and Tₙ₊₁(x) = 2x Tₙ(x) − Tₙ₋₁(x); the article also lists the rescaled family 2Tₙ(x/2), the Pₙ of Section 2.
  5. Wikipedia — Rational root theorem. If p/q in lowest terms is a root of a polynomial with integer coefficients, p divides the constant term and q the leading coefficient; for leading coefficient 1 every rational root is an integer.
  6. Wikipedia — Niven's theorem. For a rational number of degrees, the only rational values of the sine or cosine are 0, ±1/2 and ±1.
  7. OEIS A002531 and OEIS A002530. The numerators 1, 2, 5, 7, 19, 26, 71, 97, 265, 362, … and denominators 1, 1, 3, 4, 11, 15, 41, 56, 153, 209, … of the continued fraction convergents to √3.
  8. Wikipedia — Pythagorean triple. Euclid's formula gives a primitive triple exactly when m and n are coprime and exactly one of them is even, and every primitive triple comes from a unique such pair.
  9. Wikipedia — Eisenstein integer. The Eisenstein integers are the numbers x + yω with integers x, y and a primitive cube root of unity ω; their norm is N(x + yω) = x² − xy + y².
  10. C++ — Euler279.cpp. The linked immutable C++ revision contains the three range functions, their bounds and the checks in run_checkpoints. It is the implementation record for this article.

Problem 279 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <pthread.h>
#include <string>
#include <unistd.h>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;

struct Options {
    i64 limit = 100000000LL;
    bool run_checkpoints = true;
};

bool parse_i64_after_prefix(const std::string& arg, const std::string& prefix, i64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    i64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<i64>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(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_i64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 3;
}

u64 gcd3(const u64 a, const u64 b, const u64 c) {
    return std::gcd(a, std::gcd(b, c));
}

struct Counts {
    u64 with_60 = 0;
    u64 with_90 = 0;
    u64 with_120 = 0;

    u64 total() const {
        return with_60 + with_90 + with_120;
    }
};

u64 count_90_range(const i64 limit, const i64 m_start, const i64 m_end, const i64 m_step) {
    u64 out = 0;
    for (i64 m = m_start; m <= m_end; m += m_step) {
        for (i64 n = 1; n < m; ++n) {
            if (((m - n) & 1LL) == 0 || std::gcd(m, n) != 1) {
                continue;
            }
            const i64 primitive_perimeter = 2 * m * (m + n);
            out += static_cast<u64>(limit / primitive_perimeter);
        }
    }
    return out;
}

u64 count_120_range(const i64 limit, const i64 m_start, const i64 m_end, const i64 m_step) {
    u64 out = 0;
    for (i64 m = m_start; m <= m_end; m += m_step) {
        for (i64 n = 1; n < m; ++n) {
            if (std::gcd(m, n) != 1) {
                continue;
            }
            const i64 mm = m * m;
            const i64 nn = n * n;
            const i64 a = mm - nn;
            const i64 b = 2 * m * n + nn;
            const i64 c = mm + m * n + nn;
            if (gcd3(static_cast<u64>(a), static_cast<u64>(b), static_cast<u64>(c)) != 1) {
                continue;
            }
            const i64 primitive_perimeter = a + b + c;
            out += static_cast<u64>(limit / primitive_perimeter);
        }
    }
    return out;
}

u64 count_60_non_eq_range(const i64 limit, const i64 m_start, const i64 m_end, const i64 m_step) {
    u64 out = 0;
    for (i64 m = m_start; m <= m_end; m += m_step) {
        const i64 mm = m * m;
        const i64 n_max = (m - 1) / 2;
        for (i64 n = 1; n <= n_max; ++n) {
            if (std::gcd(m, n) != 1) {
                continue;
            }
            const i64 nn = n * n;
            const i64 a = mm - nn;
            const i64 c = mm - m * n + nn;

            const i64 b_a = 2 * m * n - nn;
            if (gcd3(static_cast<u64>(a), static_cast<u64>(b_a), static_cast<u64>(c)) == 1) {
                const i64 primitive_perimeter = a + b_a + c;
                out += static_cast<u64>(limit / primitive_perimeter);
            }

            const i64 b_b = mm - 2 * m * n;
            if (b_b > 0 && gcd3(static_cast<u64>(a), static_cast<u64>(b_b), static_cast<u64>(c)) == 1) {
                const i64 primitive_perimeter = a + b_b + c;
                out += static_cast<u64>(limit / primitive_perimeter);
            }
        }
    }
    return out;
}

i64 max_m_90(const i64 limit) {
    i64 m = 2;
    while (2 * m * (m + 1) <= limit) {
        ++m;
    }
    return m - 1;
}

i64 max_m_120(const i64 limit) {
    i64 m = 2;
    while (2 * m * m + 3 * m + 1 <= limit) {
        ++m;
    }
    return m - 1;
}

i64 max_m_60(const i64 limit) {
    i64 m = 2;
    while (true) {
        const __int128 mm = static_cast<__int128>(m) * static_cast<__int128>(m);
        const __int128 min_perimeter_a = 2 * mm + m - 1;
        const __int128 min_perimeter_b = static_cast<__int128>(3) * m * ((m + 1) / 2);
        if (min_perimeter_a > limit && min_perimeter_b > limit) {
            break;
        }
        ++m;
    }
    return m - 1;
}

struct WorkerCtx {
    i64 limit = 0;
    i64 m_max_90 = 0;
    i64 m_max_120 = 0;
    i64 m_max_60 = 0;
    i64 m_start = 0;
    i64 m_step = 0;
    u64 c90 = 0;
    u64 c120 = 0;
    u64 c60 = 0;
};

void* worker_main(void* ptr) {
    auto* ctx = static_cast<WorkerCtx*>(ptr);
    ctx->c90 = count_90_range(ctx->limit, ctx->m_start, ctx->m_max_90, ctx->m_step);
    ctx->c120 = count_120_range(ctx->limit, ctx->m_start, ctx->m_max_120, ctx->m_step);
    ctx->c60 = count_60_non_eq_range(ctx->limit, ctx->m_start, ctx->m_max_60, ctx->m_step);
    return nullptr;
}

unsigned choose_thread_count() {
    long cpu = sysconf(_SC_NPROCESSORS_ONLN);
    unsigned threads = (cpu > 0) ? static_cast<unsigned>(cpu) : 1U;
    if (threads > 8U) {
        threads = 8U;
    }
    if (threads == 0U) {
        threads = 1U;
    }
    return threads;
}

Counts count_by_parameterization(const i64 limit) {
    Counts counts;
    const i64 m_max_90 = max_m_90(limit);
    const i64 m_max_120 = max_m_120(limit);
    const i64 m_max_60 = max_m_60(limit);

    counts.with_60 = static_cast<u64>(limit / 3);
    const unsigned thread_count = choose_thread_count();
    if (thread_count <= 1U) {
        counts.with_90 = count_90_range(limit, 2, m_max_90, 1);
        counts.with_120 = count_120_range(limit, 2, m_max_120, 1);
        counts.with_60 += count_60_non_eq_range(limit, 2, m_max_60, 1);
        return counts;
    }

    std::vector<pthread_t> threads(thread_count);
    std::vector<WorkerCtx> ctx(thread_count);
    unsigned created = 0U;
    bool failed = false;
    for (unsigned t = 0; t < thread_count; ++t) {
        ctx[t] = WorkerCtx{limit,
                           m_max_90,
                           m_max_120,
                           m_max_60,
                           static_cast<i64>(2 + t),
                           static_cast<i64>(thread_count),
                           0ULL,
                           0ULL,
                           0ULL};
        if (pthread_create(&threads[t], nullptr, worker_main, &ctx[t]) != 0) {
            failed = true;
            break;
        }
        ++created;
    }
    for (unsigned t = 0; t < created; ++t) {
        pthread_join(threads[t], nullptr);
    }

    if (failed) {
        counts.with_90 = count_90_range(limit, 2, m_max_90, 1);
        counts.with_120 = count_120_range(limit, 2, m_max_120, 1);
        counts.with_60 += count_60_non_eq_range(limit, 2, m_max_60, 1);
        return counts;
    }

    for (const WorkerCtx& w : ctx) {
        counts.with_90 += w.c90;
        counts.with_120 += w.c120;
        counts.with_60 += w.c60;
    }

    return counts;
}

Counts count_by_bruteforce(const int limit) {
    Counts counts;

    for (int a = 1; a <= limit / 3; ++a) {
        for (int b = a; b <= (limit - a) / 2; ++b) {
            const int c_max = limit - a - b;
            for (int c = b; c <= c_max; ++c) {
                if (a + b <= c) {
                    continue;
                }

                const std::array<i64, 3> sides = {a, b, c};
                bool has_60 = false;
                bool has_90 = false;
                bool has_120 = false;

                for (int i = 0; i < 3; ++i) {
                    const i64 z = sides[static_cast<std::size_t>(i)];
                    const i64 x = sides[static_cast<std::size_t>((i + 1) % 3)];
                    const i64 y = sides[static_cast<std::size_t>((i + 2) % 3)];
                    const i64 zz = z * z;
                    const i64 xx_yy = x * x + y * y;

                    if (zz == xx_yy) {
                        has_90 = true;
                    }
                    if (zz == xx_yy - x * y) {
                        has_60 = true;
                    }
                    if (zz == xx_yy + x * y) {
                        has_120 = true;
                    }
                }

                counts.with_60 += static_cast<u64>(has_60);
                counts.with_90 += static_cast<u64>(has_90);
                counts.with_120 += static_cast<u64>(has_120);
            }
        }
    }

    return counts;
}

bool run_checkpoints() {
    {
        const Counts fast = count_by_parameterization(100);
        if (fast.total() != 84 || fast.with_60 != 53 || fast.with_90 != 17 || fast.with_120 != 14) {
            std::cerr << "Checkpoint failed at limit=100" << '\n';
            return false;
        }
    }

    for (const int limit : {120, 200, 300}) {
        const Counts brute = count_by_bruteforce(limit);
        const Counts fast = count_by_parameterization(limit);
        if (brute.with_60 != fast.with_60 || brute.with_90 != fast.with_90 ||
            brute.with_120 != fast.with_120 || brute.total() != fast.total()) {
            std::cerr << "Mismatch at limit=" << limit << ": brute(total=" << brute.total()
                      << ", 60=" << brute.with_60 << ", 90=" << brute.with_90
                      << ", 120=" << brute.with_120 << ") fast(total=" << fast.total()
                      << ", 60=" << fast.with_60 << ", 90=" << fast.with_90
                      << ", 120=" << fast.with_120 << ")" << '\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_checkpoints()) {
        return 2;
    }

    const Counts result = count_by_parameterization(options.limit);
    std::cout << result.total() << '\n';
    return 0;
}

Python

import math
import multiprocessing

def gcd(a, b):
    while b != 0:
        a, b = b, a % b
    return a

def gcd3(a, b, c):
    return gcd(a, gcd(b, c))

def count_90_range(limit, m_start, m_end, m_step):
    out = 0
    for m in range(m_start, m_end + 1, m_step):
        for n in range(1, m):
            if ((m - n) & 1) == 0 or gcd(m, n) != 1:
                continue
            primitive_perimeter = 2 * m * (m + n)
            out += limit // primitive_perimeter
    return out

def count_120_range(limit, m_start, m_end, m_step):
    out = 0
    for m in range(m_start, m_end + 1, m_step):
        for n in range(1, m):
            if gcd(m, n) != 1:
                continue
            mm = m * m
            nn = n * n
            a = mm - nn
            b = 2 * m * n + nn
            c = mm + m * n + nn
            if gcd3(a, b, c) != 1:
                continue
            primitive_perimeter = a + b + c
            out += limit // primitive_perimeter
    return out

def count_60_non_eq_range(limit, m_start, m_end, m_step):
    out = 0
    for m in range(m_start, m_end + 1, m_step):
        mm = m * m
        n_max = (m - 1) // 2
        for n in range(1, n_max + 1):
            if gcd(m, n) != 1:
                continue
            nn = n * n
            a = mm - nn
            c = mm - m * n + nn

            b_a = 2 * m * n - nn
            if gcd3(a, b_a, c) == 1:
                primitive_perimeter = a + b_a + c
                out += limit // primitive_perimeter

            b_b = mm - 2 * m * n
            if b_b > 0 and gcd3(a, b_b, c) == 1:
                primitive_perimeter = a + b_b + c
                out += limit // primitive_perimeter
    return out

def max_m_90(limit):
    m = 2
    while 2 * m * (m + 1) <= limit:
        m += 1
    return m - 1

def max_m_120(limit):
    m = 2
    while 2 * m * m + 3 * m + 1 <= limit:
        m += 1
    return m - 1

def max_m_60(limit):
    m = 2
    while True:
        mm = m * m
        min_perimeter_a = 2 * mm + m - 1
        min_perimeter_b = 3 * m * ((m + 1) // 2)
        if min_perimeter_a > limit and min_perimeter_b > limit:
            break
        m += 1
    return m - 1

def worker_main(args):
    limit, m_max_90, m_max_120, m_max_60, m_start, m_step = args
    c90 = count_90_range(limit, m_start, m_max_90, m_step)
    c120 = count_120_range(limit, m_start, m_max_120, m_step)
    c60 = count_60_non_eq_range(limit, m_start, m_max_60, m_step)
    return c90, c120, c60

def solve(limit=100000000):
    m_max_90 = max_m_90(limit)
    m_max_120 = max_m_120(limit)
    m_max_60 = max_m_60(limit)
    
    total_60 = limit // 3
    total_90 = 0
    total_120 = 0
    
    threads = max(1, multiprocessing.cpu_count())
    
    if threads <= 1:
        total_90 += count_90_range(limit, 2, m_max_90, 1)
        total_120 += count_120_range(limit, 2, m_max_120, 1)
        total_60 += count_60_non_eq_range(limit, 2, m_max_60, 1)
    else:
        args = []
        for t in range(threads):
            args.append((limit, m_max_90, m_max_120, m_max_60, 2 + t, threads))
            
        with multiprocessing.Pool(threads) as pool:
            results = pool.map(worker_main, args)
            
        for c90, c120, c60 in results:
            total_90 += c90
            total_120 += c120
            total_60 += c60
            
    return str(total_60 + total_90 + total_120)

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

Java

import java.util.*;
import java.util.concurrent.*;

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

    static long gcd3(long a, long b, long c) {
        return gcd(a, gcd(b, c));
    }

    static long count90Range(long limit, long mStart, long mEnd, long mStep) {
        long out = 0;
        for (long m = mStart; m <= mEnd; m += mStep) {
            for (long n = 1; n < m; ++n) {
                if (((m - n) & 1L) == 0 || gcd(m, n) != 1L) {
                    continue;
                }
                long primitivePerimeter = 2L * m * (m + n);
                out += limit / primitivePerimeter;
            }
        }
        return out;
    }

    static long count120Range(long limit, long mStart, long mEnd, long mStep) {
        long out = 0;
        for (long m = mStart; m <= mEnd; m += mStep) {
            for (long n = 1; n < m; ++n) {
                if (gcd(m, n) != 1L) {
                    continue;
                }
                long mm = m * m;
                long nn = n * n;
                long a = mm - nn;
                long b = 2L * m * n + nn;
                long c = mm + m * n + nn;
                if (gcd3(a, b, c) != 1L) {
                    continue;
                }
                long primitivePerimeter = a + b + c;
                out += limit / primitivePerimeter;
            }
        }
        return out;
    }

    static long count60NonEqRange(long limit, long mStart, long mEnd, long mStep) {
        long out = 0;
        for (long m = mStart; m <= mEnd; m += mStep) {
            long mm = m * m;
            long nMax = (m - 1) / 2;
            for (long n = 1; n <= nMax; ++n) {
                if (gcd(m, n) != 1L) {
                    continue;
                }
                long nn = n * n;
                long a = mm - nn;
                long c = mm - m * n + nn;

                long ba = 2L * m * n - nn;
                if (gcd3(a, ba, c) == 1L) {
                    long primitivePerimeter = a + ba + c;
                    out += limit / primitivePerimeter;
                }

                long bb = mm - 2L * m * n;
                if (bb > 0 && gcd3(a, bb, c) == 1L) {
                    long primitivePerimeter = a + bb + c;
                    out += limit / primitivePerimeter;
                }
            }
        }
        return out;
    }

    static long maxM90(long limit) {
        long m = 2;
        while (2L * m * (m + 1) <= limit) {
            ++m;
        }
        return m - 1;
    }

    static long maxM120(long limit) {
        long m = 2;
        while (2L * m * m + 3L * m + 1 <= limit) {
            ++m;
        }
        return m - 1;
    }

    static long maxM60(long limit) {
        long m = 2;
        while (true) {
            // Need BigInteger or __int128 equivalent to avoid overflow, but m is around
            // 10000 for limit=1e8
            // which fits fine in 64-bit int. Let's do it safely.
            long mm = m * m;
            long minPerimeterA = 2L * mm + m - 1;
            long minPerimeterB = 3L * m * ((m + 1) / 2);
            if (minPerimeterA > limit && minPerimeterB > limit) {
                break;
            }
            ++m;
        }
        return m - 1;
    }

    public static String solve() {
        long limit = 100_000_000L;
        long mMax90 = maxM90(limit);
        long mMax120 = maxM120(limit);
        long mMax60 = maxM60(limit);

        long with60 = limit / 3;
        long with90 = 0;
        long with120 = 0;

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());

        if (threads <= 1) {
            with90 += count90Range(limit, 2, mMax90, 1);
            with120 += count120Range(limit, 2, mMax120, 1);
            with60 += count60NonEqRange(limit, 2, mMax60, 1);
        } else {
            ExecutorService executor = Executors.newFixedThreadPool(threads);
            List<Future<long[]>> futures = new ArrayList<>();

            for (int t = 0; t < threads; ++t) {
                final long mStart = 2 + t;
                final long mStep = threads;

                futures.add(executor.submit(() -> {
                    long c90 = count90Range(limit, mStart, mMax90, mStep);
                    long c120 = count120Range(limit, mStart, mMax120, mStep);
                    long c60 = count60NonEqRange(limit, mStart, mMax60, mStep);
                    return new long[] { c90, c120, c60 };
                }));
            }

            for (Future<long[]> f : futures) {
                try {
                    long[] res = f.get();
                    with90 += res[0];
                    with120 += res[1];
                    with60 += res[2];
                } catch (Exception e) {
                    e.printStackTrace();
                }
            }
            executor.shutdown();
        }

        long total = with60 + with90 + with120;
        return String.valueOf(total);
    }

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