Improve isomorphism from pc-groups and fp-groups to FinGenAbGroup

src/GAP/wrappers.jl

@@ -13,6 +13,8 @@ module GAPWrap
 using GAP
 
 # the following list is intended to be sorted according to `LC_COLLATE=C sort -f`, i.e. ignoring case
+GAP.@wrap AbelianGroup(x::GapObj, y::GapObj)::GapObj
+GAP.@wrap AbelianPcpGroup(x::GAP.Obj, y::GapObj)::GapObj
 GAP.@wrap AlgebraicExtension(x::GapObj, y::GapObj)::GapObj
 GAP.@wrap AlgExtElm(x::GapObj, y::GAP.Obj)::GapObj
 GAP.@wrap AntiSymmetricParts(x::GapObj, y::GapObj, z::GapInt)::GapObj
@@ -134,6 +136,7 @@ GAP.@wrap ImagesRepresentative(x::GapObj, y::Any)::GAP.Obj
 GAP.@wrap ImagesSource(x::GapObj)::GapObj
 GAP.@wrap ImmutableMatrix(x::GapObj, y::GapObj, z::Bool)::GapObj
 GAP.@wrap IndependentGeneratorExponents(x::Any, y::Any)::GapObj
+GAP.@wrap IndependentGeneratorsOfAbelianGroup(x::GapObj)::GapObj
 GAP.@wrap Indeterminate(x::GapObj)::GapObj
 GAP.@wrap Indeterminate(x::GapObj, y::GAP.Obj)::GapObj
 GAP.@wrap IndeterminateNumberOfUnivariateRationalFunction(x::GapObj)::Int
@@ -273,6 +276,7 @@ GAP.@wrap LibInfoCharacterTable(x::GapObj)::GapObj
 GAP.@wrap LieAlgebraByStructureConstants(x::GapObj, y::GapObj)::GapObj
 GAP.@wrap LinearCharacters(x::GapObj)::GapObj
 GAP.@wrap LinearCombination(x::GapObj, y::GapObj)::GapObj
+GAP.@wrap LinearCombinationPcgs(x::GapObj, y::GapObj)::GapObj
 GAP.@wrap ListPerm(x::GapObj)::GapObj
 GAP.@wrap MarksTom(x::GapObj)::GapObj
 GAP.@wrap MatScalarProducts(x::GapObj, y::GapObj, z::GapObj)::GapObj
@@ -305,6 +309,7 @@ GAP.@wrap Order(x::Any)::GapInt
 GAP.@wrap OrthogonalComponents(x::GapObj, y::GapObj, z::GapInt)::GapObj
 GAP.@wrap PcElementByExponentsNC(x::GapObj, y::GapObj)::GapObj
 GAP.@wrap Pcgs(x::GapObj)::GapObj
+GAP.@wrap PcpElementByExponentsNC(x::GapObj, y::GapObj)::GapObj
 GAP.@wrap PcpGroupByCollectorNC(x::GapObj)::GapObj
 GAP.@wrap PCore(x::GapObj, y::GapInt)::GapObj
 GAP.@wrap PermList(x::GapObj)::GapObj
@@ -324,6 +329,7 @@ GAP.@wrap Random(x::GapObj, y::GapObj)::GAP.Obj
 GAP.@wrap Range(x::GapObj)::GapObj
 GAP.@wrap RecognizeGroup(x::GapObj)::GapObj
 GAP.@wrap ReduceCoeffs(x::GapObj, y::GapObj)
+GAP.@wrap RelativeOrders(x::GapObj)::GapObj
 GAP.@wrap RelatorsOfFpGroup(x::GapObj)::GapObj
 GAP.@wrap Representative(x::GapObj)::GAP.Obj
 GAP.@wrap RepresentativeTom(x::GapObj, y::Int)::GapObj
@@ -351,6 +357,7 @@ GAP.@wrap Stabilizer(v::GapObj, w::Any, x::GapObj, y::GapObj, z::GapObj)::GapObj
 GAP.@wrap StringViewObj(x::Any)::GapObj
 GAP.@wrap StructureConstantsTable(x::GapObj)::GapObj
 GAP.@wrap StructureDescription(x::GapObj)::GapObj
+GAP.@wrap SubgroupNC(x::GapObj, y::GapObj)::GapObj
 GAP.@wrap SubsTom(x::GapObj)::GapObj
 GAP.@wrap SylowSubgroup(x::GapObj, y::GapInt)::GapObj
 GAP.@wrap SymmetricParts(x::GapObj, y::GapObj, z::GapInt)::GapObj

src/Groups/GAPGroups.jl

@@ -2057,8 +2057,7 @@ julia> invs
 """
 function map_word(g::Union{FPGroupElem, SubFPGroupElem}, genimgs::Vector; genimgs_inv::Vector = Vector(undef, length(genimgs)), init = nothing)
   G = parent(g)
-  Ggens = gens(G)
-  if length(Ggens) == 0
+  if ngens(G) == 0
     @req init !== nothing "use '; init =...' if there are no generators"
     return init
   end
@@ -2107,18 +2106,11 @@ x*y^4
 ```
 """
 function map_word(g::Union{PcGroupElem, SubPcGroupElem}, genimgs::Vector; genimgs_inv::Vector = Vector(undef, length(genimgs)), init = nothing)
-  G = parent(g)
-  Ggens = gens(G)
-  if length(Ggens) == 0
+  if ngens(parent(g)) == 0
+    @req init !== nothing "use '; init =...' if there are no generators"
     return init
   end
-  gX = GapObj(g)
-
-  if GAPWrap.IsPcGroup(GapObj(G))
-    l = GAP.Globals.ExponentsOfPcElement(GAP.Globals.FamilyPcgs(GapObj(G)), gX)
-  else  # GAP.Globals.IsPcpGroup(GapObj(G))
-    l = GAP.Globals.Exponents(gX)
-  end
+  l = _exponent_vector(g)
   @assert length(l) == length(genimgs)
   ll = Pair{Int, Int}[i => l[i] for i in 1:length(l)]
   return map_word(ll, genimgs, genimgs_inv = genimgs_inv, init = init)
@@ -2216,26 +2208,44 @@ Return the syllables of `g` as a list of pairs `gen => exp` where
 julia> F = @free_group(:F1, :F2);
 
 julia> syllables(F1^5*F2^-3)
-2-element Vector{Pair{Int64, Int64}}:
+2-element Vector{Pair{Int64, ZZRingElem}}:
  1 => 5
  2 => -3
 
 julia> syllables(one(F))
-Pair{Int64, Int64}[]
+Pair{Int64, ZZRingElem}[]
 
 julia> G, epi = quo(F, [F1^10, F2^10]);
 
 julia> syllables(epi(F1^5*F2^-3))
-2-element Vector{Pair{Int64, Int64}}:
+2-element Vector{Pair{Int64, ZZRingElem}}:
  1 => 5
  2 => -3
 ```
 """
 function syllables(g::Union{FPGroupElem, SubFPGroupElem})
   l = GAPWrap.ExtRepOfObj(GapObj(g))
-  return Pair{Int, Int}[l[i] => l[i+1] for i in 1:2:length(l)]
+  return Pair{Int, ZZRingElem}[l[i] => l[i+1] for i in 1:2:length(l)]
+end
+
+function _exponent_vector(g::Union{PcGroupElem, SubPcGroupElem})
+  gX = GapObj(g)
+  G = parent(g)
+  GX = GapObj(G)
+  if GAPWrap.IsPcGroup(GX)
+    return Vector{ZZRingElem}(GAPWrap.ExponentsOfPcElement(GAPWrap.FamilyPcgs(GX), gX))
+  else  # GAP.Globals.IsPcpGroup(GapObj(G))
+    return Vector{ZZRingElem}(GAP.Globals.Exponents(gX)::GapObj)
+  end
 end
 
+function exponents_of_abelianization(g::Union{FPGroupElem, SubFPGroupElem})
+  v = zeros(ZZRingElem, ngens(parent(g)))
+  for (i, e) in syllables(g)
+    v[i] = v[i] + e
+  end
+  return v
+end
 
 @doc raw"""
     letters(g::FPGroupElem)
@@ -2333,15 +2343,27 @@ true
 """
 function (G::FPGroup)(pairs::AbstractVector{Pair{T, S}}) where {T <: IntegerUnion, S <: IntegerUnion}
    n = ngens(G)
-   ll = IntegerUnion[]
+   extrep = IntegerUnion[]
    for p in pairs
      @req 0 < p.first && p.first <= n "generator number is at most $n"
      if p.second != 0
-       push!(ll, p.first)
-       push!(ll, p.second)
+       push!(extrep, p.first)
+       push!(extrep, p.second)
      end
    end
-   return G(ll)
+
+   famG = GAPWrap.ElementsFamily(GAPWrap.FamilyObj(GapObj(G)))
+   if GAP.Globals.IsFreeGroup(GapObj(G))
+     w = GAPWrap.ObjByExtRep(famG, GapObj(extrep, true))
+   else
+     # For quotients of free groups, `GAPWrap.ObjByExtRep` is not defined.
+     F = GAP.getbangproperty(famG, :freeGroup)
+     famF = GAPWrap.ElementsFamily(GAPWrap.FamilyObj(F))
+     w1 = GAPWrap.ObjByExtRep(famF, GapObj(extrep, true))
+     w = GAPWrap.ElementOfFpGroup(famG, w1)
+   end
+
+   return FPGroupElem(G, w)
 end
 
 # This format is used in the serialization of `FPGroupElem`.

src/Groups/homomorphisms.jl

@@ -636,7 +636,7 @@ function isomorphism(T::Type{PcGroup}, G::GAPGroup; on_gens::Bool=false)
          Ggens = GapObj(gens(G); recursive = true)::GapObj
          Cpcgs = GAP.Globals.PcgsByPcSequence(fam, Ggens)::GapObj
          CC = GAP.Globals.PcGroupWithPcgs(Cpcgs)::GapObj
-         CCpcgs = GAP.Globals.FamilyPcgs(CC)::GapObj
+         CCpcgs = GAPWrap.FamilyPcgs(CC)
          f = GAP.Globals.GroupHomomorphismByImages(GapObj(G), CC, Cpcgs, CCpcgs)::GapObj
          return GAPGroupHomomorphism(G, T(CC), f)
        else
@@ -658,7 +658,7 @@ error("do not know how to create a pcp group on given generators in GAP")
          if GAPWrap.IsPcGroup(C)::Bool
            Cpcgs = GAP.Globals.Pcgs(C)::GapObj
            CC = GAP.Globals.PcGroupWithPcgs(Cpcgs)::GapObj
-           CCpcgs = GAP.Globals.FamilyPcgs(CC)::GapObj
+           CCpcgs = GAPWrap.FamilyPcgs(CC)
          else
            Cpcgs = GAP.Globals.Pcp(C)::GapObj
            CC = GAP.Globals.PcpGroupByPcp(Cpcgs)::GapObj
@@ -717,6 +717,9 @@ function isomorphism(::Type{T}, A::FinGenAbGroup) where T <: GAPGroup
    # Known isomorphisms are cached in the attribute `:isomorphisms`.
    isos = get_attribute!(Dict{Tuple{Type, Bool}, Any}, A, :isomorphisms)::Dict{Tuple{Type, Bool}, Any}
    return get!(isos, (T, false)) do
+     @assert T != PcGroup "There should be a special method for type PcGroup"
+     @assert T != FPGroup "There should be a special method for type FPGroup"
+
      # find independent generators
      if is_diagonal(rels(A))
        exponents = diagonal(rels(A))
@@ -727,26 +730,15 @@ function isomorphism(::Type{T}, A::FinGenAbGroup) where T <: GAPGroup
        A2, A2_to_A = snf(A)
      end
      A_to_A2 = inv(A2_to_A)
+
      # Create an isomorphic GAP group whose `GAPWrap.GeneratorsOfGroup`
      # consists of independent elements of the orders in `exponents`.
-     if T == PcGroup
-       # We cannot guarantee that these generators form a pcgs in the case
-       # `T == PcGroup`, hence we cannot call `abelian_group(T, exponents)`.
-       if 0 in exponents
-         GapG = GAP.Globals.AbelianPcpGroup(length(exponents), GapObj(exponents; recursive = true))
-         G = PcGroup(GapG)
-       else
-         GapG = GAP.Globals.AbelianGroup(GAP.Globals.IsPcGroup, GapObj(exponents; recursive = true))
-         G = PcGroup(GAP.Globals.SubgroupNC(GapG, GAP.Globals.FamilyPcgs(GapG)))
-       end
-     else
-       G = abelian_group(T, exponents)
-       GapG = GapObj(G)
-     end
+     G = abelian_group(T, exponents)
+     GapG = GapObj(G)
 
      # `GAPWrap.GeneratorsOfGroup(GapG)` consists of independent elements
      # of the orders in `exponents`.
-     # `GAP.Globals.IndependentGeneratorsOfAbelianGroup(GapG)` chooses generators
+     # `GAPWrap.IndependentGeneratorsOfAbelianGroup(GapG)` chooses generators
      # that may differ from these generators,
      # and that belong to the exponent vectors returned by
      # `GAPWrap.IndependentGeneratorExponents(GapG, g)`.
@@ -771,7 +763,7 @@ function isomorphism(::Type{T}, A::FinGenAbGroup) where T <: GAPGroup
        end
        Ggens = newGgens
      end
-     gensindep = GAP.Globals.IndependentGeneratorsOfAbelianGroup(GapG)::GapObj
+     gensindep = GAPWrap.IndependentGeneratorsOfAbelianGroup(GapG)::GapObj
      orders = [GAPWrap.Order(g) for g in gensindep]
      exps = map(x -> x == GAP.Globals.infinity ? ZZRingElem(0) : ZZRingElem(x), orders)
      Aindep = abelian_group(exps)
@@ -798,6 +790,113 @@ function isomorphism(::Type{T}, A::FinGenAbGroup) where T <: GAPGroup
    end::GroupIsomorphismFromFunc{FinGenAbGroup, T}
 end
 
+################################################################################
+#
+#  special methods for `isomorphism` to abelian `PcGroup` or `FPGroup`:
+#  computing images and preimages is easier in these cases,
+#  since we need not go via `GAP.Globals.IndependentGeneratorsOfAbelianGroup`
+#  and `GAPWrap.IndependentGeneratorExponents`
+#
+################################################################################
+
+function isomorphism(::Type{PcGroup}, A::FinGenAbGroup)
+   # Known isomorphisms are cached in the attribute `:isomorphisms`.
+   isos = get_attribute!(Dict{Tuple{Type, Bool}, Any}, A, :isomorphisms)::Dict{Tuple{Type, Bool}, Any}
+   return get!(isos, (PcGroup, false)) do
+     # find independent generators
+     # trivial diagonal entries can get removed by `GAPWrap.AbelianPcpGroup`,
+     # thus we cannot simply take the diagonal
+     exponents = elementary_divisors(A)
+     A2, A2_to_A = snf(A)
+     n = length(exponents)
+     A_to_A2 = inv(A2_to_A)
+
+     if 0 in exponents
+       # Create a `PcpGroup` in GAP.
+       # Its `GeneratorsOfGroup` are the generators of the defining
+       # presentations, they correspond to the generators of `A2`.
+       GapG = GAPWrap.AbelianPcpGroup(length(exponents), GapObj(exponents; recursive = true))::GapObj
+       C = GAPWrap.Collector(GapG)::GapObj
+       G = PcGroup(GapG)
+
+       f = function(a::FinGenAbGroupElem)
+         diag = A_to_A2(a)
+         exps = GapObj([diag[i] for i in 1:n], recursive = true)
+         return group_element(G, GAPWrap.PcpElementByExponentsNC(C, exps))
+       end
+
+       finv = g -> A2_to_A(A2(_exponent_vector(g)))
+     else
+#TODO: As soon as https://github.com/gap-packages/polycyclic/issues/88 is fixed,
+#      we can change the code to create a `PcpGroup` in GAP also if the group
+#      is finite.
+       # Create a `PcGroup` in GAP.
+       # The generators of this group correspond to those of `A2`
+       # but they are in general only a subset of the defining pcgs.
+       # We implement the decomposition of an element into the generators
+       # via a matrix multiplication with the exponent vector of the element.
+       GapG = GAPWrap.AbelianGroup(GAP.Globals.IsPcGroup, GapObj(exponents; recursive = true))
+       Gpcgs = GAPWrap.FamilyPcgs(GapG)
+       G = PcGroup(GAPWrap.SubgroupNC(GapG, Gpcgs))
+
+       # Find the "intervals" in the pcgs that correspond to the generators.
+       indepgens = Vector{GapObj}(GAPWrap.GeneratorsOfGroup(GapG))
+       pcgs = Vector{GapObj}(Gpcgs)
+       m = length(pcgs)
+       relord = GAPWrap.RelativeOrders(Gpcgs)
+       starts = [findfirst(isequal(x), pcgs) for x in indepgens]
+       push!(starts, length(pcgs)+1)
+       M = zero_matrix(ZZ, m, n)
+       for j in 1:length(indepgens)
+         e = 1
+         for i in starts[j]:(starts[j+1]-1)
+           M[i, j] = e
+           e = e * relord[i]
+         end
+       end
+       starts = starts[1:n]
+
+       f = function(a::FinGenAbGroupElem)
+         diag = A_to_A2(a)
+         v = zeros(ZZ, m)
+         v[starts] = [diag[i] for i in 1:n]
+         exps = GapObj(v, recursive = true)
+         return group_element(G, GAPWrap.LinearCombinationPcgs(Gpcgs, GapObj(v, true)))
+       end
+
+       finv = g -> A2_to_A(A2(_exponent_vector(g) * M))
+     end
+
+     return GroupIsomorphismFromFunc(A, G, f, finv)
+   end::GroupIsomorphismFromFunc{FinGenAbGroup, PcGroup}
+end
+
+function isomorphism(::Type{FPGroup}, A::FinGenAbGroup)
+   # Known isomorphisms are cached in the attribute `:isomorphisms`.
+   isos = get_attribute!(Dict{Tuple{Type, Bool}, Any}, A, :isomorphisms)::Dict{Tuple{Type, Bool}, Any}
+   return get!(isos, (FPGroup, false)) do
+      # Do not call `abelian_group(FPGroup, ...)`,
+      # because then we would need an indirection via a group with
+      # diagonal relations.
+      # Instead, we create a group with the same defining relations.
+      n = ngens(A)
+      G = free_group(n; eltype = :syllable)
+      R = rels(A)
+      gG = gens(G)
+      gGi = map(inv, gG)
+      s = vcat(elem_type(G)[gG[i]*gG[j]*gGi[i]*gGi[j] for i in 1:n for j in (i+1):n],
+           elem_type(G)[prod([gen(G, i)^R[j,i] for i=1:n if !iszero(R[j,i])], init = one(G)) for j=1:nrows(R)])
+      F, mF = quo(G, s)
+      set_is_abelian(F, true)
+      set_is_finite(F, is_finite(A))
+      is_finite(A) && set_order(F, order(A))
+      return MapFromFunc(
+        A, F,
+        y->F([i => y[i] for i=1:n]),
+        x->A(exponents_of_abelianization(x)))
+   end::MapFromFunc{FinGenAbGroup, FPGroup}
+end
+
 ####
 mutable struct GroupIsomorphismFromFunc{R, T} <: Map{R, T, Hecke.HeckeMap, MapFromFunc}
     map::MapFromFunc{R, T}
@@ -879,31 +978,6 @@ end
 
 ####
 
-# We need not find independent generators in order to create
-# a presentation of a fin. gen. abelian group.
-function isomorphism(::Type{FPGroup}, A::FinGenAbGroup)
-   # Known isomorphisms are cached in the attribute `:isomorphisms`.
-   isos = get_attribute!(Dict{Tuple{Type, Bool}, Any}, A, :isomorphisms)::Dict{Tuple{Type, Bool}, Any}
-   return get!(isos, (FPGroup, false)) do
-      n = ngens(A)
-      G = free_group(n; eltype = :syllable)
-      R = rels(A)
-      gG = gens(G)
-      gGi = map(inv, gG)
-      s = vcat(elem_type(G)[gG[i]*gG[j]*gGi[i]*gGi[j] for i in 1:n for j in (i+1):n],
-           elem_type(G)[prod([gen(G, i)^R[j,i] for i=1:n if !iszero(R[j,i])], init = one(G)) for j=1:nrows(R)])
-      F, mF = quo(G, s)
-      set_is_abelian(F, true)
-      set_is_finite(F, is_finite(A))
-      is_finite(A) && set_order(F, order(A))
-      return MapFromFunc(
-        A, F,
-        y->F([i => y[i] for i=1:n]),
-        x->sum([w.second*gen(A, w.first) for w = syllables(x)], init = zero(A)))
-   end::MapFromFunc{FinGenAbGroup, FPGroup}
-end
-
-
 # We need not find independent generators in order to create
 # a presentation of a `FPModule` over a finite field.
 # Note that additively, the given module is an elementary abelian p-group