K-term Continued Fractions

Another interesting set of examples in SICP are the exercises dealing with infinite continued fractions, and their k-term truncations. [That is, Exercises 1.37, 1.38, and 1.39.]

As given in the book, an infinite continued fraction looks like this:

Of course that expression doesn't ever stop, so we approximate it with a k-term continued fraction, where we cut it off after k terms. So it looks like this:

Since we have a finite representation now, we can easily write code to work with such structures! So we'll do just that, using Scala as our JVM language.

Coding the Function

So first, we need to decide how we're actually going to represent one of these continued fractions. One nice way to look at this (and the way that the SICP book uses, coincidentally enough) is to instead suppose that we have two infinite sequences (N)_i and (D)_i , out of which we construct the fraction. At any point in the fraction, N_i is divided by its corresponding denominator D_i plus the rest of the fraction.

Then we can construct the k-term by applying the same idea, except this time we have the finite sequences N₁ , N₂ , …, N_k and D₁ , D₂ , …, D_k .

So for our function, we want to take in three parameters:

• a function n which takes an integer index i and returns the numerator N_i ,
• a function d which takes an integer index i and returns the denominator D_i ,
• an integer k that determines the end of the finite continued fraction.

Now if you look back at the structure for the k-term continued fraction, you can see that it will be easiest for us to actually start from the bottom kth term and work our way back up.

Suppose that we start with a variable result set to the initial value 0. Then at each point we divide n(i) by the addition of d(i) and whatever result was, and store the answer back into result.

So we divide n(k) by (d(k) + 0), and then divide n(k-1) by (d(k-1) + n(k)/d(k)), etc.

This translates directly into a Scala function here:


And that code will work just fine. But there's something sort of off about it, right? We have that variable declaration and everything, but all we're doing is accumulating a result as we iterate across a data structure.

Using Scala gives us access to all of the awesome techniques of functional paradigm programming. In particular, all we're doing here is applying a simple fold to a list of indices. So we could write the function more idiomatically as follows:


Aside from how they're written though, both of the functions are pretty much equivalent. They both have O(k) time-complexity (where k is the number of terms we want) and run with constant O(1) space (they only need that one extra variable to store the result at a given point.)

So let's not worry any more about what the function looks like, and instead start using it to compute some interesting continued fractions!

An Approximation for the Golden Ratio Reciprocal

The golden ratio is defined as follows:

This number pops up all over the place (for example, it's the limit of the ratio between successive Fibonacci numbers, as we saw in a previous post.)

Since we find the golden ratio everywhere, it's not too surprising that we'd find its reciprocal turning up in the simplest infinite continued fraction:

So we can use our function at the Scala REPL to get an approximation for the reciprocal as follows.

     continuedFraction(i => 1, i => 1, 10)
     //> res0: Double = 0.6179775280898876

     continuedFraction(i => 1, i => 1, 100)
     //> res1: Double = 0.6180339887498948

     2 / (1 + Math.sqrt(5))
     //> res2: Double = 0.6180339887498948

That's a pretty awesome approximation! The 100-term fraction gives us the same Double precision as working out the square root of 5 and performing the division.

An Approximation for e

The book also suggests a continued fraction representation for e-2, attributed to Euler:

Here the functions to generate the numerator and denominator sequences are a little trickier.

We have that the numerators are all 1, so that's fine. But the denominator sequence looks like this (the indices for each element are included below):

   1  2  1  1  4  1  1  6  1   1   8   1   1  10   1   1  12 ...
   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17 ...

We can notice that the only indices that don't correspond to 1 are 2, 5, 8, 11, 14, 17, etc. Let's see if we can discover a pattern when we isolate those indices:

   2  4  6   8  10  12
   2  5  8  11  14  17

Hmm. We have the multiples of 2 on the top, and we sort of have the multiples of 3 on the bottom, if we shift the indices forward by one:

   2  4  6   8  10  12
   3  6  9  12  15  18

Nice! So we see here that the indices that didn't correspond to 1's before are simply the multiples of 3 shifted back by one. Now if we divide these new index values by 3, we get the following:

   2  4  6  8  10  12
   1  2  3  4   5   6

Perfect. We now have an obvious relationship between the two sets of numbers here; we simply times each of the bottom values by 2.

So to develop a function that will output the value for the denominator at the index i, we do the following:

1.  Add one to the index i.
2.  If this new value is not divisible by 3, return 1.
3.  Otherwise, divide the new value by 3 and then multiply it by 2, then return that value.

So using this, we can write the ancillary function for the denominator values, and then compute the k-terms for the e approximation (remembering to add 2; our continued fraction gives us e-2 after all.)

     def d(i: Int) = if ((i + 1) % 3 != 0) 1 else ((i + 1) / 3) * 2
     //> d: (i: Int)Int

     2 + continuedFraction(i => 1, d, 10)
     //> res3: Double = 2.7182817182817183
 
     2 + continuedFraction(i => 1, d, 100)
     //> res4: Double = 2.7182818284590455
 
     Math.E
     //> res5: Double = 2.718281828459045

Once again, we get quite a bit of accuracy just for the 100-term!

An Approximation to the Tangent Function

Let's do something really cool to end this post on. We have Lambert's continued fraction representation for the tangent function as follows:

Here x is the angle in radians. Note also that we have subtraction of each successive term instead of addition; this is easy to get around though. Well just say that all of the x² numerators are -(x²) instead, and we can use our addition as usual.

So we can use our continuedFraction() function to write our own tangent function:

• Our tangent function will take in the value x in radians.
• Inside the function, we will define the numerator and denominator functions. The numerator function returns x if the index is 1, and -(x²) otherwise (the numerator and denominator functions are closures, so we have x available in our scope.) The denominator returns the odd number defined by 2i – 1.
• We'll go for the 100-term finite continued fraction, so we have a good approximation.

The code for our tangent function is then as follows:


Then we can use our tangent function as usual. We'll compare the results against the normal Math.tan() function:

     tan(2)
     //> res6: Double = -2.185039863261519

     Math.tan(2)
     //> res7: Double = -2.185039863261519

     tan(Math.PI)
     //> res8: Double = -1.4135798584282297E-16

     Math.tan(Math.PI)
     //> res9: Double = -1.2246467991473532E-16

     tan(Math.PI / 4)
     //> res10: Double = 1.0

     Math.tan(Math.PI / 4)
     //> res11: Double = 0.9999999999999999

So we can even use the continued fraction idea to implement a working tangent function! And it's pretty accurate, too!