MasonProtter · November 1, 2022 01:03 · Nov 1, 2022
diff --git a/PSL2Z.jl b/PSL2Z.jl
@@ -0,0 +1,76 @@
+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
No results found