diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts
index 32af371ea7..9609aaf411 100644
--- a/docs/src/.vitepress/config.mts
+++ b/docs/src/.vitepress/config.mts
@@ -95,6 +95,7 @@ export default defineConfig({
                 { text: 'Introduction', link: 'manual/algebras/intro'},
                 { text: 'Basics', link: 'manual/algebras/basics'},
                 { text: 'Structure constant algebras', link: 'manual/algebras/structureconstant'},
+                { text: 'Ideals', link: 'manual/algebras/ideals'},
               ]
             },
             {
diff --git a/docs/src/manual/algebras/basics.md b/docs/src/manual/algebras/basics.md
index 9fc7ff026f..154bc3ba04 100644
--- a/docs/src/manual/algebras/basics.md
+++ b/docs/src/manual/algebras/basics.md
@@ -9,7 +9,7 @@ end
 
 ## Creation of algebras
 
-See the corresponding sections in ?????
+See the corresponding sections on [structure constant algebras](@ref SCA).
 ```@docs
 zero_algebra(::Field)
 ```
diff --git a/docs/src/manual/algebras/structureconstant.md b/docs/src/manual/algebras/structureconstant.md
index 3cc4f003a2..f7c3117989 100644
--- a/docs/src/manual/algebras/structureconstant.md
+++ b/docs/src/manual/algebras/structureconstant.md
@@ -1,4 +1,4 @@
-# Structure constant algebras
+# [Structure constant algebras](@id SCA)
 
 ```@meta
 CurrentModule = Hecke
diff --git a/src/AlgAss/AbsAlgAss.jl b/src/AlgAss/AbsAlgAss.jl
index 00d854684b..979ab75f20 100644
--- a/src/AlgAss/AbsAlgAss.jl
+++ b/src/AlgAss/AbsAlgAss.jl
@@ -727,7 +727,7 @@ of the basis elements of $A$ and a map from $B$ to $A$.
 function regular_matrix_algebra(A::Union{ StructureConstantAlgebra, GroupAlgebra })
   K = base_ring(A)
   B = matrix_algebra(K, [ representation_matrix(A[i], :right) for i = 1:dim(A) ], isbasis = true)
-  return B, hom(B, A, identity_matrix(K, dim(A)), identity_matrix(K, dim(A)))
+  return B, hom(B, A, identity_matrix(K, dim(A)), identity_matrix(K, dim(A)); check = false)
 end
 
 ###############################################################################
@@ -1117,7 +1117,7 @@ function product_of_components_with_projection(A::AbstractAssociativeAlgebra, a:
   for b in basis(A)
     push!(imgs, sum(injstoB[i](injs[i]\b) for i in eachindex(a); init = zero(B)))
   end
-  p = hom(A, B, basis_matrix(imgs))
+  p = hom(A, B, basis_matrix(imgs); check = false)
   _transport_refined_wedderburn_decomposition_forward(p)
   return B, p
 end
diff --git a/src/AlgAss/AlgGrp.jl b/src/AlgAss/AlgGrp.jl
index 68365e56dc..1e1517f7ef 100644
--- a/src/AlgAss/AlgGrp.jl
+++ b/src/AlgAss/AlgGrp.jl
@@ -10,7 +10,7 @@ denominator_of_multiplication_table(A::GroupAlgebra{QQFieldElem}) = one(ZZ)
 
 base_ring(A::GroupAlgebra{T}) where {T} = A.base_ring::parent_type(T)
 
-base_ring_type(::Type{GroupAlgebra{T}}) where {T} = parent_type(T)
+base_ring_type(::Type{GroupAlgebra{T, S, R}}) where {T, S, R} = parent_type(T)
 
 Generic.dim(A::GroupAlgebra) = size(multiplication_table(A, copy = false), 1)
 
diff --git a/src/AlgAss/AlgMat.jl b/src/AlgAss/AlgMat.jl
index a9bc5532ca..162c072e65 100644
--- a/src/AlgAss/AlgMat.jl
+++ b/src/AlgAss/AlgMat.jl
@@ -14,7 +14,13 @@ base_ring_type(::Type{MatAlgebra{T, S}}) where {T, S} = parent_type(T)
 
 coefficient_ring(A::MatAlgebra{T, S}) where {T, S} = A.coefficient_ring::base_ring_type(S)
 
-basis(A::MatAlgebra) = A.basis::Vector{elem_type(A)}
+function basis(A::MatAlgebra; copy::Bool = true)
+  if copy
+    return Base.copy(A.basis::Vector{elem_type(A)})
+  else
+    return A.basis::Vector{elem_type(A)}
+  end
+end
 
 has_one(A::MatAlgebra) = true
 
diff --git a/src/AlgAss/Elem.jl b/src/AlgAss/Elem.jl
index 6ea54af3b3..976dbfc966 100644
--- a/src/AlgAss/Elem.jl
+++ b/src/AlgAss/Elem.jl
@@ -644,6 +644,11 @@ function (A::GroupAlgebra)(a::GroupAlgebraElem)
   return a
 end
 
+function (A::MatAlgebra)(a::MatAlgebraElem)
+  @assert parent(a) == A "Wrong parent"
+  return a
+end
+
 # For polynomial substitution
 for T in subtypes(AbstractAssociativeAlgebra)
   @eval begin
diff --git a/src/AlgAss/Ideal.jl b/src/AlgAss/Ideal.jl
index 8103443586..6b741ff0c3 100644
--- a/src/AlgAss/Ideal.jl
+++ b/src/AlgAss/Ideal.jl
@@ -7,17 +7,39 @@ algebra(a::AbsAlgAssIdl) = a.algebra
 
 iszero(a::AbsAlgAssIdl) = (a.iszero == 1)
 
+@doc raw"""
+    dim(a::AbsAlgAssIdl) -> Int
+
+Return the vector space dimension an ideal.
+"""
+dim(a::AbsAlgAssIdl) = nrows(basis_matrix(a, copy = false))
+
 ###############################################################################
 #
 #  String I/O
 #
 ###############################################################################
 
-function show(io::IO, a::AbsAlgAssIdl)
-  print(io, "Ideal of ")
-  print(io, algebra(a))
-  println(io, " with basis matrix")
-  print(io, basis_matrix(a, copy = false))
+function show(io::IO, A::AbsAlgAssIdl)
+  if is_terse(io)
+    print(io, "Ideal")
+  else
+    io = pretty(io)
+    print(io, "Ideal of dimension ", dim(A), " in ")
+    print(terse(io), Lowercase(), algebra(A))
+  end
+end
+
+function show(io::IO, mime::MIME"text/plain", a::AbsAlgAssIdl)
+  println(io, "Ideal")
+  io = pretty(io)
+  println(io, Indent(), "of dimension ", dim(a))
+  println(io, "in ", Lowercase(), algebra(a))
+  println(io, "with basis matrix")
+  print(io, Indent())
+  show(io, MIME"text/plain"(), basis_matrix(a, copy = false))
+  print(io, Dedent())
+  print(io, Dedent())
 end
 
 ################################################################################
@@ -66,9 +88,9 @@ Returns the basis of $a$.
 function basis(a::AbsAlgAssIdl; copy::Bool = true)
   assure_has_basis(a)
   if copy
-    return deepcopy(a.basis)
+    return Base.copy(a.basis)::Vector{elem_type(algebra(a))}
   else
-    return a.basis
+    return a.basis::Vector{elem_type(algebra(a))}
   end
 end
 
@@ -79,9 +101,19 @@ Returns the basis matrix of $a$ with respect to the basis of the algebra.
 """
 function basis_matrix(a::AbsAlgAssIdl; copy::Bool = true)
   if copy
-    return deepcopy(a.basis_matrix)
+    return deepcopy(a.basis_matrix)::dense_matrix_type(base_ring_type(algebra(a)))
   else
-    return a.basis_matrix
+    return a.basis_matrix::dense_matrix_type(base_ring_type(algebra(a)))
+  end
+end
+
+function basis_matrix_solve_context(a::AbsAlgAssIdl)
+  if isdefined(a, :basis_matrix_solve_ctx)
+    return a.basis_matrix_solve_ctx::Solve.solve_context_type(base_ring(algebra(a)))
+  else
+    c = Solve.solve_init(basis_matrix(a, copy = false))
+    a.basis_matrix_solve_ctx = c
+    return c
   end
 end
 
@@ -91,15 +123,13 @@ end
 #
 ################################################################################
 
-@doc raw"""
-    in(x::AbstractAssociativeAlgebraElem, a::AbsAlgAssIdl) -> Bool
+function in(x, a::AbsAlgAssIdl)
+  c = basis_matrix_solve_context(a)
+  return can_solve(c, coefficients(x, copy = false); side = :left)
+end
 
-Returns `true` if $x$ is an element of $a$ and `false` otherwise.
-"""
-function in(x::T, a::AbsAlgAssIdl{S, T, U}) where {S, T, U}
-  A = algebra(a)
-  M = matrix(base_ring(A), 1, dim(A), coefficients(x, copy = false))
-  return rank(vcat(basis_matrix(a, copy = false), M)) == nrows(basis_matrix(a, copy = false)) # so far we assume nrows(basis_matrix) == rank(basis_matrix)
+function is_subset(a::AbsAlgAssIdl, b::AbsAlgAssIdl)
+  return all(in(b), basis(a, copy = false))
 end
 
 ################################################################################
@@ -112,18 +142,21 @@ function _test_ideal_sidedness(a::AbsAlgAssIdl, side::Symbol)
   A = algebra(a)
   ba = basis(a, copy = false)
   t = A()
-  for i = 1:dim(A)
-    for j = 1:length(ba)
-      if side == :left
+  for i in 1:dim(A)
+    for j in 1:length(ba)
+      if side === :left || side === :twosided
         t = mul!(t, A[i], ba[j])
-      elseif side == :right
+        if !(t in a)
+          return false
+        end
+      elseif side === :right || side === :twosided
         t = mul!(t, ba[j], A[i])
+        if !(t in a)
+          return false
+        end
       else
         error("side must be either :left or :right")
       end
-      if !(t in a)
-        return false
-      end
     end
   end
   return true
@@ -181,16 +214,13 @@ end
 #
 ################################################################################
 
-@doc raw"""
-    +(a::AbsAlgAssIdl, b::AbsAlgAssIdl) -> AbsAlgAssIdl
+function +(a::AbsAlgAssIdl{S}, b::AbsAlgAssIdl{S}) where {S}
+  @req algebra(a) === algebra(b) "Ideals must have same algebra"
 
-Returns $a + b$.
-"""
-function +(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where {S, T, U}
-  if iszero(a)
-    return deepcopy(b)
-  elseif iszero(b)
-    return deepcopy(a)
+  if is_zero(a)
+    return b
+  elseif is_zero(b)
+    return a
   end
 
   M = vcat(basis_matrix(a), basis_matrix(b))
@@ -198,19 +228,15 @@ function +(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where {S, T, U}
   if r != nrows(M)
     M = sub(M, 1:r, 1:ncols(M))
   end
-  return ideal(algebra(a), M; M_in_rref=true)
+  return _ideal_from_matrix(algebra(a), M; M_in_rref=true)
 end
 
-@doc raw"""
-    *(a::AbsAlgAssIdl, b::AbsAlgAssIdl) -> AbsAlgAssIdl
-
-Returns $a \cdot b$.
-"""
-function *(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where {S, T, U}
+function *(a::AbsAlgAssIdl{S}, b::AbsAlgAssIdl{S}) where {S}
+  @req algebra(a) === algebra(b) "Ideals must have same algebra"
   if iszero(a)
-    return deepcopy(a)
+    return a
   elseif iszero(b)
-    return deepcopy(b)
+    return b
   end
 
   A = algebra(a)
@@ -223,20 +249,15 @@ function *(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where {S, T, U}
       elem_to_mat_row!(M, ii + j, ba[i]*bb[j])
     end
   end
-  return ideal(algebra(a), M)
+  return _ideal_from_matrix(algebra(a), M)
 end
 
-@doc raw"""
-    ^(a::AbsAlgAssIdl, e::Union{ Int, ZZRingElem }) -> AbsAlgAssIdl
-
-Returns $a^e$.
-"""
 ^(A::AbsAlgAssIdl, e::Int) = Base.power_by_squaring(A, e)
 ^(A::AbsAlgAssIdl, e::ZZRingElem) = Base.power_by_squaring(A, BigInt(e))
 
 function one(a::AbsAlgAssIdl)
   A = algebra(a)
-  return ideal(A, identity_matrix(base_ring(A), dim(A)); side=:twosided, M_in_rref=true)
+  return _ideal_from_matrix(A, identity_matrix(base_ring(A), dim(A)); side=:twosided, M_in_rref=true)
 end
 
 function Base.copy(a::AbsAlgAssIdl)
@@ -246,11 +267,11 @@ end
 function *(x::AbstractAssociativeAlgebraElem, a::AbsAlgAssIdl)
   @assert is_left_ideal(a) "Not a left ideal"
   if iszero(a)
-    return deepcopy(a)
+    return a
   end
 
   basis_a = basis(a, copy = false)
-  return ideal_from_gens(algebra(a), [ x*basis_a[i] for i = 1:length(basis_a) ])
+  return _ideal_from_kgens(algebra(a), [ x*basis_a[i] for i = 1:length(basis_a) ])
 end
 
 function *(a::AbsAlgAssIdl, x::AbstractAssociativeAlgebraElem)
@@ -260,21 +281,18 @@ function *(a::AbsAlgAssIdl, x::AbstractAssociativeAlgebraElem)
   end
 
   basis_a = basis(a, copy = false)
-  return ideal_from_gens(algebra(a), [ basis_a[i]*x for i = 1:length(basis_a) ])
+  return _ideal_from_kgens(algebra(a), [ basis_a[i]*x for i = 1:length(basis_a) ])
 end
+
 ################################################################################
 #
+#
 #  Equality
 #
 ################################################################################
 
-@doc raw"""
-    ==(a::AbsAlgAssIdl, b::AbsAlgAssIdl) -> Bool
-
-Returns `true` if $a$ and $b$ are equal and `false` otherwise.
-"""
 function ==(a::AbsAlgAssIdl, b::AbsAlgAssIdl)
-  algebra(a) != algebra(b) && return false
+  algebra(a) !== algebra(b) && return false
   return basis_matrix(a, copy = false) == basis_matrix(b, copy = false)
 end
 
@@ -284,17 +302,41 @@ end
 #
 ################################################################################
 
-@doc raw"""
-    ideal_from_gens(A::AbstractAssociativeAlgebra, b::Vector{ <: AbstractAssociativeAlgebraElem},
-                    side::Symbol = :nothing)
-      -> AbsAlgAssIdl
+# side is necessary
+function ideal(A::AbstractAssociativeAlgebra, b::Vector; side::Symbol)
+  @req side in (:left, :right, :twosided) "Side must be :left, :right or :twosided"
+  if length(b) == 0
+    M = zero_matrix(base_ring(A), 0, dim(A))
+    return _ideal_from_matrix(A, M; side, M_in_rref=true)
+  end
+  B = basis(A)
+
+  kgens = elem_type(A)[]
+  for i in 1:length(b)
+    for j in 1:dim(A)
+      el = b[i]
+      if side == :left || side == :twosided
+        Bel = B[j] * el
+        if side == :twosided
+          for k in 1:dim(A)
+            push!(kgens, Bel * B[k])
+          end
+        else
+          push!(kgens, Bel)
+        end
+      end
+      if side == :right
+        push!(kgens, el * B[j])
+      end
+    end
+  end
+  return _ideal_from_kgens(A, kgens; side = side)
+end
 
-Returns the ideal of $A$ generated by the elements of `b` as a subspace of $A$.
-"""
-function ideal_from_gens(A::AbstractAssociativeAlgebra, b::Vector{T}, side::Symbol = :nothing) where { T <: AbstractAssociativeAlgebraElem }
+function _ideal_from_kgens(A::AbstractAssociativeAlgebra, b::Vector{<:AbstractAssociativeAlgebraElem}; side = nothing)
   if length(b) == 0
     M = zero_matrix(base_ring(A), 0, dim(A))
-    return ideal(A, M; side, M_in_rref=true)
+    return _ideal_from_matrix(A, M; side, M_in_rref=true)
   end
 
   @assert parent(b[1]) == A
@@ -303,60 +345,21 @@ function ideal_from_gens(A::AbstractAssociativeAlgebra, b::Vector{T}, side::Symb
   for i = 1:length(b)
     elem_to_mat_row!(M, i, b[i])
   end
-  return ideal(A, M; side)
+  return _ideal_from_matrix(A, M; side)
 end
 
-@doc raw"""
-    ideal(A::AbstractAssociativeAlgebra, x::AbstractAssociativeAlgebraElem) -> AbsAlgAssIdl
+left_ideal(A::AbstractAssociativeAlgebra, x...; kw...) = ideal(A, x...; side = :left, kw...)
 
-Returns the twosided principal ideal of $A$ generated by $x$.
-"""
-function ideal(A::AbstractAssociativeAlgebra, x::AbstractAssociativeAlgebraElem)
-  t1 = A()
-  t2 = A()
-  M = zero_matrix(base_ring(A), dim(A)^2, dim(A))
-  for i = 1:dim(A)
-    t1 = mul!(t1, A[i], x)
-    ii = (i - 1)*dim(A)
-    for j = 1:dim(A)
-      t2 = mul!(t2, t1, A[j])
-      elem_to_mat_row!(M, ii + j, t2)
-    end
-  end
+right_ideal(A::AbstractAssociativeAlgebra, x...; kw...) = ideal(A, x...; side = :right, kw...)
 
-  return ideal(A, M; side=:twosided)
-end
-
-@doc raw"""
-    ideal(A::AbstractAssociativeAlgebra, x::AbstractAssociativeAlgebraElem, action::Symbol) -> AbsAlgAssIdl
-
-Returns the ideal $x \cdot A$ if `action == :left`, and $A \cdot x$ if
-`action == :right`.
-"""
-function ideal(A::AbstractAssociativeAlgebra, x::AbstractAssociativeAlgebraElem, action::Symbol)
-  M = representation_matrix(x, action)
-  a = ideal(A, M)
+twosided_ideal(A::AbstractAssociativeAlgebra, x...; kw...) = ideal(A, x...; side = :twosided, kw...)
 
-  if action == :left
-    a.isright = 1
-  elseif action == :right
-    a.isleft = 1
-  end
+*(A::AbstractAssociativeAlgebra, x::NCRingElement) = left_ideal(A, x)
 
-  return a
-end
+*(x::NCRingElement, A::AbstractAssociativeAlgebra) = right_ideal(A, x)
 
 @doc raw"""
-    *(A::AbstractAssociativeAlgebra, x::AbstractAssociativeAlgebraElem) -> AbsAlgAssIdl
-    *(x::AbstractAssociativeAlgebraElem, A::AbstractAssociativeAlgebra) -> AbsAlgAssIdl
-
-Returns the ideal $A \cdot x$ or $x \cdot A$ respectively.
-"""
-*(A::AbstractAssociativeAlgebra, x::AbstractAssociativeAlgebraElem) = ideal(A, x, :right)
-*(x::AbstractAssociativeAlgebraElem, A::AbstractAssociativeAlgebra) = ideal(A, x, :left)
-
-@doc raw"""
-    ideal(A::AbstractAssociativeAlgebra, M::MatElem; side::Symbol = :nothing, M_in_rref::Bool = false)
+    _ideal_from_matrix(A::AbstractAssociativeAlgebra, M::MatElem; side::Symbol = :nothing, M_in_rref::Bool = false)
       -> AbsAlgAssIdl
 
 Returns the ideal of $A$ with basis matrix $M$.
@@ -365,13 +368,13 @@ set to `:right`/`:left`/`:twosided` respectively.
 If `M_in_rref == true`, it is assumed that $M$ is already in row reduced echelon
 form.
 """
-function ideal(A::AbstractAssociativeAlgebra, M::MatElem; side::Symbol = :nothing, M_in_rref::Bool = false)
+function _ideal_from_matrix(A::AbstractAssociativeAlgebra, M::MatElem; side = nothing, M_in_rref::Bool = false)
   @assert base_ring(M) == base_ring(A)
   @assert ncols(M) == dim(A)
   if !M_in_rref
     r, N = rref(M)
     if r == 0
-      a = AbsAlgAssIdl{typeof(A), typeof(M)}(A, zero_matrix(base_ring(A), 0, dim(A)))
+      a = AbsAlgAssIdl{typeof(A)}(A, zero_matrix(base_ring(A), 0, dim(A)))
       a.iszero = 1
       return a
     end
@@ -382,12 +385,12 @@ function ideal(A::AbstractAssociativeAlgebra, M::MatElem; side::Symbol = :nothin
     end
   end
   if M_in_rref && nrows(M) == 0
-    a = AbsAlgAssIdl{typeof(A), typeof(M)}(A, M)
+    a = AbsAlgAssIdl{typeof(A)}(A, M)
     a.iszero = 1
     return a
   end
 
-  a = AbsAlgAssIdl{typeof(A), typeof(M)}(A, M)
+  a = AbsAlgAssIdl{typeof(A)}(A, M)
   _set_sidedness(a, side)
   a.iszero = 2
   return a
@@ -395,7 +398,7 @@ end
 
 # Helper function to set the side-flags
 # side can be :right, :left or :twosided
-function _set_sidedness(a::Union{ AbsAlgAssIdl, AlgAssAbsOrdIdl, AlgAssRelOrdIdl }, side::Symbol)
+function _set_sidedness(a::Union{ AbsAlgAssIdl, AlgAssAbsOrdIdl, AlgAssRelOrdIdl }, side)
   if side == :right
     a.isleft = 0
     a.isright = 1
@@ -405,9 +408,11 @@ function _set_sidedness(a::Union{ AbsAlgAssIdl, AlgAssAbsOrdIdl, AlgAssRelOrdIdl
   elseif side == :twosided
     a.isleft = 1
     a.isright = 1
-  else
+  elseif side === nothing || side === :nothing
     a.isleft = 0
     a.isright = 0
+  else
+    error("Not a valid side")
   end
   return nothing
 end
@@ -419,17 +424,19 @@ end
 ################################################################################
 
 @doc raw"""
-    quo(A::AbstractAssociativeAlgebra, a::AbsAlgAssIdl) -> AbstractAssociativeAlgebra, AbsAlgAssMor
+    quo(A::AbstractAssociativeAlgebra, a::AbsAlgAssIdl)
+                                    -> AbstractAssociativeAlgebra, AbsAlgAssMor
 
-Returns the quotient algebra $A/a$ and the projection map $A \to A/a$.
+Return the quotient algebra $A/a$ and the projection map $A \to A/a$.
 """
-function quo(A::S, a::AbsAlgAssIdl{S, T, U}) where { S, T, U }
-  @assert A == algebra(a)
+function quo(A::S, a::AbsAlgAssIdl{S}) where {S}
+  @req A === algebra(a) "Ideal not in the algebra"
   K = base_ring(A)
+  d = dim(A)
 
   # First compute the vector space quotient
   Ma = basis_matrix(a, copy = false)
-  M = hcat(deepcopy(transpose(Ma)), identity_matrix(K, dim(A)))
+  M = hcat(transpose(Ma), identity_matrix(K, d))
   r = rref!(M)
   pivot_cols = Vector{Int}()
   j = 1
@@ -446,38 +453,37 @@ function quo(A::S, a::AbsAlgAssIdl{S, T, U}) where { S, T, U }
   end
 
   # We now have the basis (basis of the quotient, basis of the ideal)
-  n = dim(A) - nrows(Ma)
-  M = vcat(zero_matrix(K, n, dim(A)), Ma)
-  oneK = K(1)
-  zeroK = K()
+  n = d - nrows(Ma)
+  M = vcat(zero_matrix(K, n, d), Ma)
+  oneK = one(K)
+  zeroK = zero(K)
   for i = 1:n
     M[i, pivot_cols[i]] = oneK
   end
   iM = inv(M)
 
-  N = sub(M, 1:n, 1:dim(A))
-  NN = sub(iM, 1:dim(A), 1:n)
+  N = sub(M, 1:n, 1:d)
+  NN = sub(iM, 1:d, 1:n)
 
   # Lift a basis of the quotient to A
   quotient_basis = Vector{elem_type(A)}(undef, n)
-  b = zero_matrix(K, 1, n)
+  b = elem_type(K)[zero(K) for i in 1:n]
   for i = 1:n
-    b[1, i] = oneK
+    b[i] = oneK
     bN = b*N
-    quotient_basis[i] = A([ bN[1, i] for i = 1:dim(A) ])
-    b[1, i] = zeroK
+    quotient_basis[i] = A(bN; copy = true)
+    b[i] = zeroK
   end
 
   # Build the multiplication table
   t = A()
-  s = zero_matrix(K, 1, dim(A))
+  s = elem_type(K)[zero(K) for i in d]
   mult_table = Array{elem_type(K), 3}(undef, n, n, n)
-  for i = 1:n
-    for j = 1:n
+  for i in 1:n
+    for j in 1:n
       t = mul!(t, quotient_basis[i], quotient_basis[j])
-      elem_to_mat_row!(s, 1, t)
-      sNN = s*NN
-      mult_table[i, j, :] = [ sNN[1, k] for k = 1:n ]
+      sNN = coefficients(t, copy = false) * NN
+      mult_table[i, j, :] = sNN
     end
   end
 
@@ -493,15 +499,19 @@ end
 Given ideals $b \subseteq a$, this function returns the quotient algebra $a/b$
 and the projection map $a \to a/b$.
 """
-function quo(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where { S, T, U }
-  @assert algebra(a) == algebra(b)
+function quo(a::AbsAlgAssIdl{S}, b::AbsAlgAssIdl{S}) where {S}
+  @req algebra(a) === algebra(b) "Ideals must have same algebra"
+  @req is_subset(b, a) "Second ideal must be a subsets of the first ideal"
+  @req _test_ideal_sidedness(b, :twosided) "Second ideal must be two-sided"
+
   A = algebra(a)
+  d = dim(A)
   K = base_ring(A)
 
   # First compute the vector space quotient
   Ma = basis_matrix(a, copy = false)
   Mb = basis_matrix(b, copy = false)
-  M = hcat(deepcopy(transpose(Mb)), deepcopy(transpose(Ma)))
+  M = hcat(transpose(Mb), transpose(Ma))
   r = rref!(M)
   pivot_cols = Vector{Int}()
   j = 1
@@ -519,21 +529,21 @@ function quo(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where { S, T, U
 
   # Build the basis matrix for the quotient
   n = nrows(Ma) - nrows(Mb)
-  M = zero_matrix(K, n, dim(A))
+  M = zero_matrix(K, n, d)
   for i = 1:n
-    for j = 1:dim(A)
+    for j = 1:d
       M[i, j] = deepcopy(Ma[pivot_cols[i], j])
     end
   end
 
   # Lift a basis of the quotient to A
   quotient_basis = Vector{elem_type(A)}(undef, n)
-  b = zero_matrix(K, 1, n)
+  b = [zero(K) for i in 1:n]
   for i = 1:n
-    b[1, i] = one(K)
+    b[i] = one(K)
     bM = b*M
-    quotient_basis[i] = A([ bM[1, j] for j in 1:dim(A) ])
-    b[1, i] = zero(K)
+    quotient_basis[i] = A(bM; copy = true)
+    b[i] = zero(K)
   end
 
   # Another basis matrix for a: basis of the quotient + basis of b
@@ -541,14 +551,13 @@ function quo(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where { S, T, U
 
   # Build the multiplication table
   t = A()
-  s = zero_matrix(K, 1, dim(A))
   mult_table = Array{elem_type(K), 3}(undef, n, n, n)
+  Nctx = solve_init(N)
   for i = 1:n
     for j = 1:n
       t = mul!(t, quotient_basis[i], quotient_basis[j])
-      elem_to_mat_row!(s, 1, t)
-      y = solve(N, s, side = :left)
-      mult_table[i, j, :] = [ y[1, k] for k = 1:n ]
+      y = solve(Nctx, coefficients(t, copy = false), side = :left)
+      mult_table[i, j, :] = view(y, 1:n)
     end
   end
 
@@ -557,21 +566,17 @@ function quo(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where { S, T, U
   AtoB = AbsAlgAssMor{typeof(A), typeof(B), typeof(M)}(A, B)
 
   function _image(x::AbstractAssociativeAlgebraElem)
-    t, y = can_solve_with_solution(N, matrix(K, 1, dim(A), coefficients(x, copy = false)), side = :left)
+    t, y = can_solve_with_solution(Nctx, coefficients(x, copy = false), side = :left)
     if t
-      return B([ y[1, i] for i in 1:dim(B) ])
+      return B(y[1:dim(B)]; copy = false)
     else
       error("Element is not in the domain")
     end
   end
 
   function _preimage(x::AbstractAssociativeAlgebraElem)
-    t = zero_matrix(K, 1, dim(B))
-    for i = 1:dim(B)
-      t[1, i] = x.coeffs[i]
-    end
-    tt = t*M
-    return A([ tt[1, i] for i in 1:dim(A) ])
+    tt = coefficients(x, copy = false) * M
+    return A(tt; copy = false)
   end
 
   AtoB.header.image = _image
@@ -586,27 +591,35 @@ end
 ################################################################################
 
 # TODO: implement for ::Type{AbsAlgAssIdl}
-Random.gentype(a::AbsAlgAssIdl) = elem_type(algebra(a))
+#Random.gentype(a::AbsAlgAssIdl) = elem_type(algebra(a))
 
-function rand(rng::AbstractRNG, a_sp::Random.SamplerTrivial{<:AbsAlgAssIdl})
-  a = a_sp[]
-  A = algebra(a)
-  x = A()
-  for b in basis(a, copy = false)
-    x += rand(rng, base_ring(A))*b
+function RandomExtensions.maketype(I::AbsAlgAssIdl, _)
+  return elem_type(algebra(I))
+end
+
+function RandomExtensions.make(I::AbsAlgAssIdl, vs...)
+  R = base_ring(algebra(I))
+  if length(vs) == 1 && elem_type(R) == Random.gentype(vs[1])
+    RandomExtensions.Make(I, vs[1]) # forward to default Make constructor
+  else
+    RandomExtensions.Make(I, make(R, vs...))
   end
-  return x
 end
 
-function rand(a::AbsAlgAssIdl, rng::AbstractUnitRange{Int})
+function rand(rng::AbstractRNG, a_sp::Random.SamplerTrivial{<:Make2{<:NCRingElem, <:AbsAlgAssIdl}})
+  a, v = a_sp[][1:end]
   A = algebra(a)
   x = A()
   for b in basis(a, copy = false)
-    x += rand(base_ring(A), rng)*b
+    x = add!(x, rand(rng, v)*b)
   end
   return x
 end
 
+rand(rng::AbstractRNG, a::AbsAlgAssIdl, v...) = rand(rng, make(a, v...))
+
+rand(a::AbsAlgAssIdl, v...) = rand(Random.GLOBAL_RNG, a, v...)
+
 ################################################################################
 #
 #  Reduction of element modulo ideal
@@ -643,8 +656,8 @@ end
 #
 ################################################################################
 
-function left_principal_generator(a::AbsAlgAssIdl{S, T, U}) where { S <: MatAlgebra, T, U }
-  @assert is_left_ideal(a) "Not a left ideal"
+function left_principal_generator(a::AbsAlgAssIdl{S}) where {S <: MatAlgebra}
+  @req is_left_ideal(a) "Not a left ideal"
   A = algebra(a)
   if dim(A) != degree(A)^2*dim_of_coefficient_ring(A)
     error("Only implemented for full matrix algebras")
@@ -674,7 +687,7 @@ function left_principal_generator(a::AbsAlgAssIdl{S, T, U}) where { S <: MatAlge
   return x
 end
 
-function right_principal_generator(a::AbsAlgAssIdl{S, T, U}) where { S <: MatAlgebra, T, U }
+function right_principal_generator(a::AbsAlgAssIdl{S}) where {S <: MatAlgebra}
   @assert is_right_ideal(a) "Not a right ideal"
   A = algebra(a)
   if dim(A) != degree(A)^2*dim_of_coefficient_ring(A)
diff --git a/src/AlgAss/Map.jl b/src/AlgAss/Map.jl
index 34cb6cb0e8..83de797722 100644
--- a/src/AlgAss/Map.jl
+++ b/src/AlgAss/Map.jl
@@ -220,7 +220,8 @@ function hom(A::R, B::S, M::T; check = true) where {R <: AbstractAssociativeAlge
   return h
 end
 
-function hom(A::R, B::S, M::T, N::T) where {R <: AbstractAssociativeAlgebra, S <: AbstractAssociativeAlgebra, T <: MatElem}
+function hom(A::R, B::S, M::T, N::T; check = true) where {R <: AbstractAssociativeAlgebra, S <: AbstractAssociativeAlgebra, T <: MatElem}
+  # TODO: add check
   return AbsAlgAssMor{R, S, T}(A, B, M, N)
 end
 
diff --git a/src/AlgAss/Types.jl b/src/AlgAss/Types.jl
index 693b937e61..d69398e353 100644
--- a/src/AlgAss/Types.jl
+++ b/src/AlgAss/Types.jl
@@ -236,20 +236,21 @@ end
 
 # S is the type of the algebra, T = elem_type(S) and U is the type of matrices
 # over the base ring of the algebra
-mutable struct AbsAlgAssIdl{S, T, U}
+mutable struct AbsAlgAssIdl{S}
   algebra::S
-  basis::Vector{T}
-  basis_matrix::U
+  basis#::Vector{elem_type(algebra)}
+  basis_matrix#::dense_matrix_type(base_ring_Type(A))
+  basis_matrix_solve_ctx#solve_context_type(...)
 
   isleft::Int                      # 0 Not known
                                    # 1 Known to be a left ideal
                                    # 2 Known not to be a left ideal
   isright::Int                     # as for isleft
 
-  iszero::Int
+  iszero::Int                      # as for isleft
 
-  function AbsAlgAssIdl{S, T, U}(A::S) where {S, T, U}
-    I = new{S, T, U}()
+  function AbsAlgAssIdl{S}(A::S) where {S}
+    I = new{S}()
     I.algebra = A
     I.isleft = 0
     I.isright = 0
@@ -257,8 +258,8 @@ mutable struct AbsAlgAssIdl{S, T, U}
     return I
   end
 
-  function AbsAlgAssIdl{S, U}(A::S, M::U) where {S, U}
-    I = new{S, elem_type(S), U}()
+  function AbsAlgAssIdl{S}(A::S, M::MatElem) where {S}
+    I = new{S}()
     I.algebra = A
     I.basis_matrix = M
     I.isleft = 0
diff --git a/src/AlgAss/radical.jl b/src/AlgAss/radical.jl
index 42704a8de0..b6cc2ec2be 100644
--- a/src/AlgAss/radical.jl
+++ b/src/AlgAss/radical.jl
@@ -32,7 +32,7 @@
 Returns the Jacobson radical of $A$.
 """
 function radical(A::AbstractAssociativeAlgebra)
-  return ideal_from_gens(A, _radical(A), :twosided)
+  return _ideal_from_kgens(A, _radical(A); side = :twosided)
 end
 
 function _radical(A::AbstractAssociativeAlgebra)
diff --git a/src/AlgAssAbsOrd/Ideal.jl b/src/AlgAssAbsOrd/Ideal.jl
index 065dc110eb..71489f799d 100644
--- a/src/AlgAssAbsOrd/Ideal.jl
+++ b/src/AlgAssAbsOrd/Ideal.jl
@@ -1090,7 +1090,7 @@ function _islocally_free_left(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl, p::Union{Int,
       push!(gensJ, toIpI(bb*c))
     end
   end
-  JinIpI = ideal_from_gens(IpI, gensJ)
+  JinIpI = _ideal_from_kgens(IpI, gensJ)
   IJ, toIJ = quo(IpI, JinIpI)
 
   a = O()
@@ -2049,7 +2049,7 @@ function maximal_integral_ideal_containing(I::AlgAssAbsOrdIdl, p::Union{ ZZRingE
   B, 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) ])
+  JinC = _ideal_from_kgens(C, [ toC(OPtoB(toOP(O(b)))) for b in absolute_basis(J) ])
   y = left_principal_generator(JinC)
   m = matrix(y)
   r = rref!(m)
diff --git a/src/AlgAssRelOrd/Eichler.jl b/src/AlgAssRelOrd/Eichler.jl
index dbc94a9839..d6451446ca 100644
--- a/src/AlgAssRelOrd/Eichler.jl
+++ b/src/AlgAssRelOrd/Eichler.jl
@@ -65,8 +65,8 @@ function _eichler_find_transforming_unit_maximal(M::T, N::T) where { T <: Union{
   OpO, toOpO = quo(O, p*O, p)
   B, toB = _as_matrix_algebra(OpO)
 
-  I = ideal_from_gens(B, [ toB(toOpO(O(b))) for b in absolute_basis(M) ])
-  J = ideal_from_gens(B, [ toB(toOpO(O(b))) for b in absolute_basis(N) ])
+  I = _ideal_from_kgens(B, [ toB(toOpO(O(b))) for b in absolute_basis(M) ])
+  J = _ideal_from_kgens(B, [ toB(toOpO(O(b))) for b in absolute_basis(N) ])
 
   # Compute the image of 1 under the canonical projections O -> O/M respectively O -> O/N
   Fq = base_ring(B)
@@ -83,7 +83,7 @@ function _eichler_find_transforming_unit_maximal(M::T, N::T) where { T <: Union{
       break
     end
   end
-  v = matrix(Fq, degree(B), 1, [ vv[i, nonZeroCol] for i = 1:degree(B) ])
+  v = matrix(Fq, degree(B), 1, elem_type(Fq)[ vv[i, nonZeroCol] for i = 1:degree(B) ])
   ww = mod(one(B), J)
   nonZeroCol = 0
   for i = 1:degree(B)
@@ -97,7 +97,7 @@ function _eichler_find_transforming_unit_maximal(M::T, N::T) where { T <: Union{
       break
     end
   end
-  w = matrix(Fq, degree(B), 1, [ ww[i, nonZeroCol] for i = 1:degree(B) ])
+  w = matrix(Fq, degree(B), 1, elem_type(Fq)[ ww[i, nonZeroCol] for i = 1:degree(B) ])
 
   b = ceil(Int, degree(base_ring(B))*dim(B)*log2(BigInt(characteristic(base_ring(B)))))
   # A minimal set of generators of around b elements should generate B^\times
diff --git a/src/AlgAssRelOrd/Ideal.jl b/src/AlgAssRelOrd/Ideal.jl
index 78eaa52f12..9cc4b5ce24 100644
--- a/src/AlgAssRelOrd/Ideal.jl
+++ b/src/AlgAssRelOrd/Ideal.jl
@@ -1096,7 +1096,7 @@ function is_locally_free(O::AlgAssRelOrd, I::AlgAssRelOrdIdl, p::Union{ AbsNumFi
       push!(gensJ, toIpI(bb*c))
     end
   end
-  JinIpI = ideal_from_gens(IpI, gensJ)
+  JinIpI = _ideal_from_kgens(IpI, gensJ)
   IJ, toIJ = quo(IpI, JinIpI)
 
   a = O()
@@ -1352,7 +1352,7 @@ function maximal_integral_ideal_containing(I::AlgAssRelOrdIdl, p::Union{ AbsNumF
   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) ])
+  JinC = _ideal_from_kgens(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)
diff --git a/src/exports.jl b/src/exports.jl
index da71484562..b5e6d6d2bd 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -613,6 +613,7 @@ export lattice_with_local_conditions
 export lcm
 export leading_coefficient
 export leech_lattice
+export left_ideal
 export left_order
 export level
 export lll
@@ -801,6 +802,7 @@ export restrict_scalars
 export restrict_scalars_with_map
 export resultant
 export riemann_roch_space
+export right_ideal
 export right_order
 export ring_of_integers
 export ring_of_multipliers
@@ -911,6 +913,7 @@ export trivial_morphism
 export trred
 export trred_matrix
 export twists
+export twosided_ideal
 export uniformizer
 export unit_group
 export unit_group_fac_elem
diff --git a/test/AlgAss/AbsAlgAss.jl b/test/AlgAss/AbsAlgAss.jl
index abbd7fa066..1051cffd86 100644
--- a/test/AlgAss/AbsAlgAss.jl
+++ b/test/AlgAss/AbsAlgAss.jl
@@ -112,7 +112,7 @@
     A = StructureConstantAlgebra(f2g3)
     fg = A(QQFieldElem[-5, 5, -2, 6, 3, 1, 0, 0, 0, 0, 0, 0, 0]) # = f*g
     J = radical(A)
-    I = ideal(A, fg)
+    I = left_ideal(A, fg)
     @test I == J
 
     f = x^2 + 1
@@ -124,7 +124,7 @@
     A = StructureConstantAlgebra(g2h3)
     gh = A(map(K, [10, -5, -28, -13, 2, 1, 0, 0, 0, 0, 0, 0])) # = g*h
     J = radical(A)
-    I = ideal(A, gh)
+    I = left_ideal(A, gh)
     @test I == J
 
     G = small_group(8, 4)
@@ -138,7 +138,7 @@
             0 0 0 0 1 0 0 1;
             0 0 0 0 0 1 0 1;
             0 0 0 0 0 0 1 1]
-    @test I == ideal(A, bI)
+    @test I == Hecke._ideal_from_matrix(A, bI)
     ge = [A(g) - A(one(G)) for g in G]
     @test all(in(I), ge)
     AS, AStoA = StructureConstantAlgebra(A)
diff --git a/test/AlgAss/Ideal.jl b/test/AlgAss/Ideal.jl
index c2bef42730..c337e4b6c9 100644
--- a/test/AlgAss/Ideal.jl
+++ b/test/AlgAss/Ideal.jl
@@ -1,11 +1,47 @@
 @testset "Ideals in algebras" begin
 
+  let
+    A = matrix_algebra(GF(3), 2)
+    I = A[2]*A
+    @test algebra(I) === A
+    @test ideal(A, [A[2]], side = :right) == I
+    I = A*A[2]
+    @test algebra(I) === A
+    @test ideal(A, [A[2]], side = :left) == I
+    I = A[2]*A + A[3]*A
+    @test algebra(I) === A
+    @test ideal(A, [A[2], A[3]], side = :right) == I
+    I = A*A[2] + A*A[3]
+    @test algebra(I) === A
+    @test ideal(A, [A[2], A[3]], side = :left) == I
+
+    sprint(show, "text/plain", I) isa String
+    sprint(show, I) isa String
+    @test !is_zero(I)
+    @test is_one(I)
+
+    @test basis(I) !== basis(I; copy = false)
+
+    I = A*zero(A)
+    @test algebra(I) === A
+    @test ideal(A, elem_type(A)[], side = :left) == I
+    @test is_zero(I)
+    @test !is_one(I)
+  end
+
+
   @testset "Left / Right" begin
     A = matrix_algebra(GF(3), 2)
     I = A[2]*A
 
-    @test is_left_ideal(I) == false
-    @test is_right_ideal(I) == true
+    @test !is_left_ideal(I)
+    @test is_right_ideal(I)
+
+    J = ideal(A, [A[1]]; side = :twosided)
+    @test J == one(A) * A
+
+    J = twosided_ideal(A, [A[2]])
+    @test J == one(A) * A
   end
 
   @testset "Quotients" begin
@@ -70,8 +106,8 @@
     # see https://github.com/thofma/Hecke.jl/issues/1399 .
     G = small_group(8, 4)
     A, _ = StructureConstantAlgebra(group_algebra(FlintQQ, G))
-    I = ideal(A, one(A))
-    J = ideal(A, sum(A[i] for i in 1:8))
+    I = left_ideal(A, one(A))
+    J = left_ideal(A, sum(A[i] for i in 1:8))
     Q, AtoQ = quo(I, J)
 
     @test dim(Q) == 7
@@ -108,11 +144,20 @@
     E = elem_type(A)
     @test rand(I) isa E
     @test rand(rng, I) isa E
-    @test rand(rng, I, 2, 3) isa Matrix{E}
+    @test rand(rng, Hecke.RandomExtensions.make(I), 2, 3) isa Matrix{E}
 
     Random.seed!(rng, rand_seed)
     a = rand(rng, I)
     Random.seed!(rng, rand_seed)
     @test a == rand(rng, I)
   end
+
+  let
+    A = matrix_algebra(GF(3), 2)
+    J = Hecke.ideal(A, [A[1]]; side = :twosided)
+    a = @inferred Hecke.left_principal_generator(J)
+    @test A * a == J
+    a = @inferred Hecke.right_principal_generator(J)
+    @test a * A == J
+  end
 end
diff --git a/test/AlgAss/radical.jl b/test/AlgAss/radical.jl
index 6961a7f647..748e820d8b 100644
--- a/test/AlgAss/radical.jl
+++ b/test/AlgAss/radical.jl
@@ -6,7 +6,7 @@
   A = StructureConstantAlgebra(f2g3)
   fg = A(QQFieldElem[-5, 5, -2, 6, 3, 1, 0, 0, 0, 0, 0, 0, 0]) # = f*g
   J = radical(A)
-  I = ideal(A, fg)
+  I = left_ideal(A, fg)
   @test I == J
 
   f = x^2 + 1
@@ -18,7 +18,7 @@
   A = StructureConstantAlgebra(g2h3)
   gh = A(map(K, [10, -5, -28, -13, 2, 1, 0, 0, 0, 0, 0, 0])) # = g*h
   J = radical(A)
-  I = ideal(A, gh)
+  I = left_ideal(A, gh)
   @test I == J
 
   G = small_group(8, 4)
@@ -32,7 +32,7 @@
            0 0 0 0 1 0 0 1;
            0 0 0 0 0 1 0 1;
            0 0 0 0 0 0 1 1]
-    @test I == ideal(A, bI)
+    @test I == Hecke._ideal_from_matrix(A, bI)
     ge = [A(g) - A(one(G)) for g in G]
     @test all(in(I), ge)
     AS, AStoA = StructureConstantAlgebra(A)
@@ -65,6 +65,6 @@
   for K in [ GF(2), GF(4), Native.GF(2), Native.GF(2, 2), QQ, rationals_as_number_field()[1]]
     A = matrix_algebra(K, [ matrix(K, 2, 2, [ 1, 0, 0, 0 ]), matrix(K, 2, 2, [ 0, 1, 0, 0 ]), matrix(K, 2, 2, [ 0, 0, 0, 1]) ]) # i. e. upper triangular matrices
     I = radical(A)
-    @test nrows(basis_matrix(I, copy = false)) == 1
+    @test dim(I) == 1
   end
 end