Sleeping Cyborg: tag

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¹, 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/d, d > 0, that is closer to x than any approximation with a smaller denominator.

The first 16 terms output are:

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

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.)

View the code 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 + 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 monster³:
$3 + \frac{1}{2 + \frac{1}{5 + \frac{1}{5 + \frac{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⁴. 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⁵ function $x=\frac{ax' + b}{cx' + d}$ . This function establishes bounds on the value of x. Assume that x' is between 0 and infinity⁶. 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 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.

seq_cmp(seqa, seqb)

This one is pretty much straight from Gosper⁷:

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⁸ to checking whether the terms encountered so far, reversed, are a greater continued fraction than the remaining terms⁹.

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. ↩
According to doi:10.2307/2589152, it’s a special case of the arctan function. ↩
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. ↩
R. W. Gosper’s unpublished paper on continued fractions ↩
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. ↩

You didn’t get a post for July 22, so here’s some Python which finds approximations for pi.

Update 20 Dec 2012: See also the revised version utilizing continued fractions.

In honor of everyone’s favourite mathematical constant, March 14, or 3/14 is Pi Day. Hooray!

Next up: Tau Day on June 28th, or 6/28. After that we get “European Pi Approximation Day” on July 22nd, or 22/7. Such fun we can look forward to!

Ever since I first ran into this proof a few years ago, I’ve always thought it rather elegant. It came up (read: I brought it up) while studying for a proofs test, so I thought I would post it here.

Claim

$\sqrt 2$ is irrational.

Proof (by contradiction)

Suppose $\sqrt 2$ is a rational number. This means that it can be expressed as an irreducible fraction $\frac{p}{q}$ where both and are integers. So $\sqrt{2} = \frac{p}{q}$ , thus $2 = (\frac{p}{q})^2$ , thus $2 = \frac{p^2}{q^2}$ , so 2q^2 = p^2 . Since p^2 has a factor of 2, p^2 is even, and therefore is also even, meaning that p = 2k where is an integer. Substituting, 2q^2 = (2k)^2 , thus 2q^2 = 4k^2 , thus q^2 = 2k^2 . Since q^2 has a factor of 2, it is also even, so is even. Therefore, and are both even, and so share a common factor of 2. But we assumed that $\frac{p}{q}$ was irreducible. This is a contradiction. Therefore $\sqrt 2$ is irrational. □

The Ackermann function is a mathematical function that returns really really large numbers. For example, A(4,2) is an integer with over 19,700 decimal digits.

Graham’s number is a holy-moly massively horrifyingly huge number which makes a googolplex look like diddly squit. Actually, it makes numbers which make numbers which make googolplex look like diddly squit look like ant droppings look vanishingly small. It it unimaginably huge. I do not have words to explain how much massively huger than pretty much any number you can think of it is.

As you can imagine, $A(g_{64}, g_{64})$ is quite a large number.

PS— I’ll answer the $e^{\pi i}$ question soon, I want to type out a proper explanation for that.

Sleeping Cyborg

Approximating Pi Redux

General Concepts

pi_noncanon_seq()

canonicalize(seq)

condense()

seq_cmp(seqa, seqb)

best_rational_approximations(seq)

Approximating Pi

Pi Day!

Proof of the Square Root of 2’s Irrationality

Claim

Proof (by contradiction)

bearjcc asked: What number do you get when you call the Ackermann function with Graham’s number as arguments?