Continued Fractions

I was inspired recently by a wonderful Mathologer video to implement continued fractions in Scala. I chose Scala for two reasons: first (as you would know if you've been here before) is that it's my favorite language and, second, it requires lazy programming.  What are continued fractions? And why are they so interesting?

They've been around, essentially, since Euclid, although they hit their heyday in the 17th century with the likes of Leibniz, Wallis, and Euler (who came a little later), and many other less well known mathematicians. The particular subject of the video (all his videos are great, BTW) is the solution to the "Strand puzzle." This was a puzzle set in the Strand magazine (the one that published the Sherlock Holmes stories) by one of the great puzzle setters of all time: Henry Dudeney.

A man lives at house number n on a long street of N houses (all on the same side), numbered from 1 thru N. He notices that if he sums all the house numbers on one side of his (L) and all of the house numbers on the other side (R), then L = R. What is his house number?

Supposedly, the brilliant mathematician Ramanujan immediately recognized that the solutions were related to the convergents of a particular continued fraction that gives solutions to the equation Y^2 -- 2 X^2 = 1 -- one of the so-called Pell equations, after John Pell (1611-1685). I will let you watch the video to get the full story.

So, the general Pell equation Y^2 -- x X^2 = 1 relates to the following continued fraction:

When we set x = 2 (for the Strand magazine problem), we get successive approximations to the square root of two.

Here's another rather elegant definition of phi, the golden ratio:

So, how do we code this in Scala? The most important class is called ConFrac and is defined thus:

class ConFrac(val b: Long, co: => Option[CF]) {

  lazy val tailOption: Option[CF] = co

...
}

You might notice that this is not a case class. That's because we need the second parameter to be lazily evaluated, i.e. "call by name". Since it's not a case class, we need to ensure that we can reference the first parameter, b, in other places, using the val keyword. This doesn't work, of course, for the second parameter because of course that is, in reality, a function that generates an Option[CF]. Hence the lazy val (or def) defined as tailOption. Obviously, we need to define CF too:

case class CF(a: Long, c: ConFrac)

This is a case class, and defines the numerator a and identifies the next ConFrac in the series as c.

How do we get the value of our continued fraction? Well, it's an infinite series so we can't evaluate the whole thing. But, we can look at the convergents which asymptotically approach the true value. Here's the method for getting the convergents:

  def convergents: LazyList[Rational] = {
    def inner(an_2: BigInt, an_1: BigInt, bn_2: BigInt, bn_1: BigInt, w: => LazyList[Pair]): LazyList[Rational] = w match {
      case LazyList() => LazyList()
      case p #:: tail =>
        val an = p.b * an_1 + p.a * an_2
        val bn = p.b * bn_1 + p.a * bn_2
        Rational(an, bn) #:: inner(an_1, an, bn_1, bn, tail)
    }

    val h #:: z = coefficients
    Rational(h.b) #:: inner(1, h.b, 0, 1, z)
  }

Note that the inner method is not tail-recursive. I don't think there's a way to make it tail-recursive. This method returns a LazyList[Rational], and we can evaluate it as far as we want. Scala fans will note that I used a pattern in the definition val h #:: z = coefficients. This is so much more elegant than any other way of getting the information I needed. Of course, I have to be sure that coefficients is not empty. It makes use of another method coefficients:

  def coefficients: LazyList[Pair] = {
    def inner(_b: Long, a: Long, co: Option[CF]): LazyList[Pair] =
      Pair(_b, a) #:: {
        co match {
          case None => LazyList()
          case Some(cf) => inner(cf.c.b, cf.a, cf.c.tailOption)
        }
      }

    inner(b, 0, tailOption)
  }

This uses another case class Pair which is is simply two Longs: b and a. The result of coefficients is a LazyList of Pairs.  The form of continued fraction described here looks as follows where the nth element of the result is (bn, an) and a0 is always zero:

So, we can get values from a ConFrac. But how do we construct them in the first place? Obviously, declaring new ConFrac(1, new ConFrac(...)) wouldn't be so convenient! Instead, there are methods in the ConFrac companion object as follows:

  def apply(ps: LazyList[Pair]): ConFrac = ps match {
    case p #:: LazyList() => new ConFrac(p.b, None)
    case p #:: tail => new ConFrac(p.b, Some(CF(p.a, ConFrac(tail))))
  }

  def simple(xs: LazyList[Long]): ConFrac = ConFrac(Pair.zip(xs, LazyList.continually(1L)))

Pair.zip is a method in Pair's companion object which zips the two lazy lists together and then maps each resulting tuple to a Pair. The purpose of simple is to make it easier to define a so-called simple continuous fraction (where all of the a terms are unity).

Here are some of the definitions of simple LazyLists:

• phi (the golden ratio): LazyList.continually(1L)
• e (Euler's number): 2L #:: LazyList.iterate(Seq(1L, 2, 1)) { case Seq(x, y, z) => Seq(x, y + 2, z) }.flatten

And here is one of the definitions of lazy lists for generalized continued fractions:

• pi (one of several): Pair.zip(0L +: LazyList.from(0).map(2L * _ + 1), 4L +: LazyList.from(1).map(x => 1L * x * x))

It's been fun working on this project. There is one method that is currently not well implemented. It is the toDouble method. The signature is: def toDouble(epsilon: Double): Option[Double] and it is required to provide an approximation that is correct to within epsilon. Such a method is required if mathematics with known precision is being performed. If you have a required precision of, say, pi to 1E-20, then you need to be able to ensure that the result is true to that precision.

One thing that has tripped me up a few times working with LazyLists is that there are many methods which LazyList inherits through Seq which are not "call by name." An example of this is +: which prepends an element to a LazyList. But, unfortunately, it doesn't do it lazily ;)

I think it will work out to be a good assignment for my Scala class!