How to draw the picture of a closed horocycle in $\text{SL}_2(\mathbb{Z}) \backslash \mathbb{H}$?

Question

I am trying to draw the image of the path $[0,1] + \frac{1}{7}i \in \mathbb{H}$ under the image of $\text{SL}_2(\mathbb{Z})$.

This is the example of a horocycle - the image of a horocycle under the $\text{SL}_2(\mathbb{Z})$ group action. Here is the code I used to generate the points (in numPy):

N = 7

x = np.arange(0,1,0.0001)
z = 1.0j/N + x


for t in range(10):
    z = ( (z.real % 1) - 0.5 ) + z.imag*1j
    z = (np.abs(z) < 1)*(-1.0/z) + (np.abs(z) > 1 )*z

This seems pretty consistent with the algorithm given by William Stein.

Do the following until $z$ is in $\mathcal{F}$:

1. Replace $z$ by $z+n$ where $n\in \mathbb{Z}$ is an integer such that $|\text{Re}(z+n)| \leq \frac{1}{2}$.

2. If $|z|<1$, replace $z$ by $-1/z$.

Answer 1

What to expect

I assume that since you write about a closed horocycle in your question title, you actually mean the horocycle $\mathbb R+\tfrac17i$ instead of the horocycle arc $[0,1]+\tfrac17i$ you use in your description. On the other hand, since translations by $1$ turn that arc into the complete horocycle, it probably doesn't make much of a difference either way.

I started with a quick experiment using Cinderella. There I declared three transformations: an inversion in the unit circle (which thanks to the symmetry of this problem works just as well as the $z\mapsto-\tfrac1z$ you have, even though it's in fact $z\mapsto\tfrac1{\bar z}$), a translation by one unit to the right and its inverse. Then I defined a symmetry group consisting of these, and applied it to the line $y=\tfrac17$. The result gave me an idea what to expect:

I've also created a better version, with more iterations, better highlighting of fundamental domain, and generally more hand-tuning (click for bigger version):

So how do you draw this? Naively speaking, you just draw those circles that intersect the fundamental domain. There are $7$ big circles and $2$ smaller ones. The big ones have radius $\tfrac72$ and center $x$ coordinates $\{-3,-2,-1,0,1,2,3\}$. The smaller ones have radius $\tfrac78$ and center $x$ coordinates $\pm\tfrac12$. The $y$ coordinate of the center is always equal to the radius, since a horocycle touches the real axis. Getting the infinitely many smaller circles in a picture like the ones above is a lot more work but not needed if you focus on one fundamental domain only.

Debugging

So where did your code go wrong? Take another look at this part:

( (z.real % 1) - 0.5 )

This should better be

(((z.real + 0.5) % 1) - 0.5)

or some such, so that it is idempotent for values already in the range $(-\tfrac12,\tfrac12)$. You could also write that line as

z = z - np.round(z.real)

I could find the error by just plugging $z=\tfrac17i$ into your algorithm, which I knew should end up at $7i$ but got some real component thanks to this bug.

How I created my image

The picture I created above came out of some Sage computation. The core idea is to represent a circle as a homogeneous integer coordinate vector, namely as $(x,y,x^2+y^2-r^2,1)^T$ or some multiple thereof. The elementary transformations would be

$$ S=\begin{pmatrix}-1&0&0&0\\0&1&0&0\\0&0&0&1\\0&0&1&0\end{pmatrix}\qquad T=\begin{pmatrix}1&0&0&1\\0&1&0&0\\2&0&1&1\\0&0&0&1\end{pmatrix}\qquad T^{-1}=\begin{pmatrix}1&0&0&-1\\0&1&0&0\\-2&0&1&1\\0&0&0&1\end{pmatrix} $$

The line $\mathbb R+\tfrac17i$ can be written as $(0,7,2,0)^T$; it's image under operation $S$ is the circle $(0,7,0,2)$, a circle with center at $(0,\tfrac72)$. So starting from this, it is possible to enumerate circles like this:

from __future__ import division
import math, collections
queue = collections.deque()
seen = set()

def see(v):
  if v not in seen:
    seen.add(v)
    queue.append(v)

see((0, 7, 2, 0))
for t in range(10000):  # Limit chosen arbitrarily
  v = queue.popleft()
  see((-v[0], v[1], v[3], v[2]))
  see((v[0]+v[3], v[1], 2*v[0]+v[2]+v[3], v[3]))
  see((v[0]-v[3], v[1], -2*v[0]+v[2]+v[3], v[3]))

def dehom(v):
  x = v[0]/v[3]
  y = v[1]/v[3]
  r = math.sqrt(x*x + y*y - v[2]/v[3])
  return (x, y, r)

circles = [dehom(v) for v in seen if v[3]]
circles.sort(key = lambda xyr: (-xyr[2], xyr[0])) # sort by r desc, then x asc
for x, y, r in circles:
  if r > 0.005 and x+r > -2 and x-r < 2:  # Limit roughly to picture scope
    print("circle center {}, {} radius {}".format(x, y, r))

You can easily replace the print statement with an SVG write, or some draw command, or some such. Using integer arithmetic avoids rounding errors when checking for presence in the set. Usually you'd want some step somewhere to normalize the representatives, in order to avoid having different multiples of the same vector in the set. But since all three operations leave the second coordinate unmodified, we can only ever get one possible representative, namely the one with second coordinate equal to $7$.