Numerical implementation of Archimedes' method of exhaustion

Archimedes' method of exhaustion is a famous way to approximate pi via inscribed and circumscribed polygons. From it you get the inequalities n sin(pi/n)<pi<n tan(pi/n) (n>=3) in modern notation. You can get these sines and tangents for n=3,6,12,24,... (or n=4,8,16,...) using half angle formulas. This relies on only the four basic operations and square roots.

I think this is a fun example for showing mathematically equivalent methods of implementing a procedure that are far from being equivalent in floating point. For example, using cos(t/2)=sqrt((1+cos(t))/2)), sin(t/2)=sqrt((1-cos(t))/2), tan(t/2)=sin(t/2)/cos(t/2) is a bad scheme because of how quickly catastrophic cancellation kicks in from the 1-cos. You can circumvent it by multiplying top and bottom by sqrt(1+cos(t)), obtaining sin(t/2)=sin(t)/sqrt(2+2cos(t))=sin(t)/(2cos(t/2)). As a bonus, now each step of the scheme only needs one square root!

You can get the same thing if you come at it from the tangent side: tan(t/2)=(sqrt(1+tan(t)^2 ) - 1)/tan(t) is a bad scheme , multiplying top and bottom by 1+sqrt(1+tan(t)^2 ) gives tan(t/2)=tan(t)/(1+sqrt(1+tan(t)^2 )) as an alternative scheme. And again, in this scheme you basically only need one square root per step, as sin(t/2)=tan(t/2)/sec(t/2)=tan(t/2)/sqrt(1+tan(t/2)^(2)), so you get the square root you will need in the next step in the course of getting the sine you need in this step.

In my undergrad numerical analysis class, we did basically this type of manipulation with the quadratic formula, to get a better scheme for the small root when b^2 >>|4ac|. But I think this example is a bit more fun, maybe.