I blogged about this topic before. The algorithm provided there is very good. But there is a better algorithm involving the same amount of construction effort. Moreover, it is much easier to understand. The error will be comparable for large angles but the order of the error is much smaller (9 instead of 5). This blog is also a lesson on how to determine the order of an error numerically.

We can use the geometry tools of Euler Math Toolbox (EMT) to do the construction of the old blog entry (or my program C.a.R. if you prefer a visual construction). Here is the code for this.

>load geometry Numerical and symbolic geometry. >O=[0,0]; c=circleWithCenter(O,1); >setPlotRange(0,1.2,0,1.2); >color(black); plotCircle(c,""); >phi=85°; P=[cos(phi),sin(phi)]; >A=[1,0]; s=lineThrough(A,P); >plotPoint(P,"P",[0,-1]); plotPoint(A,"A"); plotSegment(P,A,""); >B=A+(P-A)/3; l1=lineThrough(O,B); {P1,P2}:=lineCircleIntersections(l1,c); >l2=perpendicular(B,s); {P3,P4}:=lineCircleIntersections(l2,c); >plotPoint(B,"B",[-1,-1]); >color(green); plotPoint(P1,""); plotSegment(O,P1,""); >color(red); plotPoint(P4,""); plotSegment(B,P4,""); >color(black); plotPoint([cos(phi/3),sin(phi/3)],""):

The green construction is the easiest one. The red one is a bit better. We get a big improvement if we trisect the chord between the red and the green point again. So in fact, we have two trisections of segments to do, one for B and one for the final result. Let us plot the errors of the green and the red constructions up to 45°.

I did that with the following code in EMT.

>function map first (phi) ... $ global A,O,c; $ P=[cos(phi),sin(phi)]; $ B=A+(P-A)/3; $ return atan2(B[1],B[2]); $ endfunction >phi=linspace(0°,45°,450); appr1=first(phi); >function map second (phi) ... $ if phi~=0 then return 0; endif; $ global A,O,c; $ P=[cos(phi),sin(phi)]; $ B=A+(P-A)/3; $ l=perpendicular(B,lineThrough(A,P)); $ {P1,P2}=lineCircleIntersections(l,c); $ return atan2(P2[1],P2[2]); $ endfunction >appr2=second(phi); >plot2d(deg(phi),deg([appr1;appr2]-phi/3),color=[green,red]):

Now, the second trisection by chance gives us a better approximation. It was discovered or re-discovered by my colleague Hans Fischer.

>function map third (phi) ... $ a=first(phi); A=[cos(a),sin(a)]; $ b=second(phi); B=[cos(b),sin(b)]; $ C=A+(B-A)*2/3; $ return atan2(C[1],C[2]); $ endfunction >plot2d(deg(phi),deg(phi/3-third(phi)),color=blue):

As you see the error is less than 1/100 degree up to 45°. However, we can do better and much easier.

For this, we simply trisect the error of the first trisection using the first approximation (simple trisection of the chord). This is done by taking the approximate trisection times 3 and trisecting the error again. Then add this to the first approximation. The error of this is the following.

>function map fourth (phi) ... $ a=first(phi); $ return a+first(phi-3*a); $ endfunction >plot2d(deg(phi),deg(fourth(phi)-phi/3),color=cyan):

This looks about the same as the Fischer construction, but for small angles it is much better. To measure the order of the error of the first approximation, we assume

\(\alpha = \dfrac{\phi}{3} + c \phi^3 + \ldots\)

If that is the case, the following plot should look almost linear as it does.

>plot2d("abs(first(x)-x/3)^(1/3)",0,pi/2):

The constant c can be determined from the slope of this line.

>x=0.0001; ((first(x)-x/3)^(1/3)/x)^3 -0.0123456813891

The order of this approximation is 3 (as the angle tends to 0 and taking the absolute error).

In the same manner, we can determine that the order of approximation of the Fischer method is 5. That is already very good.

But obviously, our last method is even better with an order of 9. This is obvious since we first have an error with exponent 3 and then make this error smaller with the same exponent again. Here is the code to check this numerically.

>plot2d("abs(fourth(x)-x/3)^(1/9)",0,pi/2):

As you see, for small angles the error is so small that it becomes 0 numerically.

For this simple method, we could even analyze the Taylor series with an algebra system such as Maxima in EMT. I leave that to the reader.