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 , where is the radius of the bounding circle (so in our case)and 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 .

But the problem with that is, that with the exception of the start and end points, we don’t know what 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 ).

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.

Advertisement