A clear descent into geekdom, I started searching for rational approximations to pi. Everybody knows that 22/7ths is close to pi, but what's the next best one?
A brief mathemetical review: a rational number is one that can be expressed as an integer divided by another integer (i.e., a/b). The math question is, what a's and b's work well and how are they distributed. The computer science question is what's an easy algorithm to find the next best one given an existing one.
So I wrote a program to figure out the answers (exhaustively...I'm a computer ``scientist''). To jump to the chase, here's the:
Because you were going to ask, yes the program generalizes:
If you can't wait, here's the summaries:
#h a b r e new 3 1 3.0000000000000000 1.4159265358979312e-01 * 13 4 3.2500000000000000 1.0840734641020688e-01 * 16 5 3.2000000000000002 5.8407346410207062e-02 * 19 6 3.1666666666666665 2.5074013076873403e-02 * 22 7 3.1428571428571428 1.2644892673496777e-03 * 179 57 3.1403508771929824 1.2417763968106676e-03 * 201 64 3.1406250000000000 9.6765358979311600e-04 * 223 71 3.1408450704225350 7.4758316725809237e-04 * 245 78 3.1410256410256410 5.6701256415214729e-04 * 267 85 3.1411764705882352 4.1618300155787935e-04 * 289 92 3.1413043478260869 2.8830576370619809e-04 * 311 99 3.1414141414141414 1.7851217565167943e-04 * 333 106 3.1415094339622640 8.3219627529107498e-05 * 355 113 3.1415929203539825 2.6676418940496660e-07 * 52163 16604 3.1415923873765359 2.6621325721620792e-07 *
IMHO, only 22/7 and 355/113 are even close to ``useful''.
#h a b r e new 3 1 3.0000000000000000 1.4159265358979312e-01 * 6 2 3.0000000000000000 1.4159265358979312e-01 . 9 3 3.0000000000000000 1.4159265358979312e-01 . 13 4 3.2500000000000000 1.0840734641020688e-01 * 16 5 3.2000000000000002 5.8407346410207062e-02 * 19 6 3.1666666666666665 2.5074013076873403e-02 * 22 7 3.1428571428571428 1.2644892673496777e-03 * 25 8 3.1250000000000000 1.6592653589793116e-02 . 28 9 3.1111111111111112 3.0481542478681956e-02 . 31 10 3.1000000000000001 4.1592653589793027e-02 . 35 11 3.1818181818181817 4.0225528228388541e-02 . 38 12 3.1666666666666665 2.5074013076873403e-02 . 41 13 3.1538461538461537 1.2253500256360628e-02 . 44 14 3.1428571428571428 1.2644892673496777e-03 . 47 15 3.1333333333333333 8.2593202564598123e-03 . 50 16 3.1250000000000000 1.6592653589793116e-02 . 53 17 3.1176470588235294 2.3945594766263678e-02 . 57 18 3.1666666666666665 2.5074013076873403e-02 . 60 19 3.1578947368421053 1.6302083252312194e-02 . 63 20 3.1499999999999999 8.4073464102067952e-03 .
This table shows the basic structure: you'll find a cluster of better approximations, then a long stretch of not very good ones.
I'm not a mathematician (in fact, I had to have ispell help me spell it, but that's another story). What's here is probably obvious to someone with real training. (In fact, see below for a better way to figure rational approxmations, for some definition of better.) But perhaps you'll find it a little bit amusing.
If you want to know more about pi (much more), go to Yahoo and search for pi, or (much better still), talk to your local math teacher.
These are my comments about the program and the algorithm. For full details, see the program source code.
usage: rational_approximation [-v] [-l denominator-limit] [-t target-value]
Outputs a list of rational approximations to pi. By default, data is shown only when there is a new best approximation. The verbose flag (-v) shows data for each new denominator. The debug flag (-d) shows additional information. "-l X" sets the last denominator to try to X (by default, infinity), and "-t Y" sets the value to search for (by default pi).
The algorithm is to start with a rational approximation (initially the closest integer, 3/1 for pi) and refine that approximation by testing each successive denominator.
We make use of the fact that is we know a rational number a/b which we want to use to approximate R, then we know that the closest approximation with the next higher denominator is c/(b+1) which is one of (c'-1)/(b+1), c'/(b+1), (c'+1)/(b+1) where c' = a + int(a/b).
Output format is a table of the format
#h a b r e new 3 1 3.0000000000000000 1.4159265358979312e-01 * 13 4 3.2500000000000000 1.0840734641020688e-01 * 16 5 3.2000000000000002 5.8407346410207062e-02 * 19 6 3.1666666666666665 2.5074013076873403e-02 * 22 7 3.1428571428571428 1.2644892673496777e-03 *
where each column is a different value. The first two (a, b) are the numerator and denominator of the rational approximation. The third (r) is the decimal expansion of a/b. The fourth (e) is the error, i.e. abs(pi - a/b). The fifth (new) is a flag indicating when the approximation is a new best.
The algorithm is not written using numerically sound techniques, so its accuracy is limited (although this has not yet been a limit to the program's amusement value.) The problem is that it depends on the accuracy of subtracting two floating point numbers. Since floating point numbers are of finite precision (unlike the mathematical concept of real numbers that they approximate), this means that you cannot trust this program to approximate pi closer than the number of digits (minus one) present in your computer's implementation of floating point. For IEEE floating point present in all modern computers gives you at least 14 decimal digits, although I should look up (or get the energy to compute) exactly how accurate it is. Eventually the program should be re-written using a bignum package to get around this problem.
Another comment about the algorithm: currently we iterate through all denominators. This seems easy and works well since computers are fast. It might be better instead to use the error value reported to jump more quickly to the next plausible candidate. (I.e., 22/7ths is accurate to ~1.3e-3, and numbers expressed in units of a/7 have a resolution of ~1.4e-1. From these two we should be able to figure out that we can't possibly do better for a while. I haven't worked out the exact relationship yet.)
It turns out there's a much more clever way to compute rational approximations if your definition of better is converges fasters. In fact, it finds 104348/33215 in four steps. I'll let the person who explained it to me explain it to you (thanks, Brian). (This method takes advantage of the error as I suggested above but didn't work through.)
In the process of writing up this page, I came across a couple of other interesting web pages relating to numbers: