Programmatic adventures in the hyperbolic plane


Recursion is dangerous. So dangerous that in my professional life I’m explicitly banned from using it in code. After all who would want their aeroplane to stall, their car to crash or their fridge to defrost because of a stack overflow.

But it is also, most (but not all) programmers would agree, beautiful and more than a bit mysterious. Indeed one of the biggest selling popular science/mathematics books of all time centres on recursion as a possible partial explanation of human consciousness.

In the last week or so I’ve been attempting to use recursion in R to draw a tessellating pattern of squares on the conformal representation of the hyperbolic plane – an idea from M.C. Escher via Nobel laureate Sir Roger Penrose.

The code here is where I’ve got to.

#!/usr/bin/env Rscript

library(ggplot2)
library(ggforce)

basicLengths<-seq(from<-0, to<-1.0, by<-1/250000)
projLengths<-basicLengths * (2 /(1 + basicLengths * basicLengths))

getBestMatch<-function(numIn)
{
  repVec<-basicLengths[which.min(abs(numIn - projLengths))]
  return(repVec)
}


drawConformalLine <-function(startX, startY, endX, endY)
{
  #Length and angle of line to starting conformal point
  lengthToStart <- sqrt(startX * startX + startY * startY)
  angleToStart <- atan2(startY, startX)
  #And to end Point
  lengthToEnd <- sqrt(endX * endX + endY * endY)
  angleToEnd <- atan2(endY, endX)
  #Project
  newStartLength<-lengthToStart * (2 / (1 + lengthToStart * lengthToStart))
  newEndLength<-lengthToEnd * (2 / (1 + lengthToEnd * lengthToEnd))
  newStartX<-newStartLength * cos(angleToStart)
  newStartY<-newStartLength * sin(angleToStart)
  newEndX<-newEndLength * cos(angleToEnd)
  newEndY<-newEndLength * sin(angleToEnd)

  xpoints<-seq(from<-newStartX, to<-newEndX, by<-(newEndX - newStartX)/50)
  ypoints<-seq(from<-newStartY, to<-newEndY, by<-(newEndY - newStartY)/50)
  euclidLine<-data.frame(xpoints, ypoints)
  #Conformal line (reverse projection)
  lengths<-sqrt(euclidLine$xpoints * euclidLine$xpoints + euclidLine$ypoints * euclidLine$ypoints)
  newLengths<-vector()
  for (i in 1:length(lengths))
  {
    newLengths<-c(newLengths, getBestMatch(lengths[i]))
  }
  angles<-atan2(euclidLine$ypoints, euclidLine$xpoints)
  newPoints<-data.frame(xp = cos(angles) * newLengths, yp = sin(angles) * newLengths)
  return(newPoints)
}

recursiveSquares<-function(setOfPoints, base, radius, sign)
{
  if (base + radius > 0.75){
    return (setOfPoints)
  }
  setOfPoints<-rbind(setOfPoints, recursiveSquares(setOfPoints, base + radius, radius/2, sign))
  positiveCorner = (sign * base + sign * radius)
  negativeCorner = (sign * base)
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, positiveCorner, positiveCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, negativeCorner, negativeCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, negativeCorner, negativeCorner, positiveCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, positiveCorner, positiveCorner, positiveCorner))
  return (setOfPoints)
}


partitions<-100
distances<-seq(from=1, to=partitions, by=1)
distances<-exp(distances/partitions)/(2 * distances)
circles<-data.frame(x0=rep(0, partitions), y0=rep(0, partitions), r=distances)
t <- ggplot() + geom_circle(aes(x0=x0, y0=y0, r=1-r, color=rgb(r,1-r,1)), data=circles)
t <- t + geom_circle(aes(x0=0, y0=0, r=1), color="black", data=data.frame(x0=0, y0=0, r=1)) + coord_fixed()
distances<- 1-distances
n_distances<-(distances * 2)/(1 + distances * distances)
n_circles<-data.frame(x0=rep(0, partitions), y0=rep(0, partitions), r=n_distances)
y <- ggplot() + geom_circle(aes(x0=x0, y0=y0, r=r, color=rgb(r,1-r,1)), data=n_circles)
y <- y + geom_circle(aes(x0=0, y0=0, r=1), color="black", data=data.frame(x0=0, y0=0, r=1)) + coord_fixed()
square<-0.4
squarePoints<-data.frame()
squarePoints<-recursiveSquares(squarePoints, 0, square, 1)
nPoints<-data.frame()
nPoints<-recursiveSquares(nPoints, 0, square, -1)
squarePoints<-rbind(squarePoints, nPoints)
t<-t + geom_point(aes(color="magenta", x=squarePoints$xp, y=squarePoints$yp))
t<-t + geom_point(aes(color="magenta", x = -squarePoints$xp, y=squarePoints$yp)) 


I had several problems to solve, and recursion seemed a good fit for some of them, though in the end I’ve only been able to make partial progress.

Here are some of the problems and how I tackled them:

  1. Finding the straight lines

In the conformal representation of the hyperbolic plane ‘straight lines’ – in the sense of the shortest routes between two points – are actually represented by the segment of the circle joins the two points and that meets the boundary of the representation of the plane orthogonally (ie., at right angles).

For points that lie on the same polar angle in the circle that’s easy – the ‘straight line’ that joins them is simply the Euclidean line that passes through both of them and the origin (centre of the circular representation of the plane). The circle they are linked by in the conformal representation is of infinite radius.

But what about other, non-degenerate, cases? Perhaps there is an algorithmic method to work this out? But if there is I don’t know what it is and I couldn’t find it online either!

What I do know (again thanks to Sir Roger) is that is we substitute the projective representation for the conformal representation we can join the points by Euclidean lines. To go from the conformal to the projective representation of the hyperbolic plane by preserving the polar angle of points but magnifying the distance from the origin by \frac{2R^2}{R^2+r^2_c}, where R is the radius of the bounding circle (so R=1 in our case)and r_c is the length of the of the Euclidean line from the centre (origin) to the point.

The great advantage of the projective representation is that straight lines between points are just what we think of them as being in Euclidean geometry – in other words both the shortest length between two points and at a constant angle.

So to get the points on the line in the conformal representation the thing to do seemed to be to convert the points to the projective representation, calculate the straight line (very easy) and shrink everything back by the factor \frac{1+r^2_c}{2} .

But the problem with that is, that with the exception of the start and end points, we don’t know what r^2_c is to begin with – indeed that is more or less what we are trying to calculate.

My answer to that is to generate lookup tables (relatively quickly done using R’s vectorisation) and a function that indexes one to the other:

basicLengths<-seq(from<-0, to<-1.0, by<-1/250000)
projLengths<-basicLengths * (2 /(1 + basicLengths * basicLengths))

getBestMatch<-function(numIn)
{
  repVec<-basicLengths[which.min(abs(numIn - projLengths))]
  return(repVec)
}

2. Generate the squares

This is where the recursion comes in – I settled on the idea of four squares touching at each vertex, meaning each square as we travelled out from the centre the Euclidean lengths halve.

Here is (a simplified) version of the code:

recursiveSquares<-function(setOfPoints, base, radius)
{
  if (base + radius > 0.75){
    return (setOfPoints)
  }
  setOfPoints<-rbind(setOfPoints, recursiveSquares(setOfPoints, base + radius, radius/2, sign))
  positiveCorner = (base + radius)
  negativeCorner = (base)
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, positiveCorner, positiveCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, negativeCorner, negativeCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, negativeCorner, negativeCorner, positiveCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, positiveCorner, positiveCorner, positiveCorner))
  return (setOfPoints)
}

Those of you who have read the classic Structure and Interpretation of Computer Programs may immediately recognise a weakness here – by not using tail recursion I am putting more stress on the stack.

And it’s true that an alternative tail recursive form (this is the full code) works:

recursiveSquares<-function(setOfPoints, base, radius, sign)
{
  if (base + radius > 0.75){
    return (setOfPoints)
  }
  positiveCorner = (sign * base + sign * radius)
  negativeCorner = (sign * base)
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, positiveCorner, positiveCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(positiveCorner, negativeCorner, negativeCorner, negativeCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, negativeCorner, negativeCorner, positiveCorner))
  setOfPoints<-rbind(setOfPoints, drawConformalLine(negativeCorner, positiveCorner, positiveCorner, positiveCorner))
  recursiveSquares(setOfPoints, base + radius, radius/2, sign)
}

But it still is liable to bust the stack – indeed the whole system quite finely balanced. Make the initial squares too small and it is easy to break the stack – though the tail recursive form does seem a little more robust.

3. Getting a tessellation/exploiting the dihedral symmetry

Unfortunately the tendency of the stack to blow up makes it very difficult (as far as I can see anyway) to write some elegant code that will seamlessly deliver the tessellations sought. But some simple group theory can help us get some more patterns.

The pattern above shows rotational and reflective symmetry – in other words the shape has the symmetries of square (dihedral group D_4 ).

We use this in the code through the <code>sign</code> factor in the calls to <code>recursiveSquares</code> and then in the two lines:

t<-t + geom_point(aes(color="magenta", x=squarePoints$xp, y=squarePoints$yp))
t<-t + geom_point(aes(color="magenta", x = -squarePoints$xp, y=squarePoints$yp)) 

These effectively reflect earlier results.

Who has been given the cash changing problem?


Some class, somewhere, has obviously been given the recursive money changing problem as a piece of work, because I have had several hundred visits in the last week from people seeking to get a grip on it.

Here’s the best solution, either buy Structure and Interpretation of Computer Programs or simply read it for free online.

The recursive cash changing algorithm


This is another one from Structure and Interpretation of Computer Programs which is a great book, but not a particularly easy read (not least because Scheme, the language in which the programs in the book is written is just different from C/Java etc). This example annoyed me for some time so I am writing down an explanation in case someone else wants to look it up…

Anyway this one is about the number of ways of to change cash. The algorithm is based on the fact that you can break down the number of ways to sort change of the sum n with t types of coin as follows:

  • All ways to sort n with t-1 coins; plus
  • All ways to sort n-d with t coins, where d is the denomination of the coin excluded in the above case.

(Think about it for five pence with two pence and one penny coins. One way to sort with pennies and two ways to change three pence with a two pences and pennies.)

To solve this recursively we have to define the degenerative cases – ie the points at which the recursion stops. These are:

  • If the sum to be cashed to change is 0 then there is 1 way to do it;
  • If the sum to be cashed to change is less than 0 then there are 0 ways to change it;
  • If there are 0 coins then there are 0 ways to offer the change.

So how does it work … hopefully the graphic below shows how (changing 10 pence with fives, twos and ones):

Coin changing algorithm tree

The idea is that f(n, {c}) returns the number of ways n can be changed with the set of coins {c}. In this tree the number of leaf nodes (10) is the answer. I have cheated here – the turquoise leaves are not degenerate cases – each has to be worked down to a case where the one would apply if the algorithm is to be universal, though we can simply add an extra rule when the set of coins being used consists of only 1p coins – that is 1 way of sorting. (But if the universal rule was applied then the number of computing steps increases greatly – as each of those 1p cases has to be worked down to a degenerate case).

Update (27 March 2011): There were actually a couple of mistakes in how the terminating cases were shown in the graphic, fixed now I hope.