SmallHyperbolic/adj_psl2_eigvals.jl

using RamanujanGraphs
using LinearAlgebra
using Arblib
using ArgParse

using Logging
using Dates

import RamanujanGraphs.Primes: isprime

include(joinpath(@__DIR__, "src", "eigen_utils.jl"))

function SL2p_gens(p::Integer)
    @assert isprime(p)
    if p == 31
        a, b = let
            a = SL₂{p}([8 14; 4 11])
            b = SL₂{p}([23 0; 14 27])
            @assert isone(a^10)
            @assert isone(b^10)

            a, b
        end
    elseif p == 41
        a, b = let
            a = SL₂{p}([0 28; 19 35])
            b = SL₂{p}([38 27; 2 9])
            @assert isone(a^10)
            @assert isone(b^10)

            a, b
        end
    elseif 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

    return a, b
end

function adjacency(ϱ, a, b; prec = 256)
    order_a = findfirst(i -> isone(a^i), 1:100)
    order_b = findfirst(i -> isone(b^i), 1:100)
    @assert !isnothing(order_a) && order_a > 1
    @assert !isnothing(order_b) && order_b > 1

    k = order_a - 1 + order_b - 1

    A = AcbMatrix(ϱ(a), prec = prec)
    B = AcbMatrix(ϱ(b), prec = prec)
    res = sum(A^i for i = 1:order_a-1) + sum(B^i for i = 1:order_b-1)
    #return Arblib.scalar_div!(res, res, k)
    return res
end

function parse_our_args()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "-p"
        help = "the prime p for which to use PSL(2,p)"
        arg_type = Int
        required = true
        "-a"
        help = "generator a (optional)"
        "-b"
        help = "generator b (optional)"
        "--ab"
        help = "array of generators a and b (optional)"
        "--precision"
        help = "set the precision of computations"
        arg_type = Int
        default = 128
    end

    result = parse_args(s)
    for key in ["a", "b", "ab"]
        val = get(result, key, "")
        if val != nothing
            result[key] = eval(Meta.parse(val))
        else
            delete!(result, key)
        end
    end
    val = get(result, "ab", "")
    if val != ""
        result["a"] = val[1]
        result["b"] = val[2]
    end
    result
end

parsed_args = parse_our_args()

const p = let p = parsed_args["p"]
    isprime(p) ||
        @error "You need to provide a prime, ex: `julia adj_psl2_eigvals.jl -p 31`"
    p
end

const PRECISION = parsed_args["precision"]
const LOGFILE = joinpath("log", "SL(2,$p)_eigvals_$(now()).log")

open(LOGFILE, "w") do io
    @info "Logging into $LOGFILE"
    with_logger(SimpleLogger(io)) do

        @info "Arguments:" args = parsed_args

        a, b = SL2p_gens(p)
        a = SL₂{p}(get(parsed_args, "a", a))
        b = SL₂{p}(get(parsed_args, "b", b))
        @info "Generators" a b

        Borel_cosets = let p = p, (a, b) = (a, b)
            SL2p, sizes =
                RamanujanGraphs.generate_balls([a, b, inv(a), inv(b)], radius = 21)
            @assert sizes[end] == RamanujanGraphs.order(SL₂{p})
            RamanujanGraphs.CosetDecomposition(SL2p, Borel(SL₂{p}))
        end

        all_large_evs = Arb[]
        let α = RamanujanGraphs.generator(RamanujanGraphs.GF{p}(0))

            for j = 0:(p-1)÷4
                h = PrincipalRepr(
                    α => unit_root((p - 1) ÷ 2, j, prec = PRECISION),
                    Borel_cosets,
                )

                @time adj = adjacency(h, a, b, prec = PRECISION)

                try
                    @time evs = let evs = safe_eigvals(adj)
                        count_multiplicites(evs)
                    end
                    append!(all_large_evs, [real(first(x)) for x in evs[1:2]])

                    @info "Principal Series Representation $j" evs[1:2] evs[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(α))

            if p % 4 == 1
                ub = (p - 1) ÷ 4
                ζ = unit_root((p + 1) ÷ 2, 1, prec = PRECISION)
            else # p % 4 == 3
                ub = (p + 1) ÷ 4
                ζ = unit_root((p + 1), 1, prec = PRECISION)
            end

            for k = 1:ub

                h = DiscreteRepr(
                    RamanujanGraphs.GF{p}(1) => unit_root(p, prec = PRECISION),
                    β => ζ^k,
                )

                @time adj = adjacency(h, a, b, prec = PRECISION)

                try
                    @time evs = let evs = safe_eigvals(adj)
                        count_multiplicites(evs)
                    end
                    append!(all_large_evs, [real(first(x)) for x in evs[1:2]])

                    @info "Discrete Series Representation $k" evs[1:2] evs[end]
                catch ex
                    @error "Discrete Series Representation $k : failed" ex
                    ex isa InterruptException && rethrow(ex)
                end
            end
        end
        all_large_evs = sort(all_large_evs, rev = true)
        λ = all_large_evs[2]
        ε = (λ - 3) / 5
        α = acos(ε)
        α_deg = (α / pi) * 180
        @info "Certified values:" λ ε α α_deg
    end # with_logger
end # open(logfile)

#
# using RamanujanGraphs.LightGraphs
# using Arpack
#
# Γ, eigenvalues = let p = 109,
#     a = PSL₂{p}([  0  1; 108  11]),
#     b = PSL₂{p}([ 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(RamanujanGraphs.order(PSL₂{p}), S)
#
#     @assert all(LightGraphs.degree(Γ,i) == length(S) for i in vertices(Γ))
#     @assert LightGraphs.nv(Γ) == RamanujanGraphs.order(PSL₂{p})
#     A = adjacency_matrix(Γ)
#     @time eigenvalues, _ = eigs(A, nev=5)
#     @show Γ eigenvalues
#     Γ, eigenvalues
# end
#
# let p = 131,
#     a = PSL₂{p}([-58 -24; -58 46]),
#     b = PSL₂{p}([0 -3; 44 -12])
#
#     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(RamanujanGraphs.order(PSL₂{p}), S)
#
#     @assert all(LightGraphs.degree(Γ,i) == length(S) for i in vertices(Γ))
#     @assert LightGraphs.nv(Γ) == RamanujanGraphs.order(PSL₂{p})
#     A = adjacency_matrix(Γ)
#     @time eigenvalues, _ = eigs(A, nev=5)
#     @show Γ eigenvalues
#     Γ, eigenvalues
# end