## 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)
}

{
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))
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)
}
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))
}


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.

## An idea I saw on normblog

Norman Geras is a great man. He’s a social, not a computer, scientist, and this year I have been mainly reading computer books. Still, here’s an idea I have pinched off his website (he got it from here). Match the prompts to books you have read in the last year – in my case mainly computer books.

One time on holiday: Understanding the Linux Virtual Memory Manager

Weekends at my house are: Code Complete

My neighbour is: Programming Pearls

My boss is: The Mythical Man Month

My superhero secret identity is: The Annotated Turing

You wouldn’t like me when I’m angry because: Mussolini’s Italy

I’d win a gold medal in: The Little Schemer

I’d pay good money for: Structure and Interpretation of Computer Programs

If I were Prime Minister I would: Groovy in Action

When I don’t have good books, I: P, NP, and NP-Completeness

Loud talkers at the cinema should be: Albion’s Fatal Tree

## My first R program

Having used Groovy (which makes the scripting environment feel familiar) and some Scheme (via Structure and Interpretation of Computer Programs), R does feel completely alien, but it still feels like a steep learning curve.

But here’s my short script –

unpatched <- read.csv("~/unpatched.txt")
unpatchcons <- transform(unpatched, realm=realm*60 + reals)
attach(unpatchcons)
linelog<-lm(realm~size)
plot(size, realm, log="y")
abline(reg=linelog, untf=TRUE, col="blue",lty=3)
detach(unpatchcons)


And here’s the graph (of Linux kernel compile times) it generates – the blue line is obviously a very bad fit!

## Something more on Fibonacci numbers

My earlier blog about the Fibonacci series gets a lot of hits, so I thought I would write something more, as clearly there is interest.

Once again this is from “Structure and Interpretation of Computer Programs” (available for free here electronically):

Let $\psi = (1 - \sqrt{5})/2$ and $\phi = (1 + \sqrt{5}/2)$. Then $Fib(n) = (\phi^n - \psi^n)/\sqrt{5}$.

So $Fib(1) = (\frac{1}{2} +\frac{\sqrt{5}}{2} - \frac{1}{2} + \frac{\sqrt{5}}{2})/\sqrt{5} = 1$.

(1)$Fib(n) = (\phi^n - \psi^n)/\sqrt{5} = Fib(n - 1) + Fib(n - 2)$

(2) Then $Fib(n -1) + Fib( n - 2) = \frac{\phi^{n-1} - \psi^{n-1}}{\sqrt{5}} + \frac{\phi^{n-2} - \psi^{n-2}}{\sqrt{5}}$

(3) Simplifying (2) $Fib(n) = \frac{\phi^{n-1} - \psi^{n-1} + \phi^{n-2} - \psi^{n-2}}{\sqrt{5}}$

(4) Simplifying further – from (1): $\phi^n - \psi^n = \phi^{n-1} - \psi^{n-1} + \phi^{n-2} - \psi^{n-2}$

(5) $\phi^n - \psi^n = \phi^{n - 2}(\phi + 1) -\psi^{n- 2} - \psi^{n - 1}$

(6) We know that $\phi^2 = (\phi + 1)$ so we can restate (5) as $\phi^n - \psi^n = \phi^n - \psi^{n - 2} - \psi^{n - 1}$

So (7) $\psi^n = \psi^{n - 2} + \psi^{n - 1}$

But $\psi^2 = \psi + 1$ also and (7) can be restated as $\psi^n = \psi^{n-2}(\psi + 1) = \psi^{n-2}\psi^2 = \psi^n$

In other words, $Fib(n) = (\phi^n - \psi^n)/\sqrt{5}$.

## The lambda calculus and closures in Groovy

No sooner had I written about the lambda calculus and Structure and Interpretation of Computer Programs than I sat in a lecture on closures in Groovy and was presented with a structure like this (which multiplies two numbers, in this case 3 and 4):

Which immediately reminded me of one of Alonzo Church‘s formulations of lambdas – eg:

More to come…

## Coming up… the lambda calculus

Another thing that The Annotated Turing taught me is what all those lambdas that I have seen over the last 25 years were about, or at least it introduced me to what they were about.

So I have just ordered a copy of Structure and Interpretation of Computer Programs: and I can guess that some blog posts will eventually follow.

Oh no. Maybe I am turning into a LISP hacker – having finally reached the age that all these guys were when I first became aware of them (older, actually).