PSL2Z.jl

using StaticArrays, Plots, Dictionaries

struct PSL2{T <: Real}
    a::T
    b::T
    c::T
    d::T
    function PSL2(a::T,b::T,c::T,d::T) where {T}
        (a*d - b*c ≈ 1) || error("Bad PSL₂($T) parameters, must have ad - bc = 1, got $((;a, b, c, d))")
        a >= 0 ? new{T}(a,b,c,d) : new{T}(-a,-b,-c,-d)
    end
end
StaticArrays.SMatrix((a,b,c,d)::PSL2) = SA[a b; c d]

PSL2(M::SMatrix{2, 2}) = PSL2(M[1,1], M[1,2], M[2,1], M[2, 2])

Base.one(::PSL2{T}) where {T} = PSL2(true, false, false, true)
Base.one(::Type{PSL2{T}}) where {T} = PSL2(one(T), zero(T), zero(T), zero(T))
Base.one(::Type{PSL2}) = PSL2(true, false, false, true)

Base.eltype(::Type{PSL2{T}}) where {T} = T
Base.iterate(g::PSL2, i=1) = i ∈ 1:4 ? (getfield(g, i), i + 1) : nothing
    
((a,b,c,d)::PSL2)(z) = (a*z + b)/(c*z + d)
Base.:(*)(g::PSL2, z::Number) = g(z)

Base.inv((a,b,c,d)::PSL2) = PSL2(d, -b, -c, a)

Base.:(*)(g::PSL2, h::PSL2) = PSL2(SMatrix(g) * SMatrix(h))
Base.:(∘)(g::PSL2, h::PSL2) = g * h

function Base.:(^)(g::PSL2, n::Integer) 
    gn = one(g)
    if n < 0
        return inv(g)^-n
    end
    while n > 0
        gn = gn * g
        n  = n - 1
    end
    gn
end

let
    T = PSL2(1, 1, 0, 1)
    ∂Δ = cis.(range(0, π, length=100))
    p = plot(ylims=(0, 1.5), xlims=(-1, 1), aspect_ratio=1.0, legend=nothing)

    𝒪 = let N = 4
        v = Indices(PSL2{Int}[])
        for a ∈ 1:N, b ∈ -N:N, c ∈ -N:N
            if rem((1 + b*c), a) == 0
                set!(v, PSL2(a, b, c, (1 + b*c)÷a))
            end
        end
        set!(v, PSL2(0, 1, -1, 0))
        v
    end
    pts = [x*im for x in range(0, 1, length=20)]
    append!(pts, [cis(θ) for θ ∈ range(π/2, π/2 + π/6, length=10)])
    append!(pts, [-1 + cis(θ) for θ ∈ range((π/2 - π/6), 0, length=20)])
    reflect(z) = -real(z) + imag(z)*im

    foreach(𝒪) do h
        foreach(-2:2) do n
            g = T^n ∘ h
            plot!(Shape(reim.(reflect.(g.(pts)))), color=:black)
            plot!(g.(∂Δ), color=:black)
        end
    end
    
    xticks!([0, 0.5], ["0", "1/2"])
    yticks!([1], ["1"])
    xlabel!(""); ylabel!("")
    p
end