From 44339f91f0ff7a416386bf31c348e31e63074b4a Mon Sep 17 00:00:00 2001
From: kalmarek <kalmar@amu.edu.pl>
Date: Tue, 26 May 2020 00:10:24 +0200
Subject: [PATCH] add eigenvalues computations

---
 adj_psl2_eigvals.jl | 127 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 127 insertions(+)
 create mode 100644 adj_psl2_eigvals.jl

diff --git a/adj_psl2_eigvals.jl b/adj_psl2_eigvals.jl
new file mode 100644
index 0000000..cfd6749
--- /dev/null
+++ b/adj_psl2_eigvals.jl
@@ -0,0 +1,127 @@
+using RamanujanGraphs
+using LinearAlgebra
+using Nemo
+
+include("src/nemo_utils.jl")
+
+const p = try
+    @assert length(ARGS) == 2 && ARGS[1] == "-p"
+    p = parse(Int, ARGS[2])
+    RamanujanGraphs.Primes.isprime(p)
+    @assert p % 4 == 1
+catch ex
+    @error "You need to provide a prime `-p` which is congruent to 1 mod 4."
+    rethrow(ex)
+end
+
+const CC = AcbField(256)
+
+SL2p = let
+    if p == 109
+        a,b = let
+            a = SL₂{p}([  0  1; 108 11])
+            b = SL₂{p}([57  2;   52 42])
+            @assert isone(a^10)
+            @assert isone(b^10)
+
+            a, b
+        end
+    elseif p == 131
+        a,b = let
+            a = SL₂{p}([-58 -24; -58 46])
+            b = SL₂{p}([0 -3; 44 -12])
+            @assert isone(a^10)
+            @assert isone(b^10)
+
+            a, b
+        end
+    else
+        @warn "no special set of generators for prime $p"
+        a,b = let
+            a = SL₂{p}(1, 0, 1, 1)
+            b = SL₂{p}(1, 1, 0, 1)
+            a, b
+        end
+    end
+
+    E, sizes = RamanujanGraphs.generate_balls([a,b, inv(a), inv(b)], radius=21);
+    @assert sizes[end] == RamanujanGraphs.order(SL₂{p})
+    E
+end
+
+let Borel_cosets = Bcosets = RamanujanGraphs.CosetDecomposition(SL2p, Borel(SL₂{p})),
+    α = RamanujanGraphs.generator(RamanujanGraphs.GF{p}(0))
+
+    for j in 0:(p-1)÷4
+        try
+            h = PrincipalRepr(
+                α => root_of_unity(CC, (p-1)÷2, j),
+                Borel_cosets)
+
+            @time adjacency = let
+                A = matrix(CC, h(SL2p[2]))
+                B = matrix(CC, h(SL2p[3]))
+                sum(A^i for i in 1:4) + sum(B^i for i in 1:4)
+            end
+
+            @time ev = let evs = safe_eigvals(adjacency)
+                _count_multiplicites(evs)
+            end
+            if length(ev) == 1
+                @info "Principal Series Representation $j" ev[1]
+            else
+                @info "Principal Series Representation $j" ev[1:2] ev[end]
+            end
+        catch ex
+            @error "Principal Series Representation $j failed" ex
+            ex isa InterruptException && rethrow(ex)
+        end
+    end
+end
+
+let α = RamanujanGraphs.generator(RamanujanGraphs.GF{p}(0)),
+    β = RamanujanGraphs.generator_min(QuadraticExt(α))
+    for j = 1:(p-1)÷4
+        try
+            h = DiscreteRepr(
+                RamanujanGraphs.GF{p}(1) => root_of_unity(CC, p),
+                β => root_of_unity(CC, (p+1)÷2, j*(p-1)÷4))
+
+            @time adjacency = let
+                A = matrix(CC, h(SL2p[2]))
+                B = matrix(CC, h(SL2p[3]))
+                sum(A^i for i in 1:4) + sum(B^i for i in 1:4)
+            end
+
+            @time ev = let evs = safe_eigvals(adjacency)
+                _count_multiplicites(evs)
+            end
+
+            @info "Discrete Series Representation $j" ev[1:2] ev[end]
+        catch ex
+            @error "Discrete Series Representation $j : failed" ex
+            ex isa InterruptException && rethrow(ex)
+        end
+    end
+end
+
+
+# using RamanujanGraphs.LightGraphs
+# using Arpack
+#
+# Γ, eigenvalues = let q = 109
+#     a = RamanujanGraphs.PSL₂{q}([  0  1
+#     108 11])
+#     b = RamanujanGraphs.PSL₂{q}([57  2
+#     52 42])
+#
+#     S = unique([[a^i for i in 1:4]; [b^i for i in 1:4]])
+#
+#     @info "Generating set S of $(eltype(S))" S
+#     @time Γ, verts, vlabels, elabels = RamanujanGraphs.cayley_graph((q^3 - q)÷2, S)
+#     @assert all(LightGraphs.degree(Γ,i) == length(S) for i in vertices(Γ))
+#     A = adjacency_matrix(Γ)
+#     @time eigenvalues, _ = eigs(A, nev=5)
+#     @show Γ eigenvalues
+#     Γ, eigenvalues
+# end