Sleeping Cyborg

Jonathan David Page talks about whatever he happens to be thinking about. Sometimes other people join in.

Email · @parathetic (Twitter) · @jdpage (Github)
Subscribe to feed
 

Links

A collection of cool people and projects.

Approximating Pi Redux

by on 20 December 2012
in , ,
with some comments, maybe.

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 do, 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. I decided to do so. A detailed explanation follows the source code, 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 \over d\), \(d \gt 0\), that is closer to \(x\) than any approximation with a smaller denominator.

The first 16 terms output are:

rational decimal
\(3 \over 1\) 3.000000000000000
\(13 \over 4\) 3.250000000000000
\(16 \over 5\) 3.200000000000000
\(19 \over 6\) 3.166666666666667
\(22 \over 7\) 3.142857142857143
\(179 \over 57\) 3.140350877192982
\(201 \over 64\) 3.140625000000000
\(223 \over 71\) 3.140845070422535
\(245 \over 78\) 3.141025641025641
\(267 \over 85\) 3.141176470588235
\(289 \over 92\) 3.141304347826087
\(311 \over 99\) 3.141414141414141
\(333 \over 106\) 3.141509433962264
\(355 \over 113\) 3.141592920353983
\(52163 \over 16604\) 3.141592387376536
\(52518 \over 16717\) 3.141592390979242

This suggests that as far as rational approximations of pi go, \(355 \over 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.)

Sorry about the horizontal scrolling; the code is wrapped to 80 characters, but my blog is rather narrower than that. If it bothers you, you can view it on GitHub.

General Concepts

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 + \frac{1}{2}\). But not really 2, more like \(2 + \frac{1}{5}\). But not really 5, more like \(5 + \frac{1}{5}\). But not really 5, more like \(5 + \frac{1}{2}\). Actually, exactly like \(5 + \frac{1}{2}\). This leaves us with the following monster: 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 nicely set equations. $$3 + \cfrac{1}{2 + \cfrac{1}{5 + \cfrac{1}{5 + \cfrac{1}{2}}}}$$ 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.

pi_noncanon_seq()

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 pi. 1 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.

canonicalize(seq)

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 homographic Strictly speaking, if \(c = 0\), then it isn’t homographic. It’s a linear polynomial. function \(x=\frac{ax’ + b}{cx’ + d}\). This function establishes bounds on the value of \(x\). Assume that \(x’\) is between 0 and infinity. 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\). 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 $$x = \frac{a\left(p + \frac{q}{x''}\right) + b}{c\left(p + \frac{q}{x''}\right) + d}$$ Multiplying the numerator and denominator by \(x''\) simplifies this to $$x = \frac{\left(pa+b\right)x'' + qa}{\left(pc+d\right)x'' + qc}$$ 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. $$x = t + \left(\frac{ax’ + b}{cx’ + d} – t\right)$$ $$x = t + \frac{ax’ + b – t\left(cx’ + d\right)}{cx’ + d}$$ $$x = t + \frac{1}{\frac{cx’ + d}{\left(a-tc\right)x’ + \left(b-td\right)}}$$ 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.

condense()

This function follows the same pattern as canonicalize(), but instead of spitting out the integer 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.

seq_cmp(seqa, seqb)

This one is pretty much straight from Gosper:2

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.

best_rational_approximations(seq)

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 equivalent Because Gosper2 says so! to checking whether the terms encountered so far, reversed, are a greater continued fraction than the remaining terms. 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.

  1. a According to doi:10.2307/2589152, it’s a special case of the arctan function.
  2. a b R. W. Gosper’s unpublished paper on continued fractions