Enumerating Bianchi circles Background: Katherine Stange describes Schmidt arrangements in "Visualising the arithmetic of imaginary quadratic fields", arXiv:1410.0417. Given an imaginary quadratic field $K$, we study the Bianchi group $\mathrm{PSL}_2(\mathcal{O}_K)$, which is the group of Möbius transformations with coefficients in the ring of integers of $K$. The image of $\mathbb R$ under a group element is called a $K$-Bianchi circle, and the set of $K$-Bianchi circles is called a Schmidt arrangement.
Here's an example from Stange's image gallery, taking $K=\mathbb Q(\sqrt{-7})$ and drawing all circles with curvature up to $30\sqrt 7$ that intersect the fundamental domain of $\mathcal{O}_K$:

My question: How can I recreate such images for arbitrary $K$, starting from rational integer arithmetic? What algorithms would I need to write C++ code from scratch to enumerate a Schmidt arrangement? (I prefer not to pull a Sage package off the shelf.)
It's a straightforward task to implement the basic arithmetic operations of $K$, $\mathcal{O}_K$, and $\mathrm{GL}_2(\mathcal{O}_K)$. Then, a brute-force approach is to generate lots of matrices in $\mathrm{GL}_2(\mathcal{O}_K)$, check if each has determinant $1$, and if so, draw the appropriate circle. There are many problems with this brute-force approach:


*

*It wastes most of its time inspecting matrices without determinant $1$, especially as we search for large-norm coefficients and high-curvature circles.

*It also wastes time on circles that are outside the bounds of the illustration.

*It revisits many group elements that generate the same circle. We should quotient out the stabilizer of $\mathbb R$, namely the modular group $\mathrm{PSL}_2(\mathbb Z)$.

*It's unclear how many matrices need to be tested before we can say that we've enumerated all the circles in a diagram like Stange's.


Whenever I reach for a more clever approach, I'm slowed down by the fact that $\mathcal{O}_K$ isn't necessarily a unique factorization domain, and even when it is, it isn't necessarily Euclidean. How can we enumerate something like $\mathrm{PSL}_2(\mathcal{O}_K)$ with reasonable efficiency when $K$ is so unruly? 
Edit: To explain my last remark, we're looking for matrices $\pmatrix{a&b\\c&d}$ with coefficients in $\mathcal{O}_K$ where $ad-bc=1$. Each row and column of such a matrix consists of two coprime numbers. How does one enumerate pairs of coprime numbers in a non-Euclidean ring? (Or is this the wrong sub-problem to tackle?)
 A: For the record, I've figured out how to cobble together enough inequalities to make the brute-force approach work. Here's a really quick outline. I won't prove that it works, but it does. I would still appreciate a more elegant solution!
I'll follow Stange's notation, so a general element of the Bianchi group is $\pmatrix{\alpha&\gamma\\\beta&\delta}$.


*

*Fix a maximum curvature $M$. We will draw all Bianchi circles having curvature bounded by $M$. Note that the curvature is $i(\beta\overline\delta-\delta\overline\beta)$, so we will begin the outer loops by enumerating $\beta$ and $\delta$.

*Enumerate $\beta\in\mathcal{O}_K$ bounded by $N(\beta)\leq M^2$, and such that $\beta\geq0$ in the dictionary order. (By the projective symmetry, we can choose the sign of $\beta$ without loss of generality.)

*Enumerate $\delta\in\mathcal{O}_K$ bounded by $N(\beta)\leq N(\delta)\leq \frac{4M^2}{3N(\beta)}$, and such that $|\Re(\delta/\beta)|\leq1/2$. (We only need for $\delta/\beta$ to cover one fundamental domain of the modular group.) If the curvature is indeed under $M$, then continue to the innermost loop:

*Enumerate $\alpha\in\mathcal{O}_K$ bounded by $N(\alpha)\leq\frac{(4+|\Delta|)N(\beta)}{16}$, where $\Delta$ is the discriminant of $K$. (This eliminates some, but not all, of the redundancy of generating circles that differ only by a translation in $\mathcal{O}_K$.)

*Solve $\alpha\delta-\beta\gamma=1$ for $\gamma$. If $\gamma\in\mathcal{O}_K$, then draw the appropriate circle and any of its translates by $\mathcal{O}_K$ in the viewing area.


The three nested loops seem like they would blow up as $M$ becomes large, but in practice it only takes a few seconds to get up to $M=600$, which is plenty for most illustrations, like this one:

A: Running out of time.  In general the generators of Bianchi groups may be hard to compute, as this paper by Swan suggests [1]
However, §4 of Stange's paper looks very useful for computing things.  


*

*4.1 $K$-bianchi circles only intersect at $K$ points.  

*4.2 Two $K$-bianchi circles may only intersect at angle $\theta$ iff $e^{i\theta} \in \mathbb{O}_K$.  Mainly this applies to $K = \mathbb{Q}(i), \mathbb{Q}(e^{2\pi i /3})$

*4.3 Let $a/b \in K$ be a fraction with $a,b \in \mathcal{O}_K$ and  $(a)+(b)=(1)$ "coprime" as ideals.  Then we have neighboring "farey fractions" 
$$ \left[\begin{array}{c} a \\ b \end{array}\right] , 
u\left[\begin{array}{c} a \\ b \end{array}\right] + k\tau\left[\begin{array}{c} c \\ d \end{array}\right] \in \mathbb{P}^1_K
$$
Here $u \in \mathcal{O}^\ast_K$ be a unit, which is $\{1,- 1\}$ if $K = \mathbb{Q}(\sqrt{-d})$ with $d > 3$.  Here $k \in \mathbb{Z}$ and $\tau$ is [not quite sure].  $\tau$ must be some input from the class group of $K$. is $\sqrt{-d}$ for $d \not \equiv 3(4)$ or $\frac{1 + \sqrt{-d}}{2}$ for $d \equiv 3(4)$.
The original paper by Asmus Schmidt is interesting reading.  He did not have computer so we definitely have a leg up on him.  See also Stange's most recent paper The Appolonian Structure of Bianchi Groups.  She includes this:

In this case we are reduced to computing coprime pairs of numbers $(a)+(b)=1$ which may or may not be easier than just finding solutions to $ad-bc=1$.
