diff --git a/docs/src/quad_forms/Zgenera.md b/docs/src/quad_forms/Zgenera.md
index 0aa864d8ae..d83fa5c693 100644
--- a/docs/src/quad_forms/Zgenera.md
+++ b/docs/src/quad_forms/Zgenera.md
@@ -41,7 +41,7 @@ Zgenera(sig_pair::Tuple{Int,Int}, determinant::Union{Int,ZZRingElem})
 ```
 ### From other genus symbols
 ```@docs
-orthogonal_sum(G1::ZGenus, G2::ZGenus)
+direct_sum(G1::ZGenus, G2::ZGenus)
 ```
 
 ## Attributes of the genus
@@ -127,9 +127,9 @@ gram_matrix(S::ZpGenus)
 rescale(S::ZpGenus, a::RationalUnion)
 ```
 
-### Orthogonal sums
+### Direct sums
 ```@docs
-orthogonal_sum(S1::ZpGenus, S2::ZpGenus)
+direct_sum(S1::ZpGenus, S2::ZpGenus)
 ```
 
 ### Embeddings/Representations
diff --git a/docs/src/quad_forms/basics.md b/docs/src/quad_forms/basics.md
index 024aff8d97..3c651b3401 100644
--- a/docs/src/quad_forms/basics.md
+++ b/docs/src/quad_forms/basics.md
@@ -211,13 +211,35 @@ is_represented_by(H2, H)
 ```
 ---
 
+## Categorical constructions
+One can construct direct sums of spaces of the same kind. Since those are also
+direct products, they are called biproducts in this context. Depending on the user
+usage, one of the following three methods can be called to obtain the direct
+sum of a finite collection of spaces. Note that the corresponding copies
+of the original spaces in the direct sum are pairwise orthogonal.
+
+```@docs
+direct_sum(x::Vector{AbstractSpace})
+direct_product(x::Vector{AbstractSpace})
+biproduct(x::Vector{AbstractSpace})
+```
+
+### Example
+
+```@repl 2
+using Hecke # hide
+E, b = cyclotomix_field_as_cm_extensions(7);
+H = hermitian_space(E, 3);
+H2 = hermitian_space(E, E[-1 0 0; 0 1 0; 0 0 -1]);
+H3, inj, proj = biproduct(H, H2)
+isone(compose(inj[1], proj[1]))
+iszero(compose(inj[1], proj[2]))
+```
 ## Orthogonality operations
 
 ```@docs
 orthogonal_complement(::AbstractSpace, ::MatElem)
 orthogonal_projection(::AbstractSpace, ::MatElem)
-orthogonal_sum(::AbstractSpace, ::AbstractSpace)
-direct_sum(x::Vararg{QuadSpace})
 ```
 
 ### Example
@@ -226,15 +248,8 @@ direct_sum(x::Vararg{QuadSpace})
 using Hecke # hide
 K, a = CyclotomicRealSubfield(7);
 Kt, t = K["t"];
-E, b = number_field(t^2-a*t+1, "b");
 Q = quadratic_space(K, K[0 1; 1 0]);
-H = hermitian_space(E, 3);
-H2 = hermitian_space(E, E[-1 0 0; 0 1 0; 0 0 -1]);
 orthogonal_complement(Q, matrix(K, 1, 2, [1 0]))
-H3, map1, map2 = orthogonal_sum(H, H2);
-H3
-map1
-map2
 ```
 ---
 
diff --git a/docs/src/quad_forms/discriminant_group.md b/docs/src/quad_forms/discriminant_group.md
index a348e22a9e..694a428097 100644
--- a/docs/src/quad_forms/discriminant_group.md
+++ b/docs/src/quad_forms/discriminant_group.md
@@ -123,8 +123,9 @@ brown_invariant(T::TorQuadModule)
 is_genus(T::TorQuadModule, signature_pair::Tuple{Int, Int})
 ```
 
-### Orthogonal sums
+### Categorical constructions
 ```@docs
-orthogonal_sum(T::TorQuadModule, U::TorQuadModule)
-direct_sum(x::Vararg{TorQuadModule})
+direct_sum(x::Vector{TorQuadModule})
+direct_product(x::Vector{TorQuadModule})
+biproduct(x::Vector{TorQuadModule})
 ```
diff --git a/docs/src/quad_forms/genusherm.md b/docs/src/quad_forms/genusherm.md
index df03770985..edfc55c83a 100644
--- a/docs/src/quad_forms/genusherm.md
+++ b/docs/src/quad_forms/genusherm.md
@@ -435,8 +435,8 @@ length(representatives(G1))
 ## Sum of genera
 
 ```@docs
-orthogonal_sum(::HermLocalGenus, ::HermLocalGenus)
-orthogonal_sum(::HermGenus, ::HermGenus)
+direct_sum(::HermLocalGenus, ::HermLocalGenus)
+direct_sum(::HermGenus, ::HermGenus)
 ```
 
 ### Examples
@@ -458,8 +458,8 @@ gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [1, 0, 0]), map(E, [a, 0, 0]), ma
 L = hermitian_lattice(E, gens, gram = D);
 g2 = genus(L, p);
 G2 = genus(L);
-orthogonal_sum(g1, g2)
-orthogonal_sum(G1, G2)
+direct_sum(g1, g2)
+direct_sum(G1, G2)
 ```
 
 ---
diff --git a/docs/src/quad_forms/integer_lattices.md b/docs/src/quad_forms/integer_lattices.md
index 55c33e409b..a298b711fa 100644
--- a/docs/src/quad_forms/integer_lattices.md
+++ b/docs/src/quad_forms/integer_lattices.md
@@ -132,11 +132,16 @@ divisibility(::ZLat, ::Union{Vector, QQMatrix})
 
 ## Embeddings
 
+### Categorical constructions
+```@docs
+direct_sum(x::Vector{ZLat})
+direct_product(x::Vector{ZLat})
+biproduct(x::Vector{ZLat})
+```
+
 ### Orthogonal sublattices
 ```@docs
-orthogonal_sum(::ZLat, ::ZLat)
 orthogonal_submodule(::ZLat, ::ZLat)
-direct_sum(x::Vararg{ZLat})
 irreducible_components(::ZLat)
 ```
 
diff --git a/docs/src/quad_forms/lattices.md b/docs/src/quad_forms/lattices.md
index e277886766..ddf8425f2d 100644
--- a/docs/src/quad_forms/lattices.md
+++ b/docs/src/quad_forms/lattices.md
@@ -256,6 +256,17 @@ ambient_space(rescale(Lquad,3*a))
 pseudo_matrix(Lquad)
 ```
 
+## Categorical constructions
+Given finite collections of lattices, one can construct their direct sums, which
+are also direct products in this context. They are also sometimes called biproducts.
+Depending on the user usage, it is possible to call one of the following functions.
+
+```@docs
+direct_sum(x::Vector{AbstractLat})
+direct_product(x::Vector{AbstractLat})
+biproduct(x::Vector{AbstractLat})
+```
+
 ---
 
 ## Invariants
diff --git a/src/Deprecations.jl b/src/Deprecations.jl
index e92d36a1bd..59853de087 100644
--- a/src/Deprecations.jl
+++ b/src/Deprecations.jl
@@ -32,6 +32,10 @@
 
 @deprecate automorphisms(x::LocalField, y::Union{FlintPadicField, FlintQadicField, LocalField}) automorphism_list(x, y)
 
+# Deprecated during 0.18.*
+
+@deprecate orthogonal_sum(x::T, y::T) where T <: Union{AbstractSpace, ZGenus, ZpGenus, HermGenus, HermLocalGenus, QuadGenus, QuadLocalGenus, JorDec, LocalQuadSpaceCls, QuadSpaceCls} direct_sum(x, y)
+
 # Things that moved to Nemo
 
 # > 0.18.1
diff --git a/src/GrpAb/GrpAbFinGen.jl b/src/GrpAb/GrpAbFinGen.jl
index bd9e48b1e8..8624c0fb05 100644
--- a/src/GrpAb/GrpAbFinGen.jl
+++ b/src/GrpAb/GrpAbFinGen.jl
@@ -39,7 +39,7 @@ export abelian_group, free_abelian_group, is_snf, ngens, nrels, rels, snf, isfin
        direct_product, is_torsion, torsion_subgroup, sub, quo, is_cyclic,
        psylow_subgroup, is_subgroup, abelian_groups, flat, tensor_product,
        dual, chain_complex, is_exact, free_resolution, obj, map,
-       primary_part, is_free, is_pure, is_neat
+       primary_part, is_free, is_pure, is_neat, direct_sum, biproduct
 
 import Base.+, Nemo.snf, Nemo.parent, Base.rand, Nemo.is_snf
 
@@ -286,7 +286,7 @@ end
 @doc Markdown.doc"""
     is_snf(G::GrpAbFinGen) -> Bool
 
-Returns whether the current relation matrix of the group $G$ is in Smith
+Return whether the current relation matrix of the group $G$ is in Smith
 normal form.
 """
 is_snf(A::GrpAbFinGen) = A.is_snf
@@ -294,7 +294,7 @@ is_snf(A::GrpAbFinGen) = A.is_snf
 @doc Markdown.doc"""
     ngens(G::GrpAbFinGen) -> Int
 
-Returns the number of generators of $G$ in the current representation.
+Return the number of generators of $G$ in the current representation.
 """
 function ngens(A::GrpAbFinGen)
   if is_snf(A)
@@ -307,7 +307,7 @@ end
 @doc Markdown.doc"""
     nrels(G::GrpAbFinGen) -> Int
 
-Returns the number of relations of $G$ in the current representation.
+Return the number of relations of $G$ in the current representation.
 """
 function nrels(A::GrpAbFinGen)
   if is_snf(A)
@@ -320,7 +320,7 @@ end
 @doc Markdown.doc"""
     rels(A::GrpAbFinGen) -> ZZMatrix
 
-Returns the currently used relations of $G$ as a single matrix.
+Return the currently used relations of $G$ as a single matrix.
 """
 function rels(A::GrpAbFinGen)
   if is_snf(A)
@@ -381,7 +381,7 @@ end
 @doc Markdown.doc"""
     snf(A::GrpAbFinGen) -> GrpAbFinGen, GrpAbFinGenMap
 
-Returns a pair $(G, f)$, where $G$ is an abelian group in canonical Smith
+Return a pair $(G, f)$, where $G$ is an abelian group in canonical Smith
 normal form isomorphic to $A$ and an isomorphism $f : G \to A$.
 """
 function snf(G::GrpAbFinGen)
@@ -458,7 +458,7 @@ end
 @doc Markdown.doc"""
     isfinite(A::GrpAbFinGen) -> Bool
 
-Returns whether $A$ is finite.
+Return whether $A$ is finite.
 """
 isfinite(A::GrpAbFinGen) = is_snf(A) ? is_finite_snf(A) : is_finite_gen(A)
 
@@ -475,7 +475,7 @@ is_finite_gen(A::GrpAbFinGen) = isfinite(snf(A)[1])
 @doc Markdown.doc"""
     rank(A::GrpAbFinGen) -> Int
 
-Returns the rank of $A$, that is, the dimension of the
+Return the rank of $A$, that is, the dimension of the
 $\mathbf{Q}$-vectorspace $A \otimes_{\mathbf Z} \mathbf Q$.
 """
 rank(A::GrpAbFinGen) = is_snf(A) ? rank_snf(A) : rank_gen(A)
@@ -507,7 +507,7 @@ end
 @doc Markdown.doc"""
     order(A::GrpAbFinGen) -> ZZRingElem
 
-Returns the order of $A$. It is assumed that $A$ is finite.
+Return the order of $A$. It is assumed that $A$ is finite.
 """
 order(A::GrpAbFinGen) = is_snf(A) ? order_snf(A) : order_gen(A)
 
@@ -527,7 +527,7 @@ order_gen(A::GrpAbFinGen) = order(snf(A)[1])
 @doc Markdown.doc"""
     exponent(A::GrpAbFinGen) -> ZZRingElem
 
-Returns the exponent of $A$. It is assumed that $A$ is finite.
+Return the exponent of $A$. It is assumed that $A$ is finite.
 """
 function exponent(A::GrpAbFinGen)
   if is_snf(A)
@@ -562,7 +562,7 @@ exponent_gen(A::GrpAbFinGen) = exponent(snf(A)[1])
 @doc Markdown.doc"""
     istrivial(A::GrpAbFinGen) -> Bool
 
-Checks if $A$ is the trivial group.
+Return whether $A$ is the trivial group.
 """
 istrivial(A::GrpAbFinGen) = isfinite(A) && isone(order(A))
 
@@ -575,7 +575,7 @@ istrivial(A::GrpAbFinGen) = isfinite(A) && isone(order(A))
 @doc Markdown.doc"""
     is_isomorphic(G::GrpAbFinGen, H::GrpAbFinGen) -> Bool
 
-Checks if $G$ and $H$ are isomorphic.
+Return whether $G$ and $H$ are isomorphic.
 """
 function is_isomorphic(G::GrpAbFinGen, H::GrpAbFinGen)
   b = filter(x -> x != 1, snf(G)[1].snf) == filter(x -> x != 1, snf(H)[1].snf)
@@ -584,37 +584,85 @@ end
 
 ################################################################################
 #
-#  Direct product
+#  Direct products/direct sums/biproducts
 #
 ################################################################################
 #TODO: check the universal properties here!!!
+
 @doc Markdown.doc"""
-    direct_product(G::GrpAbFinGen...; task::Symbol = :prod) -> GrpAbFinGen, GrpAbFinGenMap, GrpAbFinGenMap
+    direct_sum(G::GrpAbFinGen...) -> GrpAbFinGen, Vector{GrpAbFinGenMap}
+
+Return the direct sum $D$ of the (finitely many) abelian groups $G_i$, together
+with the injections $G_i \to D$.
 
-Returns the direct product $D$ of the abelian groups $G_i$. `task` can be
-":sum", ":prod", ":both" or ":none" and determines which canonical maps
-are computed as well: ":sum" for the injections, ":prod" for the
-projections.
+For finite abelian groups, finite direct sums and finite direct products agree and
+they are therefore called biproducts.
+If one wants to obtain $D$ as a direct product together with the projections
+$ D \to G_i$, one should call `direct_product(G...)`.
+If one wants to obtain $D$ as a biproduct together with the projections and the
+injections, one should call `biproduct(G...)`.
+
+Otherwise, one could also call `canonical_injections(D)` or `canonical_projections(D)`
+later on.
 """
-function direct_product(G::GrpAbFinGen...  ; task::Symbol = :prod, kwargs...)
+function direct_sum(G::GrpAbFinGen...; task::Symbol = :sum, kwargs...)
+  @assert task in [:sum, :prod, :both, :none]
+  return _direct_product(:sum, G...; task = task, kwargs...)
+end
+
+@doc Markdown.doc"""
+    direct_product(G::GrpAbFinGen...) -> GrpAbFinGen, Vector{GrpAbFinGenMap}
+
+Return the direct product $D$ of the (finitely many) abelian groups $G_i$, together
+with the projections $D \to G_i$.
+
+For finite abelian groups, finite direct sums and finite direct products agree and
+they are therefore called biproducts.
+If one wants to obtain $D$ as a direct sum together with the injections $D \to G_i$,
+one should call `direct_sum(G...)`.
+If one wants to obtain $D$ as a biproduct together with the projections and the
+injections, one should call `biproduct(G...)`.
+
+Otherwise, one could also call `canonical_injections(D)` or `canonical_projections(D)`
+later on.
+"""
+function direct_product(G::GrpAbFinGen...; task::Symbol = :prod, kwargs...)
   @assert task in [:prod, :sum, :both, :none]
   return _direct_product(:prod, G...; task = task, kwargs...)
 end
 
 @doc Markdown.doc"""
-    direct_sum(G::GrpAbFinGen...; task::Symbol = :sum) -> GrpAbFinGen, GrpAbFinGenMap, GrpAbFinGenMap
+    biproduct(G::GrpAbFinGen...) -> GrpAbFinGen, Vector{GrpAbFinGenMap}, Vector{GrpAbFinGenMap}
+
+Return the direct product $D$ of the (finitely many) abelian groups $G_i$, together
+with the projections $D \to G_i$ and the injections $G_i \to D$.
+
+For finite abelian groups, finite direct sums and finite direct products agree and
+they are therefore called biproducts.
+If one wants to obtain $D$ as a direct sum together with the injections $G_i \to D$,
+one should call `direct_sum(G...)`.
+If one wants to obtain $D$ as a direct product together with the projections $D \to G_i$, 
+one should call `direct_product(G...)`.
 
-Returns the direct sum $D$ of the abelian groups $G_i$. `task` can be
-":sum", ":prod", ":both" or ":none" and determines which canonical maps
-are computed as well: ":sum" for the injections, ":prod" for the
-projections.
+Otherwise, one could also call `canonical_injections(D)` or `canonical_projections(D)`
+later on.
 """
-function direct_sum(_G::GrpAbFinGen, Gs::GrpAbFinGen...  ; task::Symbol = :sum, kwargs...)
-  G = (_G, Gs...)
+function biproduct(G::GrpAbFinGen...; task::Symbol = :both, kwargs...)
   @assert task in [:prod, :sum, :both, :none]
-  return _direct_product(:sum, G...; task = task, kwargs...)
+  return _direct_product(:prod, G...; task = task, kwargs...)
 end
 
+@doc Markdown.doc"""
+    ⊕(A::GrpAbFinGen...) -> GrpAbFinGen
+
+Return the direct sum $D$ of the (finitely many) abelian groups $G_i$.
+
+If one wants to access the injections $G_i \to D$ or the projections $D \to G_i$ later,
+one can call respectively `canonical_injections(D)` or `canonical_projections(D)`.
+"""
+⊕(A::GrpAbFinGen...) = direct_sum(A..., task = :none)
+export ⊕
+
 function _direct_product(t::Symbol, G::GrpAbFinGen...
              ; add_to_lattice::Bool = false, L::GrpAbLattice = GroupLattice, task::Symbol = :sum)
   @assert task in [:prod, :sum, :both, :none]
@@ -660,11 +708,8 @@ function _direct_product(t::Symbol, G::GrpAbFinGen...
   end
 end
 
-⊕(A::GrpAbFinGen...) = direct_sum(A..., task = :none)
-export ⊕
-
 @doc Markdown.doc"""
-    canonical_injections(G::GrpAbFinGen) -> Vector{Map}
+    canonical_injections(G::GrpAbFinGen) -> Vector{GrpAbFinGenMap}
 
 Given a group $G$ that was created as a direct product, return the
 injections from all components.
@@ -676,7 +721,7 @@ function canonical_injections(G::GrpAbFinGen)
 end
  
 @doc Markdown.doc"""
-    canonical_injection(G::GrpAbFinGen, i::Int) -> Map
+    canonical_injection(G::GrpAbFinGen, i::Int) -> GrpAbFinGenMap
 
 Given a group $G$ that was created as a direct product, return the
 injection from the $i$th component.
@@ -690,10 +735,10 @@ function canonical_injection(G::GrpAbFinGen, i::Int)
 end
 
 @doc Markdown.doc"""
-    canonical_projections(G::GrpAbFinGen) -> Vector{Map}
+    canonical_projections(G::GrpAbFinGen) -> Vector{GrpAbFinGenMap}
 
 Given a group $G$ that was created as a direct product, return the
-projections into all components.
+projections onto all components.
 """
 function canonical_projections(G::GrpAbFinGen)
   D = get_attribute(G, :direct_product)
@@ -702,7 +747,7 @@ function canonical_projections(G::GrpAbFinGen)
 end
  
 @doc Markdown.doc"""
-    canonical_projection(G::GrpAbFinGen, i::Int) -> Map
+    canonical_projection(G::GrpAbFinGen, i::Int) -> GrpAbFinGenMap
 
 Given a group $G$ that was created as a direct product, return the
 projection onto the $i$th component.
@@ -774,13 +819,12 @@ end
 
 
 @doc Markdown.doc"""
-    flat(G::GrpAbFinGen) -> GrpAbFinGen, Map
+    flat(G::GrpAbFinGen) -> GrpAbFinGenMap
 
 Given a group $G$ that is created using (iterated) direct products, or
-(iterated) tensor product,
-return a group that is a flat product: $(A \oplus B) \oplus C$
-is returned as $A \oplus B \oplus C$, (resp. $\otimes$)
-together with the  isomorphism.
+(iterated) tensor products, return a group isomorphism into a flat product:
+for $G := (A \oplus B) \oplus C$, it returns the isomorphism
+$G \to A \oplus B \oplus C$ (resp. $\otimes$).
 """
 function flat(G::GrpAbFinGen)
   s = get_attribute(G, :direct_product)
@@ -903,14 +947,14 @@ end
 @doc Markdown.doc"""
     is_torsion(G::GrpAbFinGen) -> Bool
 
-Returns true if and only if `G` is a torsion group.
+Return whether `G` is a torsion group.
 """
 is_torsion(G::GrpAbFinGen) = isfinite(G)
 
 @doc Markdown.doc"""
     torsion_subgroup(G::GrpAbFinGen) -> GrpAbFinGen, Map
 
-Returns the torsion subgroup of `G`.
+Return the torsion subgroup of `G`.
 """
 function torsion_subgroup(G::GrpAbFinGen, add_to_lattice::Bool = true,
                                           L::GrpAbLattice = GroupLattice)
@@ -933,7 +977,7 @@ end
 @doc Markdown.doc"""
     is_free(G::GrpAbFinGen) -> Bool
 
-Returns whether `G` is free or not.
+Returns whether `G` is free.
 """
 function is_free(G::GrpAbFinGen)
   T, = torsion_subgroup(G, false)
@@ -947,7 +991,7 @@ end
 ##############################################################################
 
 @doc Markdown.doc"""
-    sub(G::GrpAbFinGen, s::Vector{GrpAbFinGenElem}) -> GrpAbFinGen, Map
+    sub(G::GrpAbFinGen, s::Vector{GrpAbFinGenElem}) -> GrpAbFinGen, GrpAbFinGenMap
 
 Create the subgroup $H$ of $G$ generated by the elements in `s` together
 with the injection $\iota : H \to G$.
@@ -1012,7 +1056,7 @@ function sub(G::GrpAbFinGen, s::Vector{GrpAbFinGenElem},
 end
 
 @doc Markdown.doc"""
-    sub(s::Vector{GrpAbFinGenElem}) -> GrpAbFinGen, Map
+    sub(s::Vector{GrpAbFinGenElem}) -> GrpAbFinGen, GrpAbFinGenMap
 
 Assuming that the non-empty array `s` contains elements of an abelian group
 $G$, this functions returns the subgroup $H$ of $G$ generated by the elements
@@ -1025,7 +1069,7 @@ function sub(s::Vector{GrpAbFinGenElem},
 end
 
 @doc Markdown.doc"""
-    sub(G::GrpAbFinGen, M::ZZMatrix) -> GrpAbFinGen, Map
+    sub(G::GrpAbFinGen, M::ZZMatrix) -> GrpAbFinGen, GrpAbFinGenMap
 
 Create the subgroup $H$ of $G$ generated by the elements corresponding to the
 rows of $M$ together with the injection $\iota : H \to G$.
@@ -1115,7 +1159,7 @@ function _sub_integer_snf(G::GrpAbFinGen, n::ZZRingElem, add_to_lattice::Bool =
 end
 
 @doc Markdown.doc"""
-    sub(G::GrpAbFinGen, n::ZZRingElem) -> GrpAbFinGen, Map
+    sub(G::GrpAbFinGen, n::ZZRingElem) -> GrpAbFinGen, GrpAbFinGenMap
 
 Create the subgroup $n \cdot G$ of $G$ together
 with the injection $\iota : n\cdot G \to G$.
@@ -1153,7 +1197,7 @@ end
 ################################################################################
 
 @doc Markdown.doc"""
-    quo(G::GrpAbFinGen, s::Vector{GrpAbFinGenElem}) -> GrpAbFinGen, Map
+    quo(G::GrpAbFinGen, s::Vector{GrpAbFinGenElem}) -> GrpAbFinGen, GrpAbfinGemMap
 
 Create the quotient $H$ of $G$ by the subgroup generated by the elements in
 $s$, together with the projection $p : G \to H$.
@@ -1202,7 +1246,7 @@ function quo(G::GrpAbFinGen, s::Vector{GrpAbFinGenElem},
 end
 
 @doc Markdown.doc"""
-    quo(G::GrpAbFinGen, M::ZZMatrix) -> GrpAbFinGen, Map
+    quo(G::GrpAbFinGen, M::ZZMatrix) -> GrpAbFinGen, GrpAbFinGenMap
 
 Create the quotient $H$ of $G$ by the subgroup generated by the elements
 corresponding to the rows of $M$, together with the projection $p : G \to H$.
@@ -1438,7 +1482,7 @@ end
 @doc Markdown.doc"""
     is_cyclic(G::GrpAbFinGen) -> Bool
 
-Returns whether $G$ is cyclic.
+Return whether $G$ is cyclic.
 """
 function is_cyclic(G::GrpAbFinGen)
   if !is_snf(G)
@@ -1468,9 +1512,9 @@ function _psylow_subgroup_gens(G::GrpAbFinGen, p::IntegerUnion)
 end
 
 @doc Markdown.doc"""
-    psylow_subgroup(G::GrpAbFinGen, p::IntegerUnion) -> GrpAbFinGen, Map
+    psylow_subgroup(G::GrpAbFinGen, p::IntegerUnion) -> GrpAbFinGen, GrpAbFinGenMap
 
-Returns the $p$-Sylow subgroup of `G`.
+Return the $p$-Sylow subgroup of `G`.
 """
 function psylow_subgroup(G::GrpAbFinGen, p::IntegerUnion,
                          to_lattice::Bool = true)
@@ -1482,9 +1526,9 @@ function psylow_subgroup(G::GrpAbFinGen, p::IntegerUnion,
 end
 
 @doc Markdown.doc"""
-    primary_part(G::GrpAbFinGen, m::IntegerUnion) -> GrpAbFinGen, Map
+    primary_part(G::GrpAbFinGen, m::IntegerUnion) -> GrpAbFinGen, GrpAbFinGenMap
 
-Returns the $m$-primary part of `G`.
+Return the $m$-primary part of `G`.
 """
 function primary_part(G::GrpAbFinGen, m::IntegerUnion,
                       to_lattice::Bool = true)
@@ -1857,7 +1901,7 @@ end
 @doc Markdown.doc"""
     elementary_divisors(G::GrpAbFinGen) -> Vector{ZZRingElem}
 
-Given $G$, returns the elementary divisors of $G$, that is, the unique positive
+Given $G$, return the elementary divisors of $G$, that is, the unique positive
 integers $e_1,\dotsc,e_k$ with $e_i \mid e_{i + 1}$ and
 $G \cong \mathbf{Z}/e_1\mathbf{Z} \times \dotsb \times \mathbf{Z}/e_k\mathbf{Z}$.
 """
@@ -1878,7 +1922,7 @@ end
 @doc Markdown.doc"""
     has_quotient(G::GrpAbFinGen, invariant::Vector{Int}) -> Bool
 
-Given an abelian group $G$, returns true if it has a quotient with given elementary
+Given an abelian group $G$, return true if it has a quotient with given elementary
 divisors and false otherwise.
 """
 function has_quotient(G::GrpAbFinGen, invariants::Vector{Int})
@@ -1907,7 +1951,7 @@ end
 @doc Markdown.doc"""
     is_pure(U::GrpAbFinGen, G::GrpAbFinGen) -> Bool
 
-A subgroup `U` of `G` is called pure iff for all `n` an element in `U` 
+A subgroup `U` of `G` is called pure if for all `n` an element in `U`
 that is in the image of the multiplication by `n` map of `G` is also
 a multiple of an element in `U`.
 
@@ -1954,7 +1998,7 @@ end
 @doc Markdown.doc"""
     is_neat(U::GrpAbFinGen, G::GrpAbFinGen) -> Bool
 
-A subgroup `U` of `G` is called neat iff for all primes `p` an element in `U` 
+A subgroup `U` of `G` is called neat if for all primes `p` an element in `U`
 that is in the image of the multiplication by `p` map of `G` is also
 a multiple of an element in `U`.
 
@@ -2022,10 +2066,10 @@ end
 #TODO: a better algorithm?
 @doc Markdown.doc"""
     has_complement(f::GrpAbFinGenMap) -> Bool, GrpAbFinGenMap
-    has_complement(U::GrpAbFinGen, G::GrpAbFinGen) -> Bool, GrpAbFinGen
+    has_complement(U::GrpAbFinGen, G::GrpAbFinGen) -> Bool, GrpAbFinGenMap
 
 Given a map representing a subgroup of a group $G$, 
-or a subgroup `U` of a group `G`, returns either true and
+or a subgroup `U` of a group `G`, return either true and
 an injection of a complement in $G$, or false.
 
 See also: [`is_pure`](@ref)
@@ -2090,6 +2134,7 @@ function has_complement(U::GrpAbFinGen, G::GrpAbFinGen)
   fl && return fl, image(mp)[1]
   return fl, U
 end
+
 ################################################################################
 #
 #  Identity
diff --git a/src/QuadForm/Herm/Genus.jl b/src/QuadForm/Herm/Genus.jl
index ae7644430c..b641d12c9d 100644
--- a/src/QuadForm/Herm/Genus.jl
+++ b/src/QuadForm/Herm/Genus.jl
@@ -1,7 +1,7 @@
 export genus, representative, rank, det, uniformizer, det_representative,
        gram_matrix, hermitian_genera, hermitian_local_genera, rank,
-       orthogonal_sum, is_inert, scales, ranks, dets, is_split, is_ramified,
-       is_dyadic, norms, primes, signatures
+       is_inert, scales, ranks, dets, is_split, is_ramified, is_dyadic,
+       norms, primes, signatures
 
 ################################################################################
 #
@@ -1062,21 +1062,21 @@ end
 ################################################################################
 
 @doc Markdown.doc"""
-    orthogonal_sum(g1::HermLocalGenus, g2::HermLocalGenus) -> HermLocalGenus
+    direct_sum(g1::HermLocalGenus, g2::HermLocalGenus) -> HermLocalGenus
 
 Given two local genus symbols `g1` and `g2` for hermitian lattices over $E/K$
-at the same prime ideal $\mathfrak p$ of $\mathcal O_K$, return their orthogonal
+at the same prime ideal $\mathfrak p$ of $\mathcal O_K$, return their direct
 sum. It corresponds to the local genus symbol of the $\mathfrak p$-adic completion
-of the orthogonal sum of respective representatives of `g1` and `g2`.
+of the direct sum of respective representatives of `g1` and `g2`.
 """
-function orthogonal_sum(G1::HermLocalGenus, G2::HermLocalGenus)
+function direct_sum(G1::HermLocalGenus, G2::HermLocalGenus)
   @req prime(G1) == prime(G2) "Local genera must have the same prime ideal"
   if !G1.is_dyadic || !G2.is_ramified
     return _direct_sum_easy(G1, G2)
   else
     L1 = representative(G1)
     L2 = representative(G2)
-    L3, = orthogonal_sum(L1, L2)
+    L3, = direct_sum(L1, L2)
     return genus(L3, prime(G1))
   end
 end
@@ -1313,13 +1313,13 @@ end
 ################################################################################
 
 @doc Markdown.doc"""
-    orthogonal_sum(G1::HermGenus, G2::HermGenus) -> HermGenus
+    direct_sum(G1::HermGenus, G2::HermGenus) -> HermGenus
 
 Given two global genus symbols `G1` and `G2` for hermitian lattices over $E/K$,
-return their orthogonal sum. It corresponds to the global genus symbol of the
-orthogonal sum of respective representatives of `G1` and `G2`.
+return their direct sum. It corresponds to the global genus symbol of the
+direct sum of respective representatives of `G1` and `G2`.
 """
-function orthogonal_sum(G1::HermGenus, G2::HermGenus)
+function direct_sum(G1::HermGenus, G2::HermGenus)
   @req G1.E === G2.E "Genera must have same base field"
   E = G1.E
   LGS = local_genus_herm_type(G1.E)[]
@@ -1329,7 +1329,7 @@ function orthogonal_sum(G1::HermGenus, G2::HermGenus)
   for p in union(P1, P2)
     g1 = G1[p]
     g2 = G2[p]
-    g3 = orthogonal_sum(g1, g2)
+    g3 = direct_sum(g1, g2)
     push!(prim, p)
     push!(LGS, g3)
   end
@@ -1339,7 +1339,7 @@ function orthogonal_sum(G1::HermGenus, G2::HermGenus)
   return genus(LGS, g3)
 end
 
-Base.:(+)(G1::HermGenus, G2::HermGenus) = orthogonal_sum(G1, G2)
+Base.:(+)(G1::HermGenus, G2::HermGenus) = direct_sum(G1, G2)
 
 ################################################################################
 #
diff --git a/src/QuadForm/Lattices.jl b/src/QuadForm/Lattices.jl
index eb29cf7a15..a4a6d91ae9 100644
--- a/src/QuadForm/Lattices.jl
+++ b/src/QuadForm/Lattices.jl
@@ -1486,26 +1486,89 @@ end
 
 ################################################################################
 #
-#  Orthogonal sum
-#
-################################################################################
-
-# TODO: Make this a proper coproduct with injections?
-function orthogonal_sum(M::T, N::T) where {T <: AbstractLat}
-  @req base_ring(M) === base_ring(N) "Base rings must be equal"
-  U = ambient_space(M)
-  V = ambient_space(N)
-  W, f1, f2 = orthogonal_sum(U, V)
-  rU = rank(U)
-  rV = rank(V)
-  rW = rank(W)
-  pM = pseudo_matrix(M)
-  pN = pseudo_matrix(N)
-  MpM = matrix(pM)
-  MpN = matrix(pN)
-  H = pseudo_matrix(diagonal_matrix(MpM, MpN),
-                    vcat(coefficient_ideals(pM), coefficient_ideals(pN)))
-  return lattice(W, H), f1, f2
+#  Direct sums/direct products/biproducts
+#
+################################################################################
+
+@doc Markdown.doc"""
+    direct_sum(x::Vararg{T}) where T <: AbstractLat -> T, Vector{AbstractSpaceMor}
+    direct_sum(x::Vector{T}) where T <: AbstractLat -> T, Vector{AbstractSpaceMor}
+
+Given a collection of quadratic or hermitian lattices $L_1, \ldots, L_n$,
+return their direct sum $L := L_1 \oplus \ldots \oplus L_n$, together with
+the injections $L_i \to L$ (seen as maps between the corresponding ambient spaces).
+
+For objects of type `AbstractLat`, finite direct sums and finite direct
+products agree and they are therefore called biproducts.
+If one wants to obtain `L` as a direct product with the projections $L \to L_i$,
+one should call `direct_product(x)`.
+If one wants to obtain `L` as a biproduct with the injections $L_i \to L$ and the
+projections $L \to L_i$, one should call `biproduct(x)`.
+"""
+function direct_sum(x::Vector{T}) where T <: AbstractLat
+  @req length(x) >= 2 "Input must consist of at least two lattices"
+  W, inj = direct_sum(ambient_space.(x))
+  H = _biproduct(x)
+  return lattice(W, H), inj
+end
+
+direct_sum(x::Vararg{AbstractLat}) = direct_sum(collect(x))
+
+@doc Markdown.doc"""
+    direct_product(x::Vararg{T}) where T <: AbstractLat -> T, Vector{AbstractSpaceMor}
+    direct_product(x::Vector{T}) where T <: AbstractLat -> T, Vector{AbstractSpaceMor}
+
+Given a collection of quadratic or hermitian lattices $L_1, \ldots, L_n$,
+return their direct product $L := L_1 \times \ldots \times L_n$, together with
+the projections $L \to L_i$ (seen as maps between the corresponding ambient spaces).
+
+For objects of type `AbstractLat`, finite direct sums and finite direct
+products agree and they are therefore called biproducts.
+If one wants to obtain `L` as a direct sum with the injections $L_i \to L$,
+one should call `direct_sum(x)`.
+If one wants to obtain `L` as a biproduct with the injections $L_i \to L$ and the
+projections $L \to L_i$, one should call `biproduct(x)`.
+"""
+function direct_product(x::Vector{T}) where T <: AbstractLat
+  @req length(x) >= 2 "Input must consist of at least two lattices"
+  W, proj = direct_product(ambient_space.(x))
+  H = _biproduct(x)
+  return lattice(W, H), proj
+end
+
+direct_product(x::Vararg{AbstractLat}) = direct_product(collect(x))
+
+@doc Markdown.doc"""
+    biproduct(x::Vararg{T}) where T <: AbstractLat -> T, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
+    biproduct(x::Vector{T}) where T <: AbstractLat -> T, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
+
+Given a collection of quadratic or hermitian lattices $L_1, \ldots, L_n$,
+return their biproduct $L := L_1 \oplus \ldots \oplus L_n$, together with
+the injections $L_i \toL$ and the projections $L \to L_i$ (seen as maps
+between the corresponding ambient spaces).
+
+For objects of type `AbstractLat`, finite direct sums and finite direct
+products agree and they are therefore called biproducts.
+If one wants to obtain `L` as a direct sum with the injections $L_i \to L$,
+one should call `direct_sum(x)`.
+If one wants to obtain `L` as a direct product with the projections $L \to L_i$,
+one should call `direct_product(x)`.
+"""
+function biproduct(x::Vector{T}) where T <: AbstractLat
+  @req length(x) >= 2 "Input must consist of at least two lattices"
+  W, inj, proj = biproduct(ambient_space.(x))
+  H = _biproduct(x)
+  return lattice(W, H), inj, proj
+end
+
+biproduct(x::Vararg{AbstractLat}) = biproduct(collect(x))
+
+function _biproduct(x::Vector{T}) where T <: AbstractLat
+  px = pseudo_matrix.(x)
+  Mpx = matrix.(px)
+  H = pseudo_matrix(diagonal_matrix(Mpx),
+                    reduce(vcat, coefficient_ideals.(px)))
+  return H
 end
 
 ################################################################################
diff --git a/src/QuadForm/Quad/Genus.jl b/src/QuadForm/Quad/Genus.jl
index 03c19cc4c7..d7d6c77fb6 100644
--- a/src/QuadForm/Quad/Genus.jl
+++ b/src/QuadForm/Quad/Genus.jl
@@ -324,7 +324,7 @@ function genus(J::JorDec)
   return g
 end
 
-function orthogonal_sum(J1::JorDec{S, T, U}, J2::JorDec{S, T, U}) where {S, T, U}
+function direct_sum(J1::JorDec{S, T, U}, J2::JorDec{S, T, U}) where {S, T, U}
   @req J1.p === J2.p "Jordan decompositions must be over same prime"
   if !(J1.is_dyadic)
     i1 = 1
@@ -367,7 +367,7 @@ function orthogonal_sum(J1::JorDec{S, T, U}, J2::JorDec{S, T, U}) where {S, T, U
     return JorDec(J1.p, _sca, _rk, _dets)
   else
     # Lazy
-    return JorDec(orthogonal_sum(lattice(J1), lattice(J2))[1], J1.p)
+    return JorDec(direct_sum(lattice(J1), lattice(J2))[1], J1.p)
   end
 end
 
@@ -569,7 +569,7 @@ function witt_invariant(G::QuadLocalGenus)
   w, d, n = G.witt[1], G.dets[1], G.ranks[1]
 
   for i in 2:length(G)
-    d, w, n = _witt_of_orthogonal_sum(d, w, n, G.dets[i], G.witt[i], G.ranks[i], p)
+    d, w, n = _witt_of_direct_sum(d, w, n, G.dets[i], G.witt[i], G.ranks[i], p)
   end
 
   G.witt_inv = w
@@ -1081,7 +1081,7 @@ end
 #
 ################################################################################
 
-function orthogonal_sum(G1::QuadLocalGenus, G2::QuadLocalGenus)
+function direct_sum(G1::QuadLocalGenus, G2::QuadLocalGenus)
   @req prime(G1) === prime(G2) "Local genera must have the same prime ideal"
   if !G1.is_dyadic
     p = prime(G1)
@@ -1108,12 +1108,12 @@ function orthogonal_sum(G1::QuadLocalGenus, G2::QuadLocalGenus)
   else
     L1 = representative(G1)
     L2 = representative(G2)
-    L3, = orthogonal_sum(L1, L2)
+    L3, = direct_sum(L1, L2)
     G3 = genus(L3, prime(G1))
   end
 
   if isdefined(G1, :jordec) && isdefined(G2, :jordec)
-    G3.jordec = orthogonal_sum(G1.jordec, G2.jordec)
+    G3.jordec = direct_sum(G1.jordec, G2.jordec)
   end
 
   return G3
@@ -1161,7 +1161,7 @@ function _direct_sum_easy(G1::QuadLocalGenus, G2::QuadLocalGenus, detclassesG2 =
   return genus(QuadLat, prime(G1), uniformizer(G1), _rk, _sca, _detclass)
 end
 
-Base.:(+)(G1::QuadLocalGenus, G2::QuadLocalGenus) = orthogonal_sum(G1, G2)
+Base.:(+)(G1::QuadLocalGenus, G2::QuadLocalGenus) = direct_sum(G1, G2)
 
 ################################################################################
 #
@@ -1915,11 +1915,11 @@ end
 
 ################################################################################
 #
-#  Orthogonal sum of genus symbols
+#  Direct sum of genus symbols
 #
 ################################################################################
 
-function orthogonal_sum(G1::QuadGenus{S, T, U}, G2::QuadGenus{S, T, U}) where {S, T, U}
+function direct_sum(G1::QuadGenus{S, T, U}, G2::QuadGenus{S, T, U}) where {S, T, U}
   @req G1.K === G2.K "Global genus symbols must be defined over the same field"
   K = G1.K
   LGS = local_genus_quad_type(K)[]
@@ -1945,7 +1945,7 @@ function orthogonal_sum(G1::QuadGenus{S, T, U}, G2::QuadGenus{S, T, U}) where {S
       dcl = fl ? 1 : -1
       g2 = genus(QuadLat, p, [(0, rank(G2), dcl)])
     end
-    g3 = orthogonal_sum(g1, g2)
+    g3 = direct_sum(g1, g2)
     push!(LGS, g3)
   end
   sig1 = G1.signatures
@@ -1955,7 +1955,7 @@ function orthogonal_sum(G1::QuadGenus{S, T, U}, G2::QuadGenus{S, T, U}) where {S
   return QuadGenus(K, G1.d * G2.d, LGS, sig3)
 end
 
-Base.:(+)(G1::QuadGenus, G2::QuadGenus) = orthogonal_sum(G1, G2)
+Base.:(+)(G1::QuadGenus, G2::QuadGenus) = direct_sum(G1, G2)
 
 ################################################################################
 #
diff --git a/src/QuadForm/Quad/Spaces.jl b/src/QuadForm/Quad/Spaces.jl
index 89a768265e..030cd829e6 100644
--- a/src/QuadForm/Quad/Spaces.jl
+++ b/src/QuadForm/Quad/Spaces.jl
@@ -243,7 +243,7 @@ end
 # wi = witt invariant
 # ni = rank
 # Lam p. 117
-function _witt_of_orthogonal_sum(d1, w1, n1, d2, w2, n2, p)
+function _witt_of_direct_sum(d1, w1, n1, d2, w2, n2, p)
   _n1 = mod(n1, 4)
   if _n1 == 0 || _n1 == 1
     disc1 = d1
@@ -1192,7 +1192,7 @@ function is_isotropic_with_vector(q::QuadSpace{QQField, QQMatrix})
     return false, z
   end
   # confirm the computation
-  v = [S[1,i] for i in 1:ncols(S)]
+  v = elem_type(base_ring(q))[S[1,i] for i in 1:ncols(S)]
   @hassert :Lattice 1 inner_product(q,v,v)==0
   @hassert :Lattice 1 !all(x==0 for x in v)
   return true, v
@@ -1265,7 +1265,7 @@ function _isotropic_subspace(q::QuadSpace{QQField, QQMatrix})
     s = (0, a)
   end
   R = representative(genus(D, s))
-  LL, iM, iR = orthogonal_sum(M, R)
+  LL, inj = direct_sum(M, R)
   MM = maximal_even_lattice(LL)
   # MM is sum of hyperbolic planes -> Simon should succeed
   ok, H = _maximal_isotropic_subspace_unimodular(MM)
@@ -1274,10 +1274,11 @@ function _isotropic_subspace(q::QuadSpace{QQField, QQMatrix})
   # the  totally isotropic subspace H has large enough dimension so that its
   # intersection with L is non-trivial (and isotropic) -> we win
   VV = ambient_space(MM)
-  iso = preimage(iM, lattice(VV, H))
+  iso = preimage(inj[1], lattice(VV, H))
   @hassert :Lattice 0 rank(iso) >0
   @hassert :Lattice 0 gram_matrix(iso)==0
-  return true, basis_matrix(iso)
+  B = basis_matrix(iso)::dense_matrix_type(base_field(M))
+  return true, B
 end
 
 function _maximal_isotropic_subspace_unimodular(L)
@@ -2159,13 +2160,13 @@ function Base.:(+)(G1::LocalQuadSpaceCls, G2::LocalQuadSpaceCls)
   r2 = dim_radical(G2)
   r = r1 + r2
   d = det_nondegenerate_part(G1)*det_nondegenerate_part(G2)
-  _,w,_ = _witt_of_orthogonal_sum(G1.det, witt_invariant(G1), dim(G1)-r1,
+  _,w,_ = _witt_of_direct_sum(G1.det, witt_invariant(G1), dim(G1)-r1,
                                  G2.det, witt_invariant(G2), dim(G2)-r2, p)
   h = _witt_hasse(w, n - r, d, p)
   return local_quad_space_class(K, p, n, d, h, r)
 end
 
-orthogonal_sum(G1::LocalQuadSpaceCls, G2::LocalQuadSpaceCls) = G1 + G2
+direct_sum(G1::LocalQuadSpaceCls, G2::LocalQuadSpaceCls) = G1 + G2
 
 @doc Markdown.doc"""
     Base.:(-)(G1::LocalQuadSpaceCls, G2::LocalQuadSpaceCls)
@@ -2469,11 +2470,11 @@ end
 
 # Direct sum
 @doc Markdown.doc"""
-    orthogonal_sum(g1::QuadSpaceCls, g2::QuadSpaceCls) -> QuadSpaceCls
+    direct_sum(g1::QuadSpaceCls, g2::QuadSpaceCls) -> QuadSpaceCls
 
 Return the isometry class of the direct sum of two representatives.
 """
-function orthogonal_sum(g1::QuadSpaceCls{S,T,U},g2::QuadSpaceCls{S,T,U}) where {S,T,U}
+function direct_sum(g1::QuadSpaceCls{S,T,U},g2::QuadSpaceCls{S,T,U}) where {S,T,U}
   @req base_ring(g1) == base_ring(g2) "must be defined over the same base ring"
   K = base_ring(g1)
   g = class_quad_type(K)(K)
@@ -2503,7 +2504,7 @@ end
 Return the isometry class of the direct sum of two representatives.
 """
 function Base.:(+)(g1::QuadSpaceCls{S,T,U},g2::QuadSpaceCls{S,T,U}) where {S,T,U}
-  return orthogonal_sum(g1, g2)
+  return direct_sum(g1, g2)
 end
 
 function Base.:(-)(g1::QuadSpaceCls{S,T,U},g2::QuadSpaceCls{S,T,U}) where {S,T,U}
diff --git a/src/QuadForm/Quad/ZGenus.jl b/src/QuadForm/Quad/ZGenus.jl
index 95cd7f0100..238562b03e 100644
--- a/src/QuadForm/Quad/ZGenus.jl
+++ b/src/QuadForm/Quad/ZGenus.jl
@@ -1,6 +1,6 @@
 export genus, rank, det, dim, prime, symbol, representative, signature,
-       oddity, excess, level, Zgenera, scale, norm, mass, orthogonal_sum,
-       quadratic_space,hasse_invariant, local_symbol, local_symbols,
+       oddity, excess, level, Zgenera, scale, norm, mass,
+       quadratic_space, hasse_invariant, local_symbol, local_symbols,
        ZGenus, ZpGenus, representatives, is_elementary, is_primary, is_unimodular,
        is_primary_with_prime, is_elementary_with_prime, automorphous_numbers,
        is_automorphous, bad_primes, signature_pair, signature_tuple
@@ -463,11 +463,11 @@ function genus(A::MatElem, p)
 end
 
 @doc Markdown.doc"""
-    orthogonal_sum(S1::ZpGenus, S2::ZpGenus) -> ZpGenus
+    direct_sum(S1::ZpGenus, S2::ZpGenus) -> ZpGenus
 
-Return the local genus of the orthogonal direct sum of two representatives.
+Return the local genus of the direct sum of two representatives.
 """
-function orthogonal_sum(S1::ZpGenus, S2::ZpGenus)
+function direct_sum(S1::ZpGenus, S2::ZpGenus)
   @req prime(S1) == prime(S2) "The local genus symbols must be over the same prime"
   if rank(S1) == rank(S2) == 0
     return ZpGenus(prime(S1), Hecke.symbol(S1))
@@ -506,16 +506,14 @@ function orthogonal_sum(S1::ZpGenus, S2::ZpGenus)
   return ZpGenus(prime(S1), symbol)
 end
 
-direct_sum(S1::ZpGenus, S2::ZpGenus) = orthogonal_sum(S1, S2)
-
 @doc Markdown.doc"""
-    orthogonal_sum(G1::ZGenus, G2::ZGenus) -> ZGenus
+    direct_sum(G1::ZGenus, G2::ZGenus) -> ZGenus
 
-Return the genus of the orthogonal direct sum of `G1` and `G2`.
+Return the genus of the direct sum of `G1` and `G2`.
 
-The orthogonal direct sum is defined via representatives.
+The direct sum is defined via representatives.
 """
-function orthogonal_sum(G1::ZGenus, G2::ZGenus)
+function direct_sum(G1::ZGenus, G2::ZGenus)
   p1, n1 = signature_pair(G1)
   p2, n2 = signature_pair(G2)
   sign_pair = (p1 + p2, n1 + n2)
@@ -524,15 +522,12 @@ function orthogonal_sum(G1::ZGenus, G2::ZGenus)
   sort(primes)
   local_symbols = []
   for p in primes
-    sym_p = orthogonal_sum(local_symbol(G1, p), local_symbol(G2, p))
+    sym_p = direct_sum(local_symbol(G1, p), local_symbol(G2, p))
     push!(local_symbols, sym_p)
   end
   return ZGenus(sign_pair, local_symbols)
 end
 
-direct_sum(S1::ZGenus, S2::ZGenus) = orthogonal_sum(S1, S2)
-
-
 ###############################################################################
 #
 # Enumeration of genus symbols
@@ -2668,7 +2663,7 @@ Return a (primitive) embedding of the integral lattice `S` into some
 For now this works only for even lattices.
 
 ```jldoctest
-julia> NS = orthogonal_sum(Zlattice(:U), rescale(root_lattice(:A, 16), -1))[1];
+julia> NS = direct_sum(Zlattice(:U), rescale(root_lattice(:A, 16), -1))[1];
 
 julia> LK3, iNS, i = embed_in_unimodular(NS, 3, 19);
 
@@ -2698,7 +2693,8 @@ function embed_in_unimodular(S::ZLat, pos, neg; primitive=true, even=true)
   R = lll(R)  # make R a bit nicer
   R = Zlattice(gram=gram_matrix(R)) # clear the history of R
 
-  SR, iS, iR = orthogonal_sum(S, R)
+  SR, inj = direct_sum(S, R)
+  iS, iR = inj
   V = ambient_space(SR)
   S = lattice(V, basis_matrix(S) * iS.matrix)
   R = lattice(V, basis_matrix(R) * iR.matrix)
diff --git a/src/QuadForm/Quad/ZLattices.jl b/src/QuadForm/Quad/ZLattices.jl
index 4eeed5f3b3..b418fe34c4 100644
--- a/src/QuadForm/Quad/ZLattices.jl
+++ b/src/QuadForm/Quad/ZLattices.jl
@@ -196,52 +196,90 @@ end
 
 ################################################################################
 #
-#  Orthogonal sum
+#  Direct sums
 #
 ################################################################################
 
+function _biproduct(x::Vector{ZLat})
+  Bs = basis_matrix.(x)
+  B = diagonal_matrix(Bs)
+  return B
+end
+
 @doc Markdown.doc"""
-    orthogonal_sum(L1::ZLat, L2::ZLat)
+    direct_sum(x::Vararg{ZLat}) -> ZLat, Vector{AbstractSpaceMor}
+    direct_sum(x::Vector{ZLat}) -> ZLat, Vector{AbstractSpaceMor}
 
-Return the orthogonal direct sum of the lattices `L1` and `L2`.
+Given a collection of $\mathbb Z$-lattices $L_1, \ldots, L_n$,
+return their direct sum $L := L_1 \oplus \ldots \oplus L_n$,
+together with the injections $L_i \to L$.
+(seen as maps between the corresponding ambient spaces).
 
-It lives in the orthogonal direct sum of their ambient spaces.
+For objects of type `ZLat`, finite direct sums and finite direct products
+agree and they are therefore called biproducts.
+If one wants to obtain `L` as a direct product with the projections $L \to L_i$,
+one should call `direct_product(x)`.
+If one wants to obtain `L` as a biproduct with the injections $L_i \to L$ and
+the projections $L \to L_i$, one should call `biproduct(x)`.
 """
-function orthogonal_sum(L1::ZLat, L2::ZLat)
-  V1 = ambient_space(L1)
-  V2 = ambient_space(L2)
-  B1 = basis_matrix(L1)
-  B2 = basis_matrix(L2)
-  B = diagonal_matrix(B1, B2)
-  V, f1, f2 = orthogonal_sum(V1, V2)
-  return lattice(V, B), f1, f2
+function direct_sum(x::Vector{ZLat})
+  @req length(x) >= 2 "Input must consist of at least two lattices"
+  W, inj = direct_sum(ambient_space.(x))
+  B = _biproduct(x)
+  return lattice(W, B), inj
 end
 
-function _orthogonal_sum_with_injections_and_projections(x::Vector{ZLat})
-  @req length(x) >= 2 "Input must contain at least two lattices"
-  y = ambient_space.(x)
-  Bs = basis_matrix.(x)
-  B = diagonal_matrix(Bs)
-  V, inj, proj = Hecke._orthogonal_sum_with_injections_and_projections(y)
-  return lattice(V, B), inj, proj
+direct_sum(x::Vararg{ZLat}) = direct_sum(collect(x))
+
+@doc Markdown.doc"""
+    direct_product(x::Vararg{ZLat}) -> ZLat, Vector{AbstractSpaceMor}
+    direct_product(x::Vector{ZLat}) -> ZLat, Vector{AbstractSpaceMor}
+
+Given a collection of $\mathbb Z$-lattices $L_1, \ldots, L_n$,
+return their direct product $L := L_1 \times \ldots \times L_n$,
+together with the projections $L \to L_i$.
+(seen as maps between the corresponding ambient spaces).
+
+For objects of type `ZLat`, finite direct sums and finite direct products
+agree and they are therefore called biproducts.
+If one wants to obtain `L` as a direct sum with the injections $L_i \to L$,
+one should call `direct_sum(x)`.
+If one wants to obtain `L` as a biproduct with the injections $L_i \to L$ and
+the projections $L \to L_i$, one should call `biproduct(x)`.
+"""
+function direct_product(x::Vector{ZLat})
+  @req length(x) >= 2 "Input must consist of at least two lattices"
+  W, proj = direct_product(ambient_space.(x))
+  B = _biproduct(x)
+  return lattice(W, B), proj
 end
 
+direct_product(x::Vararg{ZLat}) = direct_product(collect(x))
+
 @doc Markdown.doc"""
-    direct_sum(x::Vararg{ZLat}) -> ZLat, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
-    direct_sum(x::Vector{ZLat}) -> ZLat, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
+    biproduct(x::Vararg{ZLat}) -> ZLat, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
+    biproduct(x::Vector{ZLat}) -> ZLat, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
 
 Given a collection of $\mathbb Z$-lattices $L_1, \ldots, L_n$,
-return their complete direct sum $L := L_1 \oplus \ldots \oplus L_n$,
-together with the injections $L_i \to L$ and the projections $L \to L_i$
+return their biproduct $L := L_1 \oplus \ldots \oplus L_n$,
+together with the injections $L_i \to L$ and the projections $L \to L_i$.
 (seen as maps between the corresponding ambient spaces).
+
+For objects of type `ZLat`, finite direct sums and finite direct products
+agree and they are therefore called biproducts.
+If one wants to obtain `L` as a direct sum with the injections $L_i \to L$,
+one should call `direct_sum(x)`.
+If one wants to obtain `L` as a direct product with the projections $L \to L_i$,
+one should call `direct_product(x)`.
 """
-function direct_sum(x::Vararg{ZLat})
-  x = collect(x)
+function biproduct(x::Vector{ZLat})
   @req length(x) >= 2 "Input must consist of at least two lattices"
-  return _orthogonal_sum_with_injections_and_projections(x)
+  W, inj, proj = biproduct(ambient_space.(x))
+  B = _biproduct(x)
+  return lattice(W, B), inj, proj
 end
 
-direct_sum(x::Vector{ZLat}) = _orthogonal_sum_with_injections_and_projections(x)
+biproduct(x::Vararg{ZLat}) = biproduct(collect(x))
 
 @doc Markdown.doc"""
     orthogonal_submodule(L::ZLat, S::ZLat) -> ZLat
@@ -261,6 +299,13 @@ function orthogonal_submodule(L::ZLat, S::ZLat)
   return lattice(V, Ks*B)
 end
 
+### Deprecated: should be removed once fixed in Oscar!!
+
+function orthogonal_sum(x::ZLat, y::ZLat)
+  z, inj = direct_sum(x, y)
+  return z, inj[1], inj[2]
+end
+
 ################################################################################
 #
 #  String I/O
diff --git a/src/QuadForm/Spaces.jl b/src/QuadForm/Spaces.jl
index 88c3938542..24b4423bba 100644
--- a/src/QuadForm/Spaces.jl
+++ b/src/QuadForm/Spaces.jl
@@ -673,58 +673,17 @@ end
 
 ################################################################################
 #
-#  Orthogonal sum
+#  Direct sums
 #
 ################################################################################
 
-function _orthogonal_sum(V::AbstractSpace, W::AbstractSpace)
-  K = base_ring(V)
-  G = diagonal_matrix(gram_matrix(V), gram_matrix(W))
-  n = dim(V) + dim(W)
-  i1 = zero_matrix(K, dim(V), n)
-  for i in 1:dim(V)
-    i1[i, i] = 1
-  end
-  i2 = zero_matrix(K, dim(W), n)
-  for i in 1:dim(W)
-    i2[i, i + dim(V)] = 1
-  end
-  return G, i1, i2
-end
-
-@doc Markdown.doc"""
-    orthogonal_sum(V::AbstractSpace, W::AbstractSpace) -> AbstractSpace, AbstractSpaceMor, AbstractSpaceMor
-
-Given two spaces `V` and `W` of the same kind (either both hermitian or both quadratic)
-and defined over the same algebra, return their orthogonal sum $V \oplus W$. It is given with
-the two natural embeddings $V \to V\oplus W$ and $W \to V\oplus W$.
-"""
-orthogonal_sum(V::AbstractSpace, W::AbstractSpace)
-
-function orthogonal_sum(V::QuadSpace, W::QuadSpace)
-  @req base_ring(V) === base_ring(W) "Base algebra must be equal"
-  G, i1, i2 = _orthogonal_sum(V, W)
-  VplusW = quadratic_space(base_ring(V), G)
-  f1 = hom(V, VplusW, i1)
-  f2 = hom(W, VplusW, i2)
-  return VplusW, f1, f2
-end
-
-function orthogonal_sum(V::HermSpace, W::HermSpace)
-  @req base_ring(V) === base_ring(W) "Base algebra must be equal"
-  G, i1, i2 = _orthogonal_sum(V, W)
-  VplusW = hermitian_space(base_ring(V), G)
-  f1 = hom(V, VplusW, i1)
-  f2 = hom(W, VplusW, i2)
-  return VplusW, f1, f2
-end
-
-function _orthogonal_sum_with_injections_and_projections(x::Vector{<:QuadSpace})
+function _biproduct(x::Vector{T}) where T <: AbstractSpace
   @req length(x) >= 2 "Input must contain at least two quadratic spaces"
   K = base_ring(x[1])
   @req all(i -> base_ring(x[i]) === K, 2:length(x)) "All spaces must be defined over the same field"
+  @req is_quadratic(x[1]) ? all(i -> is_quadratic(x[i]), 2:length(x)) : all(i -> ishermitian(x[i]), 1:length(x)) "Spaces must be all hermitian or all quadratic" 
   G = diagonal_matrix(gram_matrix.(x))
-  V = quadratic_space(K, G)
+  V = is_quadratic(x[1]) ? quadratic_space(K, G) : hermitian_space(K, G)
   n = sum(dim.(x))
   inj = AbstractSpaceMor[]
   proj = AbstractSpaceMor[]
@@ -747,20 +706,72 @@ function _orthogonal_sum_with_injections_and_projections(x::Vector{<:QuadSpace})
 end
 
 @doc Markdown.doc"""
-    direct_sum(x::Vararg{QuadSpace}) -> QuadSpace, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
-    direct_sum(x::Vector{QuadSpace}) -> QuadSpace, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
+    direct_sum(x::Vararg{T}) where T <: AbstractSpace -> T, Vector{AbstractSpaceMor}
+    direct_sum(x::Vector{T}) where T <: AbstractSpace -> T, Vector{AbstractSpaceMor}
+
+Given a collection of quadratic or hermitian spaces $V_1, \ldots, V_n$,
+return their direct sum $V := V_1 \oplus \ldots \oplus V_n$,
+together with the injections $V_i \to V$.
+
+For objects of type `AbstractSpace`, finite direct sums and finite direct
+products agree and they are therefore called biproducts.
+If one wants to obtain `V` as a direct product with the projections $V \to V_i$,
+one should call `direct_product(x)`.
+If one wants to obtain `V` as a biproduct with the injections $V_i \to V$ and
+the projections $V \to V_i$, one should call `biproduct(x)`.
+"""
+function direct_sum(x::Vector{T}) where T <: AbstractSpace
+  @req length(x) >= 2 "Input must consist of at least two spaces"
+  V, inj, = _biproduct(x)
+  return V, inj
+end
+
+direct_sum(x::Vararg{AbstractSpace}) = direct_sum(collect(x))
+
+@doc Markdown.doc"""
+    direct_product(x::Vararg{T}) where T <: AbstractSpace -> T, Vector{AbstractSpaceMor}
+    direct_product(x::Vector{T}) where T <: AbstractSpace -> T, Vector{AbstractSpaceMor}
+
+Given a collection of quadratic or hermitian spaces $V_1, \ldots, V_n$,
+return their direct product $V := V_1 \times \ldots \times V_n$,
+together with the projections $V \to V_i$.
+
+For objects of type `AbstractSpace`, finite direct sums and finite direct
+products agree and they are therefore called biproducts.
+If one wants to obtain `V` as a direct sum with the injections $V_i \to V$,
+one should call `direct_sum(x)`.
+If one wants to obtain `V` as a biproduct with the injections $V_i \to V$ and
+the projections $V \to V_i$, one should call `biproduct(x)`.
+"""
+function direct_product(x::Vector{T}) where T <: AbstractSpace
+  @req length(x) >= 2 "Input must consist of at least two spaces"
+  V, _, proj = _biproduct(x)
+  return V, proj
+end
+
+direct_product(x::Vararg{AbstractSpace}) = direct_product(collect(x))
+
+@doc Markdown.doc"""
+    biproduct(x::Vararg{T}) where T <: AbstractSpace -> T, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
+    biproduct(x::Vector{T}) where T <: AbstractSpace -> T, Vector{AbstractSpaceMor}, Vector{AbstractSpaceMor}
+
+Given a collection of quadratic or hermitian spaces $V_1, \ldots, V_n$,
+return their biproduct $V := V_1 \oplus \ldots \oplus V_n$, together
+with the injections $V_i \to V$ and the projections $V \to V_i$.
 
-Given a collection of quadratic spaces $V_1, \ldots, V_n$,
-return their complete direct sum $V := V_1 \oplus \ldots \oplus V_n$,
-together with the injections $V_i \to V$ and the projections $V \to V_i$.
+For objects of type `AbstractSpace`, finite direct sums and finite direct
+products agree and they are therefore called biproducts.
+If one wants to obtain `V` as a direct sum with the injections $V_i \to V$,
+one should call `direct_sum(x)`.
+If one wants to obtain `V` as a direct product with the projections $V \to V_i$,
+one should call `direct_product(x)`.
 """
-function direct_sum(x::Vararg{QuadSpace})
-  x = collect(x)
-  @req length(x) >= 2 "Input must consist of at least two quadratic spaces"
-  return _orthogonal_sum_with_injections_and_projections(x)
+function biproduct(x::Vector{T}) where T <: AbstractSpace
+  @req length(x) >= 2 "Input must consisy of at least two spaces"
+  return _biproduct(x)
 end
 
-direct_sum(x::Vector{<:QuadSpace}) = _orthogonal_sum_with_injections_and_projections(x)
+biproduct(x::Vararg{AbstractSpace}) = biproduct(collect(x))
 
 ################################################################################
 #
diff --git a/src/QuadForm/Torsion.jl b/src/QuadForm/Torsion.jl
index a29e720bd2..05ca9be074 100644
--- a/src/QuadForm/Torsion.jl
+++ b/src/QuadForm/Torsion.jl
@@ -1756,87 +1756,144 @@ end
 
 ###############################################################################
 #
-#  Orthogonal sum
+#  Sums
 #
 ###############################################################################
 
 @doc Markdown.doc"""
-    orthogonal_sum(T::TorQuadModule, U::TorQuadModule)
-                                    -> TorQuadModule, TorQuadModuleMor, TorQuadModuleMor
+    +(T::TorQuadModule, U::TorQuadModule) -> TorQuadModule
 
-Return the orthogonal direct sum `S` of `T` and `U` as a quotient of the direct
-sum of their respective covers.
-
-It is given with the two injections $T \to S$ and $U \to S$.
+Given two torsion quadratic modules `T` and `U` whose covers are in the same
+ambient space, return their sum `S` defined as the quotient of the sum of their
+covers by the sum of their respective relation lattices.
 
 Note that `T` and `U` must have the same moduli, both bilinear and quadratic
 ones.
 """
-function orthogonal_sum(T::TorQuadModule, U::TorQuadModule)
+function +(T::TorQuadModule, U::TorQuadModule)
   @req modulus_bilinear_form(T) == modulus_bilinear_form(U) "T and U must have the same bilinear modulus"
   @req modulus_quadratic_form(T) == modulus_quadratic_form(U) "T and U must have the same quadratic modulus"
-  cT = cover(T)
-  rT = relations(T)
-  cU = cover(U)
-  rU = relations(U)
-  cS, cTincS, cUincS = orthogonal_sum(cT, cU)
-  rS = lattice(ambient_space(cS), block_diagonal_matrix([basis_matrix(rT), basis_matrix(rU)]))
-  geneT = [cTincS(lift(a)) for a in gens(T)]
-  geneU = [cUincS(lift(b)) for b in gens(U)]
-  S = torsion_quadratic_module(cS, rS, gens = vcat(geneT, geneU), modulus = modulus_bilinear_form(T), modulus_qf = modulus_quadratic_form(T))
-  TinS = hom(T, S, S.(geneT))
-  UinS = hom(U, S, S.(geneU))
-  return S, TinS, UinS
-end
-
-function _orthogonal_sum_with_injections_and_projections(x::Vector{TorQuadModule})
-  @req length(x) >= 2 "Input must contain at least two torsion quadratic modules"
+  @req ambient_space(cover(T)) === ambient_space(cover(U)) "Covers must be in the same ambient space"
+  cS = cover(T) + cover(U)
+  rS = relations(T) + relations(U)
+  geneT = [lift(a) for a in gens(T)]
+  geneU = [lift(b) for b in gens(U)]
+  S = torsion_quadratic_module(cS, rS, gens = unique([g for g in vcat(geneT, geneU) if !is_zero(g)]), modulus = modulus_bilinear_form(T), modulus_qf = modulus_quadratic_form(T))
+  return S
+end
+
+function _biproduct(x::Vector{TorQuadModule}; proj = true)
+  @req length(x) >= 2 "Input must consist of at least two torsion quadratic modules"
   mbf = modulus_bilinear_form(x[1])
   mqf = modulus_quadratic_form(x[1])
-  @req all(i -> modulus_bilinear_form(x[i]) == mbf, 2:length(x)) "All torsion quadratic modules must have same bilinear modulus"
-  @req all(i -> modulus_quadratic_form(x[i]) == mqf, 2:length(x)) "All torsion quadratic modules must have same quadratic modulus"
+  @req all(i -> modulus_bilinear_form(x[i]) == mbf, 2:length(x)) "All torsion quadratic modules must have the same bilinear modulus"
+  @req all(i -> modulus_quadratic_form(x[i]) == mqf, 2:length(x)) "All torsion quadratic modules must have the same quadratic modulus"
   cs = cover.(x)
   rs = relations.(x)
-  C, inj, proj = _orthogonal_sum_with_injections_and_projections(cs)
+  C, injC, projC = biproduct(cs)
   R = lattice(ambient_space(C), block_diagonal_matrix(basis_matrix.(rs)))
   gensinj = Vector{Vector{QQFieldElem}}[]
   gensproj = Vector{Vector{QQFieldElem}}[]
   inj2 = TorQuadModuleMor[]
   proj2 = TorQuadModuleMor[]
   for i in 1:length(x)
-    gene = [inj[i](lift(a)) for a in gens(x[i])]
+    gene = [injC[i](lift(a)) for a in gens(x[i])]
     push!(gensinj, gene)
   end
   S = torsion_quadratic_module(C, R, gens = reduce(vcat, gensinj), modulus = mbf, modulus_qf = mqf)
-  for i in 1:length(x)
-    gene = [proj[i](lift(a)) for a in gens(S)]
-    push!(gensproj, gene)
-  end
   for i in 1:length(x)
     T = x[i]
     iT = hom(T, S, S.(gensinj[i]))
-    pT = hom(S, T, T.(gensproj[i]))
     push!(inj2, iT)
-    push!(proj2, pT)
+  end
+  if proj
+    for i in 1:length(x)
+      gene = [projC[i](lift(a)) for a in gens(S)]
+      push!(gensproj, gene)
+    end
+    for i in 1:length(x)
+      T = x[i]
+      pT = hom(S, T, T.(gensproj[i]))
+      push!(proj2, pT)
+    end
   end
   return S, inj2, proj2
 end
 
 @doc Markdown.doc"""
-    direct_sum(x::Vararg{TorQuadModule}) -> TorQuadModule, Vector{TorQuadModuleMor}, Vector{TorQuadModuleMor}
-    direct_sum(x::Vector{TorQuadModule}) -> TorQuadModule, Vector{TorQuadModuleMor}, Vector{TorQuadModuleMor}
+    direct_sum(x::Vararg{TorQuadModule}) -> TorQuadModule, Vector{TorQuadModuleMor}
+    direct_sum(x::Vector{TorQuadModule}) -> TorQuadModule, Vector{TorQuadModuleMor}
+
+Given a collection of torsion quadratic modules $T_1, \ldots, T_n$, return
+their direct sum $T := T_1\oplus \ldots \oplus T_n$, together with the
+injections $T_i \to T$.
+
+For objects of type `TorQuadModule`, finite direct sums and finite direct products
+agree and they are therefore called biproducts.
+If one wants to obtain `T` as a direct product with the projections $T \to T_i$,
+one should call `direct_product(x)`.
+If one wants to obtain `T` as a biproduct with the injections $T_i \to T$ and the
+projections $T \to T_i$, one should call `biproduct(x)`.
+"""
+function direct_sum(x::Vector{TorQuadModule})
+  T, inj, = _biproduct(x, proj=false)
+  return T, inj
+end
+
+direct_sum(x::Vararg{TorQuadModule}) = direct_sum(collect(x))
 
-Given a collection of torsion quadratic modules $T_1, \ldots, T_n$,
-return their complete direct sum $T := T_1 \oplus \ldots \oplus T_n$,
-together with the injections $T_i \to T$ and the projections $T \to T_i$.
+@doc Markdown.doc"""
+    direct_product(x::Vararg{TorQuadModule}) -> TorQuadModule, Vector{TorQuadModuleMor}
+    direct_product(x::Vector{TorQuadModule}) -> TorQuadModule, Vector{TorQuadModuleMor}
+
+Given a collection of torsion quadratic modules $T_1, \ldots, T_n$, return
+their direct product $T := T_1\times \ldots \times T_n$, together with the
+projections $T \to T_i$.
+
+For objects of type `TorQuadModule`, finite direct sums and finite direct products
+agree and they are therefore called biproducts.
+If one wants to obtain `T` as a direct sum with the inctions $T_i \to T$,
+one should call `direct_sum(x)`.
+If one wants to obtain `T` as a biproduct with the injections $T_i \to T$ and the
+projections $T \to T_i$, one should call `biproduct(x)`.
+"""
+function direct_product(x::Vector{TorQuadModule})
+  T, _, proj = _biproduct(x)
+  return T, proj
+end
+
+direct_product(x::Vararg{TorQuadModule}) = direct_product(collect(x))
+
+@doc Markdown.doc"""
+    biproduct(x::Vararg{TorQuadModule}) -> TorQuadModule, Vector{TorQuadModuleMor}, Vector{TorQuadModuleMor}
+    biproduct(x::Vector{TorQuadModule}) -> TorQuadModule, Vector{TorQuadModuleMor}, Vector{TorQuadModuleMor}
+
+Given a collection of torsion quadratic modules $T_1, \ldots, T_n$, return
+their biproduct $T := T_1\oplus \ldots \oplus T_n$, together with the
+injections $T_i \to T$ and the projections $T \to T_i$.
+
+For objects of type `TorQuadModule`, finite direct sums and finite direct products
+agree and they are therefore called biproducts.
+If one wants to obtain `T` as a direct sum with the inctions $T_i \to T$,
+one should call `direct_sum(x)`.
+If one wants to obtain `T` as a direct product with the projections $T \to T_i$,
+one should call `direct_product(x)`.
 """
-function direct_sum(x::Vararg{TorQuadModule})
-  x = collect(x)
-  @req length(x) >= 2 "Input must consist of at least two torsion quadratic module"
-  return _orthogonal_sum_with_injections_and_projections(x)
+function biproduct(x::Vector{TorQuadModule})
+  return _biproduct(x)
+end
+
+biproduct(x::Vararg{TorQuadModule}) = biproduct(collect(x))
+
+### Deprecated: should be removed once fixed in Oscar!!
+
+function orthogonal_sum(x::TorQuadModule, y::TorQuadModule)
+  z, inj = direct_sum(x, y)
+  return z, inj[1], inj[2]
 end
+export orthogonal_sum
 
-direct_sum(x::Vector{TorQuadModule}) = _orthogonal_sum_with_injections_and_projections(x)
+_orthogonal_sum_with_injections_and_projections(x::Vector{TorQuadModule}) = biproduct(x)
 
 ###############################################################################
 #
diff --git a/test/GrpAb/GrpAbFinGen.jl b/test/GrpAb/GrpAbFinGen.jl
index d43ff09f0f..f0ceb696dd 100644
--- a/test/GrpAb/GrpAbFinGen.jl
+++ b/test/GrpAb/GrpAbFinGen.jl
@@ -183,6 +183,18 @@
     H = abelian_group([4])
     K = direct_product(G, H)[1]
     @test is_isomorphic(K, abelian_group([60]))
+    K2 = direct_sum(H, G)[1]
+    @test is_isomorphic(K, K2)
+    K3, proj, inj = biproduct(G, H)
+    for i in 1:2, j in 1:2
+      if i == j
+        @test is_injective(inj[i])
+        @test is_surjective(proj[j])
+        @test isone(compose(inj[i], proj[j]).map)
+      else
+        @test is_zero(compose(inj[i], proj[j]))
+      end
+    end
   end
 
   @testset "Torsion" begin
diff --git a/test/Misc/Matrix.jl b/test/Misc/Matrix.jl
index 91daa707e2..7f8fda03be 100644
--- a/test/Misc/Matrix.jl
+++ b/test/Misc/Matrix.jl
@@ -35,7 +35,7 @@
 end
 
 @testset "Multiplicative order" begin
-  L = orthogonal_sum(root_lattice(:A, 7), root_lattice(:E, 8))[1]
+  L = direct_sum(root_lattice(:A, 7), root_lattice(:E, 8))[1]
   ao = automorphism_group_order(L)
   aag = automorphism_group_generators(L)
   for M in aag
diff --git a/test/QuadForm/Herm/Genus.jl b/test/QuadForm/Herm/Genus.jl
index f7617935ca..9e064b7739 100644
--- a/test/QuadForm/Herm/Genus.jl
+++ b/test/QuadForm/Herm/Genus.jl
@@ -394,8 +394,8 @@
   g = @inferred genus(HermLat, E1, q, [(0, 2, 1)])
   M = representative(g)
   G2 = genus(M)
-  @test genus(representative(orthogonal_sum(G, G2))) == orthogonal_sum(G, G2)
-  @test genus(representative(orthogonal_sum(G2, G))) == orthogonal_sum(G2, G)
+  @test genus(representative(direct_sum(G, G2))) == direct_sum(G, G2)
+  @test genus(representative(direct_sum(G2, G))) == direct_sum(G2, G)
 
   rlp = real_places(K)
   sig = Dict(rlp[1] => 2, rlp[2] => 2)
@@ -436,8 +436,8 @@
   for i in 1:10
     g1 = rand(G)
     g2 = rand(G)
-    g3 = @inferred orthogonal_sum(g1, g2)
-    @test genus(representative(g3), p) == genus(orthogonal_sum(representative(g1), representative(g2))[1], p)
+    g3 = @inferred direct_sum(g1, g2)
+    @test genus(representative(g3), p) == genus(direct_sum(representative(g1), representative(g2))[1], p)
   end
 
   GG = G[1]
@@ -449,8 +449,8 @@
   for i in 1:10
     g1 = rand(G)
     g2 = rand(G)
-    g3 = @inferred orthogonal_sum(g1, g2)
-    @test genus(representative(g3), p) == genus(orthogonal_sum(representative(g1), representative(g2))[1], p)
+    g3 = @inferred direct_sum(g1, g2)
+    @test genus(representative(g3), p) == genus(direct_sum(representative(g1), representative(g2))[1], p)
   end
 
   p = prime_decomposition(maximal_order(K), 17)[1][1]
@@ -494,7 +494,7 @@
   for i in 1:10
     g1 = rand(G)
     g2 = rand(G)
-    M, = @inferred orthogonal_sum(representative(g1), representative(g2))
+    M, = @inferred direct_sum(representative(g1), representative(g2))
     @test M in (g1 + g2)
   end
 
diff --git a/test/QuadForm/Lattices.jl b/test/QuadForm/Lattices.jl
index d77af773a2..5c5a3393e9 100644
--- a/test/QuadForm/Lattices.jl
+++ b/test/QuadForm/Lattices.jl
@@ -274,3 +274,29 @@ end
   L13orth = @inferred orthogonal_submodule(L1, L13)
   @test rank(intersect(L13clos1, L13orth)) == 0
 end
+
+@testset "Direct sums" begin
+  Qx, x = polynomial_ring(FlintQQ, "x")
+  f = x - 1
+  K, a = number_field(f, "a", cached = false)
+  Kt, t = polynomial_ring(K, "t")
+  g = t^2 + 1
+  E, b = number_field(g, "b", cached = false)
+  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
+  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-6, -10*b + 10, 0]), map(E, [-6*b + 7, 37//2*b + 21//2, -3//2*b + 5//2]), map(E, [-46*b + 71, 363//2*b + 145//2, -21//2*b + 49//2])]
+  gens2 = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-6, -10*b + 10, 0]), map(E, [-6*b + 7, 37//2*b + 21//2, -3//2*b + 5//2]), map(E, [1 + a + b, 1, 0])]
+  gens3 = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [-6*b + 7, 37//2*b + 21//2, -3//2*b + 5//2]), map(E, [2 + 2*a + 2*b, 2, 0])]
+  L1 = hermitian_lattice(E, gens, gram = D)
+  L2 = hermitian_lattice(E, gens2, gram = D)
+  L3 = hermitian_lattice(E, gens3, gram = D)
+  L4 = hermitian_lattice(E, gens, gram = 2*D)
+  @test genus(direct_sum(L1, L2)[1]) == direct_sum(genus(L1), genus(L2))
+  @test genus(direct_product(L3, L4)[1]) == direct_sum(genus(L3), genus(L4))
+  L5, inj, proj = @inferred biproduct(L1, L2, L3, L4)
+  for i in 1:4, j in 1:4
+    f = compose(inj[i], proj[j])
+    m = f.matrix
+    @test i == j ? isone(m) : iszero(m)
+  end
+end
+
diff --git a/test/QuadForm/Quad/Genus.jl b/test/QuadForm/Quad/Genus.jl
index 564d8c4465..6446e5043f 100644
--- a/test/QuadForm/Quad/Genus.jl
+++ b/test/QuadForm/Quad/Genus.jl
@@ -43,7 +43,7 @@
   @test G2 != G1
   @test G2 != G1a
 
-  # test representative and orthogonal_sum
+  # test representative and direct_sum
 
   G = Hecke.quadratic_local_genera(K, p2, rank = 3, det_val = 4)
   for i in 1:10
@@ -52,8 +52,8 @@
     L1 = representative(G1)
     @test L1 in G1
     L2 = representative(G2)
-    G3 = @inferred orthogonal_sum(G1, G2)
-    L3, = orthogonal_sum(L1, L2)
+    G3 = @inferred direct_sum(G1, G2)
+    L3, = direct_sum(L1, L2)
     @test G3 == genus(L3, p2)
     @test genus(L1,p2) + genus(QuadLat,p2) == G1
   end
@@ -66,8 +66,8 @@
     G2 = rand(G)
     L1 = representative(G1)
     L2 = representative(G2)
-    G3 = @inferred orthogonal_sum(G1, G2)
-    L3, = orthogonal_sum(L1, L2)
+    G3 = @inferred direct_sum(G1, G2)
+    L3, = direct_sum(L1, L2)
     @test G3 == genus(L3, p2)
   end
 
@@ -78,8 +78,8 @@
     L1 = representative(G1)
     @test L1 in G1
     L2 = representative(G2)
-    G3 = @inferred orthogonal_sum(G1, G2)
-    L3, = orthogonal_sum(L1, L2)
+    G3 = @inferred direct_sum(G1, G2)
+    L3, = direct_sum(L1, L2)
     @test G3 == genus(L3, p2)
   end
 
@@ -89,8 +89,8 @@
     G2 = rand(G)
     L1 = representative(G1)
     L2 = representative(G2)
-    G3 = @inferred orthogonal_sum(G1, G2)
-    L3, = orthogonal_sum(L1, L2)
+    G3 = @inferred direct_sum(G1, G2)
+    L3, = direct_sum(L1, L2)
     @test G3 == genus(L3, p3)
     @test genus(L1,p3) + genus(QuadLat,p3) == G1
   end
@@ -108,8 +108,8 @@
     G2 = rand(G)
     L1 = representative(G1)
     L2 = representative(G2)
-    G3 = @inferred orthogonal_sum(G1, G2)
-    L3, = orthogonal_sum(L1, L2)
+    G3 = @inferred direct_sum(G1, G2)
+    L3, = direct_sum(L1, L2)
     @test G3 == genus(L3, p5)
   end
 
@@ -119,8 +119,8 @@
     G2 = rand(G)
     L1 = representative(G1)
     L2 = representative(G2)
-    G3 = @inferred orthogonal_sum(G1, G2)
-    L3, = orthogonal_sum(L1, L2)
+    G3 = @inferred direct_sum(G1, G2)
+    L3, = direct_sum(L1, L2)
     @test L3 in G3
   end
 
diff --git a/test/QuadForm/Quad/Lattices.jl b/test/QuadForm/Quad/Lattices.jl
index 8d3cbaed25..3179662ec6 100644
--- a/test/QuadForm/Quad/Lattices.jl
+++ b/test/QuadForm/Quad/Lattices.jl
@@ -15,7 +15,8 @@
   @test pseudo_matrix(L1) == pseudo_matrix(L)
   @test ambient_space(L1) != ambient_space(L)
 
-  LL,i,j = orthogonal_sum(L1,L1)
+  LL, inj = direct_sum(L1,L1)
+  i, j = inj
   @inferred i(L1)
   @test i(L1)+j(L1) == LL
 
diff --git a/test/QuadForm/Quad/Spaces.jl b/test/QuadForm/Quad/Spaces.jl
index 83205bdc89..3ba8424f9c 100644
--- a/test/QuadForm/Quad/Spaces.jl
+++ b/test/QuadForm/Quad/Spaces.jl
@@ -114,7 +114,7 @@
     @test is_diagonal(gram_matrix(V, B))
 
     V1 = quadratic_space(K, zero_matrix(K,2,2))
-    V2, _, _ = orthogonal_sum(V, V1)
+    V2, _, _ = biproduct(V, V1)
     B2 = @inferred orthogonal_basis(V2)
     @test is_diagonal(gram_matrix(V2, B2))
 
@@ -340,12 +340,12 @@
     g = Hecke.isometry_class(q)
     gg = Hecke.isometry_class(qq)
     gg_deg = Hecke.isometry_class(q_deg)
-    ggg = Hecke.isometry_class(orthogonal_sum(q,qq)[1])
+    ggg = Hecke.isometry_class(direct_product(q,qq)[1])
     @test g + gg == ggg
     @test g + gg - g == gg
     @test g + g + gg - g == gg+ g
     @test is_locally_represented_by(q, qq, 2) == represents(local_symbol(g,2), local_symbol(gg,2))
-    @test !is_locally_represented_by(orthogonal_sum(q,q)[1],q, 2)
+    @test !is_locally_represented_by(direct_sum(q,q)[1],q, 2)
     @test represents(gg, -1)
     @test represents(gg, 3)
     @test represents(gg, 2)
@@ -371,7 +371,7 @@
     @test Hecke.isometry_class(representative(gg + gg + g)) == gg + gg + g
     @test Hecke.isometry_class(representative(g+g+gg+gg)) == g + g + gg+gg
     gdegenerte = Hecke.isometry_class(quadratic_space(QQ,zero_matrix(QQ,2,2)))
-    h = orthogonal_sum(g,gdegenerte)
+    h = direct_sum(g,gdegenerte)
     @test is_isotropic(h)
     @test is_isotropic(local_symbol(h,2))
     @test is_isotropic(local_symbol(h,3))
@@ -387,7 +387,7 @@
     @test Hecke.is_isotropic(qq, infF2)
     @test Hecke._isisotropic_with_vector(gram_matrix(h))[1]
     @test !Hecke._isisotropic_with_vector(gram_matrix(q))[1]
-    hh,_,_ = orthogonal_sum(qq,quadratic_space(F,-gram_matrix(qq)))
+    hh,_,_ = biproduct(qq,quadratic_space(F,-gram_matrix(qq)))
     i = Hecke._maximal_isotropic_subspace(gram_matrix(hh))
     @test nrows(i)==dim(qq)
     @test i*gram_matrix(hh)*transpose(i) == 0
@@ -495,16 +495,16 @@
   end
 end
 
-@testset "orthogomal sums" begin
+@testset "direct sums" begin
   K, _ = cyclotomic_field(27, cached=false)
   V = quadratic_space(K, 4)
-  S, inj, proj = @inferred direct_sum(V, rescale(V, -1//5))
+  S, inj, proj = @inferred biproduct(V, rescale(V, -1//5))
   @test dim(S) == 8
   for i in 1:2, j in 1:2
     f = compose(inj[i], proj[j])
     m = f.matrix
     @test i != j ? iszero(m) : isone(m)
   end
-  S, inj, proj = @inferred direct_sum(V, V, V)
-  @test_throws ArgumentError direct_sum(V)
+  S, inj, proj = @inferred biproduct(V, V, V)
+  @test_throws ArgumentError direct_product(V)
 end
diff --git a/test/QuadForm/Quad/ZGenus.jl b/test/QuadForm/Quad/ZGenus.jl
index a14e86a3c5..58b91dd10e 100644
--- a/test/QuadForm/Quad/ZGenus.jl
+++ b/test/QuadForm/Quad/ZGenus.jl
@@ -172,13 +172,13 @@
   g3 = genus(diagonal_matrix(map(ZZ,[1,3,27])), 3)
   n3 = genus(matrix(ZZ,0,0,[]),3)
   n5 = genus(matrix(ZZ,0,0,[]),5)
-  @test g3 == orthogonal_sum(n3, g3)
-  @test_throws ArgumentError orthogonal_sum(n3, n5)
-  @test n3 == orthogonal_sum(n3, n3)
+  @test g3 == direct_sum(n3, g3)
+  @test_throws ArgumentError direct_sum(n3, n5)
+  @test n3 == direct_sum(n3, n3)
   @test Hecke._species_list(g3) == [1, 1, 1]
   h3 = genus(diagonal_matrix(map(ZZ,[1,3,9])), 3)
   @test Hecke._standard_mass(h3) ==  9//16
-  @test orthogonal_sum(g3,h3)==direct_sum(h3,g3)
+  @test direct_sum(g3,h3)==direct_sum(h3,g3)
 
 
   # These examples are taken from Table 2 of [CS1988]_::
@@ -255,7 +255,7 @@
   @test order(q) == 1
   L2 = Zlattice(gram=2*ZZ[2 1; 1 2])
   G2 = genus(L2)
-  @test genus(orthogonal_sum(L,L2)[1]) == orthogonal_sum(G, G2)
+  @test genus(direct_sum(L,L2)[1]) == direct_sum(G, G2)
   @test length(representatives(G2)) == 1
   @test representative(G2)===representative(G2) # caching
 
@@ -507,7 +507,7 @@ end
   L = root_lattice(:A, 4)
   G = genus(L)
   G2 = rescale(G, -1//13)
-  G22 = orthogonal_sum(G2, G)
+  G22 = direct_sum(G2, G)
   L2 = representative(G22)
   @test genus(L2) == G22
   G = Zgenera((0,8), 1)[1]
diff --git a/test/QuadForm/Quad/ZLattices.jl b/test/QuadForm/Quad/ZLattices.jl
index 5bb6999faf..49cf9dd46b 100644
--- a/test/QuadForm/Quad/ZLattices.jl
+++ b/test/QuadForm/Quad/ZLattices.jl
@@ -177,11 +177,11 @@ end
     @test gram_matrix(V) == gram_matrix(L)
   end
 
-  # orthogonal sum
+  # direct sum
 
   for L in [Lr0, Lr1, Lr2]
     for LL in [Lr0, Lr1, Lr2]
-      LLL, = @inferred orthogonal_sum(L, LL)
+      LLL, = @inferred direct_product(L, LL)
       @test gram_matrix(LLL) == diagonal_matrix(gram_matrix(L), gram_matrix(LL))
     end
   end
@@ -194,7 +194,8 @@ end
   # root lattice recognition
 
   L = Zlattice(gram=ZZ[4;])
-  LL,i,j = orthogonal_sum(L,L)
+  LL, inj = direct_sum(L,L)
+  i, j = inj
   @test LL == i(L)+j(L)
   @test L == preimage(i, LL)
   R = @inferred root_sublattice(L)
@@ -203,7 +204,7 @@ end
   R = lattice(ambient_space(L),basis_matrix(L)[1,:])
   @test rank(root_sublattice(R))==1
 
-  L = orthogonal_sum(root_lattice(:A,2),root_lattice(:D,4))[1]
+  L = biproduct(root_lattice(:A,2),root_lattice(:D,4))[1]
   R = root_lattice_recognition(L)
   @test length(R[1]) == 2
   @test (:D,4) in R[1] && (:A,2) in R[1]
@@ -247,7 +248,8 @@ end
       if rank(L) > 4 # small examples suffice
         continue
       end
-      C1L, _, f = orthogonal_sum(C1, L)
+      C1L, inj = direct_sum(C1, L)
+      f = inj[2]
       V = ambient_space(C1L)
       imL = lattice(V, f.matrix)
       Ge = automorphism_group_generators(imL, ambient_representation = true)
@@ -343,7 +345,7 @@ end
   V = quadratic_space(QQ, QQ[1 0 0; 0 1 0; 0 0 1])
   L = lattice(V, ZZ[1 -1 0; 0 1 -1])
   S = lattice(V, ZZ[1 -1 0;])
-  submod = Hecke.orthogonal_submodule(L, S)
+  submod = orthogonal_submodule(L, S)
   @test  basis_matrix(submod) == matrix(QQ, 1, 3, [1 1 -2])
 
   @test is_definite(L)
@@ -362,7 +364,7 @@ end
   V = quadratic_space(QQ,gram)
   L = lattice(V, BS)
   H = lattice(V, BH)
-  R = Hecke.orthogonal_submodule(L,H)
+  R = orthogonal_submodule(L,H)
   @test is_sublattice(L,R)
 
   # local modification
@@ -660,11 +662,11 @@ end
   @test !Hecke.is_isometric(L1, L2)
 end
 
-@testset "orthogonal sums" begin
+@testset "direct sums" begin
   L1 = root_lattice(:A, 6)
   L2 = root_lattice(:E, 7)
-  L, inj, proj = @inferred direct_sum([L1, L2])
-  @test genus(L) == orthogonal_sum(genus(L1), genus(L2))
+  L, inj, proj = @inferred biproduct([L1, L2])
+  @test genus(L) == direct_sum(genus(L1), genus(L2))
   for i in 1:2, j in 1:2
     f = compose(inj[i], proj[j])
     m = f.matrix
diff --git a/test/QuadForm/Torsion.jl b/test/QuadForm/Torsion.jl
index a39b30dd1d..b9f70505ba 100644
--- a/test/QuadForm/Torsion.jl
+++ b/test/QuadForm/Torsion.jl
@@ -234,7 +234,7 @@
   rq, i = radical_quadratic(Tsub)
   bool, j = @inferred has_complement(i)
   N = domain(j)
-  T2, _, _ = orthogonal_sum(rq, N)
+  T2, _ = direct_sum(rq, N)
   @test is_degenerate(T2)
   bool, phi = @inferred is_isometric_with_isometry(Tsub, T2)
   @test bool
@@ -287,30 +287,31 @@
   @test is_bijective(phi)
   @test all(a -> Hecke.quadratic_product(a) == (-1)*Hecke.quadratic_product(phi(a)), T1sub)
 
-  # orthogonal sum
+  # direct sums
 
   B = matrix(FlintQQ, 3, 3 ,[1, 1, 0, 1, -1, 0, 0, 1, -1])
   G = matrix(FlintQQ, 3, 3 ,[1, 0, 0, 0, 1, 0, 0, 0, 1])
   L1 = Zlattice(B, gram = G)
   qL1 = discriminant_group(L1)
   Z = torsion_quadratic_module(QQ[1;])
-  @test_throws ArgumentError orthogonal_sum(qL1, Z)
-  @test_throws ArgumentError orthogonal_sum(qL1, rescale(Z, 2))
+  @test_throws ArgumentError direct_sum(qL1, Z)
+  @test_throws ArgumentError direct_product(qL1, rescale(Z, 2))
 
   B = matrix(FlintQQ, 4, 4 ,[2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 2, 0, 1, 1, 1, 1])
   G = matrix(FlintQQ, 4, 4 ,[1//2, 0, 0, 0, 0, 1//2, 0, 0, 0, 0, 1//2, 0, 0, 0, 0, 1//2])
   L2 = Zlattice(B, gram = G)
   qL2 = discriminant_group(L2)
   Z = torsion_quadratic_module(QQ[2;])
-  q, _, _ = @inferred orthogonal_sum(qL2, Z)
+  q, _ = @inferred direct_product(qL2, Z)
   @test is_isometric_with_isometry(q, qL2)[1]
   @test modulus_bilinear_form(q) == modulus_bilinear_form(qL2)
   @test modulus_quadratic_form(q) == modulus_quadratic_form(Z)
 
-  L3, _, _ = orthogonal_sum(L1, L2)
+  L3, _ = direct_product(L1, L2)
   qL3 = discriminant_group(L3)
 
-  q, qL1inq, qL2inq = @inferred orthogonal_sum(qL1, qL2)
+  q, inj = @inferred direct_sum(qL1, qL2)
+  qL1inq, qL2inq = inj
   @test is_injective(qL1inq) && is_injective(qL2inq)
   bool, _ = is_isometric_with_isometry(qL3, q)
   @test bool
@@ -331,14 +332,14 @@
     qL = discriminant_group(L)
     @test is_elementary(qL, 9-i)
   end
-  L = orthogonal_sum(root_lattice(:A, 7), root_lattice(:D, 7))[1]
+  L = biproduct(root_lattice(:A, 7), root_lattice(:D, 7))[1]
   qL = discriminant_group(L)
   @test is_primary(qL, 2) && !is_elementary(qL, 2)
 
-  # Additional orthogonal sum
+  # Additional direct sum
   list_lat = [root_lattice(:A, i) for i in 1:7]
   list_quad = discriminant_group.(list_lat)
-  S, inj, proj = direct_sum(list_quad)
+  S, inj, proj = biproduct(list_quad)
   @test order(S) == prod(order.(list_quad))
   for i in 1:7, j in 1:7
     f = compose(inj[i], proj[j])
@@ -347,8 +348,14 @@
   end
 
   T = discriminant_group(root_lattice(:D, 6))
-  S, inj, proj = @inferred direct_sum(T, T, T)
+  S, inj, proj = @inferred biproduct(T, T, T)
 
+  U2 = hyperbolic_plane_lattice(2)
+  q = discriminant_group(U2)
+  qq, qqinq = sub(q, [q[1]+q[2]])
+  ok, qqqinq = has_complement(qqinq)
+  @test ok
+  @test is_isometric_with_isometry(qq + domain(qqqinq), q)[1]
 
   # Smith normal form
   L = Zlattice(gram=matrix(ZZ, [[2, -1, 0, 0, 0, 0],[-1, 2, -1, -1, 0, 0],[0, -1, 2, 0, 0, 0],[0, -1, 0, 2, 0, 0],[0, 0, 0, 0, 6, 3],[0, 0, 0, 0, 3, 6]]))