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?)

