"""
    RiemannMatrix

Type containing information about the Riemann Matrix, for use in computing the Riemann theta function.
"""
struct RiemannMatrix
    τ::Array{<:Number} # matrix in Siegel upper-half space
    g::Integer # genus
    X::Array{<:Real} # real part of τ
    Y::Array{<:Real} # imaginary part of τ
    T::Array{<:Real} # upper triangular matrix in Cholesky factorization of Y
    ellipsoid::Array{} # points to sum over for computing theta functions

    """
        RiemannMatrix(τ, siegel=false, ϵ=1.0e-12, nderivs=4)

    Construct a RiemannMatrix type using a matrix τ in the Siegel upper-half space, for use in computing the Riemann theta function.

    # Arguments
    - `τ::Array{<:Number}`: 2-dimensional array of complex numbers
    - `siegel::Bool=false`: do a siegel transformation on τ if true, and use the input matrix otherwise
    - `ϵ::Real=1.0e-12`: absolute error in value of theta function or its derivatives
    - `nderivs::Integer=4`: highest order of the derivatives of the theta function

    # Examples
    ```julia
    julia> RiemannMatrix([1+im -1; -1 1+im], siegel=true, nderivs=2)
    ```
    """
    function RiemannMatrix(τ::Array{<:Number}; siegel::Bool=false, ϵ::Real=1.0e-12, nderivs::Integer=4)
        τ = (siegel ? siegel_transform(τ)[2] : τ);
        X = real(τ);
        Y = imag(τ);
        g = size(Y)[1];
        T = (g > 1 ? convert(Array{Float64}, cholesky(Symmetric(Y)).U) : sqrt.(Y));
        ρ = svp(T)[2]*sqrt(π); # length of shortest vector in lattice generated by √π*T
        radius = [radius_ellipsoid(ϵ, T, Y, ρ, i) for i=0:nderivs]; # radius of ellipsoids for computing theta function, and up to nderivs derivatives.
        ellipsoid = [ellipsoid_uniform(T, r) for r in radius];
        new(τ, g, X, Y, T, ellipsoid)
    end
end


"""
    random_siegel(g)

Sample a random matrix in the Siegel upper half space with genus g.

# Arguments
- `g::Integer`: genus.

# Examples
```julia
julia> random_siegel(5)
```
"""
function random_siegel(g::Integer)
    Mx = 2*rand(g, g) - ones(g,g); # entries are random between [-1,1)
    My = 2*rand(g, g) - ones(g,g); # entries are random between [-1,1)
    τ = Symmetric((1/2)*(transpose(Mx) + Mx)) + Symmetric(transpose(My)*My)*im; # symmetric real part and psd imaginary part
    return convert(Array, τ);
end
