text:
export isinvertible, contract

@doc Markdown.doc"""
    order(a::AlgAssAbsOrdIdl) -> AlgAssAbsOrd

Returns the order of which $a$ is a (possibly fractional) left or right ideal
if such an order is known.
"""
function order(a::AlgAssAbsOrdIdl)
  if isdefined(a, :order)
    return a.order
  else
    error("No order defined")
  end
end

@doc Markdown.doc"""
    algebra(a::AlgAssAbsOrdIdl) -> AbsAlgAss

Returns the algebra containing $a$.
"""
algebra(a::AlgAssAbsOrdIdl) = a.algebra

_algebra(a::AlgAssAbsOrdIdl) = algebra(a)

# The basis matrix is (should be) in lowerleft HNF, so if the upper left corner
# is not zero, then the matrix has full rank.
isfull_lattice(a::AlgAssAbsOrdIdl) = !iszero(basis_matrix(a, copy = false)[1, 1])

# Whether I is equal to order(I)
function isone(I::AlgAssAbsOrdIdl)
  if !isone(denominator(I, order(I)))
    return false
  end
  return isone(norm(I, copy = false))
end

function one(I::AlgAssAbsOrdIdl)
  return ideal(order(I), one(order(I)))
end

function Base.copy(I::AlgAssAbsOrdIdl)
  return I
end

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    ideal(A::AbsAlgAss, M::FakeFmpqMat, M_in_hnf::Bool = false) -> AlgAssAbsOrdIdl

Returns the ideal in $A$ with basis matrix $M$.
If `M_in_hnf == true`, it is assumed that $M$ is already in lower left HNF.
"""
function ideal(A::AbsAlgAss{fmpq}, M::FakeFmpqMat, M_in_hnf::Bool = false)
  !M_in_hnf ? M = hnf(M, :lowerleft) : nothing
  return AlgAssAbsOrdIdl{typeof(A), elem_type(A)}(A, M)
end

@doc Markdown.doc"""
    ideal(A::AbsAlgAss, O::AlgAssAbsOrd, M::FakeFmpqMat, side::Symbol = :nothing,
          M_in_hnf::Bool = false)
      -> AlgAssAbsOrdIdl

Returns the ideal of $O$ in $A$ with basis matrix $M$ (in the basis of $A$).
If the ideal is known to be a right/left/twosided ideal of $O$, `side` may be
set to `:right`/`:left`/`:twosided` respectively.
If `M_in_hnf == true`, it is assumed that $M$ is already in lower left HNF.
"""
function ideal(A::AbsAlgAss{fmpq}, O::AlgAssAbsOrd, M::FakeFmpqMat, side::Symbol = :nothing, M_in_hnf::Bool = false)
  a = ideal(A, M, M_in_hnf)
  a.order = O
  _set_sidedness(a, side)
  if ismaximal_known(O) && ismaximal(O)
    if side == :left || side == :twosided
      a.left_order = O
    end
    if side == :right || side == :twosided
      a.right_order = O
    end
  end
  return a
end

@doc Markdown.doc"""
    ideal(O::AlgAssAbsOrd, x::AbsAlgAssElem) -> AlgAssAbsOrdIdl
    ideal(O::AlgAssAbsOrd, x::AlgAssAbsOrdElem) -> AlgAssAbsOrdIdl

Returns the twosided principal ideal of $O$ generated by $x$.
"""
function ideal(O::AlgAssAbsOrd{S, T}, x::T) where {S, T}
  A = algebra(O)
  @assert parent(x) === A
  if iszero(x)
    return _zero_ideal(A, O)
  end

  if iscommutative(O)
    a = ideal(O, x, :left)
    a.isright = 1
    return a
  end

  t1 = A()
  t2 = A()
  M = zero_matrix(FlintQQ, dim(A)^2, dim(A))
  b = basis_alg(O, copy = false)
  for i = 1:degree(O)
    t1 = mul!(t1, b[i], x)
    ii = (i - 1)*degree(O)
    for j = 1:degree(O)
      t2 = mul!(t2, t1, b[j])
      elem_to_mat_row!(M, ii + j, t2)
    end
  end
  M = sub(hnf(FakeFmpqMat(M), :lowerleft), nrows(M) - degree(O) + 1:nrows(M), 1:ncols(M))

  return ideal(A, O, M, :twosided, true)
end

ideal(O::AlgAssAbsOrd{S, T}, x::AlgAssAbsOrdElem{S, T}) where { S, T } = ideal(O, elem_in_algebra(x, copy = false))

@doc Markdown.doc"""
    ideal(O::AlgAssAbsOrd, x::AbsAlgAssElem, side::Symbol) -> AlgAssAbsOrdIdl
    ideal(O::AlgAssAbsOrd, x::AlgAssAbsOrdElem, side::Symbol) -> AlgAssAbsOrdIdl

Returns the ideal $O \cdot x$ if `side == :left`, and $x \cdot O$ if
`side == :right`.
"""
function ideal(O::AlgAssAbsOrd{S, T}, x::T, side::Symbol) where { S, T }
  A = algebra(O)
  @assert parent(x) === A
  if iszero(x)
    return _zero_ideal(A, O)
  end
  if side == :twosided
    return ideal(O, x)
  end

  if side == :left
    M = basis_matrix(O)*FakeFmpqMat(representation_matrix(x, :right))
  elseif side == :right
    M = basis_matrix(O)*FakeFmpqMat(representation_matrix(x, :left))
  else
    error("Option :$(side) for side not implemented")
  end
  return ideal(A, O, M, side)
end

ideal(O::AlgAssAbsOrd{S, T}, x::AlgAssAbsOrdElem{S, T}, side::Symbol) where { S, T } = ideal(O, elem_in_algebra(x, copy = false), side)

@doc Markdown.doc"""
    *(O::AlgAssAbsOrd, x::AbsAlgAssElem) -> AlgAssAbsOrdIdl
    *(O::AlgAssAbsOrd, x::AlgAssAbsOrdElem) -> AlgAssAbsOrdIdl
    *(O::AlgAssAbsOrd, x::Int) -> AlgAssAbsOrdIdl
    *(O::AlgAssAbsOrd, x::fmpz) -> AlgAssAbsOrdIdl
    *(x::AbsAlgAssElem, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    *(x::AlgAssAbsOrdElem, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    *(x::Int, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    *(x::fmpz, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl

Returns the ideal $O \cdot x$ or $x \cdot O$ respectively.
"""
*(O::AlgAssAbsOrd{S, T}, x::T) where {S, T <: AlgAssElem} = ideal(O, x, :left)
*(x::T, O::AlgAssAbsOrd{S, T}) where {S, T <: AlgAssElem} = ideal(O, x, :right)
*(O::AlgAssAbsOrd{S, T}, x::AlgAssAbsOrdElem{S, T}) where {S, T} = ideal(O, x, :left)
*(x::AlgAssAbsOrdElem{S, T}, O::AlgAssAbsOrd{S, T}) where {S, T} = ideal(O, x, :right)
*(O::AlgAssAbsOrd, x::Union{ Int, fmpz }) = ideal(O, O(x), :left)
*(x::Union{ Int, fmpz }, O::AlgAssAbsOrd) = ideal(O, O(x), :right)

@doc Markdown.doc"""
    ideal_from_lattice_gens(A::AbsAlgAss, v::Vector{ <: AbsAlgAssElem })
      -> AlgAssAbsOrdIdl

Returns the ideal in $A$ generated by the elements of `v` as a
$\mathbb Z$-lattice.
"""
function ideal_from_lattice_gens(A::AbsAlgAss{fmpq}, v::Vector{ <: AbsAlgAssElem{fmpq} })
  M = zero_matrix(FlintQQ, max(length(v), dim(A)), dim(A))
  for i = 1:length(v)
    elem_to_mat_row!(M, i, v[i])
  end
  M = hnf(FakeFmpqMat(M), :lowerleft)
  if length(v) >= dim(A)
    M = sub(M, (nrows(M) - dim(A) + 1):nrows(M), 1:dim(A))
  end

  return ideal(A, M, true)
end

@doc Markdown.doc"""
     ideal_from_lattice_gens(A::AbsAlgAss, O::AlgAssAbsOrd,
                             v::Vector{ <: AbsAlgAssElem},
                             side::Symbol = :nothing)
       -> AlgAssAbsOrdIdl

Returns the ideal of $O$ in $A$ generated by the elements of `v` as a
lattice over $\mathbb Z$.
If the ideal is known to be a right/left/twosided ideal of $O$, `side` may be
set to `:right`/`:left`/`:twosided` respectively.
"""
function ideal_from_lattice_gens(A::S, O::AlgAssAbsOrd{S, T}, v::Vector{T}, side::Symbol = :nothing) where { S <: AbsAlgAss{fmpq}, T <: AbsAlgAssElem{fmpq} }
  a = ideal_from_lattice_gens(A, v)
  a.order = O
  _set_sidedness(a, side)
  return a
end

###############################################################################
#
#  Zero ideal
#
###############################################################################

function _zero_ideal(A::AbsAlgAss{fmpq})
  a = ideal(A, FakeFmpqMat(zero_matrix(FlintZZ, dim(A), dim(A)), fmpz(1)), true)
  a.iszero = 1
  return a
end

function _zero_ideal(A::AbsAlgAss{fmpq}, O::AlgAssAbsOrd)
  a = _zero_ideal(A)
  a.order = O
  _set_sidedness(a, :twosided)
  return a
end

function iszero(a::AlgAssAbsOrdIdl)
  if a.iszero == 0
    if iszero_row(basis_matrix(a, copy = false).num, dim(algebra(a)))
      a.iszero = 1
    else
      a.iszero = 2
    end
  end
  return a.iszero == 1
end

###############################################################################
#
#  String I/O
#
###############################################################################

function show(io::IO, a::AlgAssAbsOrdIdl)
  print(io, "Ideal in ")
  show(IOContext(io, :compact => true), algebra(a))
  print(io, " with basis matrix\n")
  show(IOContext(io, :compact => true), basis_matrix(a, copy = false))
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::AlgAssAbsOrdIdl, dict::IdDict)
  b = typeof(a)(algebra(a))
  for i in fieldnames(typeof(a))
    if isdefined(a, i)
      if i != :algebra && i != :order && i != :right_order && i != :left_order && i != :basis_matrix_wrt && i != :norm && i != normred
        setfield!(b, i, Base.deepcopy_internal(getfield(a, i), dict))
      end
    end
  end
  if isdefined(a, :order)
    b.order = order(a)
  end
  if isdefined(a, :right_order)
    b.right_order = right_order(a)
  end
  if isdefined(a, :left_order)
    b.left_order = left_order(a)
  end
  for (O, M) in a.basis_matrix_wrt
    b.basis_matrix_wrt[O] = deepcopy(M)
  end
  for (O, n) in a.norm
    b.norm[O] = deepcopy(n)
  end
  for (O, n) in a.normred
    b.normred[O] = deepcopy(n)
  end

  return b
end

################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function assure_has_basis(a::AlgAssAbsOrdIdl)
  if isdefined(a, :basis)
    return nothing
  end

  M = basis_matrix(a, copy = false)
  A = algebra(a)
  b = Vector{elem_type(A)}(undef, dim(A))
  for i = 1:dim(A)
    b[i] = elem_from_mat_row(A, M.num, i, M.den)
  end
  a.basis = b
  return nothing
end

function assure_has_basis_matrix_wrt(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd)
  if haskey(a.basis_matrix_wrt, O)
    return nothing
  end

  M = hnf(basis_matrix(a, copy = false)*basis_mat_inv(O, copy = false), :lowerleft)
  a.basis_matrix_wrt[O] = M
  return nothing
end

function assure_has_basis_mat_inv(a::AlgAssAbsOrdIdl)
  if isdefined(a, :basis_mat_inv)
    return nothing
  end
  @assert isfull_lattice(a) "The ideal is not a full lattice"
  a.basis_mat_inv = inv(basis_matrix(a, copy = false))
  return nothing
end

################################################################################
#
#  Basis (matrices)
#
################################################################################

@doc Markdown.doc"""
    basis(a::AlgAssAbsOrdIdl; copy::Bool = true) -> Vector{AbsAlgAssElem}

Returns the basis of $a$.
"""
function basis(a::AlgAssAbsOrdIdl; copy::Bool = true)
  assure_has_basis(a)
  if copy
    return deepcopy(a.basis)
  else
    return a.basis
  end
end

absolute_basis(a::AlgAssAbsOrdIdl) = basis(a)

@doc Markdown.doc"""
    basis_matrix(a::AlgAssAbsOrdIdl; copy::Bool = true) -> FakeFmpqMat

Returns the basis matrix of $a$ with respect to the basis of the algebra.
"""
function basis_matrix(a::AlgAssAbsOrdIdl; copy::Bool = true)
  if copy
    return deepcopy(a.basis_matrix)
  else
    return a.basis_matrix
  end
end

@doc Markdown.doc"""
    basis_mat_inv(a::AlgAssAbsOrdIdl; copy::Bool = true) -> FakeFmpqMat

Returns the inverse of the basis matrix of $a$.
"""
function basis_mat_inv(a::AlgAssAbsOrdIdl; copy::Bool = true)
  assure_has_basis_mat_inv(a)
  if copy
    return deepcopy(a.basis_mat_inv)
  else
    return a.basis_mat_inv
  end
end

@doc Markdown.doc"""
    basis_matrix_wrt(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd; copy::Bool = true)
      -> FakeFmpqMat

Returns the basis matrix of $a$ with respect to the basis of $O$.
"""
function basis_matrix_wrt(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd; copy::Bool = true)
  assure_has_basis_matrix_wrt(a, O)
  if copy
    return deepcopy(a.basis_matrix_wrt[O])
  else
    return a.basis_matrix_wrt[O]
  end
end

@doc Markdown.doc"""
    integral_basis_matrix_wrt(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd;
                              copy::Bool = true)
      -> fmpz_mat

Returns the basis matrix of $a$ with respect to the basis of $O$ as an integral
matrix, assuming the ideal is contained in $O$.
"""
function integral_basis_matrix_wrt(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd; copy::Bool = true)
  M = basis_matrix_wrt(a, O, copy = copy)
  @assert isone(denominator(M, copy = false)) "Ideal is not integral"
  return numerator(M, copy = false)
end

################################################################################
#
#  Arithmetic
#
################################################################################

@doc Markdown.doc"""
    +(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

Returns $a + b$.
"""
function +(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  @assert algebra(a) === algebra(b)
  if iszero(a)
    return deepcopy(b)
  elseif iszero(b)
    return deepcopy(a)
  end

  d = dim(algebra(a))
  M = vcat(basis_matrix(a, copy = false), basis_matrix(b, copy = false))
  if all(i -> !iszero(M[i, i]), 1:d)
    M = sub(hnf(M, :lowerleft, triangular_top = true), (d + 1):2*d, 1:d)
  else
    M = sub(hnf(M, :lowerleft), (d + 1):2*d, 1:d)
  end
  c = ideal(algebra(a), M, true)
  if isdefined(a, :order) && isdefined(b, :order) && order(a) === order(b)
    c.order = order(a)
  end
  return c
end

@doc Markdown.doc"""
    *(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

Returns $a \cdot b$.
"""
function *(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  @assert algebra(a) === algebra(b)

  if iszero(a)
    return deepcopy(a)
  elseif iszero(b)
    return deepcopy(b)
  end

  A = algebra(a)
  d = dim(A)
  ba = basis(a, copy = false)
  bb = basis(b, copy = false)
  d2 = d^2

  M = zero_matrix(FlintQQ, d2, d)
  t = one(A)
  for i = 1:d
    i1d = (i - 1)*d
    for j = 1:d
      t = mul!(t, ba[i], bb[j])
      elem_to_mat_row!(M, i1d + j, t)
    end
  end

  H = sub(hnf(FakeFmpqMat(M), :lowerleft), (d2 - d + 1):d2, 1:d)
  c = ideal(A, H, true)

  if _left_order_known_and_maximal(a)
    c.left_order = left_order(a)
  end
  if _right_order_known_and_maximal(b)
    c.right_order = right_order(b)
  end
  if isdefined(a, :order) && isdefined(b, :order) && order(a) === order(b)
    c.order = order(a)
  end

  return c
end

@doc Markdown.doc"""
    ^(a::AlgAssAbsOrdIdl, e::Int) -> AlgAssAbsOrdIdl
    ^(a::AlgAssAbsOrdIdl, e::fmpz) -> AlgAssAbsOrdIdl

Returns $a^e$.
"""
^(A::AlgAssAbsOrdIdl, e::Int) = Base.power_by_squaring(A, e)
^(A::AlgAssAbsOrdIdl, e::fmpz) = Base.power_by_squaring(A, BigInt(e))

@doc Markdown.doc"""
    intersect(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

Returns $a \cap b$.
"""
function intersect(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  d = dim(algebra(a))
  M1 = hcat(basis_matrix(a, copy = false), basis_matrix(a, copy = false))
  M2 = hcat(FakeFmpqMat(zero_matrix(FlintZZ, d, d), fmpz(1)), basis_matrix(b, copy = false))
  M = vcat(M1, M2)
  H = sub(hnf(M, :lowerleft), 1:d, 1:d)
  c = ideal(algebra(a), H, true)

  if isdefined(a, :order) && isdefined(b, :order) && order(a) === order(b)
    c.order = order(a)
  end
  return c
end

################################################################################
#
#  Ad hoc multiplication
#
################################################################################

@doc Markdown.doc"""
    *(a::AlgAssAbsOrdIdl, x::Int) -> AlgAssAbsOrdIdl
    *(x::Int, a::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl
    *(a::AlgAssAbsOrdIdl, x::fmpz) -> AlgAssAbsOrdIdl
    *(x::fmpz, a::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl
    *(a::AlgAssAbsOrdIdl, x::fmpq) -> AlgAssAbsOrdIdl
    *(x::fmpq, a::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl
    *(a::AlgAssAbsOrdIdl{S, T}, x::T) where { S, T } -> AlgAssAbsOrdIdl{S, T}
    *(x::T, a::AlgAssAbsOrdIdl{S, T}) where { S, T } -> AlgAssAbsOrdIdl{S, T}
    *(a::AlgAssAbsOrdIdl{S, T}, x::AlgAssAbsOrdElem{S, T}) where { S, T }
      -> AlgAssAbsOrdIdl{S, T}
    *(x::AlgAssAbsOrdElem{S, T}, a::AlgAssAbsOrdIdl{S, T}) where { S, T }
      -> AlgAssAbsOrdIdl{S, T}

Returns the ideal $a*x$ respectively $x*a$.
"""
function *(a::AlgAssAbsOrdIdl, x::Union{ Int, fmpz, fmpq })
  if iszero(x)
    return _zero_ideal(algebra(a))
  end
  b = ideal(algebra(a), x*basis_matrix(a, copy = false), true)
  if isdefined(a, :left_order)
    b.left_order = left_order(a)
  end
  if isdefined(a, :right_order)
    b.right_order = right_order(a)
  end
  if isdefined(a, :order)
    b.order = order(a)
  end
  return b
end

*(x::Union{ Int, fmpz, fmpq }, a::AlgAssAbsOrdIdl) = a*x

function *(a::AlgAssAbsOrdIdl{S, T}, x::T) where { S, T <: Union{AlgAssElem, AlgGrpElem} }
  if iszero(x)
    return _zero_ideal(algebra(a))
  end
  M = basis_matrix(a, copy = false)*FakeFmpqMat(representation_matrix(x, :right))
  b = ideal(algebra(a), M)
  if isdefined(a, :left_order)
    b.left_order = left_order(a)
  end
  return b
end

function *(x::T, a::AlgAssAbsOrdIdl{S, T}) where { S, T <: Union{AlgAssElem, AlgGrpElem}}
  if iszero(x)
    return _zero_ideal(algebra(a))
  end
  M = basis_matrix(a, copy = false)*FakeFmpqMat(representation_matrix(x, :left))
  b = ideal(algebra(a), M)
  if isdefined(a, :right_order)
    b.right_order = right_order(a)
  end
  return b
end

function *(a::AlgAssAbsOrdIdl{S, T}, x::AlgAssAbsOrdElem{S, T}) where { S, T }
  b = a*elem_in_algebra(x, copy = false)
  if isdefined(a, :order) && parent(x) === order(a) && isright_ideal(a)
    b.order = order(a)
    b.isright = a.isright
    b.isleft = a.isleft
  end
  return b
end

function *(x::AlgAssAbsOrdElem{S, T}, a::AlgAssAbsOrdIdl{S, T}) where { S, T }
  b = elem_in_algebra(x, copy = false)*a
  if isdefined(a, :order) && parent(x) === order(a) && isleft_ideal(a)
    b.order = order(a)
    b.isright = a.isright
    b.isleft = a.isleft
  end
  return b
end

################################################################################
#
#  Inclusion of elements in ideals
#
################################################################################

@doc Markdown.doc"""
    in(x::AbsAlgAssElem, a::AlgAssAbsOrdIdl) -> Bool
    in(x::AlgAssAbsOrdElem, a::AlgAssAbsOrdIdl) -> Bool

Returns `true` if $x$ is in $a$ and `false` otherwise.
"""
function in(x::T, a::AlgAssAbsOrdIdl{S, T}) where { S, T }
  parent(x) !== algebra(a) && error("Algebra of element and ideal must be equal")
  A = algebra(a)
  t = FakeFmpqMat(matrix(FlintQQ, 1, dim(A), coeffs(x, copy = false)))
  t = t*basis_mat_inv(a, copy = false)
  return denominator(t, copy = false) == 1
end

in(x::AlgAssAbsOrdElem{S, T}, a::AlgAssAbsOrdIdl{S, T}) where { S, T } = in(elem_in_algebra(x, copy = false), a)

################################################################################
#
#  Equality
#
################################################################################

@doc Markdown.doc"""
    ==(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> Bool

Returns `true` if $a$ and $b$ are equal and `false` otherwise.
"""
function ==(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where { S, T }
  algebra(a) !== algebra(b) && return false
  return basis_matrix(a, copy = false) == basis_matrix(b, copy = false)
end

################################################################################
#
#  isleft/isright
#
################################################################################

# functions isright_ideal and isleft_ideal are in AlgAss/Ideal.jl

# This only works if a.order is defined, otherwise order(a) throws an error.
function _test_ideal_sidedness(a::AlgAssAbsOrdIdl, side::Symbol)
  O = order(a)
  return _test_ideal_sidedness(a, O, side)
end

function _test_ideal_sidedness(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd, side::Symbol)
  A = algebra(a)
  ba = basis(a, copy = false)
  bO = basis_alg(O, copy = false)
  t = A()
  for i = 1:dim(A)
    for j = 1:dim(A)
      if side == :left
        t = mul!(t, bO[i], ba[j])
      elseif side == :right
        t = mul!(t, ba[j], bO[i])
      else
        error("side must be either :left or :right")
      end
      if !(t in a)
        return false
      end
    end
  end
  return true
end

################################################################################
#
#  Ring of multipliers, left and right order
#
################################################################################

# Optimized version of ring_of_multipliers for the maximal order computation.
function _ring_of_multipliers_integral_ideal(I::AlgAssAbsOrdIdl, p::fmpz = fmpz(1))
  O = order(I)
  @hassert :AlgAssOrd 1 Hecke.check_associativity(algebra(O))
  @hassert :AlgAssOrd 1 Hecke.check_distributivity(algebra(O))
  if isdefined(I, :gens) && length(I.gens) < degree(O)
    B = I.gens
  else
    B = basis(I, copy = false)
  end
  m = zero_matrix(FlintZZ, degree(O)*length(B), degree(O))
  for i = 1:length(B)
    M = FakeFmpqMat(representation_matrix(B[i]))
    M = mul!(M, basis_matrix(O, copy = false), M)
    M = mul!(M, M, basis_mat_inv(I, copy = false))
    # M is now the representation matrix of B[i] as element of O multiplied
    # from the left by the inverse basis matrix of I in the basis of O.
    @hassert :AlgAssOrd 1 denominator(M, copy = false) == 1
    M = numerator(M, copy = false)
    for s = 1:degree(O)
      for t = 1:degree(O)
        m[t + (i - 1)*degree(O), s] = M[s, t]
      end
    end
  end
  # In the case of the p-radical, it is important to do this modulo p
  if p == 1
    m = hnf(m)
  else
    hnf_modular_eldiv!(m, p)
  end
  s = prod(fmpz[m[i, i] for i = 1:ncols(m)])
  if s == 1
    return O
  end
  # n is upper right HNF
  n = transpose(view(m, 1:degree(O), 1:degree(O)))
  b = FakeFmpqMat(pseudo_inv(n))
  mul!(b, b, basis_matrix(O, copy = false))
  @hassert :AlgAssOrd 1 defines_order(algebra(O), b)[1]
  O1 = Order(algebra(O), b)
  O1.disc = divexact(discriminant(O), s^2)
  return O1
end

# This computes a basis matrix for \{ x \in A | bx \subseteq a \} if
# side == :left or \{ x \in A | xb \subseteq a \} if side == :right.
function _colon_raw(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}, side::Symbol) where { S, T }
  @assert isfull_lattice(a) && isfull_lattice(b)
  A = algebra(a)
  @assert A === algebra(b)
  K = base_ring(A)
  d = dim(A)
  bb = basis(b, copy = false)
  B = fmpq_mat(basis_mat_inv(a, copy = false))
  M = zero_matrix(FlintQQ, d^2, d)
  for i = 1:d
    N = representation_matrix(bb[i], side)*B
    for s = 1:d
      for t = 1:d
        M[t + (i - 1)*d, s] = N[s, t]
      end
    end
  end
  M = sub(hnf(FakeFmpqMat(M), :upperright), 1:d, 1:d)
  N = inv(transpose(M))
  return N
end

@doc Markdown.doc"""
    ring_of_multipliers(a::AlgAssAbsOrdIdl) -> AlgAssAbsOrd

Given an ideal $a$, it returns the ring $(a : a)$.
"""
function ring_of_multipliers(a::AlgAssAbsOrdIdl, action::Symbol = :left)
  M = _colon_raw(a, a, action)
  return Order(algebra(a), hnf(M))
end

@doc Markdown.doc"""
    left_order(a::AlgAssAbsOrdIdl) -> AlgAssAbsOrd

Returns the largest order of which $a$ is a left ideal, that is
$\{ x \in A \mid xa \subseteq a\}$.
"""
function left_order(a::AlgAssAbsOrdIdl)
  if isdefined(a, :left_order)
    return a.left_order
  end

  a.left_order = ring_of_multipliers(a, :right)
  return a.left_order
end

@doc Markdown.doc"""
    right_order(a::AlgAssAbsOrdIdl) -> AlgAssAbsOrd

Returns the largest order of which $a$ is a right ideal, that is
$\{ x \in A \mid ax \subseteq a\}$.
"""
function right_order(a::AlgAssAbsOrdIdl)
  if isdefined(a, :right_order)
    return a.right_order
  end

  a.right_order = ring_of_multipliers(a, :left)
  return a.right_order
end

# Checks whether the left order is known and if it is known to be maximal
function _left_order_known_and_maximal(a::AlgAssAbsOrdIdl)
  if !isdefined(a, :left_order)
    return false
  end
  O = left_order(a)
  if ismaximal_known(O) && ismaximal(O)
    return true
  end
  return false
end

# Checks whether the right order is known and if it is known to be maximal
function _right_order_known_and_maximal(a::AlgAssAbsOrdIdl)
  if !isdefined(a, :right_order)
    return false
  end
  O = right_order(a)
  if ismaximal_known(O) && ismaximal(O)
    return true
  end
  return false
end

################################################################################
#
#  Norm
#
################################################################################

function assure_has_norm(a::AlgAssAbsOrdIdl{S, T}, O::AlgAssAbsOrd{S, T}) where { S, T }
  if haskey(a.norm, O)
    return nothing
  end

  a.norm[O] = abs(det(basis_matrix(a, copy = false))*det(basis_mat_inv(O, copy = false)))
  return nothing
end

@doc Markdown.doc"""
    norm(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd; copy::Bool = true) -> fmpq

Returns the norm of $a$ considered as an (possibly fractional) ideal of $O$.
"""
function norm(a::AlgAssAbsOrdIdl{S, T}, O::AlgAssAbsOrd{S, T}; copy::Bool = true) where { S, T }
  assure_has_norm(a, O)
  if copy
    return deepcopy(a.norm[O])
  else
    return a.norm[O]
  end
end

@doc Markdown.doc"""
    norm(a::AlgAssAbsOrdIdl; copy::Bool = true) -> fmpq

Returns the norm of $a$ considered as an (possibly fractional) ideal of
`order(a)`.
"""
norm(a::AlgAssAbsOrdIdl; copy::Bool = true) = norm(a, order(a), copy = copy)

function assure_has_normred(a::AlgAssAbsOrdIdl{S, T}, O::AlgAssAbsOrd{S, T}) where { S, T }
  if haskey(a.normred, O)
    return nothing
  end

  A = algebra(a)
  m = isqrt(dim(A))
  @assert m^2 == dim(A)
  N = norm(a, O, copy = false)
  b, n = ispower(N, m)
  @assert b "Cannot compute reduced norm. Maybe the algebra is not simple and central?"
  a.normred[O] = n
  return nothing
end

@doc Markdown.doc"""
    normred(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd; copy::Bool = true) -> fmpq

Returns the reduced norm of $a$ considered as an (possibly fractional) ideal
of $O$.
It is assumed that the algebra containing $a$ is simple and central.
"""
function normred(a::AlgAssAbsOrdIdl{S, T}, O::AlgAssAbsOrd{S, T}; copy::Bool = true) where { S, T }
  @assert issimple(algebra(a)) && iscentral(algebra(a)) "Only implemented for simple and central algebras"
  assure_has_normred(a, O)
  if copy
    return deepcopy(a.normred[O])
  else
    return a.normred[O]
  end
end

@doc Markdown.doc"""
    normred(a::AlgAssAbsOrdIdl; copy::Bool = true) -> fmpq

Returns the reduced norm of $a$ considered as an (possibly fractional) ideal
of `order(a)`.
"""
normred(a::AlgAssAbsOrdIdl; copy::Bool = true) = normred(a, order(a), copy = copy)

################################################################################
#
#  Locally free basis
#
################################################################################

@doc Markdown.doc"""
    locally_free_basis(a::AlgAssAbsOrdIdl, p::Union{ Int, fmpz })
      -> AlgAssAbsOrdElem

Returns an element $x$ of the order $O$ of $a$ such that $a_p = O_p \cdot x$
where $p$ is a prime of $\mathbb Z$.
See also `islocally_free`.
"""
locally_free_basis(I::AlgAssAbsOrdIdl, p::Union{Int, fmpz}; side::Symbol = :right) = locally_free_basis(order(I), I, p, side = side)

@doc Markdown.doc"""
    locally_free_basis(O::AlgAssAbsOrd, a::AlgAssAbsOrdIdl,
                       p::Union{ Int, fmpz }) -> AlgAssAbsOrdElem

Returns an element $x$ of $O$ such that $a_p = O_p \cdot x$ where $p$ is a
prime of $\mathbb Z$.
It is assumed that $a$ is an ideal of $O$ and $a \subseteq O$.
See also `islocally_free`.
"""
function locally_free_basis(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl, p::Union{Int, fmpz}; side::Symbol = :right)
  b, x = islocally_free(O, I, p, side = side)
  if !b
    error("The ideal is not locally free at the prime")
  end
  return x
end

@doc Markdown.doc"""
    islocally_free(a::AlgAssAbsOrdIdl, p::Union{ Int, fmpz })
      -> Bool, AlgAssAbsOrdElem

Returns a tuple `(true, x)` with an element $x$ of the order $O$ of $a$ such
that $a_p = O_p x$ if $a$ is locally free at $p$, and `(false, 0)` otherwise.
$p$ is a prime of $\mathbb Z$.
See also `locally_free_basis`.
"""
islocally_free(I::AlgAssAbsOrdIdl, p::Union{Int,fmpz}; side::Symbol = :right) = islocally_free(order(I), I, p, side = side)

# See Bley, Wilson "Computations in relative algebraic K-groups", section 4.2
@doc Markdown.doc"""
    islocally_free(O::AlgAssAbsOrd, a::AlgAssAbsOrdIdl, p::fmpz
                   side = :right) -> Bool, AlgAssAbsOrdElem

Returns a tuple `(true, x)` with an element $x$ of $O$ such that $a O_p = x
O_p$ resp. $O_p a = O_p x$ if $a$ is locally free right (resp. left) ideal at
$p$, and `(false, 0)` otherwise.  It is assumed that $a$ is an ideal of $O$ and
$a \subseteq O$.

See also `locally_free_basis`.
"""
function islocally_free(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl, p::Union{Int, fmpz}; side::Symbol = :right)
  b = _test_ideal_sidedness(I, O, side)
  d = denominator(I, O)
  !isone(d) && throw(error("Ideal must be contained in the order"))
  !b && throw(error("Ideal is not a $(side) ideal of the order"))
  if side === :left
    return _islocally_free_right(O, I, p)
  elseif side === :right
    B, mB = opposite_algebra(algebra(O))
    OB = mB(O)
    IB = mB(I)
    IB.order = mB(order(I))
    fl, x = _islocally_free_left(OB, IB, p)
    return fl, O(mB\elem_in_algebra(x))
  else
    throw(error("side (:$(side)) must be either :left or :right"))
  end
end

function _islocally_free_left(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl, p::Union{Int, fmpz})
  OpO, toOpO = AlgAss(O, p*O, p)
  J = radical(OpO)
  OJ, toOJ = quo(OpO, J)
  decOJ = decompose(OJ)
  basisIinO = [ O(b) for b in basis(I, copy = false) ]

  pI = p*I
  IpI, toIpI = AlgAss(I, p*I, p)
  gensJ = Vector{elem_type(IpI)}()
  for b in basis(J, copy = false)
    bb = toOpO\b
    for c in basisIinO
      push!(gensJ, toIpI(bb*c))
    end
  end
  JinIpI = ideal_from_gens(IpI, gensJ)
  IJ, toIJ = quo(IpI, JinIpI)

  a = O()
  for i = 1:length(decOJ)
    A, AtoOJ = decOJ[i]
    B, BtoA = _as_algebra_over_center(A)
    C, BtoC = _as_matrix_algebra(B)
    e = toOpO\(toOJ\(AtoOJ(BtoA(BtoC\C[1]))))
    basiseIJ = Vector{elem_type(IJ)}()
    for b in basisIinO
      bb = toIJ(toIpI(e*b))
      if !iszero(bb)
        push!(basiseIJ, bb)
      end
    end

    # Construct an Fq-basis for e*IJ where Fq \cong centre(A)
    Z, ZtoA = center(A)
    basisZ = elem_type(O)[ toOpO\(toOJ\(AtoOJ(ZtoA(Z[i])))) for i = 1:dim(Z) ]

    basiseIJoverZ = Vector{elem_type(O)}()
    M = zero_matrix(base_ring(IJ), dim(Z), dim(IJ))
    MM = zero_matrix(base_ring(IJ), 0, dim(IJ))
    r = 0
    for i = 1:length(basiseIJ)
      b = toIpI\(toIJ\basiseIJ[i])

      for j = 1:dim(Z)
        bb = toIJ(toIpI(basisZ[j]*b))
        elem_to_mat_row!(M, j, bb)
      end

      N = vcat(MM, M)
      s = rank(N)
      if s > r
        push!(basiseIJoverZ, b)
        MM = N
        r = s
      end
      if r == length(basiseIJ)
        break
      end
    end

    if length(basiseIJoverZ) != degree(C)
      # I is not locally free
      return false, O()
    end

    for i = 1:length(basiseIJoverZ)
      a += mod(toOpO\(toOJ\(AtoOJ(BtoA(BtoC\C[i]))))*basiseIJoverZ[i], pI)
    end
  end

  return true, mod(a, pI)
end

################################################################################
#
#  p-Radical
#
################################################################################

function pradical_meataxe(O::AlgAssAbsOrd, p::Int)
  A1, OtoA1 = quo(O, p*O, p)
  @vtime :AlgAssOrd 1 lg = gens(A1)
  lM = dense_matrix_type(base_ring(A1))[ transpose(representation_matrix(lg[i])) for i = 1:length(lg) ]
  M = ModAlgAss(lM)
  ls = minimal_submodules(M)
  l = sum(nrows(x) for x in ls)
  M1 = zero_matrix(base_ring(A1), l, degree(O))
  i = 1
  for x in ls
    for j = 1:nrows(x)
      for k = 1:degree(O)
        M1[i, k] = x[j, k]
      end
      i += 1
    end
  end
  r = rref!(M1)
  if r == degree(O)
    J = ideal(algebra(O), O, p*basis_matrix(O, copy = false), :twosided)
    J.gens = elem_type(algebra(O))[ p*one(algebra(O)) ]
    return J
  end
  M1 = view(M1, 1:r, 1:degree(O))
  dM = transpose(nullspace(M1)[2])
  g = Vector{elem_type(algebra(O))}(undef, nrows(dM) + 1)
  m = zero_matrix(FlintZZ, degree(O), degree(O))
  for i = 1:nrows(dM)
    for j = 1:ncols(dM)
      m[i, j] = lift(dM[i, j])
    end
    g[i] = elem_in_algebra(elem_from_mat_row(O, m, i), copy = false)
  end
  m = hnf_modular_eldiv!(m, fmpz(p))
  m = m*basis_matrix(O, copy = false)
  g[nrows(dM) + 1] = p*one(algebra(O))
  J = ideal(algebra(O), O, m, :twosided)
  J.gens = g
  return J
end

@doc Markdown.doc"""
    pradical(O::AlgAssAbsOrd, p::Int) -> AlgAssAbsOrdIdl

Given an order $O$ and a prime $p$, it returns the radical of the ideal
generated by $p$.
"""
function pradical(O::AlgAssAbsOrd, p::Int)
  # See the paper from Ronyai, Structure of finite algebra
  l = clog(fmpz(degree(O)), p)
  if l > 1
    return pradical_meataxe(O, p)
  end
  n = root(degree(O), 2)
  F = GF(p, cached = false)

  I = change_base_ring(F, n*trred_matrix(O))
  k, B = nullspace(I)
  # The columns of B give the coordinates of the elements in the order.
  if k == 0
    J = ideal(algebra(O), O, p*basis_matrix(O, copy = false), :twosided)
    J.gens = elem_type(algebra(O))[ p*one(algebra(O)) ]
    return J
  end
  # We have l == 1. In this case, we can output I: it is the standard p-trace
  # method.
  M = zero_matrix(FlintZZ, degree(O), degree(O))
  for i = 1:ncols(B)
    for j = 1:degree(O)
      M[i, j] = lift(B[j, i])
    end
  end
  M = hnf_modular_eldiv!(M, fmpz(p)) # This puts p in the "missing" pivot entries
  M = M*basis_matrix(O, copy = false)
  res = ideal(algebra(O), O, M, :twosided)
  B1 = lift(B')
  res.gens = Vector{elem_type(algebra(O))}(undef, k + 1)
  for i = 1:k
    res.gens[i] = elem_in_algebra(elem_from_mat_row(O, B1, i), copy = false)
  end
  res.gens[k + 1] = p*one(algebra(O))
  return res
end

################################################################################
#
#  Primes lying over a prime
#
################################################################################

# Returns the twosided maximal ideals of O containing I, where p*O \subseteq I.
# If strict_containment == true and I is already prime, we return an empty array.
function _maximal_ideals(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl, p::Union{Int, fmpz}; strict_containment::Bool = false)
  A1, OtoA1 = quo(O, I, p)
  @vtime :AlgAssOrd 1 lg = gens(A1)
  lM = dense_matrix_type(base_ring(A1))[ representation_matrix(lg[i]) for i = 1:length(lg) ]
  append!(lM, dense_matrix_type(base_ring(A1))[ representation_matrix(lg[i], :right) for i = 1:length(lg) ])
  M = ModAlgAss(lM)
  ls = maximal_submodules(M)
  if strict_containment && isone(length(ls)) && iszero(nrows(ls[1]))
    ls = typeof(ls[1])[]
  end
  return typeof(I)[_from_submodules_to_ideals(M, O, I, x, A1, OtoA1) for x in ls ]
end

function _from_submodules_to_ideals(M::ModAlgAss, O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl, x::Union{Zmodn_mat, Generic.Mat{Generic.ResF{fmpz}}}, A1::AlgAss, OtoA1::AbsOrdToAlgAssMor)
  @hassert :AlgAssOrd 1 begin r = rref(x)[1]; closure(x, M.action) == sub(rref(x)[2], 1:r, 1:ncols(x)) end
  m = zero_matrix(FlintZZ, nrows(x), degree(O))
  g = Vector{elem_type(algebra(O))}(undef, nrows(x))
  for i = 1:nrows(x)
    el = OtoA1\(elem_from_mat_row(A1, x, i))
    for j = 1:degree(O)
      m[i, j] = coordinates(el, copy = false)[j]
    end
    g[i] = elem_in_algebra(elem_from_mat_row(O, m, i), copy = false)
  end
  m = m*basis_matrix(O, copy = false)
  m = vcat(m, basis_matrix(I, copy = false))
  m = sub(hnf(m, :lowerleft), nrows(x) + 1:nrows(m), 1:degree(O))
  J = ideal(algebra(O), O, m, :twosided, true)
  if isdefined(I, :gens)
    append!(g, I.gens)
    J.gens = g
  else
    append!(g, basis(I, copy = false))
  end
  return J
end

@doc Markdown.doc"""
    prime_ideals_over(O::AlgAssAbsOrd, p::Union{ Int, fmpz })
      -> Vector{AlgAssAbsOrdIdl}

Returns all prime ideals of $O$ lying over the prime number $p$.
"""
function prime_ideals_over(O::AlgAssAbsOrd, p::Union{ Int, fmpz })
  @vtime :AlgAssOrd 1 max_id = collect(_maximal_ideals(O, p*O, Int(p)))
  return max_id
end

@doc Markdown.doc"""
    prime_ideals_over(O::AlgAssAbsOrd, P::AlgAssAbsOrdIdl)
      -> Vector{AlgAssAbsOrdIdl}

Given a prime ideal $P$ in an order contained in $O$, this function returns
the prime ideals of $O$ lying over $P$.
"""
function prime_ideals_over(O::AlgAssAbsOrd, P::AlgAssAbsOrdIdl)
  O1 = order(P)
  if O1 == O
    return ideal_type(O)[ P ]
  end
  M = maximal_order(O)
  lp = prime_ideals_over(M, minimum(P))
  p_critical_primes = Vector{ideal_type(O)}()
  for Q in lp
    c = contract(Q, O1)
    if c == P
      c1 = contract(Q, O)
      if !(c1 in p_critical_primes)
        push!(p_critical_primes, c1)
      end
    end
  end
  return p_critical_primes
end

################################################################################
#
#  Inverses
#
################################################################################

@doc Markdown.doc"""
    inv(a::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

Returns an ideal $b$ in `algebra(a)` such that `a*b == left_order(a)`
and `b*a == right_order(b)`.
"""
function inv(a::AlgAssAbsOrdIdl)
  if isdefined(a, :right_order)
    O = right_order(a)
    M = _colon_raw(O(1)*O, a, :right)
  else
    O = left_order(a)
    M = _colon_raw(O(1)*O, a, :left)
  end
  b = ideal(algebra(a), M)
  if isdefined(a, :left_order)
    b.right_order = left_order(a)
  end
  if isdefined(a, :right_order)
    b.left_order = right_order(a)
  end
  if isdefined(a, :order)
    b.order = order(a)
  end
  return b
end

# Tests whether a is invertible in order(a)
function isinvertible(a::AlgAssAbsOrdIdl)
  if iszero(a)
    return false, a
  end

  bmata = basis_matrix(a, copy = false)

  # If a is zero in a component of algebra(a), then it cannot be invertible
  if any(iszero(bmata[i, i]) for i in 1:degree(order(a)))
    return false, a
  end

  if order(a).ismaximal == 1
    return true, inv(a)
  end

  b = inv(a)
  c = a*b
  return isone(c), b
end

################################################################################
#
#  Divexact
#
################################################################################

@doc Markdown.doc"""
    divexact_left(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

Returns an ideal $c$ such that $a = b \cdot c$.
"""
function divexact_left(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where { S, T }
  @assert algebra(a) === algebra(b)
  M = _colon_raw(a, b, :left)
  c = ideal(algebra(a), M)

  if _right_order_known_and_maximal(a)
    c.right_order = right_order(a)
  end
  if _right_order_known_and_maximal(b)
    c.left_order = right_order(b)
  end

  return c
end

@doc Markdown.doc"""
    divexact_right(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

Returns an ideal $c$ such that $a = c \cdot b$.
"""
function divexact_right(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where { S, T }
  @assert algebra(a) === algebra(b)
  M = _colon_raw(a, b, :right)
  c = ideal(algebra(a), M)

  if _left_order_known_and_maximal(a)
    c.left_order = left_order(a)
  end
  if _left_order_known_and_maximal(b)
    c.right_order = left_order(b)
  end

  return c
end

colon(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where { S, T } = divexact_right(a, b)

################################################################################
#
#  Denominator / Numerator
#
################################################################################

# copy doesn't do anything (yet?)
@doc Markdown.doc"""
    denominator(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) -> fmpz

Returns the smallest positive integer $d$ such that $da$ is contained in $O$.
"""
function denominator(a::AlgAssAbsOrdIdl{S, T}, O::AlgAssAbsOrd{S, T}; copy::Bool = true) where { S, T }
  @assert algebra(a) === algebra(O)
  if iszero(a)
    return fmpz(1)
  end

  M = basis_matrix(a, copy = false)*basis_mat_inv(O, copy = false)
  return denominator(M, copy = false)
end

@doc Markdown.doc"""
    denominator(a::AlgAssAbsOrdIdl) -> fmpz

Returns the smallest positive integer $d$ such that $da$ is contained in
`order(a)`.
"""
denominator(a::AlgAssAbsOrdIdl; copy::Bool = true) = denominator(a, order(a), copy = copy)

# copy doesn't do anything (yet?)
function numerator(a::AlgAssAbsOrdIdl{S, T}, O::AlgAssAbsOrd{S, T}; copy::Bool = true) where { S, T }
  return denominator(a, O)*a
end

numerator(a::AlgAssAbsOrdIdl; copy::Bool = true) = numerator(a, order(a), copy = copy)

################################################################################
#
#  Ideal Set
#
################################################################################

function show(io::IO, a::AlgAssAbsOrdIdlSet)
  print(io, "Set of ideals of $(order(a))\n")
end

order(a::AlgAssAbsOrdIdlSet) = a.order

parent(I::AlgAssAbsOrdIdl) = AlgAssAbsOrdIdlSet(order(I))

function IdealSet(O::AlgAssAbsOrd)
   return AlgAssAbsOrdIdlSet(O)
end

FracIdealSet(O::AlgAssAbsOrd) = IdealSet(O)

elem_type(::Type{AlgAssAbsOrdIdlSet{S, T}}) where {S, T} = AlgAssAbsOrdIdl{S, T}

elem_type(::AlgAssAbsOrdIdlSet{S, T}) where {S, T} = AlgAssAbsOrdIdl{S, T}

parent_type(::Type{AlgAssAbsOrdIdl{S, T}}) where {S, T} = AlgAssAbsOrdIdlSet{S, T}

function Base.one(S::AlgAssAbsOrdIdlSet)
  return ideal(order(S), one(order(S)))
end

==(x::AlgAssAbsOrdIdlSet, y::AlgAssAbsOrdIdlSet) = x.order === y.order

################################################################################
#
#  Extend/contract
#
################################################################################

@doc Markdown.doc"""
    extend(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl

Returns $O \cdot a \cdot O$.
It is assumed that the order of $a$ is contained in $O$.
"""
function extend(A::AlgAssAbsOrdIdl{S, T}, O::AlgAssAbsOrd{S, T}, action::Symbol = :twosided) where S where T
  if action == :left || action == :twosided
    A = (O(1)*O)*A
  end
  if action == :right || action == :twosided
    A = A*(O(1)*O)
  end
  A.order = O
  return A
end

@doc Markdown.doc"""
    *(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    *(O::AlgAssAbsOrd, a::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

Returns $a \cdot O$ or $O \cdot a$ respectively.
It is assumed that the order of $a$ is contained in $O$.
"""
*(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) = extend(A, O, :right)
*(O::AlgAssAbsOrd, A::AlgAssAbsOrdIdl) = extend(A, O, :left)

@doc Markdown.doc"""
    contract(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    intersect(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    intersect(O::AlgAssAbsOrd, a::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

Returns $a \cap O$.
It is assumed that the order of $a$ contains $O$.
"""
function contract(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd)
  # Assumes O \subseteq order(A)

  d = degree(O)
  M1 = hcat(basis_matrix(A, copy = false), basis_matrix(A, copy = false))
  M2 = hcat(FakeFmpqMat(zero_matrix(FlintZZ, d, d), fmpz(1)), basis_matrix(O, copy = false))
  M = vcat(M1, M2)
  H = sub(hnf(M, :lowerleft), 1:d, 1:d)
  return ideal(algebra(A), O, H, :nothing, true)
end

intersect(O::AlgAssAbsOrd, A::AlgAssAbsOrdIdl) = contract(A, O)
intersect(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) = contract(A, O)

# Returns I as an ideal with order O. Mostly for compatibility with NfAbsOrdIdl.
function _as_fractional_ideal_of_smaller_order(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl)
  return ideal(algebra(O), O, basis_matrix(I))
end

################################################################################
#
#  Move ideals to another algebra
#
################################################################################

# Let m: A -> B be an injection and I an ideal in B. Then this computes the
# ideal O_A\cap I, where O_A is the maximal order of A.
function _as_ideal_of_smaller_algebra(m::AbsAlgAssMor, I::AlgAssAbsOrdIdl)
  A = domain(m)
  B = codomain(m)
  @assert dim(A) <= dim(B)
  @assert algebra(I) === B
  OA = maximal_order(A)
  # Transport OA to B
  M = zero_matrix(FlintQQ, dim(B), dim(B))
  for i = 1:dim(A)
    t = m(basis_alg(OA, copy = false)[i])
    elem_to_mat_row!(M, i, t)
  end
  # Compute the intersection of M and I
  M = FakeFmpqMat(M)
  M1 = hcat(M, M)
  M2 = hcat(FakeFmpqMat(zero_matrix(FlintZZ, dim(B), dim(B)), fmpz(1)), basis_matrix(I, copy = false))
  N = sub(hnf(vcat(M1, M2), :lowerleft), 1:dim(B), 1:dim(B))
  # Map the basis to A
  basis_in_A = Vector{elem_type(A)}(undef, dim(B))
  for i = 1:dim(B)
    t = elem_from_mat_row(B, N.num, i, N.den)
    b, s = haspreimage(m, t)
    @assert b
    basis_in_A[i] = s
  end
  J = ideal_from_lattice_gens(A, OA, basis_in_A)
  return J
end

# Let m: A -> B be an injection and O an order in B. Then this computes the
# order O_A\cap O, where O_A is the maximal order of A.
function _as_order_of_smaller_algebra(m::AbsAlgAssMor, O::AlgAssAbsOrd, OB::AlgAssAbsOrd)
  A = domain(m)
  B = codomain(m)
  @assert dim(A) <= dim(B)
  @assert algebra(O) == B
  OA = maximal_order(A)
  # Transport OA to B
  M = zero_matrix(FlintZZ, dim(B), dim(B))
  for i = 1:dim(A)
    t = OB(m(basis_alg(OA, copy = false)[i]))
    elem_to_mat_row!(M, i, t)
  end
  # Compute the intersection of M and O (in OB)
  N = basis_matrix(O, copy = false)*basis_mat_inv(OB, copy = false)
  @assert N.den == 1
  H = vcat(M, N.num)
  K = left_kernel(H)[2]
  N = sub(K, 1:dim(B), 1:dim(B))*M
  # Map the basis to A
  basis_in_A = Vector{elem_type(A)}(undef, dim(B))
  for i = 1:dim(B)
    t = elem_from_mat_row(OB, N, i)
    b, s = haspreimage(m, elem_in_algebra(t, copy = false))
    @assert b
    basis_in_A[i] = s
  end
  M = zero_matrix(base_ring(A), length(basis_in_A), dim(A))
  for i = 1:length(basis_in_A)
    elem_to_mat_row!(M, i, basis_in_A[i])
  end
  M = FakeFmpqMat(M)
  MM = sub(hnf!(M, :lowerleft), (dim(B) - dim(A) + 1):dim(B), 1:dim(A))
  OO = Order(A, MM)
  return OO
end

# Let m: A -> B be an injection, I an ideal in A and O an order of B. Then this
# computes I*O.
# WE ASSUME THAT m(A) LIES IN THE CENTRE OF B!
function _as_ideal_of_larger_algebra(m::AbsAlgAssMor, I::AlgAssAbsOrdIdl, O::AlgAssAbsOrd)
  A = domain(m)
  B = codomain(m)
  @assert dim(A) <= dim(B)
  @assert algebra(I) == A
  @assert algebra(O) == B

  M = zero_matrix(FlintQQ, dim(A)*dim(B), dim(B))
  X = elem_type(B)[ m(b) for b in basis(I, copy = false) ]
  for i = 1:dim(A)
    for j = 1:dim(B)
      t = X[i]*basis_alg(O, copy = false)[j]
      elem_to_mat_row!(M, (i - 1)*dim(B) + j, t)
    end
  end
  M = sub(hnf(FakeFmpqMat(M), :lowerleft), dim(B)*(dim(A) - 1) + 1:dim(A)*dim(B), 1:dim(B))
  return ideal(B, O, M, :nothing, true)
end

################################################################################
#
#  Hashing
#
################################################################################

function Base.hash(A::AlgAssAbsOrdIdl, h::UInt)
  return Base.hash(basis_matrix(A, copy = false), h)
end

################################################################################
#
#  Idempotents
#
################################################################################

@doc Markdown.doc"""
    idempotents(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl)
      -> AlgAssAbsOrdElem, AlgAssAbsOrdElem

Given coprime integral ideals $a$ and $b$, this function returns elements
$x \in a$ and $y \in b$ such that $x + y = 1$.
"""
function idempotents(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl)
  !(order(a) === order(b)) && error("Orders do not coincide")
  O = order(a)
  d = degree(O)

  V = zero_matrix(FlintZZ, 1 + 2*degree(O), 1 + 2*degree(O))
  V[1, 1] = 1
  u = coordinates(one(O))
  for i = 1:d
    V[1, i + 1] = u[i]
  end

  _copy_matrix_into_matrix(V, 2, 2, integral_basis_matrix_wrt(a, O, copy = false))
  _copy_matrix_into_matrix(V, 2 + d, 2, integral_basis_matrix_wrt(b, O, copy = false))
  for i = 2:d + 1
    V[i, d + i] = 1
  end

  H = hnf!(V)
  for i = 2:d + 1
    if H[1, i] != 0
      error("Ideals are not coprime")
    end
  end

  v = sub(H, 1:1, d + 2:2*d + 1)*integral_basis_matrix_wrt(a, O, copy = false)
  z = O([ v[1, i] for i = 1:d ])

  @assert -z in a
  @assert one(O) + z in b
  return -z, one(O) + z
end

################################################################################
#
#  From algebra to number field
#
################################################################################

function _as_ideal_of_number_field(I::AlgAssAbsOrdIdl, m::AbsAlgAssToNfAbsMor)
  @assert algebra(order(I)) == domain(m)
  K = codomain(m)
  OK = maximal_order(K)

  b = Vector{elem_type(OK)}()
  for i = 1:dim(domain(m))
    push!(b, OK(m(basis(I, copy = false)[i])))
  end
  return ideal_from_z_gens(OK, b)
end

function _as_ideal_of_number_field(I::FacElem{ <: AlgAssAbsOrdIdl, <: AlgAssAbsOrdIdlSet}, m::AbsAlgAssToNfAbsMor)
  @assert algebra(order(base_ring(parent(I)))) == domain(m)
  K = codomain(m)
  OK = maximal_order(K)

  if isempty(I)
    return FacElemMon(IdealSet(OK))()
  end

  bases = Vector{ideal_type(OK)}()
  exps = Vector{fmpz}()
  for (b, e) in I
    push!(bases, _as_ideal_of_number_field(b, m))
    push!(exps, e)
  end
  return FacElem(bases, exps)
end

function _as_ideal_of_algebra(I::NfAbsOrdIdl, i::Int, O::AlgAssAbsOrd, one_ideals::Vector{Vector{T}}) where { T <: AbsAlgAssElem }
  A = algebra(O)
  fields_and_maps = as_number_fields(A)
  b = Vector{elem_type(A)}()
  for j = 1:length(fields_and_maps)
    K, AtoK = fields_and_maps[j]
    if j == i
      @assert nf(order(I)) == K
      append!(b, [ AtoK\K(bb) for bb in basis(I, copy = false) ])
    else
      append!(b, one_ideals[j])
    end
  end
  return ideal_from_lattice_gens(algebra(O), O, b, :twosided)
end

function _as_ideal_of_algebra(I::FacElem{NfOrdIdl, NfOrdIdlSet}, i::Int, O::AlgAssAbsOrd, one_ideals::Vector{Vector{T}}) where { T <: AbsAlgAssElem }
  if isempty(I)
    return FacElemMon(IdealSet(O))()
  end

  base = Vector{ideal_type(O)}()
  exp = Vector{fmpz}()
  for (b, e) in I
    push!(base, _as_ideal_of_algebra(b, i, O, one_ideals))
    push!(exp, e)
  end
  return FacElem(base, exp)
end

# Returns an array of bases of the ideals O_i(1)*O_i lifted to O, where O_i
# are the maximal orders of the number fields in the decomposition of
# algebra(O).
function _lift_one_ideals(O::AlgAssAbsOrd)
  A = algebra(O)
  fields_and_maps = as_number_fields(A)

  one_ideals = Vector{Vector{elem_type(A)}}()
  for i = 1:length(fields_and_maps)
    K, AtoK = fields_and_maps[i]
    OK = maximal_order(K)
    a = OK(1)*OK
    b = Vector{elem_type(A)}()
    for i = 1:degree(K)
      push!(b, AtoK\K(basis(a, copy = false)[i]))
    end
    push!(one_ideals, b)
  end
  return one_ideals
end

################################################################################
#
#  Factorization
#
################################################################################

@doc Markdown.doc"""
    factor(I::AlgAssAbsOrdIdl) -> Dict{AlgAssAbsOrdIdl, Int}

Assumes that `algebra(A)` is commutative.
Returns the prime ideal factorization of $I$ as a dictionary.
"""
function factor(I::AlgAssAbsOrdIdl)
  @assert iscommutative(algebra(I))
  O = order(I)
  A = algebra(O)
  fields_and_maps = as_number_fields(A)

  one_ideals = _lift_one_ideals(O)

  fac = Dict{ideal_type(O), Int}()
  for i = 1:length(fields_and_maps)
    K, AtoK = fields_and_maps[i]
    J = _as_ideal_of_number_field(I, AtoK)
    fac2 = factor(J)
    for (p, e) in fac2
      P = _as_ideal_of_algebra(p, i, O, one_ideals)
      fac[P] = e
    end
  end
  return fac
end

################################################################################
#
#  Normal ideals
#
################################################################################

@doc Markdown.doc"""
    isnormal(a::AlgAssAbsOrdIdl) -> Bool

Returns `true` if $a$ is a normal ideal and `false` otherwise.
"""
isnormal(a::AlgAssAbsOrdIdl) = ismaximal(left_order(a))

################################################################################
#
#  Minimum
#
################################################################################

# Only for prime ideals so far...
function minimum(P::AlgAssAbsOrdIdl)
  N = norm(P, copy = false)
  @assert isone(denominator(N))
  N = numerator(N)
  f, p = ispower(N)
  @assert isprime(p)
  return p
end

################################################################################
#
#  Random elements
#
################################################################################

RandomExtensions.maketype(P::AlgAssAbsOrdIdl, ::Int) = elem_type(algebra(P))

function rand(rng::AbstractRNG,
              sp::SamplerTrivial{<:Make2{<:RingElem,<:AlgAssAbsOrdIdl,Int}})
  a, B = sp[][1:end]
  r = rand(rng, -B:B, dim(algebra(a)))
  b = basis(a, copy = false)
  z = algebra(a)()
  t = algebra(a)()
  for i = 1:dim(algebra(a))
    t = mul!(t, r[i], b[i])
    z = add!(z, z, t)
  end
  return z
end

rand(a::AlgAssAbsOrdIdl, B::Int) = rand(Random.GLOBAL_RNG, a, B)
rand(rng::AbstractRNG, a::AlgAssAbsOrdIdl, B::Int) = rand(rng, make(a, B))

################################################################################
#
#  Coprime representative
#
################################################################################

# Returns x \in A with Ix + a*O == O and Ix \subseteq O.
function integral_coprime_representative(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl, a::fmpz)
  A = algebra(O)
  d = denominator(I, O)
  I = d*I

  if one(O) in I + a*O
    return one(A)*d
  end

  fac_a = factor(a)
  primes = collect(keys(fac_a.fac))
  q = prod(primes)

  x = A()
  for i = 1:length(primes)
    p = primes[i]
    z = divexact(q, p)
    if one(O) in I + p*O
      # I is already coprime to this prime
      x += A(z)
      continue
    end
    b, g = islocally_free(O, I, p)
    @assert b "No local generator found for $p"
    ig = inv(elem_in_algebra(g, copy = false))
    Ig = I*ig
    y = denominator(Ig, O)
    @assert valuation(y, p) == 0
    x += ig*y*z
  end
  return x*d
end

################################################################################
#
#  Maximal integral ideals
#
################################################################################

# Computes any maximal integral ideal with left order O (if side = :left) or
# right order O (if side = :right) which contains p.
# Assumes (so far?) that the algebra is simple and O is maximal.
function maximal_integral_ideal(O::AlgAssAbsOrd, p::Union{ fmpz, Int }, side::Symbol)
  A = algebra(O)
  @assert issimple(A)
  @assert ismaximal(O)

  P = prime_ideals_over(O, p)[1] # if the algebra is simple, then there is exactly one prime lying over p

  # P is the Jacobson radical of O/pO, so O/P is a simple algebra
  B, OtoB = quo(O, P, p)
  C, BtoC, CtoB = _as_algebra_over_center(B)
  D, CtoD = _as_matrix_algebra(C)

  n = degree(D)
  if isone(n)
    return P
  end

  N = fmpq_mat(basis_matrix(P))
  t = zero_matrix(FlintQQ, 1, degree(O))
  # Now we only need to lift a basis for diag(1, ..., 1, 0)*D (side = :left) or
  # D*diag(1, ..., 1, 0) (side = :right) since these are maximal ideals of D.
  if side == :left
    jMax = n - 1
    iMax = n
  elseif side == :right
    jMax = n
    iMax = n - 1
  else
    error("Option :$(side) for side not implemented")
  end
  for j = 1:jMax
    jn = (j - 1)*n
    for i = 1:iMax
      b = (OtoB\(CtoB(CtoD\D[jn + i])))
      elem_to_mat_row!(t, 1, elem_in_algebra(b, copy = false))
      N = vcat(N, deepcopy(t))
    end
  end
  N = FakeFmpqMat(N)
  N = sub(hnf(N, :lowerleft), nrows(N) - degree(O) + 1:nrows(N), 1:degree(O))

  M = ideal(algebra(O), O, N, side, true)
  if side == :left
    M.left_order = O # O is maximal
  else
    M.right_order = O
  end
  return M
end

# Constructs a maximal integral ideal M of O such that M\cap R = p and I\subseteq M,
# where O is the left order (if side = :left) or right order (if side = :right)
# of I. It is assumed that I \subseteq O.
# M is a left ideal of O if side = :left and a right ideal if side = :right.
# Assumes (so far?) that the algebra is simple and O is maximal.
function maximal_integral_ideal_containing(I::AlgAssAbsOrdIdl, p::Union{ fmpz, Int }, side::Symbol)
  if side == :left
    O = left_order(I)
  elseif side == :right
    O = right_order(I)
  else
    error("Option :$(side) for side not implemented")
  end

  @assert issimple(algebra(O))
  @assert ismaximal(O)

  n = normred(I, O)
  if valuation(n, p) == 0
    error("Cannot find a maximal ideal for the given prime")
  end
  if n == p
    # The ideal is maximal iff its reduced norm is prime
    return I
  end

  P = prime_ideals_over(O, p)[1] # if the algebra is simple, then there is exactly one prime lying over p
  J = I + P
  if normred(J, O) == p
    return J
  end

  OP, toOP = quo(O, P, p)
  B, OPtoB, BtoOP = _as_algebra_over_center(OP)
  C, toC = _as_matrix_algebra(B)

  JinC = ideal_from_gens(C, [ toC(OPtoB(toOP(O(b)))) for b in absolute_basis(J) ])
  y = left_principal_generator(JinC)
  m = matrix(y)
  r = rref!(m)
  k = degree(C)
  @assert r < k - 1 # Otherwise J would be maximal, which we have checked...

  # We need to "add" pivots to m to get rank k - 1
  row = 1
  col = 1
  while r < k - 1
    if !iszero(m[row, col])
      row += 1
      continue
    end

    r += 1
    m[r, col] = 1
    col += 1
  end
  @assert rank(m) == k - 1
  if side == :left
    c = C*C(m)
  else
    c = C(m)*C
  end
  basis_c = basis(c, copy = false)

  t = zero_matrix(FlintQQ, length(basis_c), degree(O))
  for i = 1:length(basis_c)
    b = toOP\(BtoOP(toC\(basis_c[i])))
    elem_to_mat_row!(t, i, elem_in_algebra(b, copy = false))
  end
  m = hnf(vcat(basis_matrix(P), FakeFmpqMat(t)), :lowerleft)
  m = sub(m, length(basis_c) + 1:nrows(m), 1:ncols(m))

  M = ideal(algebra(O), O, m, side, true)
  @assert normred(M, O) == p
  if side == :left
    M.left_order = O # O is maximal
  else
    M.right_order = O
  end

  return M
end


text:
"""
A single field FE space with transient Dirichlet data (see Multifield below).
"""
struct TransientTrialFESpace
  space::SingleFieldFESpace
  dirichlet_t::Union{Function,Vector{<:Function}}
  Ud0::TrialFESpace

  function TransientTrialFESpace(space::SingleFieldFESpace,dirichlet_t::Union{Function,Vector{<:Function}})
    Ud0 = HomogeneousTrialFESpace(space)
    new(space,dirichlet_t,Ud0)
  end
end

function TransientTrialFESpace(space::SingleFieldFESpace)
  HomogeneousTrialFESpace(space)
end

"""
Time evaluation without allocating Dirichlet vals
"""
function evaluate!(Ut::TrialFESpace,U::TransientTrialFESpace,t::Real)
  if isa(U.dirichlet_t,Vector)
    objects_at_t = map( o->o(t), U.dirichlet_t)
  else
    objects_at_t = U.dirichlet_t(t)
  end
  TrialFESpace!(Ut,objects_at_t)
  Ut
end

"""
Allocate the space to be used as first argument in evaluate!
"""
function allocate_trial_space(U::TransientTrialFESpace)
    HomogeneousTrialFESpace(U.space)
end

"""
Time evaluation allocating Dirichlet vals
"""
function evaluate(U::TransientTrialFESpace,t::Real)
  Ut = allocate_trial_space(U)
  evaluate!(Ut,U,t)
  Ut
end

"""
We can evaluate at `nothing` when we do not care about the Dirichlet vals
"""
function evaluate(U::TransientTrialFESpace,t::Nothing)
  U.Ud0
end

evaluate(U::TrialFESpace,t::Nothing) = U

"""
Functor-like evaluation. It allocates Dirichlet vals in general.
"""
(U::TransientTrialFESpace)(t) = evaluate(U,t)

(U::TrialFESpace)(t) = U
(U::ZeroMeanFESpace)(t) = U
# (U::Union{TrialFESpace,ZeroMeanFESpace})(t) = U

"""
Time derivative of the Dirichlet functions
"""
∂t(U::TransientTrialFESpace) = TransientTrialFESpace(U.space,∂t.(U.dirichlet_t))

# ∂t(U::TrialFESpace) = TransientTrialFESpace(U.space,∂t.(U.dirichlet_t))
∂t(U::SingleFieldFESpace) = HomogeneousTrialFESpace(U)

∂t(U::MultiFieldFESpace) = MultiFieldFESpace(∂t.(U.spaces))

∂t(t::T) where T<:Number = zero(T)

# Testing the interface

function test_transient_trial_fe_space(Uh)
  UhX = evaluate(Uh,nothing)
  @test isa(UhX,FESpace)
  Uh0 = allocate_trial_space(Uh)
  Uh0 = evaluate!(Uh0,Uh,0.0)
  @test isa(Uh0,FESpace)
  Uh0 = evaluate(Uh,0.0)
  @test isa(Uh0,FESpace)
  Uh0 = Uh(0.0)
  @test isa(Uh0,FESpace)
  Uht=∂t(Uh)
  Uht0=Uht(0.0)
  @test isa(Uht0,FESpace)
  true
end

# Define the TransientTrialFESpace interface for stationary spaces

function evaluate!(Ut::FESpace,U::FESpace,t::Real)
  U
end

function allocate_trial_space(U::FESpace)
  U
end

function evaluate(U::FESpace,t::Real)
  U
end

function evaluate(U::FESpace,t::Nothing)
  U
end

@static if VERSION >= v"1.3"
  (U::FESpace)(t) = U
end

# Define the interface for MultiField

struct TransientMultiFieldTrialFESpace
  spaces::Vector
end

function TransientMultiFieldFESpace(spaces::Vector)
  TransientMultiFieldTrialFESpace(spaces)
end

function TransientMultiFieldFESpace(spaces::Vector{<:SingleFieldFESpace})
  MultiFieldFESpace(spaces)
end

function evaluate!(Ut::MultiFieldFESpace,U::TransientMultiFieldTrialFESpace,t::Real)
  spaces_at_t = [evaluate!(Ut.spaces[i],U.spaces[i],t) for i in 1:length(U.spaces)]
  MultiFieldFESpace(spaces_at_t)
end

function allocate_trial_space(U::TransientMultiFieldTrialFESpace)
  spaces = allocate_trial_space.(U.spaces)
  MultiFieldFESpace(spaces)
end

function evaluate(U::TransientMultiFieldTrialFESpace,t::Real)
  Ut = allocate_trial_space(U)
  evaluate!(Ut,U,t)
  Ut
end

function evaluate(U::TransientMultiFieldTrialFESpace,t::Nothing)
  MultiFieldFESpace([evaluate(fesp,nothing) for fesp in U.spaces])
end

(U::TransientMultiFieldTrialFESpace)(t) = evaluate(U,t)

function ∂t(U::TransientMultiFieldTrialFESpace)
  spaces = ∂t.(U.spaces)
  TransientMultiFieldFESpace(spaces)
end

ext: jl


text:
module ChainRulesTestUtils

using ChainRulesCore
using ChainRulesCore: frule, rrule
using ChainRulesCore: AbstractDifferential
using FiniteDifferences
using LinearAlgebra
using Test

const _fdm = central_fdm(5, 1)

export test_scalar, frule_test, rrule_test, isapprox, generate_well_conditioned_matrix

Base.isapprox(d_ad::DoesNotExist, d_fd; kwargs...) = error("Tried to differentiate w.r.t. a `DoesNotExist`")
Base.isapprox(d_ad::AbstractDifferential, d_fd; kwargs...) = isapprox(extern(d_ad), d_fd; kwargs...)

function _make_fdm_call(fdm, f, ȳ, xs, ignores)
    sig = Expr(:tuple)
    call = Expr(:call, f)
    newxs = Any[]
    arginds = Int[]
    i = 1
    for (x, ignore) in zip(xs, ignores)
        if ignore
            push!(call.args, x)
        else
            push!(call.args, Symbol(:x, i))
            push!(sig.args, Symbol(:x, i))
            push!(newxs, x)
            push!(arginds, i)
        end
        i += 1
    end
    fdexpr = :(j′vp($fdm, $sig -> $call, $ȳ, $(newxs...)))
    fd = eval(fdexpr)
    fd isa Tuple || (fd = (fd,))
    args = Any[nothing for _ in 1:length(xs)]
    for (dx, ind) in zip(fd, arginds)
        args[ind] = dx
    end
    return (args...,)
end

# Useful for LinearAlgebra tests
function generate_well_conditioned_matrix(rng, N)
    A = randn(rng, N, N)
    return A * A' + I
end

"""
    test_scalar(f, x; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

Given a function `f` with scalar input and scalar output, perform finite differencing checks,
at input point `x` to confirm that there are correct `frule` and `rrule`s provided.

# Arguments
- `f`: Function for which the `frule` and `rrule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).

All keyword arguments except for `fdm` is passed to `isapprox`.
"""
function test_scalar(f, x; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    ensure_not_running_on_functor(f, "test_scalar")

    r_res = rrule(f, x)
    f_res = frule(f, x, Zero(), 1)
    @test r_res !== nothing  # Check the rule was defined
    @test f_res !== nothing
    r_fx, prop_rule = r_res
    f_fx, f_∂x = f_res
    @testset "$f at $x, $(nameof(rule))" for (rule, fx, ∂x) in (
        (rrule, r_fx, prop_rule(1)),
        (frule, f_fx, f_∂x)
    )
        @test fx == f(x)  # Check we still get the normal value, right

        if rule == rrule
            ∂self, ∂x = ∂x
            @test ∂self === NO_FIELDS
        end
        @test isapprox(∂x, fdm(f, x); rtol=rtol, atol=atol, kwargs...)
    end
end

function ensure_not_running_on_functor(f, name)
    # if x itself is a Type, then it is a constructor, thus not a functor.
    # This also catchs UnionAll constructors which have a `:var` and `:body` fields
    f isa Type && return

    if fieldcount(typeof(f)) > 0
        throw(ArgumentError(
            "$name cannot be used on closures/functors (such as $f)"
        ))
    end
end

"""
    frule_test(f, (x, ẋ)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function for which the `frule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `ẋ`: differential w.r.t. `x` (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function frule_test(f, (x, ẋ); rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    return frule_test(f, ((x, ẋ),); rtol=rtol, atol=atol, fdm=fdm, kwargs...)
end

function frule_test(f, xẋs::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    ensure_not_running_on_functor(f, "frule_test")
    xs, ẋs = collect(zip(xẋs...))
    Ω, dΩ_ad = frule(f, xs..., NO_FIELDS, ẋs...)
    @test f(xs...) == Ω

    # Correctness testing via finite differencing.
    dΩ_fd = jvp(fdm, xs->f(xs...), (xs, ẋs))
    @test isapprox(
        collect(extern.(dΩ_ad)),  # Use collect so can use vector equality
        collect(dΩ_fd);
        rtol=rtol,
        atol=atol,
        kwargs...
    )
end

"""
    rrule_test(f, ȳ, (x, x̄)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function to which rule should be applied.
- `ȳ`: adjoint w.r.t. output of `f` (should generally be set randomly).
  Should be same structure as `f(x)` (so if multiple returns should be a tuple)
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `x̄`: currently accumulated adjoint (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function rrule_test(f, ȳ, (x, x̄)::Tuple{Any, Any}; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    ensure_not_running_on_functor(f, "rrule_test")

    # Check correctness of evaluation.
    fx, pullback = rrule(f, x)
    @test collect(fx) ≈ collect(f(x))  # use collect so can do vector equality
    (∂self, x̄_ad) = if fx isa Tuple
        # If the function returned multiple values,
        # then it must have multiple seeds for propagating backwards
        pullback(ȳ...)
    else
        pullback(ȳ)
    end

    @test ∂self === NO_FIELDS  # No internal fields
    # Correctness testing via finite differencing.
    x̄_fd = j′vp(fdm, f, ȳ, x)
    @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)
end

# case where `f` takes multiple arguments
function rrule_test(f, ȳ, xx̄s::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    ensure_not_running_on_functor(f, "rrule_test")

    # Check correctness of evaluation.
    xs, x̄s = collect(zip(xx̄s...))
    y, pullback = rrule(f, xs...)
    @test f(xs...) == y

    @assert !(isa(ȳ, Thunk))
    ∂s = pullback(ȳ)
    ∂self = ∂s[1]
    x̄s_ad = ∂s[2:end]
    @test ∂self === NO_FIELDS

    # Correctness testing via finite differencing.
    x̄s_fd = j′vp(fdm, f, ȳ, xs...)
    map(x̄s_ad, x̄s_fd) do x̄_ad, x̄_fd
        @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)
    end
end

end # module


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions to verify if the satellite is within the station visibility
#   circle.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export ground_station_accesses, ground_station_gaps,
       list_ground_station_accesses, list_ground_station_gaps,
       ground_station_visible

"""
    ground_station_accesses(orbp, vrs_e,     Δt, ECI, ECEF, vargs...; kwargs...)
    ground_station_accesses(orbp, [(WGS84)], Δt, ECI, ECEF, vargs...; kwargs...)

Compute the accesses of a satellite with orbit propagator `orbp` (see
`init_orbit_propagator`) to the ground stations defined in the vector `vrs_e`.
The analysis interval begins in the propagator epoch and lasts `Δt` [s].

The ground stations can be specified by an array of 3×1 vectors describing the
ground stations position in an ECEF frame `vrs_e` or by an array of tuples
containing the WGS84 position of each ground station `[(WGS84)]`:

    (latitude [rad], longitude [rad], altitude [m])

# Args

* `ECI`: Earth-Centered Inertial frame in which the state vector of the
         propagator is represented.
* `ECEF`: Earth-Centered, Earth-fixed frame to be used for the analysis. It
          must be the same frame used to compute the ground station position
          vector.
* `vargs...`: list of additional arguments to be passed to the function
              `r_eci_to_ecef` when converting the ECI frame to the ECEF.

# Keywords

* `θ`: Minimum elevation angle for communication between the satellite and the
       ground stations [rad]. (**Default** = 10ᵒ)
* `reduction`: A function that receives a boolean vector with the visibility
               between the satellite and each ground station. It must return a
               boolean value indicating if the access must be computed or not.
               This is useful to merge access time between two or more stations.
               (**Default** = `v->|(v...)` *i.e.* compute the access if at least
               one ground station is visible)

"""
ground_station_accesses(orbp, gs_wgs84::Tuple, vargs...; kwargs...) =
    ground_station_accesses(orbp, [gs_wgs84], vargs...; kwargs...)

function ground_station_accesses(orbp, vgs_wgs84::AbstractVector{T}, vargs...;
                                 kwargs...) where T<:Tuple

    vrs_e = [geodetic_to_ecef(gs_wgs84...) for gs_wgs84 in vgs_wgs84]
    return ground_station_accesses(orbp, vrs_e, vargs...; kwargs...)
end

ground_station_accesses(orbp, rs_e::AbstractVector{T}, vargs...; kwargs...) where
    T<:Number = ground_station_accesses(orbp, [rs_e], vargs...; kwargs...)

function ground_station_accesses(orbp, vrs_e::AbstractVector{T}, Δt::Number,
                                 ECI::Union{T_ECIs, T_ECIs_IAU_2006},
                                 ECEF::Union{T_ECEFs, T_ECEFs_IAU_2006}, vargs...;
                                 θ::Number = 10*pi/180,
                                 reduction::Function = v->|(v...),
                                 step::Number = 60) where T<:AbstractVector

    # Time vector of the analysis.
    t = 0:step:Δt

    # Get the epoch of the propagator.
    JD₀ = get_epoch(orbp)

    # Matrix that will contain the accesses.
    accesses = Matrix{DateTime}(undef,0,2)

    # State to help the computation.
    state = :initial

    # Lambda function to check if the ground station is visible.
    f(t)::Bool = begin
        r_i, v_i   = propagate!(orbp, t)
        r_e        = r_eci_to_ecef(DCM, ECI, ECEF, JD₀ + t/86400, vargs...)*r_i
        visibility = [ground_station_visible(r_e, rs_e, θ) for rs_e in vrs_e]

        return reduction(visibility)
    end

    access_beg = DateTime(now())
    access_end = DateTime(now())

    for k in t
        # Check if the ground station is visible.
        gs_visible = f(k)

        # Handle the initial case.
        if state == :initial
            if gs_visible
                access_beg = jd_to_date(DateTime, JD₀)
                state = :visible
            else
                state = :not_visible
            end
        # Handle transitions.
        elseif (state == :not_visible) && gs_visible
            # Refine to find the edge.
            k₀ = k - step
            k₁ = k
            kc = find_crossing(f, k₀, k₁, false, true)

            state = :visible
            access_beg = jd_to_date(DateTime, JD₀ + kc/86400)

        elseif (state == :visible) && !gs_visible
            # Refine to find the edge.
            k₀ = k - step
            k₁ = k
            kc = find_crossing(f, k₀, k₁, true, false)

            state = :not_visible
            access_end = jd_to_date(DateTime, JD₀ + kc/86400)

            accesses = vcat(accesses, [access_beg access_end])
        end
    end

    # If the analysis finished during an access, then just add the end of the
    # interval as the end of the access.
    if state == :visible
        access_end = jd_to_date(DateTime, JD₀ + Δt/86400)
        accesses   = vcat(accesses, [access_beg access_end])
    end

    return accesses
end

"""
    ground_station_gaps(args...; kwargs...)

Compute the gaps between the accesses of ground stations. The arguments and
keywords are the same as the ones used in the function
`ground_station_accesses`.

Notice that the gap analysis starts in the orbit propagator epoch and ends in
the instant defined by the argument `Δt`.

"""
function ground_station_gaps(orbp, args...; kwargs...)
    # Get the epoch of the propagator.
    JD₀ = get_epoch(orbp)
    DT₀ = jd_to_date(DateTime, JD₀)

    # Compute the list of ground station accesses.
    accesses = ground_station_accesses(orbp, args...; kwargs...)

    # Get the analysis period, which must be the second argument in `args...`.
    Δt = args[2]

    # Compute the last propagation instant.
    JD₁ = JD₀ + Δt/86400
    DT₁ = jd_to_date(DateTime, JD₁)

    # Compute the gaps between accesses.
    gaps = Matrix{DateTime}(undef,0,2)

    # Check if the simulation did not start under the visibility of a ground
    # station.
    if accesses[1,1] != DT₀
        gaps = vcat(gaps, [jd_to_date(DateTime, JD₀) accesses[1,1]])
    end

    # Check if the simulation did not end under the visibility of a ground
    # station.
    if accesses[end,2] != DT₁
        aux  = jd_to_date(DateTime, JD₁)
        gaps = vcat(gaps, [ accesses[:,2] vcat(accesses[2:end,1], aux) ])
    else
        gaps = vcat(gaps, [ accesses[1:end-1,2] accesses[2:end,1] ])
    end

    return gaps
end

"""
    list_ground_station_accesses(io, vargs...; kwargs...)

Print the ground station accesses to the io `io`. The arguments `vargs...` and
keywords `kwargs...` are those of the function `ground_station_accesses`.

Additionally, the following keywords can be used to modify the behavior of this
function:

* `format`: If `:pretty`, then a formatted table will be printed. If `:csv`,
            then the access data will be printed using the CSV format.
            (**Default** = `:pretty`)
* `time_scale`: Select the time scale of the access duration (`:s` for seconds,
                `:m` for minutes, and `:h` for hours). (**Default** = `:m`)

"""
list_ground_station_accesses(vargs...; kwargs...) =
    list_ground_station_accesses(stdout, vargs...; kwargs...)

function list_ground_station_accesses(io::IO, vargs...; format = :pretty,
                                      time_scale::Symbol = :m, kwargs...)

    # Compute the accesses.
    accesses = ground_station_accesses(vargs...; kwargs...)

    # If no access was found, then inform the user and return.
    if isempty(accesses)
        println("No access found!")
        return nothing
    end

    # Check which scale to use when printing the access duration.
    if time_scale == :s
        ts    = 1
        label = "s"
    elseif time_scale == :m
        ts    = 60
        label = "min"
    else
        ts    = 3600
        label = "h"
    end

    # Compute the accesses duration.
    access_time = map((b,e)->(e-b).value/1000/ts,
                       @view(accesses[:,1]), @view(accesses[:,2]))
    num         = length(access_time)

    # Assemble the information to be printed.
    header = ["Access #", "Beginning (UTC)", "End (UTC)", "Duration [$label]"]
    mp     = hcat(string.(collect(1:1:num)), accesses, access_time)

    # Print the information depending on the format.

    if format == :csv
        writedlm(io, vcat(header,mp), ",")
    else
        # Statistics.
        min_access, id_min = findmin(access_time)
        max_access, id_max = findmax(access_time)
        mean_access        = mean(access_time)
        total_access       = sum(access_time)

        # Mark the minimum and maximum accesses.
        mp[id_min, 1] = "MIN " * mp[id_min,1]
        mp[id_max, 1] = "MAX " * mp[id_max,1]

        # Add the statistics to the end of the table.
        mp = vcat(mp, ["" "" "         Minimum access"  min_access;
                       "" "" "            Mean access"  mean_access;
                       "" "" "         Maximum access"  max_access;
                       "" "" "           Total access"  total_access;])

        # Highlighters for the minimum and maximum accesses.
        hl_min = Highlighter((data,i,j)->i == id_min, crayon"bold red")
        hl_max = Highlighter((data,i,j)->i == id_max, crayon"bold blue")

        # Highlighter for statistics labels.
        hl = Highlighter((data,i,j)->(j == 3) && (i-num ∈ (1,2,3,4)),
                         crayon"bold")

        # Print.
        pretty_table(io, mp;
                     alignment    = [:r, :l, :l, :r],
                     body_hlines  = [num,num+3],
                     crop         = :horizontal,
                     header       = header,
                     formatters   = ft_printf("%.3f",[4]),
                     highlighters = (hl,hl_min,hl_max))
    end

    return nothing
end

"""
    list_ground_station_gaps(io, vargs...; kwargs...)

Print the ground station gaps to the io `io`. The arguments `vargs...` and
keywords `kwargs...` are those of the function `ground_station_gaps`.

Additionally, the following keywords can be used to modify the behavior of this
function:

* `format`: If `:pretty`, then a formatted table will be printed. If `:csv`,
            then the access data will be printed using the CSV format.
            (**Default** = `:pretty`)
* `time_scale`: Select the time scale of the access duration (`:s` for seconds,
                `:m` for minutes, and `:h` for hours). (**Default** = `:m`)

"""
list_ground_station_gaps(vargs...; kwargs...) =
    list_ground_station_gaps(stdout, vargs...; kwargs...)

function list_ground_station_gaps(io::IO, vargs...; format = :pretty,
                                  time_scale::Symbol = :m, kwargs...)

    # Compute the gaps.
    gaps = ground_station_gaps(vargs...; kwargs...)

    # If no gap was found, then inform the user and return.
    if isempty(gaps)
        println("No gap found!")
        return nothing
    end

    # Check which scale to use when printing the gap duration.
    if time_scale == :s
        ts    = 1
        label = "s"
    elseif time_scale == :m
        ts    = 60
        label = "min"
    else
        ts    = 3600
        label = "h"
    end

    # Compute the gaps duration.
    gap_time = map((b,e)->(e-b).value/1000/ts, @view(gaps[:,1]), @view(gaps[:,2]))
    num         = length(gap_time)

    # Assemble the information to be printed.
    header = ["   Gap #", "Beginning (UTC)", "End (UTC)", "Duration [$label]"]
    mp     = hcat(string.(collect(1:1:num)), gaps, gap_time)

    # Print the information depending on the format.

    if format == :csv
        writedlm(io, vcat(header,mp), ",")
    else
        # Statistics.
        min_gap, id_min = findmin(gap_time)
        max_gap, id_max = findmax(gap_time)
        mean_gap        = mean(gap_time)
        total_gap       = sum(gap_time)

        # Mark the minimum and maximum gaps.
        mp[id_min, 1] = "MIN " * mp[id_min,1]
        mp[id_max, 1] = "MAX " * mp[id_max,1]

        # Add the statistics to the end of the table.
        mp = vcat(mp, ["" "" "            Minimum gap"  min_gap;
                       "" "" "               Mean gap"  mean_gap;
                       "" "" "            Maximum gap"  max_gap;
                       "" "" "              Total gap"  total_gap;])

        # Highlighters for the minimum and maximum gaps.
        hl_min = Highlighter((data,i,j)->i == id_min, crayon"bold red")
        hl_max = Highlighter((data,i,j)->i == id_max, crayon"bold blue")

        # Highlighter for statistics labels.
        hl = Highlighter((data,i,j)->(j == 3) && (i-num ∈ (1,2,3,4)),
                         crayon"bold")

        # Print.
        pretty_table(io, mp;
                     alignment    = [:r, :l, :l, :r],
                     body_hlines  = [num,num+3],
                     crop         = :horizontal,
                     header       = header,
                     formatters   = ft_printf("%.3f",[4]),
                     highlighters = (hl,hl_min,hl_max))
    end

    return nothing
end

"""
    ground_station_visible(r_e::AbstractVector, rs_e::AbstractVector, θ::Number)

Check if the satellite with position vector `r_e` (ECEF) is inside the
visibility circle of a ground station with position vector `rs_e` (ECEF) and a
minimum elevation angle of `θ` [rad].

Notice that `r_e` and `rs_e` must be represented in the same ECEF frame, and
must have the same unit.

Returns `true` if the satellite is inside the visibility circle, or `false`
otherwise.

"""
function ground_station_visible(r_e::AbstractVector, rs_e::AbstractVector,
                                θ::Number)
    # Check if the satellite is within the visibility circle of the station.
    dr_e = r_e - rs_e
    cos_beta = dot( dr_e/norm(dr_e), rs_e/norm(rs_e) )

    return cos_beta > cos(π/2-θ)
end

"""
    ground_station_visible(r_e::AbstractVector, lat_s::Number, lon_s::Number, h_s::Number, θ::Number)

Check if the satellite with position vector `r_e` (ECEF) is inside the
visibility circle of a ground station with latitude `lat_s` [rad], longitude
`lon_s` [rad], altitude `h_s` (WGS-84), and a minimum elevation angle of
`θ` [rad].

Notice that the units of `r_e` and `h_s` must be the same.

Returns `true` if the satellite is inside the visibility circle, or `false`
otherwise.

"""
function ground_station_visible(r_e::AbstractVector, lat_s::Number,
                                lon_s::Number, h_s::Number, θ::Number)
    # Convert the ground station LLA to the ECEF frame.
    rs_e = geodetic_to_ecef(lat_s, lon_s, h_s)

    return ground_station_visible(r_e, rs_e, θ)
end


text:
using RNGPool
using Random

function foo!(out::Vector{Float64}, N::Int64)
  nt = Threads.nthreads()
  M::Int64 = div(N, nt)
  Threads.@threads for i in 1:nt
    rng::RNG = getRNG()
    v::Float64 = 0.0
    for j in 1:M
      v += rand(rng)
    end
    out[i] = v / M
  end
end


text:
module NURBStoolbox_demoruns

using Test;
using Plots;
using MagMechFEM_Matlab2Julia.NURBStoolbox;

#ENV[“GKSwstype”] = “nul”; # this omits those tests
Plots.default(show=false);

# TODO: find a good way to test plot output or plotting in general

@testset "Function: Demos" begin

  @test typeof(demo4surf())<:Plots.Plot;
  @test typeof(democirc())<:Plots.Plot;
  @test typeof(democoons())<:Plots.Plot;
  @test typeof(democurve())<:Plots.Plot;
  @test typeof(democylind())<:Plots.Plot;
  @test typeof(demodegelev())<:Plots.Plot;
  @test typeof(demodercrv())<:Plots.Plot;
  @test typeof(demodersrf())<:Plots.Plot;
  @test typeof(demoellip())<:Plots.Plot;
  @test typeof(demogeom())<:Plots.Plot;
  @test typeof(demohelix())<:Plots.Plot;
  @test typeof(demokntins())<:Plots.Plot;
  @test typeof(demoline())<:Plots.Plot;
  @test typeof(demorect())<:Plots.Plot;
  @test typeof(demorevolve())<:Plots.Plot;
  @test typeof(demoruled())<:Plots.Plot;
  @test typeof(demotorus())<:Plots.Plot;

end # testset

end # module


text:
mutable struct StrategyMemory{INFO<:AbstractVector, STRAT<:AbstractVector}
    I::Vector{INFO}
    t::Vector{Float32}
    σ::Vector{STRAT}
    capacity::Int
    i::Int
end

function StrategyMemory{INFO, STRAT}(sz::Int) where {INFO,STRAT}
    return StrategyMemory(
        sizehint!(INFO[], sz),
        sizehint!(Float32[], sz),
        sizehint!(STRAT[], sz),
        sz,
        0
    )
end

function StrategyMemory{INFO}(sz::Int) where {INFO}
    return StrategyMemory{INFO, Vector{Float32}}(sz)
end

function Base.push!(Mπ::StrategyMemory, I, t, σ)
    i = (Mπ.i += 1)
    k = Mπ.capacity
    if i ≤ k
        push!(Mπ.I,I)
        push!(Mπ.t,t)
        push!(Mπ.σ,σ)
    else
        j = rand(1:i)
        if j ≤ k
            Mπ.I[j] = I
            Mπ.t[j] = t
            Mπ.σ[j] = σ
        end
    end
end

function Base.empty!(Mπ::StrategyMemory)
    Mπ.i = 0
    empty!(Mπ.I)
    empty!(Mπ.t)
    empty!(Mπ.σ)
end

mutable struct AdvantageMemory{INFO<:AbstractVector, REGRET<:AbstractVector}
    I::Vector{INFO}
    t::Vector{Float32}
    r::Vector{REGRET}
    capacity::Int
    i::Int
end

function Base.length(mem::AdvantageMemory)
    return length(mem.t)
end

function Base.length(mem::StrategyMemory)
    return length(mem.t)
end

function AdvantageMemory{INFO, REGRET}(sz::Int) where {INFO,REGRET}
    return AdvantageMemory(
        sizehint!(INFO[], sz),
        sizehint!(Float32[], sz),
        sizehint!(REGRET[], sz),
        sz,
        0
    )
end

function AdvantageMemory{INFO}(sz::Int) where {INFO}
    return AdvantageMemory{INFO, Vector{Float32}}(sz)
end

function Base.push!(Mv::AdvantageMemory, I, t, r)
    i = (Mv.i += 1)
    k = Mv.capacity
    if i ≤ k
        push!(Mv.I,I)
        push!(Mv.t,t)
        push!(Mv.r,r)
    else
        j = rand(1:i)
        if j ≤ k
            Mv.I[j] = I
            Mv.t[j] = t
            Mv.r[j] = r
        end
    end
end

function Base.empty!(Mv::AdvantageMemory)
    Mv.i = 0
    empty!(Mv.I)
    empty!(Mv.t)
    empty!(Mv.r)
end