Ideal.jl

@doc raw"""
    order(a::AlgAssRelOrdIdl) -> AlgAssRelOrd

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

@doc raw"""
    algebra(a::AlgAssRelOrdIdl) -> AbstractAssociativeAlgebra

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

_algebra(a::AlgAssRelOrdIdl) = 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.
is_full_lattice(a::AlgAssRelOrdIdl) = !iszero(basis_matrix(a, copy = false)[1, 1])

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

@doc raw"""
    ideal(A::AbstractAssociativeAlgebra, M::PMat; M_in_hnf::Bool = false) -> AlgAssRelOrdIdl

Returns the ideal in $A$ with basis pseudo-matrix $M$.
If `M_in_hnf == true`, it is assumed that $M$ is already in lower left pseudo HNF.
"""
function ideal(A::AbstractAssociativeAlgebra{S}, M::PMat{S, T}; M_in_hnf::Bool = false) where { S <: NumFieldElem, T }
  !M_in_hnf ? M = pseudo_hnf(M, :lowerleft) : nothing
  return AlgAssRelOrdIdl{S, T, typeof(A)}(A, M)
end

@doc raw"""
    ideal(A::AbstractAssociativeAlgebra, O::AlgAssRelOrd, M::PMat; side::Symbol = :nothing,
          M_in_hnf::Bool = false)
      -> AlgAssRelOrdIdl

Returns the ideal of $O$ in $A$ with basis pseudo-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 pseudo HNF.
"""
function ideal(A::AbstractAssociativeAlgebra{S}, O::AlgAssRelOrd{S, T, U}, M::PMat{S, T}; side::Symbol = :nothing, M_in_hnf::Bool = false) where { S <: NumFieldElem, T, U }
  a = ideal(A, M; M_in_hnf)
  a.order = O
  _set_sidedness(a, side)
  if is_maximal_known(O) && is_maximal(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 raw"""
    ideal(O::AlgAssRelOrd, x::AbstractAssociativeAlgebraElem) -> AlgAssRelOrdIdl
    ideal(O::AlgAssRelOrd, x::AlgAssRelOrdElem) -> AlgAssRelOrdIdl

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

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

  t1 = A()
  t2 = A()
  M = zero_matrix(base_ring(A), degree(O)^2, degree(O))
  pb = pseudo_basis(O, copy = false)
  C = Vector{T}(undef, degree(O)^2)
  for i = 1:degree(O)
    i1d = (i - 1)*degree(O)
    C[i1d + i] = pb[i][2]^2
    for j = (i + 1):degree(O)
      C[i1d + j] = simplify(pb[i][2]*pb[j][2])
      C[(j - 1)*degree(O) + i] = deepcopy(C[i1d + j])
    end
  end

  for i = 1:degree(O)
    t1 = mul!(t1, pb[i][1], x)
    ii = (i - 1)*degree(O)
    for j = 1:degree(O)
      t2 = mul!(t2, t1, pb[j][1])
      elem_to_mat_row!(M, ii + j, t2)
    end
  end
  M = sub(pseudo_hnf(pseudo_matrix(M, C), :lowerleft), nrows(M) - degree(O) + 1:nrows(M), 1:ncols(M))

  return ideal(A, O, M; side=:twosided, M_in_hnf=true)
end

ideal(O::AlgAssRelOrd{S, T, U}, x::AlgAssRelOrdElem{S, T, U}) where { S, T, U } = ideal(O, elem_in_algebra(x, copy = false))

@doc raw"""
    ideal(O::AlgAssRelOrd, x::AbstractAssociativeAlgebraElem, side::Symbol) -> AlgAssRelOrdIdl
    ideal(O::AlgAssRelOrd, x::AlgAssRelOrdElem, side::Symbol) -> AlgAssRelOrdIdl

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

  M = zero_matrix(base_ring(A), degree(O), degree(O))
  pb = pseudo_basis(O, copy = false)
  for i = 1:degree(O)
    if side == :left
      elem_to_mat_row!(M, i, pb[i][1]*x)
    elseif side == :right
      elem_to_mat_row!(M, i, x*pb[i][1])
    end
  end
  M = pseudo_matrix(M, deepcopy(basis_pmatrix(O, copy = false).coeffs))

  return ideal(A, O, M; side)
end

ideal(O::AlgAssRelOrd{S, T, U}, x::AlgAssRelOrdElem{S, T, U}, action::Symbol) where { S, T, U } = ideal(O, elem_in_algebra(x, copy = false), action)

@doc raw"""
    *(O::AlgAssRelOrd, x::AbstractAssociativeAlgebraElem) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, x::AlgAssRelOrdElem) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, x::Int) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, x::ZZRingElem) -> AlgAssRelOrdIdl
    *(x::AbstractAssociativeAlgebraElem, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(x::AlgAssRelOrdElem, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(x::Int, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(x::ZZRingElem, O::AlgAssRelOrd) -> AlgAssRelOrdIdl

Returns the ideal $O \cdot x$ or $x \cdot O$ respectively.
"""
*(O::AlgAssRelOrd{S, T, U}, x::AbstractAssociativeAlgebraElem{S}) where {S, T, U} = ideal(O, x, :left)
*(x::AbstractAssociativeAlgebraElem{S}, O::AlgAssRelOrd{S, T, U}) where {S, T, U} = ideal(O, x, :right)
*(O::AlgAssRelOrd{S, T, U}, x::AlgAssRelOrdElem{S, T, U}) where {S, T, U} = ideal(O, x, :left)
*(x::AlgAssRelOrdElem{S, T, U}, O::AlgAssRelOrd{S, T, U}) where {S, T, U} = ideal(O, x, :right)
*(O::AlgAssRelOrd{S, T, U}, x::Union{ Int, ZZRingElem }) where {S, T, U} = ideal(O, O(x), :left)
*(x::Union{ Int, ZZRingElem }, O::AlgAssRelOrd{S, T, U}) where {S, T, U} = ideal(O, O(x), :right)

@doc raw"""
    ideal(O::AlgAssRelOrd, a::AbsSimpleNumFieldOrderFractionalIdeal) -> AlgAssRelOrdIdl
    ideal(O::AlgAssRelOrd, a::RelNumFieldOrderFractionalIdeal) -> AlgAssRelOrdIdl

Returns the ideal $a \cdot O$ where $a$ is a fractional ideal of `base_ring(O)`.
"""
function ideal(O::AlgAssRelOrd{S, T, U}, a::T) where {S, T, U}
  @assert order(a) === base_ring(O)
  d = degree(O)
  pb = pseudo_basis(O, copy = false)
  M = basis_matrix(O)
  PM = pseudo_matrix(M, [ a*pb[i][2] for i = 1:d ])
  return ideal(algebra(O), O, PM; side=:twosided)
end

@doc raw"""
    ideal(O::AlgAssRelOrd, a::AbsNumFieldOrderIdeal) -> AlgAssRelOrdIdl
    ideal(O::AlgAssRelOrd, a::RelNumFieldOrderIdeal) -> AlgAssRelOrdIdl

Returns the ideal $a \cdot O$ where $a$ is an ideal of `base_ring(O)`.
"""
function ideal(O::AlgAssRelOrd{AbsSimpleNumFieldElem, AbsSimpleNumFieldOrderFractionalIdeal}, a::AbsNumFieldOrderIdeal)
  @assert order(a) === base_ring(O)
  aa = fractional_ideal(order(a), a, ZZRingElem(1))
  return ideal(O, aa)
end

function ideal(O::AlgAssRelOrd, a::RelNumFieldOrderIdeal)
  @assert order(a) === base_ring(O)

  aa = fractional_ideal(order(a), basis_pmatrix(a); M_in_hnf=true)
  return ideal(O, aa)
end

@doc raw"""
    *(O::AlgAssRelOrd, a::AbsNumFieldOrderIdeal) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, a::RelNumFieldOrderIdeal) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, a::AbsSimpleNumFieldOrderFractionalIdeal) -> AlgAssRelOrdIdl
    *(O::AlgAssRelOrd, a::RelNumFieldOrderFractionalIdeal) -> AlgAssRelOrdIdl
    *(a::AbsNumFieldOrderIdeal, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(a::RelNumFieldOrderIdeal, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(a::AbsSimpleNumFieldOrderFractionalIdeal, O::AlgAssRelOrd) -> AlgAssRelOrdIdl
    *(a::RelNumFieldOrderFractionalIdeal, O::AlgAssRelOrd) -> AlgAssRelOrdIdl

Returns the ideal $a \cdot O$ where $a$ is a (fractional) ideal of `base_ring(O)`.
"""
*(O::AlgAssRelOrd{S, T, U}, a::T) where {S, T <: Union{AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal, AbsNumFieldOrderFractionalIdeal, RelNumFieldOrderFractionalIdeal}, U} = ideal(O, a)

*(a::T, O::AlgAssRelOrd{S, T, U}) where {S, T <: Union{AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal, AbsNumFieldOrderFractionalIdeal, RelNumFieldOrderFractionalIdeal}, U} = ideal(O, a)

*(O::AlgAssRelOrd, a::Union{AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal}) = ideal(O, a)

*(a::Union{AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal}, O::AlgAssRelOrd) = ideal(O, a)

@doc raw"""
    ideal_from_lattice_gens(A::AbstractAssociativeAlgebra, v::Vector{ <: AbstractAssociativeAlgebraElem })
      -> AlgAssRelOrdIdl

Returns the ideal in $A$ generated by the elements of `v` as a lattice over
`base_ring(maximal_order(base_ring(A)))`.
"""
function ideal_from_lattice_gens(A::AbstractAssociativeAlgebra{S}, v::Vector{ <: AbstractAssociativeAlgebraElem{S} }) where { S <: NumFieldElem }
  M = zero_matrix(base_ring(A), max(length(v), dim(A)), dim(A))
  for i = 1:length(v)
    elem_to_mat_row!(M, i, v[i])
  end
  PM = pseudo_hnf(pseudo_matrix(M), :lowerleft)
  if length(v) >= dim(A)
    PM = sub(PM, (nrows(PM) - dim(A) + 1):nrows(PM), 1:dim(A))
  end

  return ideal(A, PM; M_in_hnf=true)
end
@doc raw"""
    ideal_from_lattice_gens(A::AbstractAssociativeAlgebra, O::AlgAssRelOrd,
                            v::Vector{ <: AbstractAssociativeAlgebraElem },
                            side::Symbol = :nothing)
      -> AlgAssRelOrdIdl

Returns the ideal of $O$ in $A$ generated by the elements of `v` as a lattice
over `base_ring(O)`.
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::AbstractAssociativeAlgebra{S}, O::AlgAssRelOrd{S, T, U}, v::Vector{ <: AbstractAssociativeAlgebraElem{S} }, side::Symbol = :nothing) where { S <: NumFieldElem, T, U }
  a = ideal_from_lattice_gens(A, v)
  a.order = O
  _set_sidedness(a, side)
  return a
end

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

function _zero_ideal(A::AbstractAssociativeAlgebra{S}) where { S <: NumFieldElem }
  a = ideal(A, pseudo_matrix(zero_matrix(base_ring(A), dim(A), dim(A))); M_in_hnf=true)
  a.iszero = 1
  return a
end

function _zero_ideal(A::AbstractAssociativeAlgebra{S}, O::AlgAssRelOrd{S, T, U}) where { S <: NumFieldElem, T, U }
  a = _zero_ideal(A)
  a.order = O
  _set_sidedness(a, :twosided)
  return a
end

function iszero(a::AlgAssRelOrdIdl)
  if a.iszero == 0
    if is_zero_row(basis_matrix(a, copy = false), 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::AlgAssRelOrdIdl)
  print(io, "Ideal in ")
  show(IOContext(io, :compact => true), algebra(a))
  print(io, " with basis pseudo-matrix\n")
  show(IOContext(io, :compact => true), basis_pmatrix(a, copy = false))
end

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

function Base.deepcopy_internal(a::AlgAssRelOrdIdl, 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_pmatrix_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_pmatrix_wrt
    b.basis_pmatrix_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_pseudo_basis(a::AlgAssRelOrdIdl)
  if isdefined(a, :pseudo_basis)
    return nothing
  end
  if !isdefined(a, :basis_pmatrix)
    error("No pseudo_basis and no basis_pmatrix defined.")
  end
  P = basis_pmatrix(a, copy = false)
  A = algebra(a)
  K = base_ring(A)
  pb = Vector{Tuple{elem_type(A), fractional_ideal_type(order_type(K))}}()
  for i = 1:dim(A)
    t = elem_from_mat_row(A, P.matrix, i)
    push!(pb, (t, deepcopy(P.coeffs[i])))
  end
  a.pseudo_basis = pb
  return nothing
end

function assure_has_basis_pmatrix_wrt(a::AlgAssRelOrdIdl, O::AlgAssRelOrd)
  if haskey(a.basis_pmatrix_wrt, O)
    return nothing
  end

  PM = basis_pmatrix(a)
  PM.matrix = PM.matrix*basis_mat_inv(O, copy = false)
  PM = pseudo_hnf(PM, :lowerleft)
  a.basis_pmatrix_wrt[O] = PM
  return nothing
end

function assure_has_basis_matrix(a::AlgAssRelOrdIdl)
  if isdefined(a, :basis_matrix)
    return nothing
  end
  a.basis_matrix = basis_pmatrix(a).matrix
  return nothing
end

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

################################################################################
#
#  Pseudo basis / basis pseudo-matrix
#
################################################################################

@doc raw"""
    pseudo_basis(a::AlgAssRelOrdIdl; copy::Bool = true)

Returns the pseudo basis of $a$, i. e. a vector $v$ of pairs $(e_i, a_i)$ such
that $a = \bigoplus_i a_i e_i$, where $e_i$ is an element of `algebra(a)`
and $a_i$ is a fractional ideal of `base_ring(left_order(a))`.
"""
function pseudo_basis(a::AlgAssRelOrdIdl; copy::Bool = true)
  assure_has_pseudo_basis(a)
  if copy
    return deepcopy(a.pseudo_basis)
  else
    return a.pseudo_basis
  end
end

@doc raw"""
    basis_pmatrix(a::AlgAssRelOrdIdl; copy::Bool = true) -> PMat

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

@doc raw"""
    absolute_basis(a::AlgAssRelOrdIdl) -> Vector{ <: AbstractAssociativeAlgebraElem }

Returns a basis of $a$ as $\mathbb Z$-module.
"""
function absolute_basis(a::AlgAssRelOrdIdl)
  A = algebra(a)
  pb = pseudo_basis(a, copy = false)
  res = Vector{elem_type(A)}()
  for i = 1:dim(A)
    for b in basis(pb[i][2])
      push!(res, b*pb[i][1])
    end
  end
  return res
end

@doc raw"""
    basis_pmatrix_wrt(a::AlgAssRelOrdIdl, O::AlgAssRelOrd; copy::Bool = true)
      -> PMat

Returns the basis pseudo-matrix of $a$ with respect to the basis of $O$.
"""
function basis_pmatrix_wrt(a::AlgAssRelOrdIdl, O::AlgAssRelOrd; copy::Bool = true)
  assure_has_basis_pmatrix_wrt(a, O)
  if copy
    return deepcopy(a.basis_pmatrix_wrt[O])
  else
    return a.basis_pmatrix_wrt[O]
  end
end

################################################################################
#
#  (Inverse) basis matrix
#
################################################################################

@doc raw"""
    basis_matrix(a::AlgAssRelOrdIdl; copy::Bool = true) -> MatElem

Returns the basis matrix of $a$, that is the basis pseudo-matrix of $a$ without
the coefficient ideals.
"""
function basis_matrix(a::AlgAssRelOrdIdl; copy::Bool = true)
  assure_has_basis_matrix(a)
  if copy
    return deepcopy(a.basis_matrix)
  else
    return a.basis_matrix
  end
end

@doc raw"""
    basis_mat_inv(a::AlgAssRelOrdIdl; copy::Bool = true) -> MatElem

Returns the inverse of the basis matrix of $a$.
"""
function basis_mat_inv(a::AlgAssRelOrdIdl; 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

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

@doc raw"""
    +(a::AlgAssRelOrdIdl, b::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl

Returns $a + b$.
"""
function +(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where {S, T, U}
  @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_pmatrix(a, copy = false), basis_pmatrix(b, copy = false))
  M = sub(pseudo_hnf(M, :lowerleft), (d + 1):2*d, 1:d)
  c = ideal(algebra(a), M; M_in_hnf=true)
  if isdefined(a, :order) && isdefined(b, :order) && order(a) === order(b)
    c.order = order(a)
  end
  return c
end

@doc raw"""
    *(a::AlgAssRelOrdIdl, b::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl

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

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

  A = algebra(a)
  d = dim(A)
  d2 = d^2

  pba = pseudo_basis(a, copy = false)
  pbb = pseudo_basis(b, copy = false)
  M = zero_matrix(base_ring(A), d2, d)
  C = Vector{fractional_ideal_type(order_type(base_ring(A)))}(undef, d2)
  t = one(A)
  for i = 1:d
    i1d = (i - 1)*d
    for j = 1:d
      t = mul!(t, pba[i][1], pbb[j][1])
      elem_to_mat_row!(M, i1d + j, t)
      C[i1d + j] = simplify(pba[i][2]*pbb[j][2])
    end
  end

  H = sub(pseudo_hnf(pseudo_matrix(M, C), :lowerleft), (d2 - d + 1):d2, 1:d)
  c = ideal(A, H; M_in_hnf=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 raw"""
    ^(a::AlgAssRelOrdIdl, e::Int) -> AlgAssRelOrdIdl
    ^(a::AlgAssRelOrdIdl, e::ZZRingElem) -> AlgAssRelOrdIdl

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

@doc raw"""
    intersect(a::AlgAssRelOrdIdl, b::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl

Returns $a \cap b$.
"""
function intersect(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where {S, T, U}
  @assert algebra(a) === algebra(b)
  d = dim(algebra(a))
  Ma = basis_pmatrix(a)
  Mb = basis_pmatrix(b)
  M1 = hcat(Ma, deepcopy(Ma))
  z = zero_matrix(base_ring(Ma.matrix), d, d)
  M2 = hcat(pseudo_matrix(z, Mb.coeffs), Mb)
  M = vcat(M1, M2)
  H = sub(pseudo_hnf(M, :lowerleft), 1:d, 1:d)
  c = ideal(algebra(a), H; M_in_hnf=true)

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

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

@doc raw"""
    *(a::AlgAssRelOrdIdl, x::AbsNumFieldOrderElem) -> AlgAssRelOrdIdl
    *(x::AbsNumFieldOrderElem, a::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl
    *(a::AlgAssRelOrdIdl, x::RelNumFieldOrderElem) -> AlgAssRelOrdIdl
    *(x::RelNumFieldOrderElem, a::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl
    *(a::AlgAssRelOrdIdl, x::Int) -> AlgAssRelOrdIdl
    *(x::Int, a::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl
    *(a::AlgAssRelOrdIdl, x::ZZRingElem) -> AlgAssRelOrdIdl
    *(x::ZZRingElem, a::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl
    *(a::AlgAssRelOrdIdl{S, T, U}, x::S) where { S, T, U } -> AlgAssRelOrdIdl{S, T, U}
    *(x::S, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U } -> AlgAssRelOrdIdl{S, T, U}
    *(a::AlgAssRelOrdIdl{S, T, U}, x::AbstractAssociativeAlgebraElem{S}) where { S, T, U }
      -> AlgAssRelOrdIdl{S, T, U}
    *(x::AbstractAssociativeAlgebraElem{S}, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
      -> AlgAssRelOrdIdl{S, T, U}
    *(a::AlgAssRelOrdIdl{S, T, U}, x::AlgAssRelOrdElem{S, T, U}) where { S, T, U }
      -> AlgAssRelOrdIdl{S, T, U}
    *(x::AlgAssRelOrdElem{S, T, U}, a::AlgAssRelOrdIdl{S, T}) where { S, T, U }
      -> AlgAssRelOrdIdl{S, T, U}

Returns the ideal $a*x$ respectively $x*a$.
"""
function *(a::AlgAssRelOrdIdl{S, T, U}, x::Union{ Int, ZZRingElem, AbsNumFieldOrderElem, RelNumFieldOrderElem, S }) where { S <: NumFieldElem, T, U }
  if iszero(x)
    return _zero_ideal(algebra(a))
  end
  b = ideal(algebra(a), x*basis_pmatrix(a, copy = false); M_in_hnf=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

# Let's assume the base ring of algebra(a) is commutative (it should be a field)
*(x::Union{ Int, ZZRingElem, AbsNumFieldOrderElem, RelNumFieldOrderElem, S }, a::AlgAssRelOrdIdl{S, T, U}) where { S <: NumFieldElem, T, U } = a*x

function *(a::AlgAssRelOrdIdl{S, T, U}, x::AbstractAssociativeAlgebraElem{S}) where { S <: NumFieldElem, T, U }
  if iszero(x)
    return _zero_ideal(algebra(a))
  end
  A = algebra(a)
  M = zero_matrix(base_ring(A), dim(A), dim(A))
  for i = 1:dim(A)
    t = pseudo_basis(a, copy = false)[i][1]*x
    elem_to_mat_row!(M, i, t)
  end
  PM = pseudo_matrix(M, [ deepcopy(pseudo_basis(a, copy = false)[i][2]) for i = 1:dim(A) ])
  b = ideal(A, PM)
  if isdefined(a, :left_order)
    b.left_order = left_order(a)
  end
  return b
end

function *(x::AbstractAssociativeAlgebraElem{S}, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
  if iszero(x)
    return _zero_ideal(algebra(a))
  end
  A = algebra(a)
  M = zero_matrix(base_ring(A), dim(A), dim(A))
  for i = 1:dim(A)
    t = x*pseudo_basis(a, copy = false)[i][1]
    elem_to_mat_row!(M, i, t)
  end
  PM = pseudo_matrix(M, [ deepcopy(pseudo_basis(a, copy = false)[i][2]) for i = 1:dim(A) ])
  b = ideal(A, PM)
  if isdefined(a, :right_order)
    b.right_order = right_order(a)
  end
  return b
end

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

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

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

@doc raw"""
    in(x::AbstractAssociativeAlgebraElem, a::AlgAssRelOrdIdl) -> Bool
    in(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl) -> Bool

Returns `true` if $x$ is in $a$ and `false` otherwise.
"""
function in(x::AbstractAssociativeAlgebraElem{S}, a::AlgAssRelOrdIdl{S, T, U}) where { S, T, U }
  parent(x) !== algebra(a) && error("Algebra of element and ideal must be equal")
  A = algebra(a)
  b_pmat = basis_pmatrix(a, copy = false)
  t = matrix(base_ring(A), 1, dim(A), coefficients(x, copy = false))
  t = t*basis_mat_inv(a, copy = false)
  for i = 1:dim(A)
    if !(t[1, i] in b_pmat.coeffs[i])
      return false
    end
  end
  return true
end

in(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl) = in(elem_in_algebra(x, copy = false), a)

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

@doc raw"""
    ==(a::AlgAssRelOrdIdl, b::AlgAssRelOrdIdl) -> Bool

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

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

# functions is_right_ideal and is_left_ideal are in StructureConstantAlgebra/Ideal.jl

# This only works if a.order is defined, otherwise order(a) throws an error.

function _test_ideal_sidedness(a::AlgAssRelOrdIdl, O::AlgAssRelOrd, side::Symbol)
  b = ideal(O, one(O))

  if side == :right
    c = a*b
  elseif side == :left
    c = b*a
  else
    error("side must be either :left or :right")
  end

  return _spans_subset_of_pseudohnf(basis_pmatrix(c, copy = false), basis_pmatrix(a, copy = false), :lowerleft)
end

function _test_ideal_sidedness(a::AlgAssRelOrdIdl, side::Symbol)
  O = order(a)
  b = ideal(O, one(O))

  if side == :right
    c = a*b
  elseif side == :left
    c = b*a
  else
    error("side must be either :left or :right")
  end

  return _spans_subset_of_pseudohnf(basis_pmatrix(c, copy = false), basis_pmatrix(a, copy = false), :lowerleft)
end

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

# This computes a basis pseudo-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::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}, side::Symbol) where { S, T, U }
  @assert is_full_lattice(a) && is_full_lattice(b)
  A = algebra(a)
  @assert A === algebra(b)
  K = base_ring(A)
  d = dim(A)
  pba = pseudo_basis(a, copy = false)
  pbb = pseudo_basis(b, copy = false)
  M = zero_matrix(K, d^2, d)
  for i = 1:d
    N = representation_matrix(pbb[i][1], side)*basis_mat_inv(a, copy = false)
    for s = 1:d
      for t = 1:d
        M[t + (i - 1)*d, s] = N[s, t]
      end
    end
  end
  invcoeffs = [ simplify(inv(pba[i][2])) for i = 1:d ]
  C = Array{T}(undef, d^2)
  for i = 1:d
    for j = 1:d
      C[(i - 1)*d + j] = simplify(pbb[i][2]*invcoeffs[j])
    end
  end
  PM = pseudo_matrix(M, C)
  PM = sub(pseudo_hnf(PM, :upperright, true), 1:d, 1:d)
  N = inv(transpose(PM.matrix))
  PN = pseudo_matrix(N, [ simplify(inv(I)) for I in PM.coeffs ])
  return PN
end

function ring_of_multipliers(a::AlgAssRelOrdIdl, action::Symbol = :left)
  PM = _colon_raw(a, a, action)
  return Order(algebra(a), PM)
end

@doc raw"""
    left_order(a::AlgAssRelOrdIdl) -> AlgAssRelOrd

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

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

@doc raw"""
    right_order(a::AlgAssRelOrdIdl) -> AlgAssRelOrd

Returns the largest order of which $a$ is a right ideal, that is
$\{ x \in A \mid ax \subseteq a\}$.
"""
function right_order(a::AlgAssRelOrdIdl)
  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::AlgAssRelOrdIdl)
  if !isdefined(a, :left_order)
    return false
  end
  O = left_order(a)
  if is_maximal_known(O) && is_maximal(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::AlgAssRelOrdIdl)
  if !isdefined(a, :right_order)
    return false
  end
  O = right_order(a)
  if is_maximal_known(O) && is_maximal(O)
    return true
  end
  return false
end

################################################################################
#
#  Reduction of element modulo ideal
#
################################################################################

function mod!(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl)
  A = algebra(a)
  O = parent(x)
  @assert order(a) === O
  b = coordinates(x, copy = false)
  PM = basis_pmatrix_wrt(a, O, copy = false) # PM is assumed to be in lower left pseudo hnf
  t = parent(b[1])()
  t1 = parent(b[1])()
  for i = degree(O):-1:1
    t = add!(t, mod(b[i], PM.coeffs[i]), -b[i])
    for j = 1:i
      t1 = mul!(t1, t, PM.matrix[i, j])
      b[j] = add!(b[j], b[j], t1)
    end
  end

  t = A()
  B = pseudo_basis(O, copy = false)
  zero!(x.elem_in_algebra)
  for i = 1:degree(O)
    t = mul!(t, B[i][1], A(b[i]))
    x.elem_in_algebra = add!(x.elem_in_algebra, x.elem_in_algebra, t)
  end

  return x
end

@doc raw"""
    mod(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl) -> AlgAssRelOrdElem

Returns $y$ in `parent(x)` such that $x \equiv y \mod a$ and the coefficients
of $y$ are reduced modulo $a$.
Assumes `parent(x) == order(a)` and that $a$ is an integral ideal of `order(a)`.
"""
function mod(x::AlgAssRelOrdElem, a::AlgAssRelOrdIdl)
  return mod!(deepcopy(x), a)
end

function mod!(x::AlgAssRelOrdElem, Q::RelOrdQuoRing)
  return mod!(x, ideal(Q))
end

function mod(x::AlgAssRelOrdElem, Q::RelOrdQuoRing)
  return mod(x, ideal(Q))
end

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

# Assumes, that det(basis_matrix(a)) == 1
function assure_has_norm(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}) where { S, T, U }
  if haskey(a.norm, O)
    return nothing
  end
  if iszero(a)
    a.norm[O] = _algebra(O)()*O
    return nothing
  end

  # Now consider a as an ideal of O (but without computing a basis pseudo matrix
  # in the basis of O).
  c = basis_pmatrix(a, copy = false).coeffs
  d = inv_coeff_ideals(O, copy = false)
  n = c[1]*d[1]
  for i = 2:dim(algebra(a))
    n *= c[i]*d[i]
  end
  # The basis matrix of a in the basis of O would be
  # basis_matrix(a)*basis_mat_inv(O), so its determinant is
  # inv(det(basis_matrix(O))), since det(basis_matrix(a)) == 1.
  n *= inv(det(basis_matrix(O, copy = false)))
  simplify(n)
  a.norm[O] = n
  return nothing
end

@doc raw"""
    norm(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}; copy::Bool = true)
      where { S, T, U } -> T

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

@doc raw"""
    norm(a::AlgAssRelOrdIdl{S, T, U}; copy::Bool = true) where { S, T, U } -> T

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

function assure_has_normred(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}) where { S, T, U }
  if haskey(a.normred, O)
    return nothing
  end
  if iszero(a)
    a.normred[O] = norm(a, O)
    return nothing
  end

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

@doc raw"""
    normred(a::AlgAssRelOrdIdl{S, T}, O::AlgAssRelOrd{S, T}; copy::Bool = true)
      where { S, T, U } -> T

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::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}; copy::Bool = true) where { S, T, U }
  @assert is_simple(algebra(a)) && is_central(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 raw"""
    normred(a::AlgAssRelOrdIdl{S, T, U}; copy::Bool = true) where { S, T, U } -> T

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

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

@doc raw"""
    locally_free_basis(a::AlgAssRelOrdIdl, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal })
      -> AlgAssRelOrdElem

Returns an element $x$ of the order $O$ of $a$ such that $a_p = O_p \cdot x$
where $p$ is a prime ideal of `base_ring(O)`.
See also `is_locally_free`.
"""
locally_free_basis(I::AlgAssRelOrdIdl, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) = locally_free_basis(order(I), I, p)

@doc raw"""
    locally_free_basis(O::AlgAssRelOrd, a::AlgAssRelOrdIdl,
                       p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) -> AlgAssRelOrdElem

Returns an element $x$ of $O$ such that $a_p = O_p \cdot x$ where $p$ is a
prime ideal of `base_ring(O)`.
It is assumed that $a$ is an ideal of $O$ and $a \subseteq O$.
See also `is_locally_free`.
"""
function locally_free_basis(O::AlgAssRelOrd, I::AlgAssRelOrdIdl, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal })
  b, x = is_locally_free(O, I, p)
  if !b
    error("The ideal is not locally free at the prime")
  end
  return x
end

@doc raw"""
    is_locally_free(a::AlgAssRelOrdIdl, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal })
      -> Bool, AlgAssRelOrdElem

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 ideal of `base_ring(O)`.
See also `locally_free_basis`.
"""
is_locally_free(I::AlgAssRelOrdIdl, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) = is_locally_free(order(I), I, p)

# See Bley, Wilson "Computations in relative algebraic K-groups", section 4.2
@doc raw"""
    is_locally_free(O::AlgAssRelOrd, a::AlgAssRelOrdIdl,
                   p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) -> Bool, AlgAssRelOrdElem

Returns a tuple `(true, x)` with an element $x$ of $O$ such that $a_p = O_p x$
if $a$ is locally free at $p$, and `(false, 0)` otherwise.
$p$ is a prime ideal of `base_ring(O)`.
It is assumed that $a$ is an ideal of $O$ and $a \subseteq O$.
See also `locally_free_basis`.
"""
function is_locally_free(O::AlgAssRelOrd, I::AlgAssRelOrdIdl, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal })
  pO = p*O
  OpO, toOpO = StructureConstantAlgebra(O, pO, p)
  J = radical(OpO)
  OJ, toOJ = quo(OpO, J)
  decOJ = decompose(OJ)

  pI = pO*I
  IpI, toIpI = StructureConstantAlgebra(I, pI, p)
  basisIpI = elem_type(O)[ toIpI\IpI[i] for i = 1:dim(IpI) ]
  gensJ = Vector{elem_type(IpI)}()
  for b in basis(J, copy = false)
    bb = toOpO\b
    for c in basisIpI
      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)
    _T = MatAlgebra{elem_type(base_ring(B)), dense_matrix_type(base_ring(B))}
    C, BtoC = _as_matrix_algebra(B)::Tuple{_T, AbsAlgAssMor{StructureConstantAlgebra{elem_type(base_ring(B))}, _T, dense_matrix_type(base_ring(B))}}
    e = toOpO\(toOJ\(AtoOJ(BtoA(BtoC\C[1]))))
    basiseIJ = Vector{elem_type(IJ)}()
    for b in basisIpI
      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
#
################################################################################

# See Friedrichs: "Berechnung von Maximalordnungen über Dedekindringen", Algorithmus 5.1
@doc raw"""
    pradical(O::AlgAssRelOrd, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal })
      -> AlgAssRelOrdIdl

Returns the ideal $\sqrt{p \cdot O}$ where $p$ is a prime ideal of `base_ring(O)`.
"""
function pradical(O::AlgAssRelOrd, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal })
  K = base_ring(algebra(O))
  OK = maximal_order(K)
  pO = p*O
  OpO, OtoOpO = StructureConstantAlgebra(O, pO, p)
  J = radical(OpO)

  if isempty(basis(J, copy = false))
    return pO
  end

  N = basis_pmatrix(pO, copy = false)
  m = numerator(det(N), copy = false)
  t = pseudo_matrix(zero_matrix(K, 1, degree(O)))
  for b in basis(J, copy = false)
    bb = OtoOpO\b
    elem_to_mat_row!(t.matrix, 1, elem_in_algebra(bb, copy = false))
    N = vcat(N, deepcopy(t))
  end
  N = sub(pseudo_hnf_full_rank_with_modulus(N, m, :lowerleft), nrows(N) - degree(O) + 1:nrows(N), 1:degree(O))
  return ideal(algebra(O), O, N; side=:twosided, M_in_hnf=true)
end

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

# See Friedrichs: "Berechnung von Maximalordnungen über Dedekindringen", Algorithmus 5.23
@doc raw"""
    prime_ideals_over(O::AlgAssRelOrd, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal })
      -> Vector{AlgAssRelOrdIdl}

Returns all prime ideals of $O$ lying over $p$ where $p$ is a prime ideal
of `base_ring(O)`.
"""
function prime_ideals_over(O::AlgAssRelOrd, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal })
  return _prime_ideals_over(O, pradical(O, p), p)
end

# prad is expected to be pradical(O, p). If strict_containment == true and prad
# is already prime, we return an empty array.
function _prime_ideals_over(O::AlgAssRelOrd, prad::AlgAssRelOrdIdl, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }; strict_containment::Bool = false)
  K = base_ring(algebra(O))
  OK = order(p)

  A, OtoA = StructureConstantAlgebra(O, prad, p)
  decA = decompose(A)

  if length(decA) == 1
    if strict_containment
      return ideal_type(O)[]
    else
      return [ prad ]
    end
  end

  lifted_components = Vector{typeof(basis_pmatrix(prad, copy = false))}()
  for k = 1:length(decA)
    N = zero_matrix(K, dim(decA[k][1]), degree(O))
    for i = 1:dim(decA[k][1])
      b = OtoA\(decA[k][2](decA[k][1][i]))
      elem_to_mat_row!(N, i, elem_in_algebra(b, copy = false))
    end
    push!(lifted_components, pseudo_matrix(N))
  end

  primes = Vector{ideal_type(O)}()
  for i = 1:length(decA)
    N = basis_pmatrix(prad, copy = false)
    m = numerator(det(N), copy = false)
    for j = 1:length(decA)
      if i == j
        continue
      end

      N = vcat(N, lifted_components[j])
    end
    N = sub(pseudo_hnf_full_rank_with_modulus(N, m, :lowerleft), nrows(N) - degree(O) + 1:nrows(N), 1:degree(O))
    push!(primes, ideal(algebra(O), O, N; side=:twosided, M_in_hnf=true))
  end

  return primes
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::AlgAssRelOrd, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }, side::Symbol)
  A = algebra(O)
  @assert is_simple(A)
  @assert is_maximal(O)

  K = base_ring(algebra(O))
  OK = base_ring(O)

  P = pradical(O, p) # if the algebra is simple, then the pradical is the unique 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, CtoB = _as_algebra_over_center(B)
  D, CtoD = _as_matrix_algebra(C)

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

  N = basis_pmatrix(P)
  m = numerator(det(N), copy = false)
  t = zero_matrix(K, 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, pseudo_matrix(deepcopy(t), [ K(1)*OK ]))
    end
  end
  N = sub(pseudo_hnf_full_rank_with_modulus(N, m, :lowerleft), nrows(N) - degree(O) + 1:nrows(N), 1:degree(O))

  M = ideal(algebra(O), O, N; side, M_in_hnf=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::AlgAssRelOrdIdl, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }, 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 is_simple(algebra(O))
  @assert is_maximal(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 = pradical(O, p) # if the algebra is simple, then the pradical is the unique prime lying over p
  J = I + P
  if normred(J, O) == p
    return J
  end

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

  JinC = ideal_from_gens(C, elem_type(C)[ toC(BtoOP\(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(base_ring(algebra(O)), 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
  PM = vcat(basis_pmatrix(P), pseudo_matrix(t))
  n = numerator(det(basis_pmatrix(P, copy = false)), copy = false)
  PM = sub(pseudo_hnf_full_rank_with_modulus(PM, n, :lowerleft), length(basis_c) + 1:nrows(PM), 1:ncols(PM))

  M = ideal(algebra(O), O, PM; side, M_in_hnf=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

################################################################################
#
#  Factorization into maximal ideals
#
################################################################################

# Computes maximal integral ideals M_1, ..., M_k such that I = M_1\cdots M_k,
# and such that the right order of any factor is the left order of the next
# factor.
# Always considers I as an ideal of its left order.
function factor(I::AlgAssRelOrdIdl)
  @assert !iszero(I)
  O = left_order(I)
  @assert is_maximal(O)
  I != 1*O || error("I must be proper")

  factors = Vector{ideal_type(O)}()
  n = normred(I, O)
  fac_n = factor(n)
  primes = collect(keys(fac_n))
  sort!(primes, lt = (p, q) -> minimum(p, copy = false) < minimum(q, copy = false))
  fac_n[primes[end]] -= 1 # We don't need to find the "last" maximal ideal
  for p in primes
    for i = 1:fac_n[p]
      M = maximal_integral_ideal_containing(I, p, :left)
      push!(factors, M)
      I = divexact_left(I, M)
    end
  end
  push!(factors, I)

  return factors
end

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

@doc raw"""
    rand(a::AlgAssRelOrdIdl, B::Int) -> AbstractAssociativeAlgebraElem

Returns a random element of $a$ whose coefficient size is controlled by $B$.
"""
function rand(a::AlgAssRelOrdIdl, B::Int)
  pb = pseudo_basis(a, copy = false)
  z = algebra(a)()
  for i = 1:degree(order(a))
    t = rand(pb[i][2], B)
    z += t*pb[i][1]
  end
  return z
end

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

# Returns x \in A with Ix + a*O == O and Ix \subseteq O.
# a should be an ideal of base_ring(O).
function integral_coprime_representative(O::AlgAssRelOrd, I::AlgAssRelOrdIdl, a::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal })
  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))

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

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

@doc raw"""
    inv(a::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl

Returns an ideal $b$ in `algebra(a)` such that `a*b == left_order(a)`
and `b*a == right_order(b)`.
"""
function inv(a::AlgAssRelOrdIdl)
  if isdefined(a, :right_order)
    O = right_order(a)
    PM = _colon_raw(O(1)*O, a, :right)
  else
    O = left_order(a)
    PM = _colon_raw(O(1)*O, a, :left)
  end
  b = ideal(algebra(a), PM)
  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

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

@doc raw"""
    divexact_left(a::AlgAssRelOrdIdl, b::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl

Returns an ideal $c$ such that $a = b \cdot c$.
"""
function divexact_left(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}; check::Bool=true) where { S, T, U }
  @assert algebra(a) === algebra(b)
  PM = _colon_raw(a, b, :left)
  c = ideal(algebra(a), PM)

  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 raw"""
    divexact_right(a::AlgAssRelOrdIdl, b::AlgAssRelOrdIdl) -> AlgAssRelOrdIdl

Returns an ideal $c$ such that $a = c \cdot b$.
"""
function divexact_right(a::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}; check::Bool=true) where { S, T, U }
  @assert algebra(a) === algebra(b)
  PM = _colon_raw(a, b, :right)
  c = ideal(algebra(a), PM)

  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::AlgAssRelOrdIdl{S, T, U}, b::AlgAssRelOrdIdl{S, T, U}) where { S, T, U } = divexact_right(a, b)

################################################################################
#
#  Denominator
#
################################################################################

@doc raw"""
    denominator(a::AlgAssRelOrdIdl, O::AlgAssRelOrd) -> ZZRingElem

Returns the smallest positive integer $d$ such that $da$ is contained in $O$.
"""
function denominator(a::AlgAssRelOrdIdl{S, T, U}, O::AlgAssRelOrd{S, T, U}) where { S, T, U }
  @assert algebra(a) === algebra(O)
  if iszero(a)
    return ZZRingElem(1)
  end

  n = dim(algebra(a))
  PM = basis_pmatrix_wrt(a, O, copy = false)
  pb = pseudo_basis(O, copy = false)
  inv_coeffs = inv_coeff_ideals(O, copy = false)
  d = ZZRingElem(1)
  for i = 1:n
    for j = 1:i
      d = lcm(d, denominator(simplify(PM.matrix[i, j]*PM.coeffs[i]*inv_coeffs[j])))
    end
  end
  return d
end


@doc raw"""
    denominator(a::AlgAssRelOrdIdl) -> ZZRingElem

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

# Assumes that I is "locally integral at p", i. e. I_p \subseteq O_p.
# Returns x in R\setminus p such that Ix \subseteq O, where R = order(p)
function coprime_denominator(I::AlgAssRelOrdIdl{T1, T2}, O::AlgAssRelOrd{T1, T2}, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) where { T1, T2 }
  basis_O, basis_I, MO, MI = coprime_bases(O, I, p)
  OK = order(p)
  S = Dict{ideal_type(OK), Int}()
  for i = 1:degree(O)
    ai = inv(basis_O[i][2])
    for j = 1:degree(O)
      c = basis_I[j][2]*ai*MI[i, j]
      if iszero(norm(c))
        continue
      end
      d = denominator(c)
      facD = factor(d)
      for (q, e) in facD
        qdec = prime_decomposition(OK, q)
        for (Q, _) in qdec
          v = valuation(c, Q)
          if v >= 0
            continue
          end
          if haskey(S, Q)
            f = S[Q]
            S[Q] = max(f, -v)
          else
            S[Q] = -v
          end
        end
      end
    end
  end
  if haskey(S, p)
    error("The ideal is not locally integral at p")
  end
  primes = Vector{ideal_type(OK)}()
  vals = Vector{Int}()
  for (q, e) in S
    push!(primes, q)
    push!(vals, e)
  end
  push!(primes, p)
  push!(vals, 0)
  return approximate_nonnegative(vals, primes)
end

# Assumes N \subseteq M.
# Returns pseudo bases (a_i, alpha_i)_i of M and (b_i, beta_i)_i of N, a basis
# matrix of M_p and a basis matrix (b_{ij})_{i,j} of N_p in the basis of M_p,
# such that M_p = \bigoplus_i R_p alpha_i and N_p = \bigoplus_i R_p beta_i,
# where R = order(p) and M = \bigoplus a_i alpha_i and N = \bigoplus b_i beta_i,
# beta_i = sum_j b_{ij} alpha_j.
# Top level functions to avoid "type mix-ups" (like AlgAssRelOrd together with RelNumFieldOrderIdeal).
coprime_bases(O::RelNumFieldOrder{S, T, U}, I::RelNumFieldOrderIdeal{S, T, U}, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) where { S, T, U } = _coprime_bases(O, I, p)

coprime_bases(I::RelNumFieldOrderIdeal{S, T, U}, J::RelNumFieldOrderIdeal{S, T, U}, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) where { S, T, U } = _coprime_bases(I, J, p)

coprime_bases(O::AlgAssRelOrd{S, T, U}, I::AlgAssRelOrdIdl{S, T, U}, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) where { S, T, U } = _coprime_bases(O, I, p)

coprime_bases(I::AlgAssRelOrdIdl{S, T, U}, J::AlgAssRelOrdIdl{S, T, U}, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) where { S, T, U } = _coprime_bases(I, J, p)

function _coprime_bases(M::Union{ RelNumFieldOrder{S, T, U}, AlgAssRelOrd{S, T, U}, RelNumFieldOrderIdeal{S, T, U}, AlgAssRelOrdIdl{S, T, U} }, N::Union{ RelNumFieldOrderIdeal{S, T, U}, AlgAssRelOrdIdl{S, T, U} }, p::Union{ AbsNumFieldOrderIdeal, RelNumFieldOrderIdeal }) where { S, T, U }
  A = _algebra(M)
  @assert A === _algebra(N)
  OK = order(p)
  u = elem_in_nf(uniformizer(p), copy = false)
  iu = inv(u)
  basis_M = Vector{Tuple{elem_type(A), fractional_ideal_type(OK)}}()
  for (b, c) in pseudo_basis(M, copy = false)
    v = valuation(c, p)
    if v == 0
      push!(basis_M, (deepcopy(b), deepcopy(c)))
    elseif v < 0
      v = -v
      b = b*iu^v
      c = c*u^v
      push!(basis_M, (b, c))
    else
      b = b*u^v
      c = c*iu^v
      push!(basis_M, (b, c))
    end
  end

  basis_N = Vector{Tuple{elem_type(A), fractional_ideal_type(OK)}}()
  for (b, c) in pseudo_basis(N, copy = false)
    v = valuation(c, p)
    if v == 0
      push!(basis_N, (deepcopy(b), deepcopy(c)))
    elseif v < 0
      v = -v
      b = b*iu^v
      c = c*u^v
      push!(basis_N, (b, c))
    else
      b = b*u^v
      c = c*iu^v
      push!(basis_N, (b, c))
    end
  end

  mat_M = zero_matrix(_base_ring(A), dim(A), dim(A))
  for i = 1:dim(A)
    elem_to_mat_row!(mat_M, i, basis_M[i][1])
  end

  mat_N = zero_matrix(_base_ring(A), dim(A), dim(A))
  for i = 1:dim(A)
    elem_to_mat_row!(mat_N, i, basis_N[i][1])
  end
  mat_N = mat_N*inv(mat_M)

  return basis_M, basis_N, mat_M, mat_N
end

################################################################################
#
#  Preimage
#
################################################################################

function _as_ideal_of_smaller_algebra(m::AbsAlgAssMor, I::AlgAssRelOrdIdl)
  A = domain(m)
  B = codomain(m)
  K = base_ring(A)
  @assert K === base_ring(B)
  @assert dim(A) <= dim(B)
  @assert algebra(I) === B
  OA = maximal_order(A)
  # Transport OA to B
  M = zero_matrix(base_ring(B), dim(B), dim(B))
  c = AbsSimpleNumFieldOrderFractionalIdeal[]
  PB = pseudo_basis(OA, copy = false)
  for i = 1:dim(A)
    t = m(PB[i][1])
    elem_to_mat_row!(M, i, t)
    push!(c, PB[i][2])
  end
  for i in (dim(A)+1):dim(B)
    push!(c, one(K) * maximal_order(K))
  end
  J = _intersect_modules(basis_pmatrix(I), pseudo_matrix(M, c), false)
  # Compute the intersection of M and I
  # Map the basis to A
  gen = elem_type(A)[]
  for i = 1:nrows(J)
    t = elem_from_mat_row(B, J.matrix, i)
    b, s = has_preimage_with_preimage(m, t)
    @assert b
    for el in absolute_basis(J.coeffs[i]) # replace with generators
      push!(gen, A(el) * s)
    end
  end
  J = ideal_from_lattice_gens(A, OA, gen)
  return J
end