Speed up orbits of permutation groups on integers (#4307)

Co-authored-by: ThomasBreuer <[email protected]>

examples/GSets_timing.jl

@@ -0,0 +1,42 @@
+# "G-sets of permutation groups"
+
+
+  # natural constructions (determined by the types of the seeds)
+  G = symmetric_group(10)
+
+  Omega = gset(G)
+  @benchmark order(stabilizer(Omega, 1)[1])
+  @benchmark order(stabilizer(Omega, Set([1, 2]))[1])
+  @benchmark order(stabilizer(Omega, [1, 2])[1])
+  @benchmark order(stabilizer(Omega, (1, 2))[1])
+
+  Omega = gset(G, [Set([1, 2])])  # action on unordered pairs
+  @benchmark order(stabilizer(Omega, Set([1, 2]))[1])
+  @benchmark order(stabilizer(Omega, Set([Set([1, 2]), Set([1, 3])]))[1])
+  @benchmark order(stabilizer(Omega, [Set([1, 2]), Set([1, 3])])[1])
+  @benchmark order(stabilizer(Omega, (Set([1, 2]), Set([1, 3])))[1])
+
+
+  Omega = gset(G, [[1, 2]])  # action on ordered pairs
+  @benchmark order(stabilizer(Omega, [1, 2])[1])
+  @benchmark order(stabilizer(Omega, Set([[1, 2], [1, 3]]))[1])
+  @benchmark order(stabilizer(Omega, [[1, 2], [1, 3]])[1])
+  @benchmark order(stabilizer(Omega, ([1, 2], [1, 3]))[1])
+
+
+  Omega = gset(G, [(1, 2)])  # action on ordered pairs (repres. by tuples)
+  @benchmark order(stabilizer(Omega, (1, 2))[1])
+  @benchmark order(stabilizer(Omega, Set([(1, 2), (1, 3)]))[1])
+  @benchmark order(stabilizer(Omega, [(1, 2), (1, 3)])[1])
+  @benchmark order(stabilizer(Omega, ((1, 2), (1, 3)))[1])
+
+
+  # constructions by explicit action functions
+  G = symmetric_group(6)
+  omega = [0,1,0,1,0,1]
+  Omega = gset(G, permuted, [omega, [1,2,3,4,5,6]])
+
+  @benchmark order(stabilizer(Omega, omega)[1])
+  @benchmark order(stabilizer(Omega, Set([omega, [1,0,0,1,0,1]]))[1])
+  @benchmark order(stabilizer(Omega, [omega, [1,0,0,1,0,1]])[1])
+  @benchmark order(stabilizer(Omega, (omega, [1,0,0,1,0,1]))[1])

src/Groups/action.jl

@@ -526,7 +526,9 @@ julia> S = stabilizer(G, [1, 1, 2, 2, 3], permuted);  order(S[1])
 4
 ```
 """
-function stabilizer(G::GAPGroup, pnt::Any, actfun::Function)
+stabilizer(G::GAPGroup, pnt::Any, actfun::Function) = _stabilizer_generic(G, pnt, actfun)
+
+function _stabilizer_generic(G::GAPGroup, pnt::Any, actfun::Function)
     return Oscar._as_subgroup(G, GAPWrap.Stabilizer(GapObj(G), pnt,
         GapObj(gens(G), recursive = true), GapObj(gens(G)),
         GapObj(actfun)))
@@ -535,7 +537,10 @@ end
 # natural stabilizers in permutation groups
 # Construct the arguments on the GAP side such that GAP's method selection
 # can choose the special method.
-function stabilizer(G::PermGroup, pnt::T) where T <: Oscar.IntegerUnion
+# - stabilizer in a perm. group of an integer via `^`
+# - stabilizer in a perm. group of a vector of integers via `on_tuples`
+# - stabilizer in a perm. group of a set of integers via `on_sets`
+function stabilizer(G::PermGroup, pnt::T) where T <: IntegerUnion
     return Oscar._as_subgroup(G, GAPWrap.Stabilizer(GapObj(G),
         GapObj(pnt),
         GAP.Globals.OnPoints))  # Do not use GAPWrap.OnPoints!
@@ -553,6 +558,20 @@ function stabilizer(G::PermGroup, pnt::AbstractSet{T}) where T <: Oscar.IntegerU
         GAP.Globals.OnSets))  # Do not use GAPWrap.OnSets!
 end
 
+# now the same with given action function,
+# these calls may come from delegations from G-sets
+function stabilizer(G::PermGroup, pnt::T, actfun::Function) where T <: IntegerUnion
+    return (actfun == ^) ? stabilizer(G, pnt) : _stabilizer_generic(G, pnt, actfun)
+end
+
+function stabilizer(G::PermGroup, pnt::Vector{T}, actfun::Function) where T <: IntegerUnion
+    return actfun == on_tuples ? stabilizer(G, pnt) : _stabilizer_generic(G, pnt, actfun)
+end
+
+function stabilizer(G::PermGroup, pnt::AbstractSet{T}, actfun::Function) where T <: IntegerUnion
+    return actfun == on_sets ? stabilizer(G, pnt) : _stabilizer_generic(G, pnt, actfun)
+end
+
 # natural stabilizers in matrix groups
 stabilizer(G::MatrixGroup{ET,MT}, pnt::AbstractAlgebra.Generic.FreeModuleElem{ET}) where {ET,MT} = stabilizer(G, pnt, *)
 

src/Groups/gsets.jl

@@ -332,6 +332,17 @@ orbit(G::GAPGroup, omega) = gset_by_type(G, [omega], typeof(omega))
 
 orbit(G::Union{GAPGroup, FinGenAbGroup}, fun::Function, omega) = GSetByElements(G, fun, [omega])
 
+
+function gap_action_function(Omega::GSet)
+  f = action_function(Omega)
+  (f == ^) && return GAP.Globals.OnPoints
+  f == on_tuples && return GAP.Globals.OnTuples
+  f == on_sets && return GAP.Globals.OnSets
+  # etc.
+  return GapObj(f) # generic fallback
+end
+
+
 """
     orbit(Omega::GSet, omega)
 
@@ -349,7 +360,11 @@ julia> length(orbit(Omega, 1))
 4
 ```
 """
-function orbit(Omega::GSetByElements{<:GAPGroup, S}, omega::S) where S
+orbit(Omega::GSetByElements{<:GAPGroup, S}, omega::S) where S = _orbit_generic(Omega, omega)
+
+function _orbit_generic(Omega::GSetByElements{<:GAPGroup, S}, omega::S) where S
+    # In this generic function, we delegate the loop to GAP, but we act
+    # with Julia group elements on Julia objects via Julia functions.
     G = acting_group(Omega)
     acts = GapObj(gens(G))
     gfun = GapObj(action_function(Omega))
@@ -366,6 +381,38 @@ function orbit(Omega::GSetByElements{<:GAPGroup, S}, omega::S) where S
 end
 #T check whether omega lies in Omega?
 
+# special cases where we convert the objects to GAP
+# (the group elements as well as the objects they act on),
+# in order to use better methods on the GAP side:
+# - orbit of a perm. group on integers via `^`
+# - orbit of a perm. group on vectors of integers via `on_tuples`
+# - orbit of a perm. group on sets of integers via `on_sets`
+function orbit(Omega::GSetByElements{PermGroup, S}, omega::S) where S <: IntegerUnion
+    (action_function(Omega) == ^) || return _orbit_generic(Omega, omega)
+    return _orbit_special_GAP(Omega, omega)
+end
+
+function orbit(Omega::GSetByElements{PermGroup, S}, omega::S) where S <: Vector{<: IntegerUnion}
+    action_function(Omega) == on_tuples || return _orbit_generic(Omega, omega)
+    return _orbit_special_GAP(Omega, omega)
+end
+
+function orbit(Omega::GSetByElements{PermGroup, S}, omega::S) where S <: Set{<: IntegerUnion}
+    action_function(Omega) == on_sets || return _orbit_generic(Omega, omega)
+    return _orbit_special_GAP(Omega, omega)
+end
+
+function _orbit_special_GAP(Omega::GSetByElements{<:GAPGroup, S}, omega::S) where S
+    G = acting_group(Omega)
+    gfun = gap_action_function(Omega)
+    orb = Vector{S}(GAP.Globals.Orbit(GapObj(G), GapObj(omega), gfun)::GapObj)
+
+    res = as_gset(acting_group(Omega), action_function(Omega), orb)
+    # We know that this G-set is transitive.
+    set_attribute!(res, :orbits => [orb])
+    return res
+end
+
 function orbit(Omega::GSetByElements{FinGenAbGroup}, omega::T) where T
     return orbit_via_Julia(Omega, omega)
 end

test/Groups/gsets.jl

@@ -365,6 +365,17 @@ end
   f = x^2 + y
   orb = orbit(G, f)
   @test length(orb) == 3
+
+  F = QQBarField()
+  e = one(F)
+  s, c = sincospi(2 * e / 3)
+  mat_rot = matrix([c -s; s c])
+  G = matrix_group(mat_rot)
+  p = F.([1, 0])
+  orb = orbit(G, *, p)
+  @test length(orb) == 3
+
+
 end
 
 @testset "G-sets by right transversals" begin