implement the "first" ramanujan graph

[Krastanov]

Paper: http://math1.math.huji.ac.il/~alexlub/PAPERS/ramanujan%20graphs/ramanujanGraphs.pdf

Lubotzky, Alexander, Ralph Phillips, and Peter Sarnak. "Ramanujan graphs." Combinatorica 8.3 (1988): 261-277.

[Krastanov]

a partial implementation (for the generators) by Liam McCarthy @ UMass

using Nemo
using Oscar
using LinearAlgebra
using Random
using Graphs
using Multigraphs
using ProgressMeter
using QuantumClifford
using Oscar:coefficients
using Graphs:add_edge!,nv,ne,neighbors,Graphs

p = 5
i = 1
𝔽q, punit = FiniteField(p,i)
GL₂q = general_linear_group(2,𝔽q)
# CGL₂q, Cₘₒᵣₚₕ = Oscar.center(GL₂q) # seems to take time that scales with the size of GL₂q even though it is either 1 or 2 element group.
# PGL₂q, Pₘₒᵣₚₕ = quo(GL₂q,CGL₂q) #G is  PGL₂q
# elements = collect(PGL₂q)
# Pₘₒᵣₚₕ(elements[1])
println("DONE")

function GetValidPrime(minVal, maxVal)
    minK = floor((minVal - 1)/4)
    maxK =  ceil((maxVal - 1)/4)
    while true
        k = rand(minK : maxK)
        if is_prime(trunc(Int, 4*k + 1))
            return 4*k+1
        end
    end
end

function Find_i_and_q(minVal, maxVal, p)  
    i = -1
    q = -1
    while i == -1
        q = GetValidPrime(minVal, maxVal)
        while q == p   
            q = GetValidPrime(minVal, maxVal)
        end
        i = Find_i(q)
    end
    return (i, q)
end

function SolveFourSquares(p)
    solutions = []
    for w in 0 : floor(sqrt(p))
        for x in 0 : floor(sqrt(p))
            for y in 0 : floor(sqrt(p))
                for z in 0 : floor(sqrt(p))
                    if (w * w + x * x + y * y + z * z) == p
                                                                push!(solutions, ( w,  x,  y,  z)) 
                        w != 0 && x != 0 && y != 0 && z != 0 && push!(solutions, (-w, -x, -y, -z)) 
                                  x != 0 && y != 0 && z != 0 && push!(solutions, ( w, -x, -y, -z)) 
                        w != 0           && y != 0 && z != 0 && push!(solutions, (-w,  x, -y, -z)) 
                        w != 0 && x != 0           && z != 0 && push!(solutions, (-w, -x,  y, -z)) 
                        w != 0 && x != 0 && y != 0           && push!(solutions, (-w, -x, -y,  z)) 
                                            y != 0 && z != 0 && push!(solutions, ( w,  x, -y, -z))
                        w != 0                     && z != 0 && push!(solutions, (-w,  x,  y, -z))
                        w != 0 && x != 0                     && push!(solutions, (-w, -x,  y,  z))
                                  x != 0 && y != 0           && push!(solutions, ( w, -x, -y,  z))
                                  x != 0           && z != 0 && push!(solutions, ( w, -x,  y, -z))
                        w != 0           && y != 0           && push!(solutions, (-w,  x, -y,  z))
                        w != 0                               && push!(solutions, (-w,  x,  y,  z)) 
                                  x != 0                     && push!(solutions, ( w, -x,  y,  z)) 
                                            y != 0           && push!(solutions, ( w,  x, -y,  z)) 
                                                      z != 0 && push!(solutions, ( w,  x,  y, -z)) 
                    end
                end
            end
        end
    end
    return solutions
end

function ProcessSolutions(solutions)
    final = []
    for solution in solutions
        if isodd(solution[1]) && solution[1] > 0 && iseven(solution[2]) && iseven(solution[3]) && iseven(solution[4])
            push!(final, solution)
        end
    end
    return final
end

function Create_Generator(solutions, i)
    matrices = []
    for soln in solutions
        a =  soln[1] + i * soln[2]
        b =  soln[3] + i * soln[4]
        c = -soln[3] + i * soln[4]
        d =  soln[1] - i * soln[2]
        push!(matrices, [[a b]; [c d]])
    end
    return matrices
end

p = GetValidPrime(0, range_max)
tup = Find_i_and_q(0, range_max, p)
i = tup[1]
q = tup[2]
sfs = SolveFourSquares(p)
solutions = ProcessSolutions(sfs) 
group = Create_Generator(solutions, i)

function Find_i(q)
    r = 1

    while !isinteger(sqrt(r * q - 1))
        print(r * q)
        r = r + 1
        if r > 100
            return -1
        end
    end
    return sqrt(r * q - 1)
end

[Krastanov]

suggested readings for improvements "Randomized Algorithms in Number Theory" by Rabin and Shallit