another way to estimate π

The Atlantic City Method for computing Pi, from Marginalia at mathpages.com, via today's installment of Anarchaia.

The method is as follows:

Here's the slowest, most inefficient method known for computing p: Choose an arbitrary initial value x (real) and iterate the mapping x -> (2-x)/(3-4x) for a total of M times. Count the number of iterates that fall within a small interval around zero. For example, let N denote the number of iterates that fall between -u/2 and +u/2 for some small u. Then the limiting value of the ratio uM/N is [π] (for sufficiently large M as u goes to zero).

I'm reminded of Buffon's needle, another experimental method for computing π. If we have a floor with parallel lines drawn on it at regular intervals, and we toss a needle on the floor at random, the average number of lines on the floor which are crossed by the needle is related to π. In particular, it's 2l/(tπ) where l is the length of the needle and t is the spacing of the lines. (For those who have seen this problem before, there's often a restriction that the needle has to be shorter than the distance between the lines; then the average number of lines crossed is just the probability that a line is crossed on any given toss of the needle.) It's possible to solve this for π. It's also possible to cheat and say you did the experiment and got a very good approximation to π, by taking into account the well-known approximation 355/113 for π and setting up the probability appropriately. This is described at the Wikipedia article. I'm not sure why anyone would bother to do this.

Another way to approximate π is to pick points uniformly at random in a square of side 2, in which is inscribed a circle of radius 1; the proportion of the points which are in the circle is π/4. I doubt that anybody can pick points uniformly at random from a square with sufficient accuracy for this to work. A computer isn't bad at picking points at random, though, and this leads to the following method of approximating π, which I've heard of before -- pick X and Y uniformly at random, and independently, from [-1, 1]. (This corresponds to picking a point in the square.) If X²+Y² ≤ 1 (and thus the point is in the circle) increment a counter. Doing twenty such trials ten thousand times I get the following numbers of hits:

7788, 7792, 7806, 7806, 7814, 7825, 7833, 7834, 7841, 7842, 7846, 7850, 7851, 7851, 7870, 7877, 7897, 7914, 7929, 7934

which correspond to estimates of π/4 obtained by just placing a decimal point in front of each of these; we really have π/4 = .7853981634... It's not that surprising that these estimates are a bit spread out. The number of hits is given by a binomial distribution with n = 10000, p = π/4; thus its mean is 2500π (about 7854) and its standard deviation (10000 (π/4) (1-π/4))^1/2, or about 41.

However, I can try out the "Atlantic City Method" in the comfort of my living room. (I could try the Buffon's needle method here, since I have the right sort of floor, but I'm not that bored.) Picking a number uniformly at random from [0,1], and carrying out the "Atlantic City Method" as described above with u=.01, N=10000, I get the following values for M:

29 (four times), 30 (seven times), 32 (four times), 33 (two times), 36 (three times)

These each correspond to approximations of 1/π -- just put a decimal point in front. In reality we have 1/π = .3183.

If instead I use u = .1, the situation seems a bit better; I now get

315 (once), 316 (twice), 317 (twice), 318 (three times), 319 (six times), 320 (five times)

and again these correspond to approximations of 1/π if you put a decimal point in front; not bad! It seems surprising that increasing u should do this -- we're told that we should let u approach zero -- but we're able to reduce the ratio of the standard deviation of the number of hits to the mean, and thus make the approximations better, by increasing u.

However, if we increase u even more this doesn't hold up; for example if u = 1 then M = 3524 or 3525. The random aspect of the problem seems to go away but the constant we get is not the one we were looking for! (And if I didn't know the value of π already, then I wouldn't be able to assess this.)

Now, why should this method work? The map f: x -> (2-x)/(3-4x) is an example of a Möbius transformation, also known as a fractional linear transformation. This insight is important, because Mobius transformations can be viewed as isometries of the Riemann sphere; it's clear that this particular transformation takes real numbers to real numbers. My first instinct was that it had something to do with the magnitude of the rotation represented by this particular Mobius transformation, but that doesn't seem to be anything nice. However, near zero the projection taking the complex plane to the Riemann sphere very nearly doubles lengths. (See the picture at the Wikipedia article on the Riemann sphere.) So the length of interval u around zero is mapped to an arc of interval just slightly smaller than 2u on the "real circle" -- that is, the projection of the real axis, plus a point at infinity, onto the Riemann sphere. The sequence of iterates x, f(x), f(f(x)), ... is what one might call a "pseudo-random" sequence of points (it's certainly not random, but I'll think of it as if it were) and it lies in the arc in question a proportion N/M of the time; however it also lies in the arc a proportion just slightly less than 2u/(2π) of the time if you believe the psuedo-randomness, since the arc has length slightly less than u and the circle is a circle of radius 1. As u approaches 0 the error in this approximation goes to zero. So we get N/M = u/π, or π = uM/N. (Of course, this is a heuristic argument, and certainly not a proof, but I just wanted to convince myself that it made sense.) What's more, the actual identity of the linear transformation appears to be a red herring; I suspect this works for a rather large class of Mobius transformations. (Not all, though; for example, for some choices of f the sequence x, f(x), f(f(x)), ... approaches one of the fixed points of f.) I'm having trouble identifying other Mobius transformations that actually have this property, though, so I'm not totally confident in saying this.