Problem 177: Integer Angled Quadrilaterals

View on Project Euler

Project Euler Problem 177 Solution

EulerSolve provides an optimized solution for Project Euler Problem 177, Integer Angled Quadrilaterals, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Take a convex quadrilateral \(ABCD\) and let its diagonals meet at \(X\). The diagonals split the four corner angles into eight smaller angles. Problem 177 asks for the number of such quadrilaterals for which all eight of those angles are integers measured in degrees, with shapes identified up to rotation and reflection. The brute-force version of that question is very redundant. An arbitrary 8-tuple of positive integers does not automatically correspond to a real quadrilateral, and even when it does, the same geometric figure appears many times under different starting vertices or reversed orientations. The implementations solve the problem by finding the right five-parameter description, deriving one trigonometric closure condition, and then canonicalizing the surviving angle data under the dihedral symmetries of the quadrilateral. Mathematical Approach The useful viewpoint is to work with the diagonals and the four small triangles that they create. Once the angle data is organized correctly, the geometry separates into a linear part and a multiplicative sine-law part....

Detailed mathematical approach

Problem Summary

Let \(ABCD\) be a convex quadrilateral with the diagonals \(AC\) and \(BD\). At each vertex the diagonal through that vertex makes an angle with each of the two sides there, which gives eight corner angles; at \(A\) they are \(\angle CAD\) and \(\angle CAB\). The quadrilateral is integer angled if all eight corner angles are whole numbers of degrees. A square is integer angled, with eight angles of \(45^\circ\), and so is the quadrilateral with \(\angle DAC=20^\circ\), \(\angle BAC=60^\circ\), \(\angle ABD=50^\circ\), \(\angle CBD=30^\circ\), \(\angle BCA=40^\circ\), \(\angle DCA=30^\circ\), \(\angle CDB=80^\circ\) and \(\angle ADB=50^\circ\). We must count the integer angled quadrilaterals, with similar ones counted only once. The statement allows a computed angle to be treated as an integer if it lies within \(10^{-9}\) of one. [1]

The whole solution in one sentence: the four triangles cut out by the diagonals leave five free angles and one equation between two products of sines; sorted tables of one side of that equation and binary search find all its solutions, and each shape is counted once through the smallest of its eight labellings.

Mathematical Approach

All angles are measured in degrees, and inside formulas the degree sign is dropped.

1. Four triangles around the crossing point

Since \(ABCD\) is convex, its diagonals cross at a point \(X\) inside both of them. We name the corner angles in their order around the quadrilateral:

  • at \(A\): \(a=\angle DAC\) and \(p=\angle CAB\);
  • at \(B\): \(q=\angle ABD\) and \(r=\angle DBC\);
  • at \(C\): \(s=\angle BCA\) and \(c=\angle ACD\);
  • at \(D\): \(t=\angle CDB\) and \(u=\angle BDA\).

The triangle \(ABX\) has the angles \(p\) at \(A\) and \(q\) at \(B\), the triangle \(BCX\) has \(r\) and \(s\), \(CDX\) has \(c\) and \(t\), and \(DAX\) has \(u\) and \(a\). So the four angles at \(X\) are \(180-p-q\), \(180-r-s\), \(180-c-t\) and \(180-u-a\). Opposite angles at \(X\) are equal, and two neighbouring ones add up to \(180\), because \(A\), \(X\) and \(C\) lie on one line. With \(U=p+q\) and \(V=180-U\) this gives the linear conditions

$$p+q=c+t=U,\qquad r+s=u+a=V.$$

In the example of the statement \(U=60+50=110\) and \(V=30+40=70\); indeed \(c+t=30+80\) and \(u+a=50+20\), and the diagonals cross at angles of \(70\) and \(110\).

2. Five free angles

By Section 1 the angles \(p,q,r,s\) are positive integers with \(p+q+r+s=180\). Once they are chosen, \(c\) can be any of \(1,\dots,U-1\) and \(a\) any of \(1,\dots,V-1\), which fixes \(t=U-c\) and \(u=V-a\). Choosing \(p,q,r,s\) amounts to cutting a row of \(180\) unit steps at three of its \(179\) inner points, so there are \(\binom{179}{3}=939\,929\) choices. For a given \(U\) there are \(U-1\) choices each for \(p\) and \(c\) and \(179-U\) each for \(r\) and \(a\), so altogether

$$\sum_{U=2}^{178}(U-1)^2(179-U)^2=5\,956\,330\,073$$

lists of eight positive integers satisfy the linear conditions. Most of them belong to no quadrilateral: the four triangles must also fit together.

3. The sine condition

The law of sines says that the sides of a triangle are proportional to the sines of the opposite angles. [2] In \(ABX\) the side \(XA\) lies opposite the angle \(q\) and the side \(XB\) opposite \(p\), so the four triangles give

$$\begin{aligned}\frac{XA}{XB}&=\frac{\sin q}{\sin p},&\frac{XB}{XC}&=\frac{\sin s}{\sin r},\\ \frac{XC}{XD}&=\frac{\sin t}{\sin c},&\frac{XD}{XA}&=\frac{\sin a}{\sin u},\end{aligned}$$

and the product of the four left-hand sides is \(1\). Hence every integer angled quadrilateral satisfies the sine condition

$$\sin q\,\sin s\,\sin t\,\sin a=\sin p\,\sin r\,\sin c\,\sin u.$$

4. The conditions are also sufficient

Claim. Eight positive numbers that satisfy the linear conditions and the sine condition are the corner angles of a convex quadrilateral, and they determine it up to similarity.

The proof uses one fact. In a triangle \(XYZ\) with the angles \(\varphi\) at \(Y\) and \(W-\varphi\) at \(Z\), where \(0\lt\varphi\lt W\lt180\), the law of sines gives

$$\frac{XY}{XZ}=\frac{\sin(W-\varphi)}{\sin\varphi}=\sin W\cot\varphi-\cos W.$$

As \(\cot\varphi\) decreases strictly and \(\sin W\gt0\), the ratio \(XY/XZ\) determines \(\varphi\).

Proof of the claim. Draw two lines that cross at a point \(X\), and put \(A\) and \(B\) on two neighbouring rays with \(\angle AXB=180-U\) and \(XA/XB=\sin q/\sin p\). The angles of \(ABX\) at \(A\) and \(B\) add up to \(U\), so by the fact they are \(p\) and \(q\). Next put \(C\) on the ray opposite to \(XA\) with \(XB/XC=\sin s/\sin r\), and \(D\) on the ray opposite to \(XB\) with \(XC/XD=\sin t/\sin c\); in the same way \(BCX\) gets the angles \(r\) and \(s\), and \(CDX\) gets \(c\) and \(t\). Multiplying the three ratios and using the sine condition gives

$$\frac{XA}{XD}=\frac{\sin q\,\sin s\,\sin t}{\sin p\,\sin r\,\sin c}=\frac{\sin u}{\sin a}.$$

The angles of \(DAX\) at \(D\) and \(A\) add up to \(V=u+a\), so the fact makes them \(u\) and \(a\). The four triangles fill the quadrilateral \(ABCD\); its angles \(a+p\), \(q+r\), \(s+c\) and \(t+u\) are all smaller than \(U+V=180\), so it is convex, and its corner angles are the eight given numbers. Conversely, in any quadrilateral with these corner angles the four triangles around \(X\) have these angles, which fixes the angle between the diagonals and the ratios \(XA:XB:XC:XD\), and hence the quadrilateral up to similarity.

5. The example of the statement

For the example, \((a,p,q,r,s,c,t,u)\)\({}=(20,60,50,30,40,30,80,50)\), and the sine condition reads \(\sin50\sin40\sin80\sin20\)\({}=\sin60\sin30\sin30\sin50\). The product-to-sum rule \(\sin x\sin y\)\({}=\tfrac12\bigl(\cos(x-y)-\cos(x+y)\bigr)\) gives \(\sin20\sin40=\tfrac12(\cos20-\tfrac12)\), and the similar rule \(\cos20\sin80=\tfrac12(\sin100+\sin60)\) together with \(\sin100=\sin80\) then gives

$$\sin20\,\sin40\,\sin80=\tfrac14\sin60=\tfrac{\sqrt3}{8},$$

the sine companion of Morrie's law \(\cos20\cos40\cos80=\tfrac18\). [5] As \(\sin60\sin30\sin30=\tfrac{\sqrt3}{8}\) too, both sides equal \(\tfrac{\sqrt3}{8}\sin50\approx0.1658535\).

The angles at \(A\) and \(B\) are \(a+p=80\) and \(q+r=80\), so the lines \(AD\) and \(BC\) meet at a point \(O\) with \(\angle AOB=20\), and the diagonals make the angles \(20\) and \(30\) with the equal sides \(OA\) and \(OB\). This is the figure of Langley's adventitious angles, a puzzle posed in 1922 that asks for a proof of \(\angle ACD=30\). [4] Quadrilaterals whose angles between sides and diagonals are all rational in degrees are called adventitious quadrangles; Gerrit Bol classified them in 1936, through the equivalent problem of triple intersections of diagonals in regular polygons, and Bjorn Poonen and Michael Rubinstein checked his work by computer. [4] [6]

6. Cyclic quadrilaterals always satisfy the sine condition

A convex quadrilateral is cyclic, that is, its vertices lie on a circle, exactly when an angle between a side and a diagonal equals the angle between the opposite side and the other diagonal. [3] In a cyclic quadrilateral \(ABCD\), therefore, \(q=\angle ABD=\angle ACD=c\) and \(r=\angle DBC=\angle DAC=a\). Section 1 then gives \(t=U-q=p\) and \(u=V-r=s\), and both sides of the sine condition become \(\sin p\sin q\sin r\sin s\). Conversely, for every choice of \(p,q,r,s\) the list with \(c=q\) and \(a=r\) satisfies both conditions, so by Section 4 it belongs to a quadrilateral, which is cyclic because \(q=c\). Each of the \(939\,929\) choices of Section 2 therefore gives exactly one labelled cyclic integer angled quadrilateral.

7. Shapes, not labels

Two figures are similar if one has the same shape as the other or as its mirror image. [7] A quadrilateral can be labelled \(ABCD\) in eight ways: any vertex can be called \(A\), and the labels can run around in either direction. They are the eight symmetries of a square, the dihedral group of order \(8\), acting on its corners. [8] Starting one vertex later moves the list \((a,p,q,r,s,c,t,u)\) two places, to \((q,r,s,c,t,u,a,p)\); running the labels the other way, \(A\), \(D\), \(C\), \(B\), turns it into \((p,a,u,t,c,s,r,q)\). By Section 4, two lists belong to similar quadrilaterals exactly when a relabelling turns one into the other. So a shape is a class of lists, and we represent it by its lexicographically smallest list, the one that is smaller at the first place where two lists differ. The square has the single list \((45,\dots,45)\); the example of the statement has eight different lists, and its own list is the smallest.

8. Splitting the condition in two

Testing all lists of Section 2 is possible, but slow. Dividing both sides of the sine condition by \(\sin q\sin s\sin c\sin u\) and using \(t=U-c\) and \(u=V-a\) turns it into \(K(p,q,r)=R_U(a,c)\), where

$$K(p,q,r)=\frac{\sin p\,\sin r}{\sin q\,\sin s},$$

$$R_U(a,c)=\frac{\sin(U-c)\,\sin a}{\sin c\,\sin(V-a)}.$$

The left side depends only on \(p,q,r\), which fix \(U\), \(V\) and \(s=V-r\), and the right side only on \(U\) and the pair \((a,c)\). So for each \(U\) we compute the \((U-1)(V-1)\) values \(R_U(a,c)\) once and sort them. Then for each triple \((p,q,r)\) with \(p+q=U\) a binary search, which halves the range of possible positions at every step, finds where \(K(p,q,r)\) belongs in the sorted table. [9] There are \(939\,929\) triples and, over all \(U\), \(939\,929\) table entries, the same number \(\binom{179}{3}\) as in Section 2; the largest table, for \(U=90\), has \(89^2=7921\) entries. Instead of about \(6\cdot10^9\) tests we need one lookup per triple, made of two binary searches of at most \(13\) steps.

The values are bounded: by the fact of Section 4, \(\sin p/\sin q\) grows with \(p\) for fixed \(U\), and \(\sin r/\sin s\) with \(r\), so \(K\) is largest for \(p=U-1\), \(r=V-1\). As \(\sin(V-1)=\sin(U+1)\) and \(\sin(U-1)\sin(U+1)=\sin^2U-\sin^21\) by the product-to-sum rule,

$$\begin{aligned}\max K&=\frac{\sin(U-1)\,\sin(V-1)}{\sin^21}\\&=\frac{\sin^2U-\sin^21}{\sin^21}\le\cot^21\approx3282,\end{aligned}$$

and the smallest value is the reciprocal, at least \(\tan^21\approx3.05\cdot10^{-4}\). The values \(R_U(a,c)\) have the same range.

9. One search worked out

For the example, \(p=60\), \(q=50\) and \(r=30\) give \(U=110\), \(V=70\), \(s=40\) and

$$K=\frac{\sin60\,\sin30}{\sin50\,\sin40}\approx\frac{0.4330127}{0.4924039}\approx0.8793852.$$

The table for \(U=110\) has \(109\cdot69=7521\) entries, and exactly four of them are equal to this \(K\): \((a,c)=(20,30)\), the example itself; \((30,50)\), the cyclic quadrilateral with the same \(p,q,r,s\); \((60,100)\), where \(t=u=10\) and the sine condition reduces to \(\sin50\sin40=\tfrac12\sin100\), the double-angle formula, because \(\sin40=\cos50\); and \((40,70)\), which follows from the identity of Section 5 and the double-angle formula. The nearest other entries, \(0.8779904\) and \(0.8807565\), are more than \(10^{-3}\) away.

10. A complete small case

The whole method can be followed by hand if every angle must be a multiple of \(30\). Then \(p+q+r+s=180\) has \(\binom53=10\) solutions, and Section 2 gives \(9+16+9=34\) lists. Only \(\sin30=\tfrac12\), \(\sin60=\tfrac{\sqrt3}2\) and \(\sin90=1\) occur. For \(U=60\) the three triples give \(K=\tfrac12,1,2\) and the three pairs give \(R_U=\tfrac12,1,2\), so \(3\) lists pass; \(U=120\) behaves in the same way. For \(U=90\) both \(K\) and \(R_U\) take the values \(\tfrac13,1,1,3\), so \(1+2\cdot2+1=6\) lists pass. These \(12\) lists are the \(10\) cyclic lists of Section 6, \((30,30,60,60,30,30,60,60)\) and \((60,60,30,30,60,60,30,30)\), and under the relabellings they form four shapes:

  • the rectangle whose diagonals cross at \(60\), \((30,60,60,30,30,60,60,30)\), with \(2\) lists;
  • the rhombus with the angles \(60\) and \(120\), \((30,30,60,60,30,30,60,60)\), with \(2\) lists, the only shape that is not cyclic;
  • the kite with the angles \(60\), \(90\), \(120\), \(90\), \((30,30,60,30,60,60,30,60)\), with \(4\) lists;
  • half of a regular hexagon, \((30,30,30,30,90,30,30,90)\), with \(4\) lists.

The same exhaustive search with an exact sine test (Section 11) finds \(1\) shape when the angles are multiples of \(45\) (the square), \(4\) for \(30\), \(10\) for \(20\), \(24\) for \(18\), \(51\) for \(15\), \(56\) for \(12\), \(207\) for \(10\), \(209\) for \(9\), \(1040\) for \(6\) and \(1321\) for \(5\). For the units \(20\) and \(12\) every shape is cyclic, while for the unit \(10\) of the example \(103\) of the \(207\) shapes are not.

11. Rounding errors and exact equality

Rounding can make equal values differ slightly and, in principle, different values look equal. A double has a 53-bit significand, so each arithmetic operation is accurate to a relative \(2^{-53}\approx1.1\cdot10^{-16}\), [10] and the computed sines carry rounding errors of their own. The code accepts an entry when \(|R_U(a,c)-K|\le10^{-10}\) and the two products of the sine condition, computed again, differ by at most \(10^{-12}\). For the true solutions the ratios never differ by more than \(9.1\cdot10^{-13}\) and the products by more than \(1.7\cdot10^{-16}\), while every list that is not a solution has \(|R_U(a,c)-K|\gt1.16\cdot10^{-8}\). The closest call in relative terms is the list \((82,23,46,69,42,52,17,29)\), where \(\sin46\sin42\sin17\sin82\)\({}\approx0.1393583837\) and \(\sin23\sin69\sin52\sin29\)\({}\approx0.1393583907\) agree to seven digits. For every list that is not a solution, building three of its four triangles as in Section 4 leaves the fourth triangle with angles that differ from the prescribed ones by at least \(6.6\cdot10^{-7}\) degrees, far more than the tolerance \(10^{-9}\) of the statement.

Equality can also be decided without rounding. Put \(\zeta=e^{i\pi/180}\), so that \(\sin x=(\zeta^x-\zeta^{-x})/(2i)\) for an angle of \(x\) degrees. Multiplying out four such factors gives

$$16\sin x_1\sin x_2\sin x_3\sin x_4=\sum\pm\zeta^{\pm x_1\pm x_2\pm x_3\pm x_4},$$

a sum of \(16\) terms whose signs are the products of the four signs in their exponents; Poonen and Rubinstein treat products of three sines in the same way. [6] So \(16\) times the difference of the two sides of the sine condition is \(P(\zeta)\) for a polynomial \(P\) made of \(32\) terms \(\pm x^k\), where \(k\) can be taken modulo \(360\) because \(\zeta^{360}=1\). The minimal polynomial of the primitive \(360\)th root of unity \(\zeta\) over the rationals is the cyclotomic polynomial [11]

$$\begin{aligned}\Phi_{360}(x)&=\Phi_{30}(x^{12})\\&=x^{96}+x^{84}-x^{60}-x^{48}-x^{36}+x^{12}+1.\end{aligned}$$

The remainder of \(P\) after division by \(\Phi_{360}\) vanishes at \(\zeta\) exactly when \(P\) does, and since its degree is below \(96\), the degree of the minimal polynomial, only if it is zero. So the sine condition holds exactly when this integer remainder is zero.

How the Code Works

The C++ program [12] does all its work in solve. It stores \(\sin d\) for \(d=0,\dots,180\) in sinv and the quotients g[n][k] \({}=\sin(n-k)/\sin k\). For every \(U\) from \(2\) to \(178\) it fills ratio_list[U] with the \((U-1)(V-1)\) entries value, a, c, where value is g[U][c] / g[V][a], that is \(R_U(a,c)\) of Section 8, and sorts them by value.

The main loop runs over \(p=1,\dots,178\), \(q=1,\dots,178-p\) and \(r=1,\dots,V-1\), so that all eight angles are positive (Section 2), sets \(s=V-r\) and computes K. lower_bound and upper_bound find the entries with values in \([K-10^{-10},K+10^{-10}]\). For each of them the program recomputes left \({}=\sin q\sin s\sin(U-c)\sin a\) and right \({}=\sin p\sin r\sin c\sin(V-a)\), the two sides of the sine condition, and skips the entry if they differ by more than \(10^{-12}\) times max(1, |left|); a product of sines never exceeds \(1\), so this is an absolute tolerance. An accepted entry becomes the list seq \({}=(a,p,q,r,s,c,U-c,V-a)\), and canonical replaces it by its smallest relabelling (Section 7): it splits the list into four vertex pairs and compares the four rotations of the pairs and the four rotations of the reversed pairs \((p,a),(u,t),(c,s),(r,q)\) lexicographically, as std::array values. The results go into an unordered_set, whose size is printed. Every list of Section 2 belongs to exactly one triple \((p,q,r)\) and one table entry, so each labelled solution is found once and each shape is counted once.

The C++ program uses long double; on the ARM64 Mac used for the checks below it has the same 53-bit significand as double, the type of the other two versions. The Python version sorts tuples \((R_U(a,c),a,c)\), searches with bisect_left and bisect_right and keeps a set of tuples; it rebuilds the list vals of sorted values for every pair \((p,q)\) instead of once per \(U\), \(83\,653\,681\) element copies in all, about one of its four seconds. The Java version stores each entry as an array of three double values and has its own lowerBound and upperBound. Its canonical encodes each relabelling as an eight-digit number in base \(181\), below \(181^8\approx1.15\cdot10^{18}\lt2^{63}\), which fits in a long and, as all digits lie between \(1\) and \(179\), is ordered like the lists; the smallest goes into a HashSet<Long>.

Complexity and Verification

Building the tables costs one division per entry, and sorting a table of \(M\) entries about \(M\log_2M\) comparisons, about \(1.2\cdot10^7\) in all. The main loop makes two binary searches of at most \(13\) steps for each of the \(939\,929\) triples, as \(\log_27921\approx12.95\), plus a recheck of each entry in the window and, for each accepted one, a canonical form and a set insertion. The tables dominate the memory: the C++ entries take \(16\) bytes each on the test machine, about \(15\) MB. With an angle unit of \(180/n\) degrees there are about \(n^3/6\) entries and triples, so the method needs \(O(n^3\log n)\) time plus the work on the matches, and \(O(n^3)\) memory plus the set of shapes; the copies of vals add about \(n^4/12\) steps in Python, and testing every list of Section 2 would take \(\Theta(n^5)\) time. On a desktop computer the C++, Java and Python programs run in \(0.15\), \(0.4\) and \(4.2\) seconds, with peak memory use of \(27\), \(220\) (including the Java virtual machine) and \(136\) MB.

The programs contain no self-checks; they only print the count. For this article we checked the following, outside the published solutions.

  • The three programs print the same number.
  • Every pair of a triple and a table entry with \(|R_U(a,c)/K-1|\le10^{-4}\), \(155\,168\) of them non-solutions, was decided with the exact test of Section 11; the acceptance rule of the programs, evaluated in double precision, agrees with it on all of them, and the exact solutions give the same number of shapes. Outside this window \(|R_U(a,c)-K|\)\({}\gt10^{-4}\tan^21\approx3\cdot10^{-8}\) (Section 8), so the programs rightly reject those pairs, and reconstructed angles move by more than \(2\cdot10^{-5}\) degrees. The margins of Section 11 come from these pairs.
  • For the units of Section 10, an exhaustive search over all lists with the exact test (\(1\,514\,513\) lists for \(5\) degrees) gives the counts quoted there, and so do the solutions of the full problem whose angles are multiples of the unit.
  • All \(939\,929\) cyclic lists of Section 6, the square and the example are among the solutions, and the table for \(U=110\) has exactly the four matches of Section 9.

Footnotes and References

The references below supply the background facts used above. The five free angles, the proof that the sine condition suffices, the small cases and the precision margins are derived or computed in this article.

  1. Project Euler 177 — Integer Angled Quadrilaterals. The official statement defines the eight corner angles of a convex quadrilateral ABCD, gives the square and a second example, allows a tolerance of 10⁻⁹ for computed angles and asks for the number of non-similar integer angled quadrilaterals.
  2. Wikipedia — Law of sines. For a triangle with the sides a, b, c opposite the angles α, β, γ, a/sin α = b/sin β = c/sin γ = 2R, with R the radius of the circumcircle.
  3. Wikipedia — Cyclic quadrilateral. A convex quadrilateral ABCD is cyclic if and only if an angle between a side and a diagonal equals the angle between the opposite side and the other diagonal.
  4. Wikipedia — Langley's Adventitious Angles. Edward Mann Langley posed the puzzle in The Mathematical Gazette in 1922: in a triangle ABC with angles of 80° at B and C, lines from C at 30° to CA and from B at 20° to BA meet AB at F and AC at E, and ∠BEF = 30°. Quadrilaterals such as BCEF, with all angles between diagonals and sides rational, are adventitious quadrangles; they form infinite families and sporadic examples, and classifying them is equivalent to classifying triple intersections of diagonals in regular polygons, done by Gerrit Bol in 1936 and confirmed by Poonen and Rubinstein in 1998.
  5. Wikipedia — Morrie's law. Besides cos 20° · cos 40° · cos 80° = 1/8, the article gives the sine identity sin 20° · sin 40° · sin 80° = √3/8.
  6. B. Poonen and M. Rubinstein — The Number of Intersection Points Made by the Diagonals of a Regular Polygon (arXiv:math/9508209). SIAM J. Discrete Math. 11 (1998). Section 2 writes the concurrency of three chords of a circle as an equality of two products of three sines and expands the sines into exponentials, giving a vanishing sum of roots of unity; the introduction notes that this also enumerates the adventitious quadrilaterals.
  7. Wikipedia — Similarity (geometry). Two objects are similar if they have the same shape, or if one has the same shape as the mirror image of the other.
  8. Wikipedia — Dihedral group of order 8. It is the symmetry group of a square; every symmetry permutes the four corners.
  9. cppreference — std::lower_bound. In a sorted range it returns the first element not ordered before the given value, with at most log₂ N + O(1) comparisons.
  10. Wikipedia — Double-precision floating-point format. The significand has 53 bits, and the maximum relative rounding error, the machine epsilon, is 2⁻⁵³ ≈ 1.11 × 10⁻¹⁶.
  11. Wikipedia — Cyclotomic polynomial. Φₙ is the minimal polynomial over the rationals of the primitive n-th roots of unity, and Φₙ(x) = Φ_q(x^(n/q)) when q is the product of the distinct primes dividing n; Φ₃₀(x) = x⁸ + x⁷ − x⁵ − x⁴ − x³ + x + 1.
  12. C++ — Euler177.cpp. The linked immutable C++ revision contains solve and canonical; it is the implementation record for this article.

Problem 177 source code

C++

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <unordered_set>
#include <vector>

using namespace std;

namespace {

struct Entry {
    long double value;
    uint16_t a;
    uint16_t c;
};

struct ArrayHash {
    size_t operator()(const array<int, 8>& v) const noexcept {
        size_t h = 0;
        for (int x : v) {
            h = h * 181u + static_cast<unsigned>(x);
        }
        return h;
    }
};

array<int, 8> canonical(const array<int, 8>& seq) {
    array<pair<int, int>, 4> pairs = {
        make_pair(seq[0], seq[1]),
        make_pair(seq[2], seq[3]),
        make_pair(seq[4], seq[5]),
        make_pair(seq[6], seq[7]),
    };

    array<int, 8> best{};
    bool first = true;

    auto consider = [&](const array<pair<int, int>, 4>& ps, int rot) {
        array<int, 8> s;
        for (int i = 0; i < 4; ++i) {
            const auto& p = ps[(rot + i) & 3];
            s[2 * i] = p.first;
            s[2 * i + 1] = p.second;
        }
        if (first || s < best) {
            best = s;
            first = false;
        }
    };

    for (int r = 0; r < 4; ++r) consider(pairs, r);

    array<pair<int, int>, 4> ref = {
        make_pair(pairs[0].second, pairs[0].first),
        make_pair(pairs[3].second, pairs[3].first),
        make_pair(pairs[2].second, pairs[2].first),
        make_pair(pairs[1].second, pairs[1].first),
    };
    for (int r = 0; r < 4; ++r) consider(ref, r);

    return best;
}

uint64_t solve() {
    const long double pi = acosl(-1.0L);
    array<long double, 181> sinv{};
    for (int d = 0; d <= 180; ++d) {
        sinv[d] = sinl(pi * static_cast<long double>(d) / 180.0L);
    }

    vector<vector<long double>> g(181);
    for (int n = 2; n <= 178; ++n) {
        g[n].resize(n);
        for (int k = 1; k <= n - 1; ++k) {
            g[n][k] = sinv[n - k] / sinv[k];
        }
    }

    vector<vector<Entry>> ratio_list(181);
    for (int U = 2; U <= 178; ++U) {
        int V = 180 - U;
        auto& vec = ratio_list[U];
        vec.reserve(static_cast<size_t>(U - 1) * static_cast<size_t>(V - 1));
        for (int c = 1; c <= U - 1; ++c) {
            long double gU = g[U][c];
            for (int a = 1; a <= V - 1; ++a) {
                long double ratio = gU / g[V][a];
                vec.push_back({ratio, static_cast<uint16_t>(a), static_cast<uint16_t>(c)});
            }
        }
        sort(vec.begin(), vec.end(), [](const Entry& lhs, const Entry& rhs) {
            return lhs.value < rhs.value;
        });
    }

    unordered_set<array<int, 8>, ArrayHash> unique;
    unique.reserve(200000);

    const long double search_eps = 1e-10L;
    const long double verify_eps = 1e-12L;

    for (int p = 1; p <= 178; ++p) {
        for (int q = 1; q <= 178 - p; ++q) {
            int U = p + q;
            int V = 180 - U;
            auto& vec = ratio_list[U];
            long double sinp = sinv[p];
            long double sinq = sinv[q];
            for (int r = 1; r <= V - 1; ++r) {
                int s = V - r;
                long double K = (sinp * sinv[r]) / (sinq * sinv[s]);
                long double low = K - search_eps;
                long double high = K + search_eps;

                auto it1 = lower_bound(vec.begin(), vec.end(), low,
                                       [](const Entry& e, long double v) {
                                           return e.value < v;
                                       });
                auto it2 = upper_bound(vec.begin(), vec.end(), high,
                                       [](long double v, const Entry& e) {
                                           return v < e.value;
                                       });
                for (auto it = it1; it != it2; ++it) {
                    int a = it->a;
                    int c = it->c;
                    long double left = sinq * sinv[s] * sinv[U - c] * sinv[a];
                    long double right = sinp * sinv[r] * sinv[c] * sinv[V - a];
                    long double diff = fabsl(left - right);
                    long double scale = max(1.0L, fabsl(left));
                    if (diff > verify_eps * scale) continue;

                    array<int, 8> seq = {a, p, q, r, s, c, U - c, V - a};
                    unique.insert(canonical(seq));
                }
            }
        }
    }

    return unique.size();
}

}  // namespace

int main() {
    cout << solve() << "\n";
    return 0;
}

Python

import math
import bisect

def solve():
    pi = math.pi
    sinv = [math.sin(pi * d / 180.0) for d in range(181)]

    g = [[0.0] * 181 for _ in range(181)]
    for n in range(2, 179):
        for k in range(1, n):
            g[n][k] = sinv[n - k] / sinv[k]

    ratio_list = [[] for _ in range(181)]
    for U in range(2, 179):
        V = 180 - U
        vec = []
        for c in range(1, U):
            gU = g[U][c]
            for a in range(1, V):
                ratio = gU / g[V][a]
                vec.append((ratio, a, c))
        vec.sort()
        ratio_list[U] = vec

    unique = set()
    search_eps = 1e-10
    verify_eps = 1e-12

    def canonical(seq):
        pairs = [(seq[0], seq[1]), (seq[2], seq[3]), (seq[4], seq[5]), (seq[6], seq[7])]
        best = None
        def consider(ps, rot):
            nonlocal best
            s = []
            for i in range(4):
                p = ps[(rot + i) & 3]
                s.extend(p)
            s = tuple(s)
            if best is None or s < best:
                best = s
        for r in range(4):
            consider(pairs, r)
        ref = [(pairs[0][1], pairs[0][0]), (pairs[3][1], pairs[3][0]),
               (pairs[2][1], pairs[2][0]), (pairs[1][1], pairs[1][0])]
        for r in range(4):
            consider(ref, r)
        return best

    for p_val in range(1, 179):
        for q_val in range(1, 179 - p_val):
            U = p_val + q_val
            V = 180 - U
            vec = ratio_list[U]
            if not vec:
                continue
            vals = [v[0] for v in vec]
            sinp = sinv[p_val]
            sinq = sinv[q_val]
            for r in range(1, V):
                s = V - r
                K = (sinp * sinv[r]) / (sinq * sinv[s])
                lo = K - search_eps
                hi = K + search_eps
                i1 = bisect.bisect_left(vals, lo)
                i2 = bisect.bisect_right(vals, hi)
                for idx in range(i1, i2):
                    ratio, a, c_val = vec[idx]
                    left = sinq * sinv[s] * sinv[U - c_val] * sinv[a]
                    right = sinp * sinv[r] * sinv[c_val] * sinv[V - a]
                    diff = abs(left - right)
                    scale = max(1.0, abs(left))
                    if diff > verify_eps * scale:
                        continue
                    seq = (a, p_val, q_val, r, s, c_val, U - c_val, V - a)
                    unique.add(canonical(seq))

    return str(len(unique))

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

Java

import java.util.*;

public class Euler177 {
    public static String solve() {
        double pi = Math.PI;
        double[] sinv = new double[181];
        for (int d = 0; d <= 180; d++)
            sinv[d] = Math.sin(pi * d / 180.0);

        double[][] g = new double[181][181];
        for (int n = 2; n < 179; n++)
            for (int k = 1; k < n; k++)
                g[n][k] = sinv[n - k] / sinv[k];

        double[][][] ratioList = new double[181][][];
        for (int U = 2; U < 179; U++) {
            int V = 180 - U;
            List<double[]> vec = new ArrayList<>();
            for (int c = 1; c < U; c++) {
                double gU = g[U][c];
                for (int a = 1; a < V; a++)
                    vec.add(new double[] { gU / g[V][a], a, c });
            }
            vec.sort(Comparator.comparingDouble(x -> x[0]));
            ratioList[U] = vec.toArray(new double[0][]);
        }

        Set<Long> unique = new HashSet<>();
        double searchEps = 1e-10, verifyEps = 1e-12;

        for (int p = 1; p < 179; p++) {
            for (int q = 1; q < 179 - p; q++) {
                int U = p + q, V = 180 - U;
                double[][] vec = ratioList[U];
                if (vec == null || vec.length == 0)
                    continue;
                for (int r = 1; r < V; r++) {
                    int s = V - r;
                    double K = (sinv[p] * sinv[r]) / (sinv[q] * sinv[s]);
                    int i1 = lowerBound(vec, K - searchEps);
                    int i2 = upperBound(vec, K + searchEps);
                    for (int idx = i1; idx < i2; idx++) {
                        int a = (int) vec[idx][1], cv = (int) vec[idx][2];
                        double left = sinv[q] * sinv[s] * sinv[U - cv] * sinv[a];
                        double right = sinv[p] * sinv[r] * sinv[cv] * sinv[V - a];
                        double diff = Math.abs(left - right);
                        if (diff > verifyEps * Math.max(1.0, Math.abs(left)))
                            continue;
                        unique.add(canonical(a, p, q, r, s, cv, U - cv, V - a));
                    }
                }
            }
        }
        return String.valueOf(unique.size());
    }

    static int lowerBound(double[][] arr, double val) {
        int lo = 0, hi = arr.length;
        while (lo < hi) {
            int m = (lo + hi) / 2;
            if (arr[m][0] < val)
                lo = m + 1;
            else
                hi = m;
        }
        return lo;
    }

    static int upperBound(double[][] arr, double val) {
        int lo = 0, hi = arr.length;
        while (lo < hi) {
            int m = (lo + hi) / 2;
            if (arr[m][0] <= val)
                lo = m + 1;
            else
                hi = m;
        }
        return lo;
    }

    static long canonical(int... seq) {
        int[][] pairs = { { seq[0], seq[1] }, { seq[2], seq[3] }, { seq[4], seq[5] }, { seq[6], seq[7] } };
        long best = Long.MAX_VALUE;
        for (int r = 0; r < 4; r++) {
            long v = 0;
            for (int i = 0; i < 4; i++) {
                int[] pp = pairs[(r + i) & 3];
                v = v * 181 * 181 + pp[0] * 181 + pp[1];
            }
            if (v < best)
                best = v;
        }
        int[][] ref = { { pairs[0][1], pairs[0][0] }, { pairs[3][1], pairs[3][0] }, { pairs[2][1], pairs[2][0] },
                { pairs[1][1], pairs[1][0] } };
        for (int r = 0; r < 4; r++) {
            long v = 0;
            for (int i = 0; i < 4; i++) {
                int[] pp = ref[(r + i) & 3];
                v = v * 181 * 181 + pp[0] * 181 + pp[1];
            }
            if (v < best)
                best = v;
        }
        return best;
    }

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