<|start_filename|>extras/projective.jl<|end_filename|>
using Mods, Primes, SimpleGraphs, LinearAlgebraX

"""
`generate_points(p::Int)` generates the homogeneous coordinates of the 
points in a finite projective plane of prime order `p`.
"""
function generate_points(p::Integer)::Vector{HVector}
    if ~isprime(p)
        error("Modulus $p is not prime")
    end

    n = p * p + p + 1
    A = Array{HVector}(undef, n)
    count = 1

    # case c=1
    for a = 0:p-1
        for b = 0:p-1
            x = Mod{p}.([a, b, 1])
            A[count] = HVector(x)
            count += 1
        end
    end

    # case c=0, b=1
    for a = 0:p-1
        x = Mod{p}.([a, 1, 0])
        A[count] = HVector(x)
        count += 1
    end

    A[count] = HVector(Mod{p}(1), 0, 0)

    return A
end

"""
`incidence_matrix(p::Int)` creates the point/line incidence matrix 
of a finite projective plane of prime order `p`.
"""
function incidence_matrix(p::Int)::Matrix
    pts = generate_points(p)
    n = p*p + p + 1
    A = zeros(Int,n,n)
    for u = 1:n
        P = pts[u]
        for v = 1:n
            Q = pts[v]
            if dot(P,Q)==0
                A[u,v] = 1
            end
        end
    end
    return A
end

"""
`incidence_graph(p::Int)` creates the bipartite point-line 
incidence graph of a finite projective plane of prime order `p`.
"""
function incidence_graph(p::Int)::SimpleGraph{Int}
    M = incidence_matrix(p)
    r,c = size(M)
    A = [ zeros(r,c)  M; copy(M')  zeros(r,c) ]
    G = SimpleGraph(A)
end 

function incidence_hypergraph(p::Int)::SimpleHypergraph
    VV = generate_points(p)
    H = SimpleHypergraph{HVector{Mod{p}}}()
    for v in VV
        add!(H,v)
    end 

    for v in VV 
        S = [ w for w in VV if dot(v,w)==0 ]
        add!(H,Set(S))
    end 

    return H 
end 

<|start_filename|>src/rrefx.jl<|end_filename|>
export rrefx

"""
`rref(A)`  returns the row-reduced echelon form of the exact 
matrix `A`.
"""
function rrefx(A::AbstractMatrix{T}) where {T<:IntegerX}
    AA = big.(A // 1)
    return rrefx(AA)
end

function rrefx(A::AbstractMatrix{T}) where {T}
    AA = copy(A)
    rrefx!(AA)
    return AA
end


function _pivot!(A::AbstractMatrix{T}, i::Int, j::Int) where {T}
    r, c = size(A)
    s = A[i, j]
    row_scale!(A, i, invx(s))

    for k = 1:r
        if k != i
            s = -A[k, j]
            row_add_mult!(A, i, s, k)
        end
    end
end



function rrefx!(A::AbstractMatrix{T}) where {T}
    rn, nc = size(A)

    r = 1
    for j = 1:nc
        # see if column j's first one is at level r or below
        col = A[r:end, j]
        if all(col .== 0)
            continue
        end
        k = findfirst(col .!= 0)
        row_swap!(A, r, k + r - 1)
        _pivot!(A, r, j)
        r += 1
    end


end


<|start_filename|>src/eye.jl<|end_filename|>
export eye

"""
`eye(T,n)` returns an `n`-by-`n` identity matrix whose entries
are numbers of type `T` (`Float64` by default).
"""
function eye(T, n)
    A = zeros(T, n, n)
    for j = 1:n
        @inbounds A[j, j] = one(T)
    end
    return A
end

eye(n) = eye(Float64, n)


<|start_filename|>src/invx.jl<|end_filename|>
# matrix (and other) inversion

export invx

"""
`invx(A)` for a matrix `A` gives an exact matrix inverse.
"""
function invx(A::AbstractMatrix{T}) where {T}
    r, c = size(A)
    @assert r == c "Matrix must be square"

    if r == 0
        return A
    end

    B = eye(T, r)

    X = [A B]
    X = rrefx(X)

    d = [X[i, i] for i = 1:r]
    if all(d .== 1)
        return X[:, r+1:end]
    end
    @error "Matrix is not invertible"
end


function invx(x::T) where {T<:IntegerX}
    return 1 // x
end

invx(x) = 1 / x


<|start_filename|>extras/hilbert.jl<|end_filename|>
"""
`hilbert(n)` creates an `n`-by-`n` Hilbert matrix.

## Example
```
julia> hilbert(4)
4×4 Array{Rational{Int64},2}:
 1//1  1//2  1//3  1//4
 1//2  1//3  1//4  1//5
 1//3  1//4  1//5  1//6
 1//4  1//5  1//6  1//7
```
"""
function hilbert(n::Int)
    A = zeros(Rational{Int}, n, n)
    for i = 1:n
        for j = 1:n
            A[i, j] = 1 // (i + j - 1)
        end
    end
    return A
end



<|start_filename|>test/runtests.jl<|end_filename|>
using Test
using LinearAlgebra, LinearAlgebraX, Mods


@testset "Basics" begin
    A = big.(rand(Int, 10, 10) .% 100)

    include("../extras/hilbert.jl")

    H = hilbert(12)
    @test rankx(H) == 12

    @test detx(A) == det(Rational.(A))
    @test detx(A) == detx(A')
    @test A * invx(A) == eye(BigInt, 10)

    A = ones(Int, 3, 5)
    N = nullspacex(A)
    @test all(0 .== A * N)
end

@testset "Polynomials" begin
    using SimplePolynomials

    A = triu(ones(Int, 5, 5))
    p = char_poly(A)
    x = getx()
    @test p == (x - 1)^5

    A = ones(Int, 5, 5)
    p = char_poly(A)
    @test p == x^4 * (x - 5)
end

@testset "Modular" begin
    T = Mod{17}
    A = rand(T, 6, 6)
    p = char_poly(A)
    @test detx(A) == p(0)


    A = rand(Int, 20, 20) .% 100
    @test rankx(A) <= 20
end


@testset "Homogeneous" begin
    v = HVector(2, 3, 2)
    @test v[2] == 3 // 2
    @test dot(v, v) == 1
    w = HVector(20, 30, 20)
    @test v == w

    A = rand(Int, 3, 3) .% 20
    w = A * v
    @test w == HVector(A * v.data)
end

@testset "Permanent" begin
    A = ones(Int, 6, 6)
    @test permanent(A) == factorial(6)

    A = A - eye(Int, 6)
    @test permanent(A) == 265  # number of derangements of [6]

    A = triu(ones(5, 5))
    @test permanent(A) == 1
end


<|start_filename|>src/perm.jl<|end_filename|>
export permanent
"""
    permanent(A)

Compute the permanent of the square matrix `A`. Works best if `A`
is relatively sparse. For a dense matrix, this will be slow.
"""
function permanent(A::Matrix{T}) where {T<:Number}
    r, c = size(A)
    if r != c
        error("Matrix must be square")
    end
    if r == 0
        return zero(T)
    end

    allow = [findall((x) -> x != 0, A[i, :]) for i = 1:r]
    PG = PermGen(allow)

    try
        first(PG)  # make sure PG is not empty
        return sum(prod(A[i, p[i]] for i = 1:r) for p in PG)
    catch
        return zero(T)
    end
end


<|start_filename|>src/char_poly.jl<|end_filename|>
export char_poly

"""
`char_poly(A)` returns the characteristic polynomial of
the exact matrix `A`.
"""
function char_poly(A::AbstractMatrix{T}) where {T<:TypeX}
    r, c = size(A)
    @assert r == c "Matrix must be square"

    xI = zeros(SimplePolynomial, r, r)
    x = getx(T)
    for i = 1:r
        xI[i, i] = x
    end

    f = detx(xI - A)
    p = integerize(numerator(f))
end


<|start_filename|>src/homogeneous.jl<|end_filename|>
import Base: show, eltype, length, size, getindex, ==, (*), iszero
import Base: Vector, Matrix, adjoint, hash

import LinearAlgebra: dot
export HVector, HMatrix, dot, hash


"""
An `HVector` is a homogeneous vector; two `HVector`s are equal 
iff they are nonzero multiples of each other.

`HVector(a,b,c,...)` or `HVector([a,b,c,...])` creates a new `HVector `
with coordinates `(a,b,c,...)`.
"""
struct HVector{T}
    data::Vector{T}
    function HVector(dat::Vector{T}) where {T}
        data = dat .// one(T)
        S = typeof(sum(data))
        new{S}(_canonical(data))
    end

end
HVector(x...) = HVector(collect(x))
HVector(x::HVector) = HVector(x.data)

"""
An `HMatrix` is a homogenous matrix; two of these are equal if they are
nonzero multiples of each other. 

Given a two-dimensional array `A`, `HMatrix(A)` creates a new homogeneous 
matrix. 
"""
struct HMatrix{T}
    data::Matrix{T}
    function HMatrix(dat::AbstractMatrix{T}) where {T}
        data = dat .// one(T)
        S = typeof(sum(data))
        M = new{S}(_canonical(data))
        return M
    end
end
HMatrix(x::HMatrix) = HMatrix(x.data)


function _canonical(x)
    if length(x) == 0 || iszero(x)
        return x
    end
    idx = findlast(x .!= 0)
    a = x[idx]
    return x .// a
end

HObject = Union{HVector,HMatrix}


eltype(x::HObject) = eltype(x.data)
length(x::HObject) = length(x.data)
size(x::HObject) = size(x.data)
getindex(x::HObject, k...) = x.data[k...]
(==)(x::HObject, y::HObject) = x.data == y.data
iszero(x::HObject) = iszero(x.data)


"""
The `dot` produce of two `HVector`s is either `0` 
(if they are orthogonal) or `1` (otherwise).
"""
function dot(x::HVector, y::HVector)::Int
    d = dot(x.data, y.data)
    if d == 0
        return 0
    end
    return 1
end


function show(io::IO, x::HVector)
    print(io, "HVector(")
    n = length(x.data)
    for j = 1:n-1
        print(io, "$(x.data[j]), ")
    end
    print(io, x.data[end], ")")
end

function show(io::IO, x::HMatrix)
    print(io, "HMatrix: ")
    print(io, x.data)
end

function (*)(A::AbstractArray{T,2}, v::HVector) where {T}
    Av = A * v.data
    return HVector(Av)
end

function (*)(A::AbstractArray{T,2}, B::HMatrix) where {T}
    AB = A * B.data
    return HMatrix(AB)
end

(*)(A::HMatrix, v::HVector) = HVector(A.data * v.data)
(*)(A::HMatrix, B::HMatrix) = HMatrix(A.data * B.data)

adjoint(A::HMatrix) = HMatrix(adjoint(A.data))

hash(x::HObject, h::UInt = UInt(0)) = hash(x.data, h)

Vector(x::HVector) = copy(x.data)
Matrix(x::HMatrix) = copy(x.data)


<|start_filename|>src/detx.jl<|end_filename|>
export detx

"""
`detx(A::AbstractMatrix{T})` is an exact determinant of the matrix.
Here `T` can be any kind of integer, rational, `Mod`, or `GF2`.

I hope to expand this to `Polynomial`s.
"""
function detx(A::AbstractMatrix{T}) where {T<:IntegerX}
    # @info "Using IntegerX detx{$T}"
    try
        return T(det(A // 1))
    catch
        A = big.(A)
        return detx(A)
    end
end

function detx(A::AbstractMatrix{T}) where {T<:RationalX}
    # @info "Using RationalX detx{$T}"
    try
        return det(A)
    catch
        return det(big.(A))
    end
end

function detx(A::AbstractMatrix{T}) where {T}
    # @info "Using detx! on matrix of type $T"
    r, c = size(A)
    B = Matrix{Any}(undef, r, c)
    for i = 1:r
        for j = 1:c
            B[i, j] = A[i, j]
        end
    end
    # B = copy(A)
    try
        return detx!(B)
    catch
        @warn "Using cofactor expansion to calculate determinant; may be very slow."
        return cofactor_det(B)  # if all else fails!
    end
end


# detx! is the work behind det when we can't use Julia's det.
# it can modify the matrix so this is not exposed to the general public.

function detx!(A::AbstractMatrix{T}) where {T}
    r, c = size(A)
    @assert r == c "Matrix must be square"


    if r == 0
        return 1
    end

    if r == 1
        return A[1, 1]
    end

    if r == 2
        return A[1, 1] * A[2, 2] - A[1, 2] * A[2, 1]
    end

    col = A[:, 1]  # first column
    if all(col .== 0)
        return 0
    end

    k = findfirst(col .!= 0)  # we pivot here

    sign_factor = 1
    if k > 1
        sign_factor = -sign_factor
        row_swap!(A, 1, k)
    end

    factor = A[1, 1] * sign_factor  # multiply by this at the end

    row_scale!(A, 1, 1 // A[1, 1])

    for i = 2:r
        row_add_mult!(A, 1, -A[i, 1], i)
    end


    return detx!(A[2:end, 2:end]) * factor

end


<|start_filename|>src/row_ops.jl<|end_filename|>
# Elementary row operations

export row_swap!, row_scale!, row_add_mult!


"""
`row_swap!(A,i,j)` swaps rows `i` and `j` in the matrix `A`.
"""
function row_swap!(A::AbstractMatrix{T}, i::Int, j::Int) where {T}
    r, c = size(A)
    @assert (1 <= i <= r) && (1 <= j <= r) "Row index out of bounds"
    if i == j
        return
    end

    save_row = [A[i, k] for k = 1:c]

    for k = 1:c
        @inbounds A[i, k] = A[j, k]
    end

    for k = 1:c
        @inbounds A[j, k] = save_row[k]
    end
end

"""
`row_scale!(A,i,s)` multiplies all entries in row `i`
of `A` by `s`.
"""
function row_scale!(A::AbstractMatrix{T}, i::Int, s) where {T}
    r, c = size(A)
    @assert 1 <= i <= r "Row index out of bounds"

    for k = 1:c
        @inbounds A[i, k] = s * A[i, k]
    end
end


"""
`row_add_mult!(A,i,s,j)` adds `s` times row `i` to row `j`
in the matrix `A`.
"""
function row_add_mult!(A::AbstractMatrix{T}, i::Int, s, j::Int) where {T}
    r, c = size(A)
    @assert (1 <= i <= r) && (1 <= j <= r) "Row index out of bounds"

    save_row = [s * A[i, k] for k = 1:c]
    for k = 1:c
        @inbounds A[j, k] = A[j, k] + save_row[k]
    end
end


<|start_filename|>src/rankx.jl<|end_filename|>
export rankx

"""
`rankx(A)` computes the rank of the exact matrix `A`
"""
function rankx(A::AbstractMatrix{T})::Int where {T}
    r, c = size(A)

    AA = deepcopy(A)
    try
        AA //= 1  # convert to rational if need be
        try
            rrefx!(AA)
        catch
            AA = big.(AA)
            rrefx!(AA)
        end
    catch
        rrefx!(AA)
    end

    count = 0
    for i = 1:r
        if any(AA[i, :] .!= 0)
            count += 1
        end
    end
    return count
end


<|start_filename|>src/nullspacex.jl<|end_filename|>
export nullspacex

"""
`nullspacex(A)` returns an exact basis for the matrix `A`
"""
function nullspacex(A::AbstractMatrix{T}) where {T}
    r, c = size(A)
    B = rrefx(A)

    leads = Int[]
    # in each row, find first 1
    for i = 1:r
        row = B[i, :]
        if all(row .== 0)
            continue
        end
        k = findfirst(row .!= 0)
        append!(leads, k)
    end

    frees = setdiff(collect(1:c), leads)


    result = Matrix{T}(undef, c, 0)
    for f in frees
        v = zeros(T, c)
        v[f] = T(1)
        for i = 1:length(leads)
            l = leads[i]  # (i,l) is a leading 1
            if l < f
                v[l] = -B[i, f]
            end
        end
        result = [result v]
    end
    return result
end


function nullspacex(A::AbstractMatrix{T}) where {T<:IntegerX}
    return nullspacex(big.(A) // 1)
end


<|start_filename|>src/LinearAlgebraX.jl<|end_filename|>
module LinearAlgebraX
using LinearAlgebra, SimplePolynomials, Mods, Permutations


# IntegerX is any sort of real or Gaussian integer
IntegerX = Union{S,Complex{S}} where {S<:Integer}

# RationalX is a Rational or Complex Rational based on integers
RationalX = Union{Rational{S},Complex{Rational{S}}} where {S<:Integer}

TypeX = Union{IntegerX,RationalX,AbstractMod}


function _recip(x::T) where {T<:IntegerX}
    return 1 // x
end
_recip(x) = inv(x)

include("row_ops.jl")
include("detx.jl")
include("cofactor_det.jl")
include("eye.jl")
include("rrefx.jl")
include("invx.jl")
include("rankx.jl")
include("nullspacex.jl")
include("char_poly.jl")
include("homogeneous.jl")
include("perm.jl")


end # module
