Fairly trisect a pizza with two cuts. Warning : Sneaky hard!
Find x, the length of the blue chord cutting off 1/3 of a circle' s area.
How can this be hard? We just write the area of the segment as a function of the chord length, equate it to π/3, and solve for x:
Oops! How do you solve for x when it's both inside and outside an arcsine? You don't! Not with recognized functions, anyway. You can only approximate. But as well as you wish, which is the same as knowing the answer. In fact, that's how most computers take square root. The only difference is that square root is so important that somebody already told the computer how to approximate it out the tubes. For a direct attack on this equation, see www.tweedledum.com/rwg/oldpizza.htm, but notice how much simpler things get if we replace x by twice the sine of something. Looking at the diagram, "something" happens to be α/2, where α is the central angle subtended by x:
Cramming all the trig (along with a foresighted factor of 2) into one step, the problem simplifies to
(The triginverses and radexpand switches say to use arcsin(sin y) = y = 
.) The three main ways to solve (d3) are successive approximations, series reversion, and special functions. First, Newton's Iteration: Subtract the sides of the equation, making something you wish to vanish. Divide by the derivative, then subtract from the unknown:
This will usually be a good deal closer than α to the solution. From the diagram, the two radii delimiting α are roughly 2:30 and 9:30 o'clock. Clock numbers are 2π/12 = π/6 apart, so a good first approximation would be π - π/12 - π/12:
Substitute this for α in (d4), the Newton step:
Substitute this for α in (d4), the Newton step:
The number of correct digits usually doubles each time, so we need at most one more step:
So the "final" answer to our original question is
a little less than the diameter 2, which we can finally check by substituting into the original equation (d1):
That's π/3 to way better than one particle of pizza flour dust.
The second approach is Series Reversion. Expand the unknown side of the (d3) equation at α = 5π/6, our first approximation, and replace the 2π/3 by a symbol, say 2p3:
Now we can "solve" for α as a series in 2p3:
Restoring 2p3 = 2π/3 and approximating with double precision,
in agreement with (d8). Why did we stop with eight terms (7th order)? Because that sufficed in the rehearsal! Notice that (d12) is our initial approximation plus a series of powers of
This is pleasantly small, but its 
power portends only twelve digits of accuracy. Those rapidly growing surds must be rapidly shrinking, numerically. To see this, and to rewrite the series in more numerically hygienic terms, make the strange looking substitution
which is valid if 2p3 = 2π/3 and pi = π:
Substituting into our reverted series and expanding at pi = 3,
so the term ratio is actually five times smaller than (d14). This series and (d12) should match term for term, but watch what happens when we evaluate the last terms numerically:
The lhs, unlike the right, suffers from severe "subtractive cancellation". We got the right answer anyway because the smaller terms need less accuracy.
Third method, Special Functions: If we generalize the problem to cut off area A (instead of the special case π/3), then this problem becomes the very famous Kepler's equation, M = E - ε sin E, for a planet in an elliptical orbit! In Kepler's terms, the "eccentric anomaly" E (putting the planet at (a cos E, b sin E) relative to the center of the ellipse) = α, the "mean anomaly" M (= 2π/period times the time since perihelion) = 2A, and the eccentricity ε (= separation of foci/major axis a) of the ellipse = 1! That is, instead of a planet, we have a comet bouncing against the sun along a line segment. Kepler's problem was solved by Lagrange in 1770 in terms of an infinite series of Bessel functions:
Unfortunately, when ε = 1, the convergence of this series is horrifically slow--you need thousands of terms to get 2.605325, even though every third term = 0 when A = π/3. But the convergence is fairly smooth if we group the terms three at a time. Summing through n = 3, 6, 12, 24, 48, 96, 192, and 384,
Because of this smoothness, we can squeeze out three or four more digits by heuristically extrapolating the sequence, pretending that the nth term is of the form
Then if a, b, and c are values of this sequence for three consecutive n, as n goes to infinity the asymptotic value is
Using this formula on each consecutive group of three values in (d20), we get a new sequence shorter by two:
Repeating twice more (where % means the previous result),
The problem is, without an independent esimate we don't know how good this last value is, but it's usually slightly better than the previous, so 2.60532567 is fairly plausible. Caution: Pushing this heuristic too far leads to craziness, as demonstrated on page 81 of MIT AIM 304 at http://dspace.mit.edu/handle/1721.1/6088?show=full . (The bottom half of page 83 will save you reading the previous 80.)
But we can get useful convergence from Lagrange's formula by replacing the Bessel function by its infinite powerseries and then interchanging the order of summation. The result is messy looking,
where Im is the imaginary part, and W(t) is the Lambert W function, magically defined so that W(t) exp(W(t)) = t. The mess is worth it. For our ε =1, A = π/3 case, rotating ∞ by 90 degrees (i.e., replacing it by 8(!)) gives
(c25) block([ratprint:false],bfloat(eval(subst([epsilon=1.b0,\im=imagpart],apply_nouns(apply_nouns(subst(["-2 k"=-2*k,\a=%pi/3.b0,inf=8,\w=lambert_w],d24)))))))
This transformation might be interesting to astronomers for ε near 1. Unfortunately the derivation assumed ε < 2/e, and it is pure theology why it works for greater ε; it just seems to! Here's an efficiency, assuming we cache a particular sequence of polynomials 
The polynomials 
a Lambert's W function:
