diff --git a/Project.toml b/Project.toml
index 638e6ccc8a5e..455384a868df 100644
--- a/Project.toml
+++ b/Project.toml
@@ -39,7 +39,7 @@ CodecZlib = "0.7.8"
 Compat = "4.13.0"
 Distributed = "1.6"
 GAP = "0.16.3"
-Hecke = "0.39.14"
+Hecke = "0.39.15"
 JSON = "1.0.1"
 JSON3 = "1.13.2"
 LazyArtifacts = "1.6"
diff --git a/src/Groups/abelian_aut.jl b/src/Groups/abelian_aut.jl
index 7d928373abab..ad49b26d9a59 100644
--- a/src/Groups/abelian_aut.jl
+++ b/src/Groups/abelian_aut.jl
@@ -355,8 +355,12 @@ julia> order(OT)
 ```
 """
 @attr AutomorphismGroup{TorQuadModule} function orthogonal_group(T::TorQuadModule)
+  return _orthogonal_group(T, _orthogonal_group_gens(T); check=false)
+end
+
+function _orthogonal_group_gens(T::TorQuadModule)
   if is_trivial(abelian_group(T))
-    return _orthogonal_group(T, ZZMatrix[identity_matrix(ZZ, ngens(T))], check = false)
+    return ZZMatrix[identity_matrix(ZZ, ngens(T))]
   elseif is_semi_regular(T)
     # if T is semi-regular, it is isometric to its normal form for which
     # we know how to compute the isometries.
@@ -364,20 +368,20 @@ julia> order(OT)
     j = inv(i)
     gensON_mat = unique(_compute_gens(N))
     gensON = TorQuadModuleMap[hom(N, N, g) for g in gensON_mat]
-    gensOT = TorQuadModuleMap[i * g * j for g in gensON]
-    length(gensOT) > 1 ? filter!(!isone∘matrix, gensOT) : nothing
-    return _orthogonal_group(T, gensOT; check=false)
+    gensOT_mat = ZZMatrix[matrix(i * g * j) for g in gensON]
+    length(gensOT_mat) > 1 ? filter!(!isone, gensOT_mat) : nothing
+    return gensOT_mat
   elseif iszero(gram_matrix_quadratic(T))
     # in that case, we don't have any conditions regarding the
     # quadratic form, so we have all automorphisms coming
     # from the underlying abelian group
-    return _orthogonal_group(T, hom.(gens(automorphism_group(abelian_group(T)))); check=false)
+    return matrix.(gens(automorphism_group(abelian_group(T))))
   else
     # if T is not semi-regular, we distinghuish the cases whether or not
     # it splits its radical quadratic
     i = radical_quadratic(T)[2]
     gensOT_mat = has_complement(i)[1] ? _compute_gens_split_degenerate(T) : _compute_gens_non_split_degenerate(T)
-    return _orthogonal_group(T, gensOT_mat; check=false)
+    return gensOT_mat
   end
 end
 
diff --git a/src/NumberTheory/QuadFormAndIsom.jl b/src/NumberTheory/QuadFormAndIsom.jl
index 384509340367..a09498d569f4 100644
--- a/src/NumberTheory/QuadFormAndIsom.jl
+++ b/src/NumberTheory/QuadFormAndIsom.jl
@@ -11,4 +11,5 @@ include("QuadFormAndIsom/finite_group_actions.jl")
 include("QuadFormAndIsom/hermitian_miranda_morrison.jl")
 include("QuadFormAndIsom/enumeration.jl")
 include("QuadFormAndIsom/embeddings.jl")
+include("QuadFormAndIsom/gluing_factory.jl")
 include("QuadFormAndIsom/printings.jl")
diff --git a/src/NumberTheory/QuadFormAndIsom/embeddings.jl b/src/NumberTheory/QuadFormAndIsom/embeddings.jl
index eccb9f0b335a..f8d51cf0b4ea 100644
--- a/src/NumberTheory/QuadFormAndIsom/embeddings.jl
+++ b/src/NumberTheory/QuadFormAndIsom/embeddings.jl
@@ -7,18 +7,22 @@
 function __direct_sum(
     A::TorQuadModule,
     B::TorQuadModule,
-    same_ambient::Bool = false,
+    same_ambient::Bool = false;
+    as_bilinear_module::Bool=false,
   )
   if !same_ambient
-    return direct_sum(A, B; cached=false)
+    return direct_sum(A, B; cached=false, as_bilinear_module)
   end
 
   # Test that A and B have the same moduli and same ambient quadratic space
   @assert modulus_bilinear_form(A) == modulus_bilinear_form(B)
-  @assert modulus_quadratic_form(A) == modulus_quadratic_form(B)
+  @assert as_bilinear_module || modulus_quadratic_form(A) == modulus_quadratic_form(B)
   @assert ambient_space(cover(A)) === ambient_space(cover(B))
   V = ambient_space(cover(A))
 
+  mbf = modulus_bilinear_form(A)
+  mqf = as_bilinear_module ? mbf : modulus_quadratic_form(A)
+
   # Test that the modules are indeed in orthogonal direct sum
   # It should be used internally so we do not want to test that all the time
   @hassert :ZZLatWithIsom 1 iszero(inner_product(V, basis_matrix(cover(A)), basis_matrix(cover(B))))
@@ -32,7 +36,7 @@ function __direct_sum(
   gensDB = Vector{QQFieldElem}[lift(b) for b in gens(B)]
   # D is A\oplus B, and we fix a set of generators given first the generators
   # of A and then the ones of B
-  D = torsion_quadratic_module(covD, relD; gens=union!(gensDA, gensDB), modulus=modulus_bilinear_form(A), modulus_qf=modulus_quadratic_form(A))
+  D = torsion_quadratic_module(covD, relD; gens=union!(gensDA, gensDB), modulus=mbf, modulus_qf=mqf)
   IA = identity_matrix(ZZ, ngens(A))
   IB = identity_matrix(ZZ, ngens(B))
   AinD = hom(A, D, reduce(hcat, ZZMatrix[IA, zero_matrix(ZZ, ngens(A), ngens(B))]))
@@ -76,13 +80,14 @@ the respective covers of ``A`` and ``B``.
 function _direct_sum_with_embeddings_orthogonal_groups(
     A::TorQuadModule,
     B::TorQuadModule;
-    same_ambient::Bool=(ambient_space(cover(A)) === ambient_space(cover(B)))
+    as_bilinear_module::Bool=false,
+    OA::AutomorphismGroup{TorQuadModule}=__orthogonal_group(A; as_bilinear_module),
+    OB::AutomorphismGroup{TorQuadModule}=__orthogonal_group(B; as_bilinear_module),
+    same_ambient::Bool=(ambient_space(cover(A)) === ambient_space(cover(B))),
   )
-  D, inj = __direct_sum(A, B, same_ambient)
+  D, inj = __direct_sum(A, B, same_ambient; as_bilinear_module)
   AinD, BinD = inj
   OD = orthogonal_group(D)
-  OA = orthogonal_group(A)
-  OB = orthogonal_group(B)
   IA = identity_matrix(ZZ, ngens(A))
   IB = identity_matrix(ZZ, ngens(B))
 
@@ -118,7 +123,7 @@ function _get_V(
     p::IntegerUnion,
   )
   q = domain(f)
-  V, Vinq = kernel(hom(q, q, p*identity_matrix(ZZ, ngens(q)))) # TODO: remove once Hecke function implemented
+  V, Vinq = torsion_subgroup(q, p)
   fpV = restrict_endomorphism(f, Vinq; check=false)
   fpV = evaluate(mu, fpV)
   V, _ = kernel(fpV)
@@ -202,31 +207,19 @@ end
 # discriminant groups of $A$ and $B$ respectively. The third output of the
 # algorithm is the graph of $\gamma$ in $D$.
 #
-# The fourth and fifth inputs are optional isometries $f_A$ and $f_B$ of $A$
-# and $B$, respectively, to be extended along $\gamma$, assuming it to be an
-# $(f_A, f_B)$-equivariant gluing. The second output of the algorithm is
-# the extension $f_A\oplus f_B \in O(C)$ --- it is the identity by default.
-#
-# !!! note
-#     The isometries `fA` and `fB` are given in terms of their matrix
-#     representation on `A` and `B` respectively, not by their ambient
-#     representation. Similarly for the second output of the algorithm.
-#
 # !!! warning
 #     No sanity checks are performed: in particular, if $D_A$ and $D_B$ are not
-#     in orthogonal direct sum, if $\gamma$ is not a gluing, or if it is not
-#     $(f_A, f_B)$-equivariant, then the outputs might not satisfy the expected
-#     conditions.
+#     in orthogonal direct sum, or if $\gamma$ is not a gluing, then the outputs
+#     might not satisfy the expected conditions.
 #
-# If `same_ambient` is set to `true`, then the two lattices $A$ and $B$ are
-# considered as lying in the same ambient quadratic space.
-function _overlattice(
-    gamma::TorQuadModuleMap,
-    HAinD::TorQuadModuleMap,
-    HBinD::TorQuadModuleMap,
-    fA::QQMatrix = identity_matrix(QQ, rank(relations(domain(HAinD)))),
-    fB::QQMatrix = identity_matrix(QQ, rank(relations(domain(HBinD)))),
-  )
+# If we glue along the trivial subgroups of `D_A` and `D_B`, one can drop the
+# glue map gamma as first input
+
+function _overlattice_with_graph(
+  gamma::TorQuadModuleMap,
+  HAinD::TorQuadModuleMap,
+  HBinD::TorQuadModuleMap,
+)
   HA = domain(HAinD)
   HB = domain(HBinD)
   A = relations(HA)
@@ -241,25 +234,14 @@ function _overlattice(
   _FakeB = hnf(Fakeglue)
   _B = QQ(1, denominator(Fakeglue))*change_base_ring(QQ, numerator(_FakeB))
   C = lattice(ambient_space(L), _B[end-rank(A)-rank(B)+1:end, :])
-  fC = block_diagonal_matrix(QQMatrix[fA, fB])
-  _B = coordinates(basis_matrix(C), L)
-  fC = _B*fC*inv(_B)
-  @hassert :ZZLatWithIsom 1 fC*gram_matrix(C)*transpose(fC) == gram_matrix(C)
   _, graph = sub(D, D.(_glue))
-  return C, fC, graph
+  return C, graph
 end
 
-# Same as above where we glue along the trivial subgroups of `D_A` and `D_B`.
-# In that particular case, `D_A` and `D_B` are the discriminant groups of the
-# lattices considered (so HA = L^{\vee}/L for some lattice integral lattice L,
-# and same for HB), and D is the orthogonal direct sum of HA and HB in an
-# appropriate quadratic space.
-function _overlattice(
-    DAinD::TorQuadModuleMap,
-    DBinD::TorQuadModuleMap,
-    fA::QQMatrix = identity_matrix(QQ, rank(relations(domain(HAinD)))),
-    fB::QQMatrix = identity_matrix(QQ, rank(relations(domain(HBinD)))),
-  )
+function _overlattice_with_graph(
+  DAinD::TorQuadModuleMap,
+  DBinD::TorQuadModuleMap,
+)
   DA = domain(DAinD)
   DB = domain(DBinD)
   zA, zAinDA = sub(DA, TorQuadModuleElem[])
@@ -267,54 +249,61 @@ function _overlattice(
   zAinD = compose(zAinDA, DAinD)
   zBinD = compose(zBinDB, DBinD)
   gamma = hom(zA, zB, zero_matrix(ZZ, 0, 0))
-  return _overlattice(gamma, zAinD, zBinD, fA, fB)
+  return _overlattice_with_graph(gamma, zAinD, zBinD)
 end
 
-###############################################################################
+# Same as before, with the fourth and fifth inputs being optional isometries
+# $f_A$ and $f_B$ of $A$ and $B$, respectively, to be extended along $\gamma$,
+# assuming it to be an $(f_A, f_B)$-equivariant gluing. The second output of
+# the algorithm is the extension $f_A\oplus f_B \in O(C)$ --- it is the
+# identity by default.
 #
-#  Generic primitive extensions method
+# !!! note
+#     The isometries `fA` and `fB` are given in terms of their matrix
+#     representation on `A` and `B` respectively, not by their ambient
+#     representation. Similarly for the second output of the algorithm.
 #
-###############################################################################
+# !!! warning
+#     Again, no input checks are performed so if $\gamma$ is not equivariant
+#     with respect to $fA$ and $fB$, then the outputs might not satisfy the
+#     expected conditions.
+
+function _equivariant_overlattice_with_graph(
+  gamma::TorQuadModuleMap,
+  HAinD::TorQuadModuleMap,
+  HBinD::TorQuadModuleMap,
+  fA::QQMatrix = identity_matrix(QQ, rank(relations(domain(HAinD)))),
+  fB::QQMatrix = identity_matrix(QQ, rank(relations(domain(HBinD)))),
+)
+  C, graph = _overlattice_with_graph(gamma, HAinD, HBinD)
+  L = relations(codomain(HAinD))
+  fC = block_diagonal_matrix(QQMatrix[fA, fB])
+  _B = coordinates(basis_matrix(C), L)
+  fC = _B*fC*inv(_B)
+  @hassert :ZZLatWithIsom 1 fC*gram_matrix(C)*transpose(fC) == gram_matrix(C)
+  return C, fC, graph
+end
 
-# Given a torsion quadratic module $q_M$, return the associated torsion
-# bilinear module. Here $G_M$ and $f_{q_M}$ are respectively a group of
-# isometries and an isometry of $q_M$, for the quadratic form. The function
-# also returns the associated group of isometries and isometry for the torsion
-# bilinear module associated to $q_M$.
-function _change_to_bilinear_module(
-    qM::TorQuadModule,
-    GM::AutomorphismGroup{TorQuadModule},
-    fqM::TorQuadModuleMap,
-  )
-  qM = Hecke._as_finite_bilinear_module(qM) # TODO: to be changed
-  OqM = orthogonal_group(qM)
-  GM, _ = sub(OqM, elem_type(OqM)[OqM(matrix(g); check=false) for g in gens(GM)])
-  fqM = hom(qM, qM, matrix(fqM))
-  return qM, OqM, GM, fqM
+function _equivariant_overlattice_with_graph(
+  DAinD::TorQuadModuleMap,
+  DBinD::TorQuadModuleMap,
+  fA::QQMatrix = identity_matrix(QQ, rank(relations(domain(HAinD)))),
+  fB::QQMatrix = identity_matrix(QQ, rank(relations(domain(HBinD)))),
+)
+  C, graph = _overlattice_with_graph(DAinD, DBinD)
+  L = relations(codomain(HAinD))
+  fC = block_diagonal_matrix(QQMatrix[fA, fB])
+  _B = coordinates(basis_matrix(C), L)
+  fC = _B*fC*inv(_B)
+  @hassert :ZZLatWithIsom 1 fC*gram_matrix(C)*transpose(fC) == gram_matrix(C)
+  return C, fC, graph
 end
 
-# Given two torsion quadratic/bilinear modules $q_M$ and $q_N$, return the
-# list of possible orders for anti-isometric submodules of $q_M$ and $q_N$
-# respectively.
+###############################################################################
 #
-# One can choose a default value `glue_order`, useful in functions with keyword
-# arguments.
+#  Generic primitive extensions method
 #
-# TODO: To be removed once appropriate Hecke functions implemented
-function _possible_glue_orders(
-    qM::TorQuadModule,
-    qN::TorQuadModule,
-  )
-  _gcd = ZZ(1)
-  snM = reverse!(elementary_divisors(qM))
-  snN = reverse!(elementary_divisors(qN))
-  k = min(length(snM), length(snN))
-  for i in 1:k
-    mul!(_gcd, _gcd, gcd(snM[i], snN[i]))
-  end
-  pos_ord = divisors(_gcd)
-  return pos_ord
-end
+###############################################################################
 
 # Given an abelian group isomorphism $\phi$ between two torsion
 # quadratic/bilinear modules $H_M$ and $H_N$, and given two isometries
@@ -415,6 +404,17 @@ function _fitting_isometries(
   return reporb
 end
 
+function _as_sublattices(
+  L::ZZLat,
+  M::ZZLat,
+  N::ZZLat,
+)
+  V = ambient_space(L)
+  M2 = lattice(V, hcat(basis_matrix(M), zero_matrix(QQ, rank(M), degree(L) - degree(M))))
+  N2 = lattice(V, hcat(zero_matrix(QQ, rank(N), degree(L) - degree(N)), basis_matrix(N)))
+  return M2, N2
+end
+
 # We have a primitive extension `M\oplus N \to L`: we want to see M and N now
 # as sublattices of L.
 # If they are in the same ambient space, there is nothing to do.
@@ -807,7 +807,7 @@ function _primitive_extensions_generic(
           end
 
           for b in reporb
-            L, fL, graph = _overlattice(phig, HMinD, HNinD, fM, b)
+            L, fL, graph = _equivariant_overlattice_with_graph(phig, HMinD, HNinD, fM, b)
 
             if !isempty(Gs)
               any(isequal(genus(L)), Gs) || continue
@@ -837,20 +837,69 @@ end
 #
 ###############################################################################
 
+function _subgroups_orbit_representatives_and_stabilizers_primary_subtype(
+  Vinq::TorQuadModuleMap,
+  O::AutomorphismGroup{TorQuadModule},
+  p::ZZRingElem,
+  subtype::Vector{Int},
+  f::TorQuadModuleMap = id_hom(codomain(Vinq)),
+)
+  r1 = __subgroups_orbit_representatives_and_stabilizers_primary_subtype(Vinq, O, p, subtype, f)
+  r = Tuple{TorQuadModuleMap, AutomorphismGroup{TorQuadModule}}[(x[1], domain(x[2])) for x in r1]
+  return r
+end
+
+function __subgroups_orbit_representatives_and_stabilizers_primary_subtype(
+  Vinq::TorQuadModuleMap,
+  O::AutomorphismGroup{TorQuadModule},
+  p::ZZRingElem,
+  subtype::Vector{Int},
+  f::TorQuadModuleMap = id_hom(codomain(Vinq)),
+)
+  res = Tuple{TorQuadModuleMap, GAPGroupHomomorphism}[]
+
+  V = domain(Vinq)
+  q = codomain(Vinq)
+  A = abelian_group(V)
+  @assert is_snf(A)
+  flag = isone(matrix(f))
+
+  subs_it = Hecke._psubgroups(A, p; subtype=[subtype])
+
+  to_gap = get_attribute(O, :to_gap)
+  to_oscar = get_attribute(O, :to_oscar)
+  qgap = codomain(to_gap)
+  sgap = typeof(qgap)[]
+
+  for (H, i) in subs_it
+    Hgap, _ = sub(qgap, elem_type(qgap)[to_gap(Vinq(V(i(a)))) for a in gens(H)])
+    push!(sgap, Hgap)
+  end
+
+  m = gset(O, on_subgroups, sgap)
+  orbs = orbits(m)
+  for orb in orbs
+    _repgap = representative(orb)
+    _, rep = sub(q, TorQuadModuleElem[to_oscar(qgap(a)) for a in gens(_repgap)])
+    flag || is_invariant(f, rep) || continue
+    stab, jj = stabilizer(O, rep)
+    push!(res, (rep, jj))
+  end
+  return res
+end
+
 # Given the embedding of an `(O, f)`-stable finite quadratic submodule `V` of
 # `q`, compute representatives of `O`-orbits of `f`-stable submodules of `V` of
 # order `ord`. The stabilizers in `O` is also computed.
 #
 # Note that any torsion quadratic module `H` in output is given by an embedding
 # of `H` in `q`.
-#
-# TODO: Replace by dedicated function for p-group and subgroup types
 function _subgroups_orbit_representatives_and_stabilizers(
-    Vinq::TorQuadModuleMap,
-    O::AutomorphismGroup{TorQuadModule},
-    ord::IntegerUnion = -1,
-    f::Union{TorQuadModuleMap, AutomorphismGroupElem{TorQuadModule}} = id_hom(codomain(Vinq)),
-  )
+  Vinq::TorQuadModuleMap,
+  O::AutomorphismGroup{TorQuadModule},
+  ord::IntegerUnion = -1,
+  f::Union{TorQuadModuleMap, AutomorphismGroupElem{TorQuadModule}} = id_hom(codomain(Vinq)),
+)
   res = Tuple{TorQuadModuleMap, AutomorphismGroup{TorQuadModule}}[]
 
   V = domain(Vinq)
@@ -953,7 +1002,21 @@ function _subgroups_orbit_representatives_and_stabilizers_elementary(
     l::IntegerUnion = -1;
     algorithm::Symbol=:PermGroup,
   )
-  res = Tuple{TorQuadModuleMap, AutomorphismGroup{TorQuadModule}}[]
+  r1 = __subgroups_orbit_representatives_and_stabilizers_elementary(Vinq, G, ord, _p, f, l; algorithm)
+  r = Tuple{TorQuadModuleMap, AutomorphismGroup{TorQuadModule}}[(x[1], domain(x[2])) for x in r1]
+  return r
+end
+
+function __subgroups_orbit_representatives_and_stabilizers_elementary(
+    Vinq::TorQuadModuleMap,
+    G::AutomorphismGroup{TorQuadModule},
+    ord::IntegerUnion,
+    _p::IntegerUnion,
+    f::Union{TorQuadModuleMap, AutomorphismGroupElem{TorQuadModule}} = id_hom(codomain(Vinq)),
+    l::IntegerUnion = -1;
+    algorithm::Symbol=:PermGroup,
+  )
+  res = Tuple{TorQuadModuleMap, GAPGroupHomomorphism}[]
   p = ZZ(_p)
 
   V = domain(Vinq)
@@ -976,14 +1039,14 @@ function _subgroups_orbit_representatives_and_stabilizers_elementary(
   # and we return it, or if order(H0) > ord, there are no subgroups as wanted
   if order(H0) >= ord
     order(H0) > ord && return res
-    push!(res, (H0inq, G))
+    push!(res, (H0inq, id_hom(G)))
     return res
   end
   # Now the groups we look for should strictly contain H0.
   # If ord == order(V), then there is only V satisfying the given
   # conditions, and V is stabilized by the all G
   if ord == order(V)
-    push!(res, (Vinq, G))
+    push!(res, (Vinq, id_hom(G)))
     return res
   end
 
@@ -1018,11 +1081,11 @@ function _subgroups_orbit_representatives_and_stabilizers_elementary(
       # We keep only f-stable subspaces
       is_invariant(f, orbqinq) || continue
       stabq_gen = elem_type(G)[GtoOVmodH0\(inc_stab(s)) for s in gens(stab)]
-      stabq, _ = sub(G, union!(stabq_gen, gens(satV)))
+      stabq, jj = sub(G, union!(stabq_gen, gens(satV)))
       # Stabilizers should preserve the actual subspaces, by definition. so if we
       # have lifted everything properly, this should hold..
       @hassert :ZZLatWithIsom 1 is_invariant(stabq, orbqinq)
-      push!(res, (orbqinq, stabq))
+      push!(res, (orbqinq, jj))
     end
     @vprint :ZZLatWithIsom 5 " done \r                                     \r"
     return res
@@ -1072,147 +1135,16 @@ function _subgroups_orbit_representatives_and_stabilizers_elementary(
     is_invariant(f, orbqinq) || continue
 
     stabq_gen = elem_type(G)[GtoMGp\(s) for s in gens(stab)]
-    stabq, _ = sub(G, union!(stabq_gen, gens(satV)))
+    stabq, jj = sub(G, union!(stabq_gen, gens(satV)))
     # Stabilizers should preserve the actual subspaces, by definition. so if we
     # have lifted everything properly, this should hold..
     @hassert :ZZLatWithIsom 1 is_invariant(stabq, orbqinq)
-    push!(res, (orbqinq, stabq))
+    push!(res, (orbqinq, jj))
   end
   @vprint :ZZLatWithIsom 5 "   done \n"
   return res
 end
 
-# Compute `O`-orbits of `f`-stable submodules of `ker(mu(f))` which are
-# isometric, as torsion quadratic modules, to `H`. It also computes the
-# stabilizers in `O` of such subgroups. If `H` is not given, then return
-# orbits of stable submodules of order `ordH`.
-#
-# The outputs are given by embeddings of such submodules in `q`.
-#
-# The code splits the computation into primary parts since they are orthogonal
-# to each others.
-#
-# TODO: Remove this old bad code within the new infrastructure
-function _classes_isomorphic_subgroups(
-    q::TorQuadModule,
-    O::AutomorphismGroup{TorQuadModule},
-    f::Union{TorQuadModuleMap, AutomorphismGroupElem{TorQuadModule}} = id_hom(domain(O));
-    H::Union{Nothing, TorQuadModule}=nothing,
-    ordH::Union{Nothing, IntegerUnion}=nothing,
-    mu::PolyRingElem=zero(Hecke.Globals.Qx),
-  )
-  res = Tuple{TorQuadModuleMap, AutomorphismGroup{TorQuadModule}}[]
-
-  if isnothing(H)
-    @assert !isnothing(ordH)
-    @assert ordH > 0
-  else
-    ordH = order(H)
-  end
-
-  !is_divisible_by(order(q), ordH) && return res
-
-  if ordH == 1
-    _, j = sub(q, elem_type(q)[])
-    push!(res, (j, O))
-    return res
-  end
-
-  # Trivial case: we look for subgroups in a given primary part of q
-  ok, e, p = is_prime_power_with_data(ordH)
-  if ok
-    if (e == 1) || (!isnothing(H) && is_elementary(H, p))
-      _, Vinq = _get_V(f, mu, p)
-      sors = _subgroups_orbit_representatives_and_stabilizers_elementary(Vinq, O, ordH, p, f)
-    else
-      T, Tinq = kernel(evaluate(mu, f))
-      _, VinT = primary_part(T, p)
-      Vinq = compose(VinT, Tinq)
-      sors = _subgroups_orbit_representatives_and_stabilizers(Vinq, O, ordH, f)
-    end
-    if !isnothing(H)
-      filter!(d -> is_isometric_with_isometry(domain(d[1]), H)[1], sors)
-    end
-    return sors
-  end
-
-  # We inspect each primary part of q and look for orbit representatives and
-  # stabilizers of isomorphic subgroups which will be isometric to the given
-  # primary part of H.
-  #
-  # First, we cut q as an orthogonal direct sum of its primary parts
-  pds = sort!(prime_divisors(order(q)))
-  blocks = TorQuadModuleMap[primary_part(q, pds[1])[2]]
-  ni = Int[ngens(domain(blocks[1]))]
-  for i in 2:length(pds)
-    _f = blocks[end]
-    _, j = has_complement(_f)
-    _T = domain(j)
-    __f = primary_part(_T, pds[i])[2]
-    push!(blocks, compose(__f, j))
-    push!(ni, ngens(domain(__f)))
-  end
-  D, inj, proj = Hecke._biproduct(domain.(blocks))
-  phi = hom(D, q, TorQuadModuleElem[sum([blocks[i](proj[i](a)) for i in 1:length(pds)]) for a in gens(D)])
-  @hassert :ZZLatWithIsom 1 is_isometry(phi)
-
-  list_can = Vector{Tuple{TorQuadModuleMap, AutomorphismGroup{TorQuadModule}}}[]
-  # We collect the possible subgroups for each primary part, with the stabilizer
-  for i in 1:length(pds)
-    p = pds[i]
-    qpinq = blocks[i]
-    qp = domain(qpinq)
-    ordHp = p^valuation(ordH, p)
-    if !isnothing(H)
-      T, _ = primary_part(H, p)
-    end
-    Oqp, _ = restrict_automorphism_group(O, qpinq; check=false)
-    fqp = restrict_endomorphism(f, qpinq; check=false)
-    if (ordHp == p) || (is_elementary(qp, p)) || (!isnothing(H) && is_elementary(T, p))
-      _, j = _get_V(fqp, mu, p)
-      sors = _subgroups_orbit_representatives_and_stabilizers_elementary(j, Oqp, ordHp, p, fqp)
-    else
-      _, Tpinqp = kernel(evaluate(mu, fqp))
-      sors = _subgroups_orbit_representatives_and_stabilizers(Tpinqp, Oqp, ordHp, fqp)
-    end
-    if !isnothing(H)
-      filter!(d -> is_isometric_with_isometry(domain(d[1]), T)[1], sors)
-    end
-    is_empty(sors) && return res
-    push!(list_can, sors)
-  end
-
-  # We gather together: we do a big cartesian product, and we remember to
-  # reconstruct the stabilizer. Since primary parts do not talk to each other,
-  # we concatenate generators on an orthogonal direct sum of q into its primary
-  # parts (as we do for computations of orthogonal groups in the non split
-  # degenerate case)
-  for lis in Hecke.cartesian_product_iterator(list_can)
-    embs = TorQuadModuleMap[l[1] for l in lis]
-    embs = TorQuadModuleMap[hom(domain(embs[i]), q, TorQuadModuleElem[blocks[i](domain(blocks[i])(lift(embs[i](a)))) for a in gens(domain(embs[i]))]) for i in 1:length(lis)]
-    H2, _, _proj = Hecke._biproduct(domain.(embs))
-    _, H2inq = sub(q, elem_type(q)[sum([embs[i](_proj[i](g)) for i in 1:length(lis)]) for g in gens(H2)])
-    stabs = AutomorphismGroup{TorQuadModule}[l[2] for l in lis]
-    genestab = ZZMatrix[]
-
-    for i in 1:length(ni)
-      nb = sum(ni[1:i-1])
-      na = sum(ni[(i+1):end])
-      Inb = identity_matrix(ZZ, nb)
-      Ina = identity_matrix(ZZ, na)
-      append!(genestab, ZZMatrix[block_diagonal_matrix([Inb, matrix(f), Ina]) for f in gens(stabs[i])])
-    end
-
-    genestabD = TorQuadModuleMap[hom(D, D, g) for g in genestab]
-    genestas = ZZMatrix[matrix(compose(compose(inv(phi), g), phi)) for g in genestabD]
-    stab, _ = intersect(O, Oscar._orthogonal_group(q, unique(genestas); check=false))
-    @hassert :ZZLatWithIsom is_invariant(stab, H2inq)
-    push!(res, (H2inq, stab))
-  end
-
-  return res
-end
-
 ###############################################################################
 #
 # Primitive extensions
@@ -2025,7 +1957,7 @@ function admissible_equivariant_primitive_extensions(
     union!(geneA, geneB)
 
     # We compute the overlattice in this context
-    C2, fC2, _ = _overlattice(qAinD, qBinD, isometry(A), isometry(B))
+    C2, fC2, _ = _equivariant_overlattice_with_graph(qAinD, qBinD, isometry(A), isometry(B))
     C2fC2 = integer_lattice_with_isometry(C2, fC2; ambient_representation=false, check)
 
     # If not of the good type, we discard it
@@ -2164,7 +2096,7 @@ function admissible_equivariant_primitive_extensions(
 
       # We compute the overlattice in this context, keeping track whether we
       # work in a fixed ambient quadratic space
-      C2, fC2, extinD = _overlattice(phig, SAinD, SBinD, isometry(A), isometry(B))
+      C2, fC2, extinD = _equivariant_overlattice_with_graph(phig, SAinD, SBinD, isometry(A), isometry(B))
       C2fC2 = integer_lattice_with_isometry(C2, fC2; ambient_representation=false, check)
 
       # This is the type requirement: somehow, we want `(C2, fC2)` to be a "q-th root" of `(C, fC)`.
@@ -2436,7 +2368,7 @@ function _glue_stabilizers(
   actM = hom(stabM, OHM, elem_type(OHM)[OHM(restrict_automorphism(x, HMinqM; check=false)) for x in gens(stabM)])
   actN = hom(stabN, OHN, elem_type(OHN)[OHN(restrict_automorphism(x, HNinqN; check=false)) for x in gens(stabN)])
 
-  _, _, graph = _overlattice(phi, HMinD, HNinD, isometry(M), isometry(N))
+  _, _, graph = _equivariant_overlattice_with_graph(phi, HMinD, HNinD, isometry(M), isometry(N))
   disc, _stab = _glue_stabilizers(phi, actM, actN, OqMinOD, OqNinOD, graph)
   qL, fqL = discriminant_group(L)
   OqL = orthogonal_group(qL)
@@ -2462,3 +2394,167 @@ function _glue_stabilizers(
   Nh = integer_lattice_with_isometry(N)
   return _glue_stabilizers(Lf, Mg, Nh; kwargs...)
 end
+
+###############################################################################
+#
+#  To be removed at some point
+#
+###############################################################################
+
+# Compute `O`-orbits of `f`-stable submodules of `ker(mu(f))` which are
+# isometric, as torsion quadratic modules, to `H`. It also computes the
+# stabilizers in `O` of such subgroups. If `H` is not given, then return
+# orbits of stable submodules of order `ordH`.
+#
+# The outputs are given by embeddings of such submodules in `q`.
+#
+# The code splits the computation into primary parts since they are orthogonal
+# to each others.
+function _classes_isomorphic_subgroups(
+    q::TorQuadModule,
+    O::AutomorphismGroup{TorQuadModule},
+    f::Union{TorQuadModuleMap, AutomorphismGroupElem{TorQuadModule}} = id_hom(domain(O));
+    H::Union{Nothing, TorQuadModule}=nothing,
+    ordH::Union{Nothing, IntegerUnion}=nothing,
+    mu::PolyRingElem=zero(Hecke.Globals.Qx),
+  )
+  res = Tuple{TorQuadModuleMap, AutomorphismGroup{TorQuadModule}}[]
+
+  if isnothing(H)
+    @assert !isnothing(ordH)
+    @assert ordH > 0
+  else
+    ordH = order(H)
+  end
+
+  !is_divisible_by(order(q), ordH) && return res
+
+  if ordH == 1
+    _, j = sub(q, elem_type(q)[])
+    push!(res, (j, O))
+    return res
+  end
+
+  # Trivial case: we look for subgroups in a given primary part of q
+  ok, e, p = is_prime_power_with_data(ordH)
+  if ok
+    if (e == 1) || (!isnothing(H) && is_elementary(H, p))
+      _, Vinq = _get_V(f, mu, p)
+      sors = _subgroups_orbit_representatives_and_stabilizers_elementary(Vinq, O, ordH, p, f)
+    else
+      T, Tinq = kernel(evaluate(mu, f))
+      _, VinT = primary_part(T, p)
+      Vinq = compose(VinT, Tinq)
+      sors = _subgroups_orbit_representatives_and_stabilizers(Vinq, O, ordH, f)
+    end
+    if !isnothing(H)
+      filter!(d -> is_isometric_with_isometry(domain(d[1]), H)[1], sors)
+    end
+    return sors
+  end
+
+  # We inspect each primary part of q and look for orbit representatives and
+  # stabilizers of isomorphic subgroups which will be isometric to the given
+  # primary part of H.
+  #
+  # First, we cut q as an orthogonal direct sum of its primary parts
+  pds = sort!(prime_divisors(order(q)))
+  blocks = TorQuadModuleMap[primary_part(q, pds[1])[2]]
+  ni = Int[ngens(domain(blocks[1]))]
+  for i in 2:length(pds)
+    _f = blocks[end]
+    _, j = has_complement(_f)
+    _T = domain(j)
+    __f = primary_part(_T, pds[i])[2]
+    push!(blocks, compose(__f, j))
+    push!(ni, ngens(domain(__f)))
+  end
+  D, inj, proj = Hecke._biproduct(domain.(blocks))
+  phi = hom(D, q, TorQuadModuleElem[sum([blocks[i](proj[i](a)) for i in 1:length(pds)]) for a in gens(D)])
+  @hassert :ZZLatWithIsom 1 is_isometry(phi)
+
+  list_can = Vector{Tuple{TorQuadModuleMap, AutomorphismGroup{TorQuadModule}}}[]
+  # We collect the possible subgroups for each primary part, with the stabilizer
+  for i in 1:length(pds)
+    p = pds[i]
+    qpinq = blocks[i]
+    qp = domain(qpinq)
+    ordHp = p^valuation(ordH, p)
+    if !isnothing(H)
+      T, _ = primary_part(H, p)
+    end
+    Oqp, _ = restrict_automorphism_group(O, qpinq; check=false)
+    fqp = restrict_endomorphism(f, qpinq; check=false)
+    if (ordHp == p) || (is_elementary(qp, p)) || (!isnothing(H) && is_elementary(T, p))
+      _, j = _get_V(fqp, mu, p)
+      sors = _subgroups_orbit_representatives_and_stabilizers_elementary(j, Oqp, ordHp, p, fqp)
+    else
+      _, Tpinqp = kernel(evaluate(mu, fqp))
+      sors = _subgroups_orbit_representatives_and_stabilizers(Tpinqp, Oqp, ordHp, fqp)
+    end
+    if !isnothing(H)
+      filter!(d -> is_isometric_with_isometry(domain(d[1]), T)[1], sors)
+    end
+    is_empty(sors) && return res
+    push!(list_can, sors)
+  end
+
+  # We gather together: we do a big cartesian product, and we remember to
+  # reconstruct the stabilizer. Since primary parts do not talk to each other,
+  # we concatenate generators on an orthogonal direct sum of q into its primary
+  # parts (as we do for computations of orthogonal groups in the non split
+  # degenerate case)
+  for lis in Hecke.cartesian_product_iterator(list_can)
+    embs = TorQuadModuleMap[l[1] for l in lis]
+    embs = TorQuadModuleMap[hom(domain(embs[i]), q, TorQuadModuleElem[blocks[i](domain(blocks[i])(lift(embs[i](a)))) for a in gens(domain(embs[i]))]) for i in 1:length(lis)]
+    H2, _, _proj = Hecke._biproduct(domain.(embs))
+    _, H2inq = sub(q, elem_type(q)[sum([embs[i](_proj[i](g)) for i in 1:length(lis)]) for g in gens(H2)])
+    stabs = AutomorphismGroup{TorQuadModule}[l[2] for l in lis]
+    genestab = ZZMatrix[]
+
+    for i in 1:length(ni)
+      nb = sum(ni[1:i-1])
+      na = sum(ni[(i+1):end])
+      Inb = identity_matrix(ZZ, nb)
+      Ina = identity_matrix(ZZ, na)
+      append!(genestab, ZZMatrix[block_diagonal_matrix([Inb, matrix(f), Ina]) for f in gens(stabs[i])])
+    end
+
+    genestabD = TorQuadModuleMap[hom(D, D, g) for g in genestab]
+    genestas = ZZMatrix[matrix(compose(compose(inv(phi), g), phi)) for g in genestabD]
+    stab, _ = intersect(O, Oscar._orthogonal_group(q, unique(genestas); check=false))
+    @hassert :ZZLatWithIsom is_invariant(stab, H2inq)
+    push!(res, (H2inq, stab))
+  end
+
+  return res
+end
+
+
+function _change_to_bilinear_module(
+    qM::TorQuadModule,
+    GM::AutomorphismGroup{TorQuadModule},
+    fqM::TorQuadModuleMap,
+  )
+  qM = Hecke._as_finite_bilinear_module(qM) # TODO: to be changed
+  OqM = orthogonal_group(qM)
+  GM, _ = sub(OqM, elem_type(OqM)[OqM(matrix(g); check=false) for g in gens(GM)])
+  fqM = hom(qM, qM, matrix(fqM))
+  return qM, OqM, GM, fqM
+end
+
+
+function _possible_glue_orders(
+    qM::TorQuadModule,
+    qN::TorQuadModule,
+  )
+  _gcd = ZZ(1)
+  snM = reverse!(elementary_divisors(qM))
+  snN = reverse!(elementary_divisors(qN))
+  k = min(length(snM), length(snN))
+  for i in 1:k
+    mul!(_gcd, _gcd, gcd(snM[i], snN[i]))
+  end
+  pos_ord = divisors(_gcd)
+  return pos_ord
+end
diff --git a/src/NumberTheory/QuadFormAndIsom/gluing_factory.jl b/src/NumberTheory/QuadFormAndIsom/gluing_factory.jl
new file mode 100644
index 000000000000..3f7fadeb8c55
--- /dev/null
+++ b/src/NumberTheory/QuadFormAndIsom/gluing_factory.jl
@@ -0,0 +1,1637 @@
+###############################################################################
+#
+#  Helpers
+#
+###############################################################################
+
+### Gluing factory
+
+@doc raw"""
+    ambient_modules(Fac::ZZLatGluingFactory) -> TorQuadModule, TorQuadModule
+
+Return the two initial torsion quadratic modules from the factory `Fac`.
+"""
+ambient_modules(Fac::ZZLatGluingFactory) = Fac.ambient_modules
+
+@doc raw"""
+    local_classifying_groups(
+      Fac::ZZLatGluingFactory,
+    ) -> AutomorphismGroup{TorQuadModule}, AutomorphismGroup{TorQuadModule}
+
+Return the groups used to do the computations of local gluings in the factory
+`Fac`.
+"""
+local_classifying_groups(Fac::ZZLatGluingFactory) = Fac.local_classifying_groups
+
+@doc raw"""
+    conditions_modules(
+      Fac::ZZLatGluingFactory,
+    ) -> TorQuadModuleMap, TorQuadModuleMap
+
+Return the embedding of the conditions modules in the ambient modules of
+the factory `Fac`, computed via the conditions initiliazing `Fac.
+"""
+conditions_modules(Fac::ZZLatGluingFactory) = Fac.conditions_modules
+
+@doc raw"""
+    overlattice_parity(Fac::ZZLatGluingFactory) -> Symbol
+
+Return the parity type of a primitive extension as setup in the construction
+of the factory `Fac`.
+
+It can be :even for even extensions, :odd for odd extensions and :both for
+any integral extensions.
+"""
+overlattice_parity(Fac::ZZLatGluingFactory) = Fac.par
+
+@doc raw"""
+    possible_glue_order(Fac::ZZLatGluingFactor) -> Set{ZZRingElem}
+
+Return the set of possible order for the glue groups, computed in the
+initialization of the factory `Fac`.
+"""
+possible_glue_order(Fac::ZZLatGluingFactory) = Fac.glue_order
+
+@doc raw"""
+    test_overlattice(Fac::ZZLatGluingFactory, L::ZZLat) -> Bool
+
+Return whether the given lattice satisfies the output conditions of
+the factory `Fac`.
+
+This means that:
+- the parity of `L` is good,
+- the genus of `L` is among the allowed ones (if any),
+- the discriminant group of `L` is among the allowed ones (if any).
+"""
+function test_overlattice(Fac::ZZLatGluingFactory, L::ZZLat)
+  par = overlattice_parity(Fac)
+  if par === :even && !is_even(L)
+    return false
+  elseif par === :odd && is_even(L)
+    return false
+  end
+  if isdefined(Fac, :genus_over)
+    genus(L) in Fac.genus_over && return true
+  end
+  if isdefined(Fac, :form_over)
+    M = gram_matrix_quadratic(first(normal_form(discriminant_group(L))))
+    return M in Fac.form_over
+  end
+  return true
+end
+
+@doc raw"""
+    is_trivial(Fac::ZZLatGluingFactory) -> Bool
+
+Return whether the output conditions of the factory `Fac` are possible.
+"""
+function is_trivial(Fac::ZZLatGluingFactory)
+  isempty(possible_glue_order(Fac)) && return true
+  if isdefined(Fac, :genus_over)
+    return isempty(Fac.genus_over)
+  elseif isdefined(Fac, :form_over)
+    return isempty(Fac.form_over)
+  elseif isdefined(Fac, :glue_elementary_divisors)
+    return isempty(Fac.glue_elementary_divisors)
+  end
+  return false
+end
+
+### Gluing
+
+function _data(x::ZZLatGluing)
+  return x.glue_map, x.inv_glue_map, x.glue_group_left, x.stabilizer_left, x.glue_group_right, x.stabilizer_right
+end
+
+# The type ZZLatGluing is symmetrical, we just reverse the order of the fields
+# to see it from right to left
+function Oscar.inv(x::ZZLatGluing)
+  phi, iphi, i1, j1, i2, j2 = _data(x)
+  return ZZLatGluing(iphi, phi, i2, j2, i1, j1)
+end
+
+function glue_map(x::ZZLatGluing, i::Int)
+  if isone(i)
+    return x.glue_map
+  elseif i == 2
+    return x.inv_glue_map
+  end
+end
+
+glue_maps(x::ZZLatGluing) = (glue_map(x, 1), glue_map(x, 2))
+
+function glue_group(x::ZZLatGluing, i::Int)
+  if isone(i)
+    return x.glue_group_left
+  elseif i == 2
+    return x.glue_group_right
+  end
+end
+
+glue_groups(x::ZZLatGluing) = (glue_group(x, 1), glue_group(x, 2))
+
+function stabilizer_glue_group(x::ZZLatGluing, i::Int)
+  if isone(i)
+    return x.stabilizer_left
+  elseif i == 2
+    return x.stabilizer_right
+  end
+end
+
+### Gluing ambient
+
+function _gluing_ambient(
+    q1::TorQuadModule,
+    q2::TorQuadModule,
+    par::Symbol,
+  )
+  as_bilinear_module = par === :even ? false : true
+  x = _direct_sum_with_embeddings_orthogonal_groups(q1, q2; as_bilinear_module)
+  return ZZLatGluingAmbient(x...)
+end
+
+
+### More on genera
+
+# Return all genera with given form, signature and parity conditions
+# There are at most 2 of them
+# par is either :even, or :odd, or :both (actually anything other then
+# :odd or :even would do)
+function _integer_genera(
+  q::TorQuadModule,
+  sign::Tuple{Int, Int},
+  par::Symbol,
+)
+  p, n = sign
+  GKs = Set{ZZGenus}()
+  if par != :odd
+    ok, _G = is_genus_with_genus(q, (p, n); parity=2)
+    ok && push!(GKs, _G)
+  end
+
+  if par != :even
+    ok, _G = is_genus_with_genus(q, (p, n); parity=1)
+    ok && push!(GKs, _G)
+  end
+  return GKs
+end
+
+function __representatives(g::ZZGenus)
+  return get_attribute(g, :representatives, representatives(g))
+end
+
+### Handle parity
+
+@doc raw"""
+    orthogonal_group_bilinear(T::TorQuadModule) -> AutomorphismGroup{TorQuadModule}
+
+Return the orthogonal group for the bilinear form on ``T``, i.e. seeing ``T``
+as a torsion bilinear form. If the modulus of the quadratic form on ``T`` is
+equal to the modulus of its bilinear form, then this is the same as calling
+`orthogonal_group(T)`.
+"""
+@attr AutomorphismGroup{TorQuadModule} function orthogonal_group_bilinear(T::TorQuadModule)
+  return __orthogonal_group(T; as_bilinear_module=true)
+end
+
+function __as_finite_bilinear_module(
+  q::TorQuadModule,
+)
+  n = modulus_bilinear_form(q)
+  if n == modulus_quadratic_form(q)
+    return q
+  end
+  
+  qb =  torsion_quadratic_module(cover(q), relations(q); modulus=n, modulus_qf=n, gens=lift.(gens(q)))
+  return qb
+end
+
+function __orthogonal_group(
+  q::TorQuadModule;
+  as_bilinear_module::Bool=false,
+)
+  if !as_bilinear_module || modulus_bilinear_form(q) == modulus_quadratic_form(q)
+    return orthogonal_group(q)
+  end
+  qb = __as_finite_bilinear_module(q)
+  return _orthogonal_group(q, _orthogonal_group_gens(qb); check=false)
+end
+
+### To be moved to Hecke eventually
+
+# Return the elementary divisors of the largest finite abelian group
+# which embeds in both G1 and G2
+function _maximal_common_subgroup_snf(
+  G1::FinGenAbGroup,
+  G2::FinGenAbGroup,
+)
+  @assert is_finite(G1) && is_finite(G2)
+  s1 = elementary_divisors(G1)
+  s2 = elementary_divisors(G2)
+  s = Array{ZZRingElem}(undef, min(length(s1), length(s2)))
+  for i in 0:length(s)-1
+    s[end-i] = gcd(s1[end-i], s2[end-i])
+  end
+  return s
+end
+
+# Given the elementary divisors elG of a finite abelian group G, return the
+# set of all possible lists of valuations for the elementary divisors of a
+# subgroup of G whose p-Sylow subgroup has order p^v. The lists in output
+# contain only positive valuations and they are sorted in decreasing order.
+function _psubgroups_types(
+  elG::Vector{ZZRingElem},
+  p::ZZRingElem,
+  v::Int,
+)
+  @assert v > 0
+  testv = isequal(v)∘sum
+  p_subtypes = Set{Vector{Int}}()
+  if isempty(elG)
+    return p_subtypes
+  end
+  Gtype = Int[valuation(a, p) for a in elG]
+  reverse!(Gtype)
+  filter!(!=(0), Gtype)
+  types = Hecke._subpartitions(Gtype)
+  for t in types
+    testv(t) || continue
+    filter!(!=(0), t)
+    push!(p_subtypes, t)
+  end
+  return p_subtypes
+end
+
+@doc raw"""
+    torsion_subgroup(
+      G::FinGenAbGroup,
+      n::IntegerUnion,
+      add_to_lattice:Bool = true,
+    )
+
+Return the subgroup of ``G`` consisting of elements of order dividing ``n``.
+"""
+function torsion_subgroup(
+  G::FinGenAbGroup,
+  n::IntegerUnion,
+  add_to_lattice::Bool = true,
+)
+  f = FinGenAbGroupHom(G, G, ZZ(n)*identity_matrix(ZZ, ngens(G)))
+  return kernel(f, add_to_lattice)
+end
+
+@doc raw"""
+    torsion_subgroup(
+      T::TorQuadModule,
+      n::IntegerUnion,
+    )
+
+Return the subgroup of ``T`` consisting of elements of order dividing ``n``.
+"""
+function torsion_subgroup(
+  T::TorQuadModule,
+  n::IntegerUnion,
+)
+  Sab, j = torsion_subgroup(abelian_group(T), n)
+  return sub(T, TorQuadModuleElem[T(j(s)) for s in gens(Sab)])
+end
+
+###############################################################################
+#
+#  Init functions
+#
+###############################################################################
+
+# Initialize the gluing factory with the conditions from the problem.
+# By construction `Fac` should already know about the ambient modules and the
+# wanted parity for the primitive extensions. If the local classifying groups
+# have not been set, the function initialize them to be the orthogonal groups
+# (maybe as finite bilinear modules, depending on the value of Fac.par) of the
+# ambient modules.
+function init_gluing_factory!(
+  Fac::ZZLatGluingFactory;
+  left_glue_annihilator::Union{Nothing, TorQuadModuleMap}=nothing,
+  right_glue_annihilator::Union{Nothing, TorQuadModuleMap}=nothing,
+  glue_exponent::Union{Nothing, IntegerUnion}=nothing,
+  glue_order::AbstractVector{T}=Int[],
+  glue_elementary_divisors::Vector{Vector{ZZRingElem}}=Vector{ZZRingElem}[],
+  form_over::Vector{TorQuadModule}=TorQuadModule[],
+  genus_over::Vector{ZZGenus}=ZZGenus[],
+) where T <: IntegerUnion
+  init_gluing_conditions!(
+    Fac,
+    left_glue_annihilator,
+    right_glue_annihilator,
+    glue_exponent,
+    glue_order,
+    glue_elementary_divisors,
+    form_over,
+  )
+  assert_has_local_classifying_groups!(Fac)
+  return nothing
+end
+
+# Using the glue annihilators and glue_exponent, the function sets up the
+# "conditions modules", i.e. where to look for glue groups. From this, the
+# fonction filters through the remaining initial conditions to determine the
+# actual restrictions for the gluing computed in the factory. This list of
+# actual conditions consists of:
+# - a set of genera for the primitive extensions (optional: only if
+#   genus_over is not empty),
+# - a set of Gram matrix for the normal forms of the discriminant groups
+#   of the primitive extensions (optional: only if form_over is not empty)
+# - a set of elementary divisors for the glue groups (optional: only if
+#   glue_elementary_divisors is not empty)
+# - a set of orders for the glue groups.
+#
+# If genus_over (resp. form_over, glue_elementary_divisors) in input is not
+# empty but it is after filtering (depending on the conditions modules and the
+# input list glue_order), Fac.genus_over (resp. Fac.form_over,
+# Fac.glue_elementary_divisors) is set to be the empty set and the following
+# functions will abort (because of impossible conditions, Fac is "trivial").
+#
+# For Fac.glue_order, several things can happen:
+# - if the initial conditions genus_over, form_over or glue_elementary_divisors
+#   are not empty, then the function infers lists of possible orders from the
+#   non-empty ones and intersect them. If moreover glue_order is not empty,
+#   the previous list of infered orders is intersected with glue_order.
+# - otherwise, if glue_order is non-empty, the functions filters which orders
+#   could actually work in the context set up by the conditions modules.
+# - otherwise, Fac.glue_order consists of all the divisors of the order of the
+#   maximal abelian group embedding in both of the conditions modules.
+# In the two first cases, if the final list of orders infered from the initial
+# conditions is empty, Fac.glue_order is the empty set and the rest of the
+# algorithm will abort as well.
+function init_gluing_conditions!(
+  Fac::ZZLatGluingFactory,
+  left_glue_annihilator::Union{Nothing, TorQuadModuleMap}=nothing,
+  right_glue_annihilator::Union{Nothing, TorQuadModuleMap}=nothing,
+  glue_exponent::Union{Nothing, IntegerUnion}=nothing,
+  glue_order::AbstractVector{T}=ZZRingElem[],
+  glue_elementary_divisors::Vector{Vector{ZZRingElem}}=Vector{ZZRingElem}[],
+  form_over::Vector{TorQuadModule}=TorQuadModule[],
+  genus_over::Vector{ZZGenus}=ZZGenus[],
+) where T <: IntegerUnion
+
+  # Remember if some conditions were initially set
+  ignore_genus_over = isempty(genus_over)
+  ignore_form_over = isempty(form_over)
+  ignore_elem_divs = isempty(glue_elementary_divisors)
+  empty_order_cond = isempty(glue_order)
+  ignore_glue_order = ignore_genus_over && ignore_form_over && ignore_elem_divs && empty_order_cond
+
+  qM, qN = ambient_modules(Fac)
+
+  # Left glue group annihilated by the endomorphism left_glue_annihilator
+  if !isnothing(left_glue_annihilator)
+    @assert domain(left_glue_annihilator) === codomain(left_glue_annihilator) === qM
+    VM, VMinqM = kernel(left_glue_annihilator)
+  else
+    VMinqM = id_hom(qM)
+    VM = qM
+  end
+
+  # Right glue group annihilated by the endomorphism right_glue_annihilator
+  if !isnothing(right_glue_annihilator)
+    @assert domain(right_glue_annihilator) === codomain(right_glue_annihilator) === qN
+    VN, VNinqN = kernel(right_glue_annihilator)
+  else
+    VNinqN = id_hom(qN)
+    VN = qN
+  end
+
+  # If glue_exponent is a positive integer, the exponent of the glue groups
+  # should divide it
+  if !isnothing(glue_exponent) && glue_exponent > 0
+    filter!(Base.Fix1(is_divisible_by, glue_exponent)∘last, elementary_divisors)
+    VM, _VMinVM = torsion_subgroup(domain(VMinqM), e)
+    VMinqM = compose(_VMinVM, VMinqM)
+    VN, _VNinVN = torsion_subgroup(domain(VNinqN), e)
+    VNinqN = compose(_VNinVN, VNinqN)
+  end
+
+  # This fixes the conditions modules where we do everything
+  Fac.conditions_modules = (VMinqM, VNinqN)
+
+  elG = _maximal_common_subgroup_snf(abelian_group(VM), abelian_group(VN))
+  pds = isempty(elG) ? ZZRingElem[] : prime_divisors(last(elG))
+  ged = Dict{ZZRingElem, Dict{Vector{Int}, Vector{ZZLatGluing}}}(p => Dict{Vector{Int}, Vector{ZZLatGluing}}() for p in pds)
+  Fac.glue_group_parent_snf = elG
+  Fac.primes_of_interest = Set(pds)
+  Fac.local_gluings_primary = ged
+
+  _genus_over = Set(genus_over)
+
+  k1 = prod(elG; init=ZZ(1))
+  _glue_order = Set(filter!(>(0), glue_order))
+  for o in _glue_order
+    !is_divisible_by(k1, o) && delete!(_glue_order, o)
+  end
+
+  _glue_elem_divs = Set{Vector{ZZRingElem}}()
+  for v in glue_elementary_divisors
+    if isempty(v)
+      push!(_glue_elem_divs, v)
+      continue
+    end
+    any(<=(0), v) && continue
+    length(v) > length(elG) && continue
+    flag = is_divisible_by(elG[end-length(v)+1], first(v))
+    !flag && continue
+    for i in 1:length(v)-1
+      flag &= is_divisible_by(v[i+1], v[i])
+      flag &= is_divisible_by(elG[end-i+1], v[end-i+1])
+      !flag && break
+    end
+    if flag
+      push!(_glue_elem_divs, v)
+    end
+  end
+
+  _form_over = Set{QQMatrix}()
+  for q in form_over
+    !isone(modulus_bilinear_form(q)) && continue
+    if Fac.par === :even
+      modulus_quadratic_form(q) == 2 || continue
+    elseif Fac.par === :odd
+      modulus_quadratic_form(q) == 1 || continue
+    else
+      modulus_quadratic_form(q) in [1, 2] || continue
+    end
+    push!(_form_over, gram_matrix_quadratic(first(normal_form(q))))
+  end
+
+
+  k2 = order(qM)*order(qN)
+  _tmp = Set{ZZRingElem}()
+  for G in _genus_over
+    bool, _k = divides(k2, numerator(abs(det(G))))
+    if !bool
+      delete!(_genus_over, G)
+      continue
+    end
+    bool, _k = is_square_with_sqrt(_k)
+    if !bool
+      delete!(_genus_over, G)
+      continue
+    end
+    if empty_order_cond || (_k in _glue_order)
+      push!(_tmp, _k)
+    else
+      delete!(_genus_over, G)
+      continue
+    end
+  end
+  for q in _form_over
+    bool, _k = divides(k2, order(q))
+    if !bool
+      delete!(_form_over, q)
+      continue
+    end
+    bool, _k = is_square_with_sqrt(_k)
+    if !bool
+      delete!(_form_over, q)
+      continue
+    end
+    if empty_order_cond || (_k in _glue_order)
+      push!(_tmp, _k)
+    else
+      delete!(_genus_over, G)
+      continue
+    end
+  end
+
+  if !ignore_genus_over || !ignore_form_over
+    if empty_order_cond
+      union!(_glue_order, _tmp)
+    else
+      intersect!(_glue_order, _tmp)
+    end
+  end
+
+  if !ignore_elem_divs
+    _tmp = Set{ZZRingElem}()
+    for v in _glue_elem_divs
+      _k = prod(v; init=ZZ(1))
+      if (ignore_genus_over && ignore_form_over && empty_order_cond) || (_k in _glue_order)
+        push!(_tmp, _k)
+      else
+        delete!(_glue_elem_divs, v)
+        continue
+      end
+    end
+    if ignore_genus_over && ignore_form_over && empty_order_cond
+      union!(_glue_order, _tmp)
+    else
+      intersect!(_glue_order, _tmp)
+    end
+  end
+
+  if !ignore_glue_order
+    Fac.glue_order = _glue_order
+  else
+    Fac.glue_order = Set(prime_divisors(k1))
+  end
+
+  if !ignore_genus_over
+    Fac.genus_over = _genus_over
+  end
+
+  if !ignore_form_over
+    Fac.form_over = _form_over
+  end
+
+  if !ignore_elem_divs
+    Fac.glue_elementary_divisors = _glue_elem_divs
+  end
+  return nothing
+end
+
+function assert_has_local_classifying_groups!(Fac::ZZLatGluingFactory)
+  isdefined(Fac, :local_classifying_groups) && return nothing
+  q1, q2 = ambient_modules(Fac)
+  as_bilinear_module = (Fac.par !== :even)
+  Oq1 = __orthogonal_group(q1; as_bilinear_module)
+  Oq2 = __orthogonal_group(q2; as_bilinear_module)
+  Fac.local_classifying_groups = (Oq1, Oq2)
+  return nothing
+end
+
+###############################################################################
+#
+#  Local gluings
+#
+###############################################################################
+
+function __subgroups_orbit_representatives_and_stabilizers_primary_subtype(
+  Fac::ZZLatGluingFactory,
+  Winq::TorQuadModuleMap,
+  O::AutomorphismGroup{TorQuadModule},
+  p::ZZRingElem,
+  subtype::Vector{Int},
+  i::Int = -1,
+)
+  W = domain(Winq)
+  flag = isdefined(Fac, :Ctx) && (0 < i <= length(Fac.Ctx.modules)) && (Fac.Ctx.modules[i] === codomain(Winq))
+  if flag && haskey(Fac.Ctx.orb_and_stab, (i, p, subtype))
+    subs = Fac.Ctx.orb_and_stab[(i, p, subtype)]
+  elseif first(subtype) == 1
+    W = domain(Winq)
+    _, WWinW = torsion_subgroup(W, p)
+    WWinq = compose(WWinW, Winq)
+    subs = Oscar.__subgroups_orbit_representatives_and_stabilizers_elementary(WWinq, O, p^(length(subtype)), p)
+    if flag
+      Fac.Ctx.orb_and_stab[(i, p, subtype)] = subs
+    end
+  else
+    subs = __subgroups_orbit_representatives_and_stabilizers_primary_subtype(Winq, O, p, subtype)
+    if flag
+      Fac.Ctx.orb_and_stab[(i, p, subtype)] = subs
+    end
+  end
+  return subs
+end
+
+function _local_gluings_primary!(
+  Fac::ZZLatGluingFactory,
+  p::ZZRingElem,
+  subtype::Vector{Int},
+)
+  as_bilinear_module = (Fac.par !== :even)
+  ged = Fac.local_gluings_primary
+  @assert haskey(ged, p)
+  if haskey(ged[p], subtype)
+    return ged[p][subtype]
+  end
+ 
+  loc_glue_p = ZZLatGluing[]
+  flag1, flag2 = false, false
+  i1, i2 = -1, -1
+  if isdefined(Fac, :Ctx)
+    vi = Fac.vertex_identification
+    if vi[1] > 0
+      flag1 = true
+      i1 = vi[1]
+    end
+    if vi[2] > 0
+      flag2 = true
+      i2 = vi[2]
+    end
+  end
+
+  O1, O2 = local_classifying_groups(Fac)
+  V1inq1, V2inq2 = conditions_modules(Fac)
+  V1 = domain(V1inq1)
+  V2 = domain(V2inq2)
+
+  W1, W1inV1 = primary_part(V1, p)
+  W1inq1 = compose(W1inV1, V1inq1)
+  W2, W2inV2 = primary_part(V2, p)
+  W2inq2 = compose(W2inV2, V2inq2)
+
+  subs1 = __subgroups_orbit_representatives_and_stabilizers_primary_subtype(Fac, W1inq1, O1, p, subtype, i1)
+  subs2 = __subgroups_orbit_representatives_and_stabilizers_primary_subtype(Fac, W2inq2, O2, p, subtype, i2)
+
+  for (H1inq1, stab1) in subs1
+    H1 = domain(H1inq1)
+    for (H2inq2, stab2) in subs2
+      H2 = domain(H2inq2)
+      ok, phi = is_anti_isometric_with_anti_isometry(H1, H2; as_bilinear_module)
+      !ok && continue
+      x = ZZLatGluing(phi, inv(phi), H1inq1, stab1, H2inq2, stab2)
+      push!(loc_glue_p, x)
+    end
+  end
+  ged[p][subtype] = loc_glue_p
+  return loc_glue_p
+end
+
+function _trivial_gluing(
+  Fac::ZZLatGluingFactory,
+)
+  q1, q2 = ambient_modules(Fac)
+  s1, s2 = id_hom.(local_classifying_groups(Fac))
+  z1, z1inq1 = sub(q1, TorQuadModuleElem[])
+  z2, z2inq2 = sub(q2, TorQuadModuleElem[])
+  phi = hom(z1, z2, zero_matrix(ZZ, 0, 0))
+  iphi = hom(z2, z1, zero_matrix(ZZ, 0, 0))
+  return ZZLatGluing(phi, iphi, z1inq1, s1, z2inq2, s2)
+end
+
+function _local_glue_maps_ord(Fac::ZZLatGluingFactory, o::ZZRingElem)
+  isone(o) && return ZZLatGluing[_trivial_gluing(Fac)]
+
+  elG = Fac.glue_group_parent_snf
+  _loc = Vector{ZZLatGluing}[]
+  for (p, v) in factor(o)
+    ptypes = _psubgroups_types(elG, p, v)
+    loc_p = ZZLatGluing[]
+    for subtype in ptypes
+      append!(loc_p, _local_gluings_primary!(Fac, p, subtype))
+    end
+    isempty(loc_p) && return ZZLatGluing[]
+    push!(_loc, loc_p)
+  end
+
+  return merge_glue_maps(Fac, _loc)
+end
+
+function _local_glue_maps_eldiv(Fac::ZZLatGluingFactory, v::Vector{ZZRingElem})
+  isempty(v) && return ZZLatGluing[_trivial_gluing(Fac)]
+
+  pds = Fac.primes_of_interest
+  _loc = Vector{ZZLatGluing}[]
+  for p in pds
+    subtype = reverse!(Int[valuation(a, p) for a in v])
+    filter!(!=(0), subtype)
+    isempty(subtype) && continue
+    loc_p = _local_gluings_primary!(Fac, p, subtype)
+    isempty(loc_p) && return ZZLatGluing[]
+    push!(_loc, loc_p)
+  end
+
+  return merge_glue_maps(Fac, _loc)
+end
+
+function merge_glue_maps(
+  Fac::ZZLatGluingFactory,
+  _loc::Vector{Vector{ZZLatGluing}},
+)
+  if isempty(_loc)
+    return ZZLatGluing[]
+  elseif isone(length(_loc))
+    return only(_loc)
+  end
+
+  res = ZZLatGluing[]
+  q1, q2 = ambient_modules(Fac)
+  for x in Hecke.cartesian_product_iterator(_loc; inplace=true)
+    gens1 = TorQuadModuleElem[]
+    gens2 = TorQuadModuleElem[]
+    for y in x
+      H1 = domain(glue_group(y, 1))
+      H2 = domain(glue_group(y, 2))
+      for a in gens(H1)
+        push!(gens1, q1(lift(a)))
+      end
+      for a in gens(H2)
+        push!(gens2, q2(lift(a)))
+      end
+    end
+    H1, H1inq1 = sub(q1, gens1)
+    H2, H2inq2 = sub(q2, gens2)
+    j1 = stabilizer_glue_group(first(x), 1)
+    j2 = stabilizer_glue_group(first(x), 2)
+    for i in 2:length(_loc)
+      _j1 = stabilizer_glue_group(x[i], 1)
+      _j2 = stabilizer_glue_group(x[i], 2)
+      stab1, jj, _ = intersect(domain(j1), domain(_j1))
+      j1 = compose(jj, j1)
+      stab2, jj, _ = intersect(domain(j2), domain(_j2))
+      j2 = compose(jj, j2)
+    end
+    @hassert :ZZLatWithIsom 1 domain(j1) == first(stabilizer(codomain(j1), H1inq1))
+    @hassert :ZZLatWithIsom 1 domain(j2) == first(stabilizer(codomain(j2), H2inq2))
+    phi = hom(H1, H2, gens(H1), gens(H2))
+    @assert is_anti_isometry(phi; as_bilinear_module=(Fac.par !== :even))
+    push!(res, ZZLatGluing(phi, inv(phi), H1inq1, j1, H2inq2, j2))
+  end
+  return res
+end
+
+function local_glue_maps(
+  Fac::ZZLatGluingFactory,
+)
+  res = ZZLatGluing[]
+  if isdefined(Fac, :glue_elementary_divisors)
+    eldiv = Fac.glue_elementary_divisors
+    for v in eldiv
+      append!(res, _local_glue_maps_eldiv(Fac, v))
+    end
+  else
+    l = possible_glue_order(Fac)
+    for o in l
+      append!(res, _local_glue_maps_ord(Fac, o))
+    end
+  end
+  return res
+end
+
+###############################################################################
+#
+#  Operations on gluings
+#
+###############################################################################
+
+function _split_orbit_left_group(
+  x::ZZLatGluing,
+  O1::AutomorphismGroup{TorQuadModule},
+)
+  res = ZZLatGluing[]
+  phi, _, H1inq1, s1, H2inq2, s2 = _data(x)
+  @assert domain(O1) === codomain(H1inq1)
+
+  H1, H2 = domain(H1inq1), domain(H2inq2)
+  q1 = codomain(H1inq1)
+  stab1 = domain(s1)
+  Oq1 = codomain(s1)
+
+  elq1 = elementary_divisors(q1)
+  if isone(order(q1)) || elq1[1] == elq1[end]
+    iso1 = isomorphism(PermGroup, Oq1)
+  else
+    iso1 = id_hom(Oq1)
+  end
+
+  @vprintln :ZZLatWithIsom 1 "Split orbit of left glue group"
+  splits = double_cosets(codomain(iso1), first(iso1(stab1)), first(iso1(O1)))
+  @vprintln :ZZLatWithIsom 1 "    done: $(length(splits))"
+  for _h in splits
+    h = hom(iso1\(representative(_h)))
+    ih = inv(h)
+    _H1, _H1inq1 = sub(q1, elem_type(q1)[h(H1inq1(a)) for a in gens(H1)])
+    _phi = hom(_H1, H2, elem_type(H2)[phi(H1(lift(ih(_H1inq1(a))))) for a in gens(_H1)])
+    _iphi = inv(_phi)
+    _stab1, _ = Oscar._as_subgroup(Oq1, Oscar.GAPWrap.ConjugateSubgroup(GapObj(stab1), GapObj(Oq1(h))))
+    _stab1, _, j1 = intersect(_stab1, O1)
+    push!(res, ZZLatGluing(_phi, _iphi, _H1inq1, j1, H2inq2, s2))
+  end
+  return res
+end
+
+function _split_orbit_right_group(
+  x::ZZLatGluing,
+  O2::AutomorphismGroup{TorQuadModule},
+)
+  return inv.(_split_orbit_left_group(inv(x), O2))
+end
+
+function _pullback_left(
+  x::ZZLatGluing,
+  f::TorQuadModuleMap,
+)
+  phi, iphi, i1, s1, i2, s2 = _data(x)
+  @assert codomain(f) === codomain(i1)
+  invf = inv(f)
+  q0 = domain(f)
+  Oq0 = _orthogonal_group(q0, TorQuadModuleMap[f * hom(g) * invf for g in gens(codomain(s1))]; check=false)
+  _, s0 = sub(Oq0, elem_type(Oq0)[Oq0(f * hom(s1(g)) * invf; check=false) for g in gens(domain(s1))])
+  H0, i0 = sub(q0, TorQuadModuleElem[f\(i1(a)) for a in gens(domain(i1))])
+  H0toH1 = hom(H0, domain(i1), TorQuadModuleElem[i1\(f(i0(a))) for a in gens(H0)])
+  psi = H0toH1 * phi
+  ipsi = inv(psi)
+  return ZZLatGluing(psi, ipsi, i0, s0, i2, s2)
+end
+
+function _pullback_right(
+  x::ZZLatGluing,
+  f::TorQuadModuleMap,
+)
+  return inv(_pullback_left(inv(x), f))
+end
+
+function _all_glue_maps(
+  x::ZZLatGluing,
+)
+  res = ZZLatGluing[]
+  phi, iphi, H1inq1, s1, H2inq2, s2 = _data(x)
+  q1 = codomain(H1inq1)
+  q2 = codomain(H2inq2)
+  H1 = domain(H1inq1)
+  H2 = domain(H2inq2)
+  stab1 = domain(s1)
+  stab2 = domain(s2)
+
+  OH1 = orthogonal_group(H1)
+  imOH1 = elem_type(OH1)[OH1(restrict_automorphism(x, H1inq1); check=false) for x in gens(stab1)]
+  act1 = hom(stab1, OH1, imOH1)
+  im1, _ = image(act1)
+  OH2 = orthogonal_group(H2)
+  imOH2 = elem_type(OH2)[OH2(restrict_automorphism(x, H2inq2); check=false) for x in gens(stab2)]
+  act2 = hom(stab2, OH2, imOH2)
+  im2, _ = image(act2)
+
+  elH2 = elementary_divisors(H2)
+  if isone(order(H2)) || elH2[1] == elH2[end]
+    iso2 = isomorphism(PermGroup, OH2)
+  else
+    iso2 = id_hom(OH2)
+  end
+  OH22 = first(iso2(OH2))
+  im22 = first(iso2(im2))
+  @vprintln :ZZLatWithIsom 1 "Compute right transversals"
+  orb2 = collect(right_transversal(OH22, im22))
+  @vprintln :ZZLatWithIsom 1 "    done: $(length(orb2))"
+  
+  stab1gamma = elem_type(OH2)[OH2(iphi * hom(g) * phi; check=false) for g in gens(im1)]
+  S1, _ = sub(OH2, stab1gamma)
+  S2 = first(iso2(S1))
+  omega = gset(S2, (x,g) -> x*OH22(g), orb2)
+  for _g in orbits(omega)
+    g = hom(iso2\representative(_g))
+    phig = compose(phi, g)
+    iphig = compose(inv(g), iphi)
+    push!(res, ZZLatGluing(phig, iphig, H1inq1, s1, H2inq2, s2))
+  end
+  return res
+end
+
+function _form_over(
+  x::ZZLatGluing,
+)
+  phi = glue_map(x, 1)
+  i1, i2 = glue_groups(x)
+  H1, H2 = domain(i1), domain(i2)
+  q1, q2 = codomain(i1), codomain(i2)
+  D, (j1, j2) = direct_sum(q1, q2; cached=false, as_bilinear_module=true)
+  H1inD = i1 * j1
+  H2inD = i2 * j2
+  _glue = Vector{QQFieldElem}[lift(H1inD(a)) + lift(H2inD(phi(a))) for a in gens(H1)]
+  ext, _ = sub(D, D.(_glue))
+  perp, _ = orthogonal_submodule(D, ext)
+  disc = torsion_quadratic_module(cover(perp), cover(ext))
+  return disc
+end
+
+function _overlattice(
+  x::ZZLatGluing,
+  D::ZZLatGluingAmbient,
+)
+  phi, iphi, i1, s1, i2, s2 = _data(x)
+  H1, H2 = domain(i1), domain(i2)
+  q1, q2 = codomain(i1), codomain(i2)
+  H1inD = i1 * D.j1
+  H2inD = i2 * D.j2
+  O1, O2 = orthogonal_group(H1), orthogonal_group(H2)
+  im1 = elem_type(O1)[O1(restrict_automorphism(g, i1); check=false) for g in gens(domain(s1))]
+  act1 = hom(domain(s1), O1, im1)
+  im2 = elem_type(O2)[O2(restrict_automorphism(g, i2); check=false) for g in gens(domain(s2))]
+  act2 = hom(domain(s2), O2, im2)
+  V, extinD = _overlattice_with_graph(phi, H1inD, H2inD)
+  _disc, _stab = _glue_stabilizers(phi, act1, act2, D.k1, D.k2, extinD)
+  qV = discriminant_group(V)
+  OqV = orthogonal_group(qV)
+  phi2 = hom(qV, _disc, elem_type(_disc)[_disc(lift(x)) for x in gens(qV)])
+  iphi2 = inv(phi2)
+  GV, _ = sub(OqV, elem_type(OqV)[OqV(phi2 * g * iphi2; check=false) for g in _stab])
+  M, T = _as_sublattices(V, relations(q1), relations(q2))
+  return V, M, T, GV
+end
+
+###############################################################################
+#
+#  Primitive extensions
+#
+###############################################################################
+
+### Local glue maps
+
+function _local_glue_maps(
+  q1::TorQuadModule,
+  q2::TorQuadModule,
+  parity::Symbol;
+  Ctx=nothing,
+  vi::Tuple{Int, Int}=(-1, -1),
+  kwargs...,
+)
+  Fac = ZZLatGluingFactory(q1, q2, parity)
+  if !isnothing(Ctx)
+    Fac.Ctx = Ctx
+    Fac.vertex_identification = vi
+  end
+  init_gluing_factory!(Fac; kwargs...)
+  is_trivial(Fac) && return ZZLatGluing[]
+  res = local_glue_maps(Fac)
+  return Fac, res
+end
+
+### Unimodular case
+
+function unimodular_primitive_extensions(
+  M::ZZLat,
+  N::ZZLat;
+  right_action::Union{MatGroup{QQFieldElem, QQMatrix}, Nothing}=nothing,
+  right_discriminant_action::Union{AutomorphismGroup{TorQuadModule}, Nothing}=nothing,
+  left_action::Union{MatGroup{QQFieldElem, QQMatrix}, Nothing}=nothing,
+  left_discriminant_action::Union{AutomorphismGroup{TorQuadModule}, Nothing}=nothing,
+  first::Bool=false,
+  exist_only::Bool=false,
+  even::Bool=(is_even(M) && is_even(N)), #low priority
+  parity::Symbol=(even ? :even : :both),
+)
+  @req is_integral(M) && is_integral(N) "Only available for integral lattices"
+
+  if !isnothing(right_discriminant_action)
+    GMbar = right_discriminant_action
+  elseif !isnothing(right_action)
+    GMbar, _ = image(discriminant_representation(M, right_action; full=false))
+  elseif first || exist_only
+    qM = discriminant_group(M)
+    GMbar = Oscar._orthogonal_group(qM, TorQuadModuleMap[id_hom(qM)]; check=false)
+  else
+    GMbar, _ = image_in_Oq(M)
+  end
+
+  if !isnothing(left_discriminant_action)
+    GNbar = left_discriminant_action
+  elseif !isnothing(left_action)
+    GNbar, _ = image(discriminant_representation(N, left_action; full=false))
+  elseif first || exist_only
+    qN = discriminant_group(N)
+    GNbar = Oscar._orthogonal_group(qN, TorQuadModuleMap[id_hom(qN)]; check=false)
+  else
+    GNbar, _ = image_in_Oq(N)
+  end
+
+  if !is_even(M) || !is_even(N)
+    (parity == :even) && Tuple{ZZLat, ZZLat, ZZLat}[]
+  end
+  return _unimodular_primitive_extensions(M, N, parity, GMbar, GNbar; first, exist_only)
+end
+
+function _unimodular_primitive_extensions(
+  M::ZZLat,
+  N::ZZLat,
+  parity::Symbol,
+  GM::AutomorphismGroup{TorQuadModule},
+  GN::AutomorphismGroup{TorQuadModule};
+  first::Bool=false,
+  exist_only::Bool=false,
+)
+  results = Tuple{ZZLat, ZZLat, ZZLat}[]
+  qM = discriminant_group(M)
+  @assert qM === domain(GM)
+  qN = discriminant_group(N)
+  @assert qN === domain(GN)
+
+  parity == :even && (!is_even(M) || !is_even(N)) && return false, results
+
+  as_bilinear_module = (parity !== :even)
+  ok, phi = is_anti_isometric_with_anti_isometry(qM, qN; as_bilinear_module)
+  !ok && return false, results
+  if exist_only
+    if parity != :odd || !is_even(M) || !is_even(N) || !is_anti_isometry(phi; as_bilinear_module=false)
+      return true, results
+    end
+  end
+
+  D, inj = direct_sum(qM, qN; cached=false, as_bilinear_module)
+  qMinD, qNinD = inj
+  if first
+    if parity != :odd || !is_even(M) || !is_even(N) || !is_anti_isometry(phi; as_bilinear_module=false)
+      L, _ = _overlattice_with_graph(phi, qMinD, qNinD)
+      M2, N2 = _as_sublattices(L, M, N)
+      push!(results, (L, M2, N2))
+      return true, results
+    end
+  end
+
+  iphi = inv(phi)
+  if parity == :even
+    OqN = orthogonal_group(qN)
+  else
+    OqN = orthogonal_group_bilinear(qN)
+  end
+
+  genC = elem_type(OqN)[OqN(iphi * hom(g) * phi; check=false) for g in gens(GM)]
+  GMphi, _ = sub(OqN, genC)
+  elqN = elementary_divisors(qN)
+  if isone(order(qN)) || elqN[1] == elqN[end]
+    iso = isomorphism(PermGroup, OqN)
+  else
+    iso = id_hom(OqN)
+  end
+  reps = double_cosets(codomain(iso), Base.first(iso(GN)), Base.first(iso(GMphi)))
+  for _g in reps
+    g = iso\(representative(_g))
+    phig = compose(phi, hom(g))
+    L, _ = _overlattice_with_graph(phig, qMinD, qNinD)
+    if parity == :even
+      is_even(L) || continue
+    elseif parity == :odd
+      !is_even(L) || continue
+    end
+    exist_only && return true, results
+    M2, N2 = _as_sublattices(L, M, N)
+    push!(results, (L, M2, N2))
+    first && return true, results
+  end
+  return length(results) > 0, results
+end
+
+function _primitive_extensions_coprime_left(
+  M::ZZLat,
+  N::ZZLat,
+  parity::Symbol,
+  GM::AutomorphismGroup{TorQuadModule},
+  GN::AutomorphismGroup{TorQuadModule};
+  first::Bool=false,
+  exist_only::Bool=false,
+)
+  results = Tuple{ZZLat, ZZLat, ZZLat}[]
+  qM = discriminant_group(M)
+  @assert qM === domain(GM)
+  qN = discriminant_group(N)
+  @assert qN === domain(GN)
+
+  pds = prime_divisors(order(qM))
+  gensHN = TorQuadModuleElem[]
+  for p in pds
+    _H, j = primary_part(qN, p)
+    append!(gensHN, j.(gens(_H)))
+  end
+  HN, HNinqN = sub(qN, gensHN)
+
+  as_bilinear_module = (parity !== :even)
+  ok, phi = is_anti_isometric_with_anti_isometry(qM, HN; as_bilinear_module)
+  !ok && return false, results
+  if exist_only
+    if parity != :odd || !is_even(M) || !is_even(N) || !is_anti_isometry(phi; as_bilinear_module=false)
+      return true, results
+    end
+  end
+
+  D, inj = direct_sum(qM, qN; cached=false, as_bilinear_module)
+  qMinD, qNinD = inj
+  HNinD = compose(HNinqN, qNinD)
+  if first
+    if parity != :odd || !is_even(M) || !is_even(N) || !is_anti_isometry(phi; as_bilinear_module=false)
+      L, _ = _overlattice_with_graph(phi, qMinD, HNinD)
+      M2, N2 = _as_sublattices(L, M, N)
+      push!(results, (L, M2, N2))
+      return true, results
+    end
+  end
+
+  iphi = inv(phi)
+  if parity == :even
+    is_even(M) && is_even(N) || return results
+    OHN = orthogonal_group(HN)
+  else
+    OHN = orthogonal_group_bilinear(HN)
+  end
+  
+  imHN, resHN = restrict_automorphism_group(GN, HNinqN)
+  genC = elem_type(OHN)[OHN(iphi * hom(g) * phi; check=false) for g in gens(GM)]
+  GMphi, _ = sub(OHN, genC)
+  elHN = elementary_divisors(HN)
+  if isone(order(HN)) || elHN[1] == elHN[end]
+    iso = isomorphism(PermGroup, OHN)
+  else
+    iso = id_hom(OHN)
+  end
+  reps = double_cosets(codomain(iso), Base.first(iso(imHN)), Base.first(iso(GMphi)))
+  for _g in reps
+    g = iso\(representative(_g))
+    phig = compose(phi, hom(g))
+    L, _ = _overlattice_with_graph(phig, qMinD, HNinD)
+    if parity == :even
+      is_even(L) || continue
+    elseif parity == :odd
+      !is_even(L) || continue
+    end
+    exist_only && return true, results
+    M2, N2 = _as_sublattices(L, M, N)
+    push!(results, (L, M2, N2))
+    first && return true, results
+  end
+  return length(results) > 0, results
+end
+
+function _primitive_extensions_coprime_right(
+  M::ZZLat,
+  N::ZZLat,
+  parity::Symbol,
+  GM::AutomorphismGroup{TorQuadModule},
+  GN::AutomorphismGroup{TorQuadModule};
+  kwargs...
+)
+  r = _primitive_extensions_coprime_left(N, M, parity, GN, GM; kwargs...)
+  t = [1,3,2]
+  return eltype(r)[x[t] for x in r]
+end
+
+###############################################################################
+#
+#  Primitive embeddings
+#
+###############################################################################
+
+### Unimodular top lattice
+
+function _primitive_embeddings_in_unimodular(
+  G1::ZZGenus,
+  G2s::Vector{ZZGenus};
+  check::Bool=true,
+)
+  @assert is_unimodular(G1)
+  results = Vector{ZZLat, ZZLat, ZZLat}[]
+
+  G2s = filter(<(rank(G1))∘rank, G2s)
+  filter!(<(signature_pair(G1))∘signature_pair, G2s)
+  if is_even(G1)
+    filter!(is_even, G2s)
+  end
+
+  for G2 in G2s
+    append!(results, _primitive_embeddings_in_unimodular_safe(G1, G2))
+  end
+  return results
+end
+
+function _primitive_embeddings_in_unimodular(
+  G1::ZZGenus,
+  G2::ZZGenus,
+)
+  @assert is_unimodular(G1)
+  if is_even(G1) && !is_even(G2)
+    return results
+  elseif rank(G1) <= rank(G2)
+    return results
+  elseif !(signature(G2) < signature(G1))
+    return results
+  end
+  return _primitive_embeddings_in_unimodular_safe(G1, G2)
+end
+
+function _primitive_embeddings_in_unimodular_safe(
+  G1::ZZGenus,
+  G2::ZZGenus,
+)
+  results = Tuple{ZZLat, ZZLat, ZZLat}[]
+  if is_even(G1)
+    parity = :even
+    par = :even
+  else
+    parity = :odd
+    par = :both
+  end
+  sign = signature_pair(G1) .- signature_pair(G2)
+  q = rescale(discriminant_group(G2), -1; cached=false)
+  GKs = _integer_genera(q, sign, par)
+  isempty(GKs) && return results
+ 
+  Ns = ZZLat[]
+  for GK in GKs
+    append!(Ns, __representatives(GK))
+  end
+  Ms = __representatives(G2)
+  for M in Ms, N in Ns
+    append!(results, unimodular_primitive_extensions(M, N; parity))
+  end
+  return results
+end
+
+### Coprime det
+
+function _primitive_embeddings_coprime_det_safe(
+  G1::ZZGenus,
+  G2::ZZGenus,
+)
+  results = Tuple{ZZLat, ZZLat, ZZLat}[]
+  if is_even(G1)
+    parity = :even
+    par = :even
+  else
+    parity = :odd
+    par = :both
+  end
+  sign = signature_pair(G1) .- signature_pair(G2)
+  q1 = discriminant_group(G1)
+  q2 = discriminant_group(G2)
+  q, _ = direct_sum(q1, rescale(q2, -1; cached=false); cached=false, as_bilinear_module=(par !== :even))
+  GKs = _integer_genera(q, sign, par)
+  isempty(GKs) && return results
+
+  Ns = ZZLat[]
+  for GK in GKs
+    append!(Ns, __representatives(GK))
+  end
+  Ms = __representatives(G2)
+  for M in Ms
+    GM, _ = image_in_Oq(M)
+    for N in Ns
+      GN, _ = image_in_Oq(N)
+      tmp = _primitive_extensions_coprime_left(M, N, parity, GM, GN)
+      append!(results, tmp)
+    end
+  end
+  return results
+end
+
+### General case
+
+function _primitive_embeddings(
+  G1s::Vector{ZZGenus},
+  G2s::Vector{ZZGenus},
+)
+  res = Tuple{ZZLat, ZZLat, ZZLat}[]
+  for G2 in G2s
+    append!(res, _primitive_embeddings(G1s, G2))
+  end
+  return res
+end
+
+function _primitive_embeddings(
+  G1s::Vector{ZZGenus},
+  G2::ZZGenus,
+)
+  results = Tuple{ZZLat, ZZLat, ZZLat}[]
+  G1s = filter(>(rank(G2))∘rank, G1s)
+  filter!(>(signature_pair(G2))∘signature_pair, G1s)
+  if !is_even(G2)
+    filter!(!is_even, G1s)
+  end
+  q2 = discriminant_group(G2)
+  Ctx = ZZLatGluingCtx()
+  push!(Ctx.modules, q2)
+
+  for G1 in G1s
+    if is_unimodular(G1)
+      append!(results, _primitive_embeddings_in_unimodular_safe(G1, G2))
+    elseif isone(gcd(numerator(det(G1)), numerator(det(G2))))
+      append!(results, _primitive_embeddings_coprime_det_safe(G1, G2))
+    else
+      append!(results, _primitive_embeddings_generic_safe(G1, G2; Ctx, vi=(-1, 1), q2))
+    end
+  end
+  return results
+end
+
+function _primitive_embeddings(
+  G1::ZZGenus,
+  G2::ZZGenus,
+)
+  if is_even(G1) && !is_even(G2)
+    return results
+  elseif rank(G1) <= rank(G2)
+    return results
+  elseif !(signature(G2) < signature(G1))
+    return results
+  end
+  if is_unimodular(G1)
+    return _primitive_embeddings_in_unimodular_safe(G1, G2)
+  elseif isone(gcd(numerator(det(G1)), numerator(det(G2))))
+    return _primitive_embeddings_coprime_det_safe(G1, G2)
+  else
+    return _primitive_embeddings_generic_safe(G1, G2)
+  end
+end
+
+function _primitive_embeddings_generic_safe(
+  G1::ZZGenus,
+  G2::ZZGenus;
+  Ctx=nothing,
+  vi::Tuple{Int, Int}=(-1, -1),
+  q2::TorQuadModule=discriminant_group(G2),
+  Ms::Vector{ZZLat}=ZZLat[],
+)
+  results = Tuple{ZZLat, ZZLat, ZZLat}[]
+  R = rescale(representative(G1), -1; cached=false)
+  U = hyperbolic_plane_lattice()
+  if is_even(G1)
+    parity = :even
+    T, _ = direct_sum(R, U; cached=false)
+  else
+    parity = :both
+    T, _ = direct_sum(R, U, U; cached=false)
+  end
+  q1n = discriminant_group(T)
+  signK = signature_pair(G1) .- signature_pair(G2)
+
+  Fac, lgm = _local_glue_maps(q2, q1n, parity; Ctx, vi=reverse(vi))
+  isempty(lgm) && return results
+  DKs = Dict{ZZGenus, Vector{ZZLat}}()
+  
+  for x in lgm
+    qx = _form_over(x)
+    qK = rescale(qx, -1; cached=false)
+    GKs = _integer_genera(qK, signK, parity)
+    isempty(GKs) && continue
+    # At that point, we know that we have an extension
+    for GK in GKs
+      haskey(DKs, GK) && continue
+      DKs[GK] = representatives(GK)
+    end
+    isempty(Ms) && append!(Ms, __representatives(G2))
+    for M in Ms
+      qM = discriminant_group(M)
+      D = _gluing_ambient(qM, q1n, parity)
+      GM, _ = image_in_Oq(M)
+      _ok, phiM = is_isometric_with_isometry(qM, q2)
+      @assert _ok
+      xM = _pullback_left(x, phiM)
+      xMs = _split_orbit_left_group(xM, GM)
+      for _x in xMs, y in _all_glue_maps(_x)
+        V, M2, T2, GV = _overlattice(y, D)
+        resV = Tuple{ZZLat, ZZLat, ZZLat}[]
+        qV = domain(GV)
+        for GK in GKs, K in DKs[GK]
+          _, pe = unimodular_primitive_extensions(V, K; right_discriminant_action=GV, parity)
+          append!(resV, pe)
+        end
+        for (S, V2, W2) in resV
+          T3 = lattice_in_same_ambient_space(S, hcat(basis_matrix(T2), zero_matrix(QQ, rank(T2), degree(W2)-degree(T2))))
+          L = lll(orthogonal_submodule(S, T3))
+          genus(L) == G1 || continue
+          M3 = lattice_in_same_ambient_space(S, hcat(basis_matrix(M2), zero_matrix(QQ, rank(M2), degree(W2)-degree(M2))))
+          @assert is_sublattice(L, M3)
+          @assert is_primitive(L, M3)
+          N = orthogonal_submodule(L, M3)
+          bM = coordinates(basis_matrix(M3), L)
+          bN = coordinates(basis_matrix(N), L)
+          L = integer_lattice(; gram=gram_matrix(L))
+          M3 = lattice_in_same_ambient_space(L, bM)
+          N = lattice_in_same_ambient_space(L, bN)
+          push!(results, (L, M3, N))
+        end
+      end
+    end
+  end
+  return results
+end
+
+###############################################################################
+#
+#  Temporary
+#
+###############################################################################
+
+# To connect to the LMFDB database:
+# - fork https://github.com/thofma/LMFDB.jl
+# - run `using LMFDB`
+# - fetch the database: `db = LMFDB.LMFDBLite.LMFDBConnection()`
+# - load some genera: z.B. `l6=LMFDB.genera(db; rank = 6, nplus=6, det = >=(1) & <=(div(300, 8, RoundUp)));`
+
+function _primitive_embeddings_lmfdb(
+  G1s::Vector{ZZGenus},
+  G2s::Vector{ZZGenus},
+)
+  results = Tuple{String, String, ZZMatrix}[]
+  for G2 in G2s
+    append!(results, _primitive_embeddings_lmfdb(G1s, G2))
+  end
+  return results
+end
+
+function _primitive_embeddings_lmfdb(
+  G1s::Vector{ZZGenus},
+  G2::ZZGenus,
+)
+  results = Tuple{String, String, ZZMatrix}[]
+  G1s = filter(>(rank(G2))∘rank, G1s)
+  filter!(>(signature_pair(G2))∘signature_pair, G1s)
+  if !is_even(G2)
+    filter!(!is_even, G1s)
+  end
+  q2 = discriminant_group(G2)
+  Ctx = ZZLatGluingCtx()
+  push!(Ctx.modules, q2)
+
+  for G1 in G1s
+    if is_unimodular(G1)
+      append!(results, _primitive_embeddings_in_unimodular_safe_lmfdb(G1, G2))
+    elseif isone(gcd(numerator(det(G1)), numerator(det(G2))))
+      append!(results, _primitive_embeddings_coprime_det_safe_lmfdb(G1, G2))
+    else
+      append!(results, _primitive_embeddings_generic_safe_lmfdb(G1, G2; Ctx, vi=(-1, 1), q2))
+    end
+  end
+  return results
+end
+
+function _primitive_embeddings_in_unimodular_safe_lmfdb(
+  G1::ZZGenus,
+  G2::ZZGenus,
+)
+  results = Tuple{String, String, ZZMatrix}[]
+  if is_even(G1)
+    parity = :even
+    par = :even
+  else
+    parity = :odd
+    par = :both
+  end
+  sign = signature_pair(G1) .- signature_pair(G2)
+  q = rescale(discriminant_group(G2), -1; cached=false)
+  GKs = _integer_genera(q, sign, par)
+  isempty(GKs) && return results
+
+  Ns = ZZLat[]
+  for GK in GKs
+    append!(Ns, __representatives(GK))
+  end
+  Ms = __representatives(G2)
+  for M in Ms
+    label_bottom = get_attribute(M, :lmfdb_label)
+    for N in Ns
+      tmp = unimodular_primitive_extensions(M, N; parity)
+      for (T, M2, N2) in tmp
+        @assert gram_matrix(M) == gram_matrix(M2)
+        genus(T) != G1 && continue #Should never happen but still
+        for T2 in get_attribute(G1, :representatives)
+          if is_definite(G1)
+            ok, f_top = is_isometric_with_isometry(T, T2)
+          else
+            ok = is_isometric(T, T2)
+          end
+          !ok && continue
+          label_top = get_attribute(T2, :lmfdb_label)
+          if !is_definite(G1)
+            push!(results, (label_top, label_bottom, zero_matrix(ZZ, 0, 0)))
+          else
+            v = map_entries(ZZ, coordinates(basis_matrix(N2), T)*f_top)
+            push!(results, (label_top, label_bottom, v))
+          end
+        end
+      end
+    end
+  end
+  return results
+end
+
+### Coprime det
+
+function _primitive_embeddings_coprime_det_safe_lmfdb(
+  G1::ZZGenus,
+  G2::ZZGenus,
+)
+  results = Tuple{String, String, ZZMatrix}[]
+  if is_even(G1)
+    parity = :even
+    par = :even
+  else
+    parity = :odd
+    par = :both
+  end
+  sign = signature_pair(G1) .- signature_pair(G2)
+  q1 = discriminant_group(G1)
+  q2 = discriminant_group(G2)
+  q, _ = direct_sum(q1, rescale(q2, -1; cached=false); cached=false, as_bilinear_module=(par !== :even))
+  GKs = _integer_genera(q, sign, par)
+  isempty(GKs) && return results
+
+  Ns = ZZLat[]
+  for GK in GKs
+    append!(Ns, __representatives(GK))
+  end
+  Ms = __representatives(G2)
+  for M in Ms
+    GM, _ = image_in_Oq(M)
+    label_bottom = get_attribute(M, :lmfdb_label)
+    for N in Ns
+      GN, _ = image_in_Oq(N)
+      _, tmp = _primitive_extensions_coprime_left(M, N, parity, GM, GN)
+      for (T, M2, N2) in tmp
+        @assert gram_matrix(M) == gram_matrix(M2)
+        genus(T) != G1 && continue #Should never happen but still
+        for T2 in get_attribute(G1, :representatives)
+          if is_definite(G1)
+            ok, f_top = is_isometric_with_isometry(T, T2)
+          else
+            ok = is_isometric(T, T2)
+          end
+          !ok && continue
+          label_top = get_attribute(T2, :lmfdb_label)
+          if !is_definite(G1)
+            push!(results, (label_top, label_bottom, zero_matrix(ZZ, 0, 0)))
+          else
+            v = map_entries(ZZ, coordinates(basis_matrix(N2), T)*f_top)
+            push!(results, (label_top, label_bottom, v))
+          end
+        end
+      end
+    end
+  end
+  return results
+end
+
+
+function _primitive_embeddings_generic_safe_lmfdb(
+  G1::ZZGenus,
+  G2::ZZGenus;
+  Ctx=nothing,
+  vi::Tuple{Int, Int}=(-1, -1),
+  q2::TorQuadModule=discriminant_group(G2),
+  Ms::Vector{ZZLat}=ZZLat[],
+)
+  results = Tuple{String, String, ZZMatrix}[]
+  R = rescale(representative(G1), -1; cached=false)
+  U = hyperbolic_plane_lattice()
+  if is_even(G1)
+    parity = :even
+    T, _ = direct_sum(R, U; cached=false)
+  else
+    parity = :both
+    T, _ = direct_sum(R, U, U; cached=false)
+  end
+  q1n = discriminant_group(T)
+  signK = signature_pair(G1) .- signature_pair(G2)
+
+  Fac, lgm = _local_glue_maps(q2, q1n, parity; Ctx, vi=reverse(vi))
+  isempty(lgm) && return results
+  DKs = Dict{ZZGenus, Vector{ZZLat}}()
+
+  for x in lgm
+    qx = _form_over(x)
+    qK = rescale(qx, -1; cached=false)
+    GKs = _integer_genera(qK, signK, parity)
+    isempty(GKs) && continue
+    # At that point, we know that we have an extension
+    for GK in GKs
+      haskey(DKs, GK) && continue
+      DKs[GK] = representatives(GK)
+    end
+    isempty(Ms) && append!(Ms, __representatives(G2))
+    for M in Ms
+      label_bottom = get_attribute(M, :lmfdb_label)
+      qM = discriminant_group(M)
+      D = _gluing_ambient(qM, q1n, parity)
+      GM, _ = image_in_Oq(M)
+      _ok, phiM = is_isometric_with_isometry(qM, q2)
+      @assert _ok
+      xM = _pullback_left(x, phiM)
+      xMs = _split_orbit_left_group(xM, GM)
+      for _x in xMs, y in _all_glue_maps(_x)
+        V, M2, T2, GV = _overlattice(y, D)
+        resV = Tuple{ZZLat, ZZLat, ZZLat}[]
+        qV = domain(GV)
+        for GK in GKs, K in DKs[GK]
+          _, pe = unimodular_primitive_extensions(V, K; right_discriminant_action=GV, parity)
+          append!(resV, pe)
+        end
+        for (S, V2, W2) in resV
+          T3 = lattice_in_same_ambient_space(S, hcat(basis_matrix(T2), zero_matrix(QQ, rank(T2), degree(W2)-degree(T2))))
+          L = lll(orthogonal_submodule(S, T3))
+          genus(L) == G1 || continue
+          M3 = lattice_in_same_ambient_space(S, hcat(basis_matrix(M2), zero_matrix(QQ, rank(M2), degree(W2)-degree(M2))))
+          @assert is_sublattice(L, M3)
+          @assert is_primitive(L, M3)
+          N = orthogonal_submodule(L, M3)
+          bM = coordinates(basis_matrix(M3), L)
+          bN = coordinates(basis_matrix(N), L)
+          L = integer_lattice(; gram=gram_matrix(L))
+          M3 = lattice_in_same_ambient_space(L, bM)
+          N = lattice_in_same_ambient_space(L, bN)
+          for L2 in __representatives(G1)
+            if is_definite(G1)
+              ok, f_top = is_isometric_with_isometry(L, L2)
+            else
+              ok = is_isometric(L, L2)
+            end
+            !ok && continue
+            label_top = get_attribute(L2, :lmfdb_label)
+            if is_definite(G1)
+              v = map_entries(ZZ, coordinates(bN*f_top, L2))
+              push!(results, (label_top, label_bottom, v))
+            else
+              push!(results, (label_top, label_bottom, zero_matrix(ZZ, 0, 0)))
+            end
+          end
+        end
+      end
+    end
+  end
+  return results
+end
+
diff --git a/src/NumberTheory/QuadFormAndIsom/types.jl b/src/NumberTheory/QuadFormAndIsom/types.jl
index 75b7ef5288d1..d5cb00ed23e8 100644
--- a/src/NumberTheory/QuadFormAndIsom/types.jl
+++ b/src/NumberTheory/QuadFormAndIsom/types.jl
@@ -1,3 +1,9 @@
+###############################################################################
+#
+#  Basic types
+#
+###############################################################################
+
 @doc raw"""
     QuadSpaceWithIsom
 
@@ -265,6 +271,12 @@ Finite quadratic module of order 3
   end
 end
 
+#####################################################################
+#
+#  Enumeration context
+#
+#####################################################################
+
 mutable struct ZZLatWithIsomEnumCtX
   # what we want now
   power::Int
@@ -291,3 +303,107 @@ mutable struct ZZLatWithIsomEnumCtX
     return new()
   end
 end
+
+###############################################################################
+#
+#  Gluing factory
+#
+###############################################################################
+
+@doc raw"""
+    ZZLatGluingCtx
+
+A context object which stores some `TorQuadModule` together with the results of
+some orbits and stabilizers computations for some ``p``-subgroups of a given
+type.
+"""
+struct ZZLatGluingCtx
+  modules::Vector{TorQuadModule}
+  orb_and_stab::Dict{Tuple{Int, ZZRingElem, Vector{Int}}, Vector{Tuple{TorQuadModuleMap, GAPGroupHomomorphism}}}
+
+  function ZZLatGluingCtx()
+    modules = TorQuadModule[]
+    orb_and_stab = Dict{Tuple{Int, ZZRingElem, Vector{Int}}, Vector{Tuple{TorQuadModuleMap, GAPGroupHomomorphism}}}()
+    return new(modules, orb_and_stab)
+  end
+end
+
+@doc raw"""
+    ZZLatGluing
+
+Type for representatives for a double coset of glue maps between two
+discriminant forms ``q_1`` and ``q_2``, under the action some groups
+``G_1 \subset GL(q_1)`` and ``G_2 \subset GL(q_2)``. It consists:
+- a glue map ``\gamma\colon H_1\to H_2``,
+- the inverse glue map of ``H_2\to H_1``,
+- the embedding ``H_1 \to q_1``,
+- the embedding of the ``G_1``-stabilizer of ``H_1`` inside ``G_1``,
+- the embedding ``H_2 \to q_2``,
+- the embedding of the ``G_2``-stabilizer of ``H_2`` inside ``G_2``.
+"""
+struct ZZLatGluing
+  glue_map::TorQuadModuleMap
+  inv_glue_map::TorQuadModuleMap
+  glue_group_left::TorQuadModuleMap
+  stabilizer_left::GAPGroupHomomorphism
+  glue_group_right::TorQuadModuleMap
+  stabilizer_right::GAPGroupHomomorphism
+
+  function ZZLatGluing(x...)
+    return new(x...)
+  end
+end
+
+@doc raw"""
+    ZZLatGluingFactory
+
+A factory object for local gluings computations. It stores two `TorQuadModule`
+and a subgroup of their orthogonal group (as finite bilinear module). A context
+object `ZZLatGluingCtx` can be added to it in the context of gluings where one
+of the two modules is fixed and the second varies, for instance.
+"""
+mutable struct ZZLatGluingFactory
+  # Necessary input
+  ambient_modules::Tuple{TorQuadModule, TorQuadModule} # Discriminant forms of the lattices to glue
+  local_classifying_groups::Tuple{AutomorphismGroup{TorQuadModule}, AutomorphismGroup{TorQuadModule}} # By default the orthogonal groups of the ambient modules
+  par::Symbol # :even, :odd or :both
+
+  # Context
+  Ctx::ZZLatGluingCtx
+  vertex_identification::Tuple{Int, Int}
+
+  # Edge conditions
+  glue_order::Set{ZZRingElem} # Possible orders of a glue group
+  glue_elementary_divisors::Set{Vector{ZZRingElem}} # Possible elementary divisors of a glue groups
+  genus_over::Set{ZZGenus} # Possible genus of an overlattice
+  form_over::Set{QQMatrix} # Possible Gram matrix normal form of discriminant group overlattice
+
+  # Internal preparation
+  conditions_modules::Tuple{TorQuadModuleMap, TorQuadModuleMap} # Where the glue groups should be taken
+  glue_group_parent_snf::Vector{ZZRingElem} # Glue groups are isomorphic to a subgroup of an abelian group with such elementary divisors
+  primes_of_interest::Set{ZZRingElem} # Only primes which could divide the order of a glue group
+  local_gluings_primary::Dict{ZZRingElem, Dict{Vector{Int}, Vector{ZZLatGluing}}} # Intermediate storage for local gluings at some prime, and glue group of certain type
+
+  function ZZLatGluingFactory(
+    module_left::TorQuadModule,
+    module_right::TorQuadModule,
+    par::Symbol,
+  )
+    z = new((module_left, module_right))
+    z.par = par
+    return z
+  end
+end
+
+struct ZZLatGluingAmbient
+  D::TorQuadModule
+  j1::TorQuadModuleMap
+  j2::TorQuadModuleMap
+  OD::AutomorphismGroup{TorQuadModule}
+  k1::GAPGroupHomomorphism
+  k2::GAPGroupHomomorphism
+
+  function ZZLatGluingAmbient(x...)
+    return new(x...)
+  end
+end