finally correct eigenvalues!

PSL.jl

@@ -1,5 +1,6 @@
 using Nemo
 using DelimitedFiles
+using LinearAlgebra
 
 include("src/nemo_utils.jl")
 
@@ -7,7 +8,7 @@ const PRECISION = 256
 
 function parse_eval(expr_str, value, var_name)
     ex = Meta.parse(expr_str)
-	svar = :($var_name)
+    svar = :($var_name)
     return @eval begin
         let $svar = $value
             $ex
@@ -16,9 +17,9 @@ function parse_eval(expr_str, value, var_name)
 end
 
 function read_eval(fname, var_name, value)
-   a = readdlm(fname, ',', String)
-   a .= replace.(a, '/'=>"//")
-   return parse_eval.(a, value, var_name)
+    a = readdlm(fname, ',', String)
+    a .= replace.(a, '/' => "//")
+    return parse_eval.(a, value, var_name)
 end
 
 function load_discrete_repr(i, q = 109; CC = AcbField(PRECISION))
@@ -27,13 +28,19 @@ function load_discrete_repr(i, q = 109; CC = AcbField(PRECISION))
 
     ra = read_eval(
         "data/Discrete reps PSL(2, $q)/discrete_rep_$(i)_a.txt",
-        :Z, ζ)
+        :Z,
+        ζ,
+    )
     a = matrix(CC, [CC(s) for s in ra[1:degree, 1:degree]])
 
     rb = read_eval(
         "data/Discrete reps PSL(2, $q)/discrete_rep_$(i)_b.txt",
-        :Z, ζ)
+        :Z,
+        ζ,
+    )
     b = matrix(CC, [CC(s) for s in rb[1:degree, 1:degree]])
+    @assert contains(det(a), 1)
+    @assert contains(det(b), 1)
 
     return a, b
 end
@@ -44,46 +51,61 @@ function load_principal_repr(i, q = 109; CC = AcbField(PRECISION))
 
     ra = read_eval(
         "data/Principal reps PSL(2, $q)/principal_rep_$(i)_a.txt",
-        :zz, ζ)
+        :zz,
+        ζ,
+    )
     a = matrix(CC, [CC(z) for z in ra[1:degree, 1:degree]])
 
     rb = read_eval(
         "data/Principal reps PSL(2, $q)/principal_rep_$(i)_b.txt",
-        :zz, ζ)
+        :zz,
+        ζ,
+    )
+    b = matrix(CC, [CC(z) for z in rb[1:degree, 1:degree]])
+    @assert contains(det(a), 1)
+    @assert contains(det(b), 1)
 
     return a, b
 end
 
-
-
 function safe_eigvals(m::acb_mat)
+    evs = eigvals(m)
+    all(isfinite.(evs)) && return evs
     CC = base_ring(m)
     X = matrix(CC, rand(CC, size(m)))
-    return eigvals(X * m * inv(X))
+    evs = eigvals(X * m * inv(X))
+    all(isfinite.(evs)) && return evs
+    throw(ArgumentError("Could not compute eigenvalues"))
 end
 
-for i in 0:27
-	try
-		a,b = load_principal_repr(i)
-		adjacency = sum([[a^i for i in 1:4]; [b^i for i in 1:4]])
-		@time evc = safe_eigvals(adjacency)
-		ev = sort(real.(first.(evc)), lt=<, rev=true)
-		@info "Principal Series Representation $i" ev[1:2]
-	catch ex
-		@error "Principal Series Representation $i failed"
-		ex isa InterruptException && throw(ex)
-	end
-end
+if !isinteractive()
+    for i = 0:27
+        try
+            a, b = load_principal_repr(i)
+            adjacency = sum(a^i for i = 1:4) + sum(b^i for i = 1:4)
+            @time ev = let evs = safe_eigvals(adjacency)
+                _count_multiplicites(evs)
+            end
 
-for i in 1:27
-	try
-		a,b = load_discrete_repr(i)
-		adjacency = sum([[a^i for i in 1:4]; [b^i for i in 1:4]])
-		@time evc = eigvals(adjacency)
-		ev = sort(real.(first.(evc)), lt=<, rev=true)
-		@info "Discrete Series Representation $i" ev[1:2]
-	catch ex
-		@error "Discrete Series Representation $i : failed"
-		ex isa InterruptException && rethrow(ex)
-	end
+            @info "Principal Series Representation $i" ev[1:2] ev[end]
+        catch ex
+            @error "Principal Series Representation $i failed" ex
+            ex isa InterruptException && throw(ex)
+        end
+    end
+
+    for i = 1:27
+        try
+            a, b = load_discrete_repr(i)
+            adjacency = sum(a^i for i = 1:4) + sum(b^i for i = 1:4)
+            @time ev = let evs = safe_eigvals(adjacency)
+                _count_multiplicites(evs)
+            end
+
+            @info "Discrete Series Representation $i" ev[1:2] ev[end]
+        catch ex
+            @error "Discrete Series Representation $i : failed" ex
+            ex isa InterruptException && rethrow(ex)
+        end
+    end
 end

src/nemo_utils.jl

@@ -1,3 +1,5 @@
+Base.hash(a::acb, h::UInt) = h
+
 import Base.reim
 reim(x::Nemo.acb) = reim(convert(ComplexF64, x))
 
@@ -153,5 +155,5 @@ function _count_multiplicites(evs)
         push!(λ_m, (evs[i], m))
         i += m
     end
-    return λ_m
+    return sort(λ_m, lt=(a,b)->(real(first(a))<real(first(b))), rev=true)
 end