No, it isn’t Pi Day, or anything resembling it. I was editing some stuff, and I noticed the original “Approximating Pi” article. I’d been meaning to rewrite the code since I learned about continued fractions, and since I didn’t have anything better to do1, I decided to do so. A detailed explanation follows the source code2, since this is less mathematically facile than the previous version.
This version is slightly different from the previous one, in that it generates all of the best rational approximations of pi, in order. According to Wikipedia:
A best rational approximation to a real number x is a rational number n/d, d > 0, that is closer to x than any approximation with a smaller denominator.
The first 16 terms output are:
This suggests that as far as rational approximations of pi go, 355/113 looks like it might be the best bang for your buck. (The jump in size without much improvement corresponds to an abnormally large term in the continued fraction expansion.)
This code relies pretty heavily on continued fractions, which take a little bit of explaining. They’re basically a slightly arcane way of writing real numbers which has some nice properties. Let’s take 3.456 — or 432/125 — as an example. It’s kind of like 3, but a bit more. More like 3 + 1/2. But not really 2, more like 2 + 1/5. But not really 5, more like 5 + 1/5. But not really 5, more like 5 + 1/2. Actually, exactly like 5 + 1/2. This leaves us with the following monster3:
But instead of writing out all of those 1s, we could just write [3; 2, 5, 5, 2]. That’s a continued fraction. If you want to make things more complicated, you could have something other than 1s in the numerators, in which case you’d have a generalized continued fraction. (The normal kind are also known as simple continued fractions, to distinguish them from the generalized kind.)
Why would you want such a horrible thing? There are a variety of reasons (some of which I don’t even pretend to understand), but in this case, we’re going to take advantage of their approximation superpowers.
You know how I said that you can make things more complicated? Well, that’s exactly what this function does. It spits out a generalized continued fraction which is equal to pi4. We have to use a generalized continued fraction, because the simple continued fraction for pi has no known pattern, and appears to be more or less random.
Because generalized continued fractions tend to not only be slightly inconvenient to work with, but also lack the property that we’re after (more on that later), we convert it into a simple continued fraction.
Given a continued fraction, define x to be the expansion of the whole sequence, and x' to be the part that we don’t know. Thus, before we receive any terms, we have x = x'.
In this function, the variables a through d represent the homographic5 function . This function establishes bounds on the value of x. Assume that x' is between 0 and infinity6. At one end, x' = 0, so that x = b/d. On the other end, if we take the limit as x' approaches infinity, x = a/c. These represent our lower and upper and bounds of x, not respectively.
Now, if we take a term (p, q), indicating that the unknown part x' = p + q/x'', where x'' is the part that is still unknown, we can substitute into our homographic form and get
Multiplying the numerator and denominator by x'' simplifies this to
This is equivalent to the expression
a, b, c, d = p*a+b, q*a, p*c+d, q*c in the code.
We then evaluate the integer parts of the bounds, and when they are the same, we know for sure what the integer part of the number is, so we can yield that, and remove it from the total.
The reciprocation is there, of course, because we’re trying to do a simple continued fraction, and we want it to match general form above. This is equivalent to the expression
a, b, c, d = c, d, a-t*c, b-t*d.
The result is that the function takes in terms of the generalized continued fraction, until it can produce a term of the simple continued fraction, at which point it does so, removing it from the running total.
This function follows the same pattern as
canonicalize(), but instead of spitting out the integral part as soon as it knows it, waits until it has taken all of the terms (note that it takes a finite sequence), and outputs a rational number equivalent to the simple continued fraction we started with. Note that the upper bound is given, since the last term of a terminating continued fraction is always infinity.
This one is pretty much straight from Gosper7:
If we call the integer part of a continued fraction the 0th term, then we can say that the a (regular) continued fraction is an increasing function of its even numbered terms, and a decreasing function of its odd numbered terms. Thus, to compare two continued fractions, compare the terms at the first position where they differ, then reverse the decision if this position is odd numbered. If one continued fraction terminates before differing from the other, regard the shorter one to be suffixed by an infinite term.
Here, I use
None to represent infinity.
This is the function which performs the main purpose of the script. It essentially takes advantage of the fact that you can use the continued fraction to generate a best rational approximation of a number by truncating the continued fraction at some point, then converting it to a regular fraction. Furthermore, the last term of the truncated sequence can be reduced by up to half of its value to produce yet more best rational approximations. There’s a borderline case for even trailing terms, since they are divisible by two exactly; in that case, the result is a best rational approximation only if it’s better than the previous one in the sequence. This is equivalent8 to checking whether the terms encountered so far, reversed, are a greater continued fraction than the remaining terms9.
This is a lie. I do, in fact, have a number of things to do that various people might define as “better”, but I didn’t feel like doing any of them. I felt like sitting around in my bathrobe and writing Python. ↩
Sorry about the horizontal scrolling; the code is wrapped to 80 characters, but my blog is rather narrower than that. If it bothers you, there’s a link at the bottom of the code so that you can view it on GitHub. ↩
If you’re reading this in an RSS reader or something, you may want to consider clicking over to the actual site. There’s a script which converts all of the LaTeX to images. ↩
Strictly speaking, if c = 0, then it isn’t homographic. It’s a linear polynomial. ↩
Since x' is a continued fraction composed of positive integers, it’s trivial to demonstrate that x' is itself positive. Note that it’s a toroidal space, so the interval could be inside or outside the bounds, depending on a variety of things, including the sign of d. ↩
Because Gosper says so. ↩
The nonsense with
tee()is due to the fact that Python generators are so narrow-minded as to go one-way, and they can’t be cloned. Haskell lists are actually perfect for continued fractions in pretty much every way, but the Haskell version of this code is, unsurprisingly, even more incomprehensible than this, so I decided not to complicate the issue. ↩