From 80b43bca10d7acc0aea52df8586e3b91c82f9369 Mon Sep 17 00:00:00 2001 From: stertooy <5571903+stertooy@users.noreply.github.com> Date: Thu, 22 Oct 2026 13:27:05 +0200 Subject: [PATCH 1/3] Cache inverse powers in AddToIgs(Parallel) --- gap/basic/pcppara.gi | 34 +++++++++++++++++------ gap/basic/pcppcps.gi | 65 ++++++++++++++++++++++++++++++++++++++------ 2 files changed, 82 insertions(+), 17 deletions(-) diff --git a/gap/basic/pcppara.gi b/gap/basic/pcppara.gi index 5c6c48c..17dcf2b 100644 --- a/gap/basic/pcppara.gi +++ b/gap/basic/pcppara.gi @@ -76,10 +76,11 @@ end ); #F ReduceExpoPara( ind, gen, indd, pgen, rel ) ## ## Parallel version of ReduceExpo. -## Returns the positions changed in . +## Returns the positions changed in and in . ## BindGlobal( "ReduceExpoPara", function ( ind, gen, indd, pgen, rel ) - local i, j, a, b, q, cg; + local i, j, a, b, q, ci, cg; + ci := []; cg := []; for i in [ 1 .. Length( ind ) ] do if not IsBool( ind[i] ) and rel[i] = 0 then @@ -91,6 +92,7 @@ BindGlobal( "ReduceExpoPara", function ( ind, gen, indd, pgen, rel ) if q <> 0 then ind[j] := ind[j] * ind[i] ^ (- q); indd[j] := indd[j] * indd[i] ^ (- q); + AddSet( ci, j ); fi; fi; od; @@ -105,7 +107,7 @@ BindGlobal( "ReduceExpoPara", function ( ind, gen, indd, pgen, rel ) od; fi; od; - return cg; + return [ci, cg]; end ); ############################################################################# @@ -117,8 +119,9 @@ end ); ## InstallGlobalFunction( AddToIgsParallel, function( pcs, gens, ppcs, pgens ) - local coll, rels, n, inf, todo, tododo, ind, indd, g, gg, d, h, hh, a, b, - q, k, c, i, r, sub, keep, val, j, f, nrmd, pairs, added, oldc, chg; + local coll, rels, n, inf, todo, tododo, ind, indd, pows, ppows, g, gg, + d, h, hh, a, b, q, k, c, i, r, sub, keep, val, j, f, nrmd, pairs, + added, oldc, chg, t; if Length( gens ) = 0 then return [pcs, ppcs]; fi; @@ -129,7 +132,7 @@ function( pcs, gens, ppcs, pgens ) inf := 0 in rels; c := n+1; - # set up + # set up; pows[d] and ppows[d] cache negative powers of ind[d], indd[d] ind := ListWithIdenticalEntries(n, false); indd := ListWithIdenticalEntries(n, false); for i in [1..Length(pcs)] do @@ -137,6 +140,8 @@ function( pcs, gens, ppcs, pgens ) ind[d] := pcs[i]; indd[d] := ppcs[i]; od; + pows := []; + ppows := []; # initial gens still need to be reduced added := true; @@ -173,6 +178,8 @@ function( pcs, gens, ppcs, pgens ) nrmd := NormedPcpElementPara( g, gg ); ind[d] := nrmd[1]; indd[d] := nrmd[2]; + Unbind(pows[d]); + Unbind(ppows[d]); AddSet(f,d); h := ind[d]; hh := indd[d]; @@ -186,8 +193,8 @@ function( pcs, gens, ppcs, pgens ) if a > 0 and b > 0 and a mod b = 0 then # GcdPcpPara would leave h and hh unchanged. q := SmallestResidue(a/b, rels[d]); - g := g * h^-q; - gg := gg * hh^-q; + g := g * IgsNegativePower(ind, pows, d, q); + gg := gg * IgsNegativePower(indd, ppows, d, q); else pairs := GcdPcpPara(g, h, gg, hh); h := pairs[1]; @@ -198,6 +205,8 @@ function( pcs, gens, ppcs, pgens ) nrmd := NormedPcpElementPara( h, hh ); ind[d] := nrmd[1]; indd[d] := nrmd[2]; + Unbind(pows[d]); + Unbind(ppows[d]); AddSet(f, d); fi; fi; @@ -213,7 +222,14 @@ function( pcs, gens, ppcs, pgens ) oldc := c; c := TailLimit(ind, c); chg := []; - if inf then chg := ReduceExpoPara(ind, todo, indd, tododo, rels); fi; + if inf then + t := ReduceExpoPara(ind, todo, indd, tododo, rels); + for i in t[1] do + Unbind(pows[i]); + Unbind(ppows[i]); + od; + chg := t[2]; + fi; for i in chg do if Depth(todo[i]) < c then val[i] := IGSValFun(todo[i]); fi; od; diff --git a/gap/basic/pcppcps.gi b/gap/basic/pcppcps.gi index b283a2f..5bd4c4c 100644 --- a/gap/basic/pcppcps.gi +++ b/gap/basic/pcppcps.gi @@ -27,11 +27,12 @@ end ); ## #F ReduceExpo ## -## Returns the positions changed in . +## Returns the positions changed in and in . ## BindGlobal( "ReduceExpo", function( ind, gen, rel ) - local i, j, a, b, q, cg; + local i, j, a, b, q, ci, cg; + ci := []; cg := []; for i in [1..Length(ind)] do if not IsBool(ind[i]) and rel[i]=0 then @@ -42,6 +43,7 @@ BindGlobal( "ReduceExpo", function( ind, gen, rel ) q := QuoInt(a,b); if q <> 0 then ind[j] := ind[j]*ind[i]^-q; + AddSet(ci, j); fi; fi; od; @@ -55,7 +57,7 @@ BindGlobal( "ReduceExpo", function( ind, gen, rel ) od; fi; od; - return cg; + return [ci, cg]; end ); ############################################################################# @@ -161,13 +163,52 @@ BindGlobal( "GcdPcp", function(g, h) return [x, y]; end ); +############################################################################# +## +#F IgsNegativePower( , , , ) . . . . . . . . . . ind[d]^-q +## +## Caches the result in pows[d], which must be reset when ind[d] changes. +## +BindGlobal( "IgsNegativePower", function( ind, pows, d, q ) + local b, p, r, s, value; + + # cache bound, avoids making cache consist of large sparse lists + b := 16; + + # for positive powers we do nothing special + if q < 0 then return ind[d]^-q; fi; + if not IsBound(pows[d]) then pows[d] := []; fi; + p := pows[d]; + + # store the inverse of ind[d] + if not IsBound(p[1]) then p[1] := ind[d]^-1; fi; + + # attempt to retrieve from cache + if q <= b and IsBound(p[q]) then return p[q]; fi; + + # reuse a cached power whose exponent divides q + r := Minimum(Length(p), QuoInt(q, 2)); + s := q; + while r > 1 do + if IsBound(p[r]) and q mod r = 0 then + s := QuoInt(q, r); + break; + fi; + r := r - 1; + od; + + value := p[r]^s; + if q <= b then p[q] := value; fi; + return value; +end ); + ############################################################################# ## #F AddToIgs( , ) ## InstallGlobalFunction(AddToIgs, function(igs, gens) - local coll, rels, n, inf, c, ind, g, d, todo, val, j, f, h, a, b, k, - pair, t, added, oldc, chg, i; + local coll, rels, n, inf, c, ind, pows, g, d, todo, val, i, j, f, h, + a, b, q, k, pair, t, added, oldc, chg; if Length(gens) = 0 then return igs; fi; @@ -178,9 +219,10 @@ InstallGlobalFunction(AddToIgs, function(igs, gens) inf := 0 in rels; c := n+1; - # set up + # set up; pows[d] caches negative powers of ind[d] ind := ListWithIdenticalEntries(n, false); for g in igs do ind[Depth(g)] := g; od; + pows := []; # initial gens still need to be reduced added := true; @@ -204,6 +246,7 @@ InstallGlobalFunction(AddToIgs, function(igs, gens) h := ind[d]; if IsBool(h) then ind[d] := NormedPcpElement(g); + Unbind(pows[d]); AddSet(f, d); h := ind[d]; fi; @@ -215,13 +258,15 @@ InstallGlobalFunction(AddToIgs, function(igs, gens) b := LeadingExponent(h); if a > 0 and b > 0 and a mod b = 0 then # GcdPcp would leave h unchanged. - g := g * h^-SmallestResidue(a/b, rels[d]); + q := SmallestResidue(a/b, rels[d]); + g := g * IgsNegativePower(ind, pows, d, q); else pair := GcdPcp(g, h); h := pair[1]; g := pair[2]; if h <> ind[d] then ind[d] := NormedPcpElement( h ); + Unbind(pows[d]); AddSet(f, d); fi; fi; @@ -237,7 +282,11 @@ InstallGlobalFunction(AddToIgs, function(igs, gens) oldc := c; c := TailLimit(ind, c); chg := []; - if inf then chg := ReduceExpo(ind, todo, rels); fi; + if inf then + t := ReduceExpo(ind, todo, rels); + for i in t[1] do Unbind(pows[i]); od; + chg := t[2]; + fi; for i in chg do if Depth(todo[i]) < c then val[i] := IGSValFun(todo[i]); fi; od; From b77a4caee27e84dcbc7b8d59a1d9d62cfd862732 Mon Sep 17 00:00:00 2001 From: stertooy <5571903+stertooy@users.noreply.github.com> Date: Thu, 22 Oct 2026 17:02:01 +0200 Subject: [PATCH 2/3] Use IgsNegativePower in more places --- gap/basic/pcppara.gi | 23 ++++++++----- gap/basic/pcppcps.gi | 80 ++++++++++++++++++++++---------------------- 2 files changed, 54 insertions(+), 49 deletions(-) diff --git a/gap/basic/pcppara.gi b/gap/basic/pcppara.gi index 17dcf2b..4503573 100644 --- a/gap/basic/pcppara.gi +++ b/gap/basic/pcppara.gi @@ -73,12 +73,13 @@ end ); ############################################################################# ## -#F ReduceExpoPara( ind, gen, indd, pgen, rel ) +#F ReduceExpoPara( ind, gen, indd, pgen, rel, pows, ppows ) ## ## Parallel version of ReduceExpo. ## Returns the positions changed in and in . ## -BindGlobal( "ReduceExpoPara", function ( ind, gen, indd, pgen, rel ) +BindGlobal( "ReduceExpoPara", +function ( ind, gen, indd, pgen, rel, pows, ppows ) local i, j, a, b, q, ci, cg; ci := []; cg := []; @@ -90,8 +91,9 @@ BindGlobal( "ReduceExpoPara", function ( ind, gen, indd, pgen, rel ) a := Exponents( ind[j] )[i]; q := QuoInt( a, b ); if q <> 0 then - ind[j] := ind[j] * ind[i] ^ (- q); - indd[j] := indd[j] * indd[i] ^ (- q); + ind[j] := ind[j] * IgsNegativePower(ind, pows, i, q); + indd[j] := indd[j] * + IgsNegativePower(indd, ppows, i, q); AddSet( ci, j ); fi; fi; @@ -100,8 +102,9 @@ BindGlobal( "ReduceExpoPara", function ( ind, gen, indd, pgen, rel ) a := Exponents( gen[j] )[i]; q := QuoInt( a, b ); if q <> 0 then - gen[j] := gen[j] * ind[i] ^ (- q); - pgen[j] := pgen[j] * indd[i] ^ (- q); + gen[j] := gen[j] * IgsNegativePower(ind, pows, i, q); + pgen[j] := pgen[j] * + IgsNegativePower(indd, ppows, i, q); AddSet( cg, j ); fi; od; @@ -223,7 +226,7 @@ function( pcs, gens, ppcs, pgens ) c := TailLimit(ind, c); chg := []; if inf then - t := ReduceExpoPara(ind, todo, indd, tododo, rels); + t := ReduceExpoPara(ind, todo, indd, tododo, rels, pows, ppows); for i in t[1] do Unbind(pows[i]); Unbind(ppows[i]); @@ -235,6 +238,7 @@ function( pcs, gens, ppcs, pgens ) od; # add powers and commutators + # no need to clean up pows/ppows here added := false; for d in f do g := ind[d]; @@ -262,10 +266,11 @@ function( pcs, gens, ppcs, pgens ) added := true; fi; if rels[j] = 0 then - k := Comm(g, ind[j]^-1); + k := Comm(g, IgsNegativePower(ind, pows, j, 1)); if Depth(k) < c then Add(todo, k); - Add(tododo, Comm(gg, indd[j]^-1)); + Add(tododo, Comm(gg, + IgsNegativePower(indd, ppows, j, 1))); Add(val, IGSValFun(k)); added := true; fi; diff --git a/gap/basic/pcppcps.gi b/gap/basic/pcppcps.gi index 5bd4c4c..1623cc0 100644 --- a/gap/basic/pcppcps.gi +++ b/gap/basic/pcppcps.gi @@ -23,43 +23,6 @@ BindGlobal( "TailLimit", function( ind, c ) return i+1; end ); -############################################################################# -## -#F ReduceExpo -## -## Returns the positions changed in and in . -## -BindGlobal( "ReduceExpo", function( ind, gen, rel ) - local i, j, a, b, q, ci, cg; - - ci := []; - cg := []; - for i in [1..Length(ind)] do - if not IsBool(ind[i]) and rel[i]=0 then - b := LeadingExponent(ind[i]); - for j in [1..i-1] do - if not IsBool( ind[j] ) then - a := Exponents(ind[j])[i]; - q := QuoInt(a,b); - if q <> 0 then - ind[j] := ind[j]*ind[i]^-q; - AddSet(ci, j); - fi; - fi; - od; - for j in [1..Length(gen)] do - a := Exponents(gen[j])[i]; - q := QuoInt(a,b); - if q <> 0 then - gen[j] := gen[j]*ind[i]^-q; - AddSet(cg, j); - fi; - od; - fi; - od; - return [ci, cg]; -end ); - ############################################################################# ## #F CheckIgs( [, ] ) @@ -202,6 +165,43 @@ BindGlobal( "IgsNegativePower", function( ind, pows, d, q ) return value; end ); +############################################################################# +## +#F ReduceExpo( , , , ) +## +## Returns the positions changed in and in . +## +BindGlobal( "ReduceExpo", function( ind, gen, rel, pows ) + local i, j, a, b, q, ci, cg; + + ci := []; + cg := []; + for i in [1..Length(ind)] do + if not IsBool(ind[i]) and rel[i]=0 then + b := LeadingExponent(ind[i]); + for j in [1..i-1] do + if not IsBool( ind[j] ) then + a := Exponents(ind[j])[i]; + q := QuoInt(a,b); + if q <> 0 then + ind[j] := ind[j]*IgsNegativePower(ind, pows, i, q); + AddSet(ci, j); + fi; + fi; + od; + for j in [1..Length(gen)] do + a := Exponents(gen[j])[i]; + q := QuoInt(a,b); + if q <> 0 then + gen[j] := gen[j]*IgsNegativePower(ind, pows, i, q); + AddSet(cg, j); + fi; + od; + fi; + od; + return [ci, cg]; +end ); + ############################################################################# ## #F AddToIgs( , ) @@ -283,7 +283,7 @@ InstallGlobalFunction(AddToIgs, function(igs, gens) c := TailLimit(ind, c); chg := []; if inf then - t := ReduceExpo(ind, todo, rels); + t := ReduceExpo(ind, todo, rels, pows); for i in t[1] do Unbind(pows[i]); od; chg := t[2]; fi; @@ -317,7 +317,7 @@ InstallGlobalFunction(AddToIgs, function(igs, gens) added := true; fi; if rels[j] = 0 then - k := Comm(g, ind[j]^-1); + k := Comm(g, IgsNegativePower(ind, pows, j, 1)); if Depth(k) < c then Add(todo, k); Add(val, IGSValFun(k)); @@ -461,7 +461,7 @@ BindGlobal( "AddIgsToIgs", function( pcs1, pcs2 ) d := Depth( g ); od; c := TailLimit(ind, c); - if inf then ReduceExpo(ind, todo, rels); fi; + if inf then ReduceExpo(ind, todo, rels, []); fi; todo := Filtered(todo, x -> Depth(x) Date: Fri, 23 Oct 2026 01:17:39 +0200 Subject: [PATCH 3/3] Speed up AddToIgs and AddToIgsParallel Most of the time went into commutators that were computed only to be discarded because they already lay in the tail of the igs: for UT(20,5) with 1000 random elements, about 15,000 of 16,000. A table of lower bounds for commutator depths, derived from the presentation once per collector and cached in it, lets the loop skip those without computing them. The generators are now taken one at a time, each closed under powers and commutators before the next is sifted, so the tail limit drops as early as possible. In groups with infinite relative orders a generator is reduced only just before it is sifted instead of every time the igs changes; with many generators that reduction dominated the running time. In infinite groups with few generators the running time depends strongly on processing order, and the new order is slower on some inputs. Setup: gap> ut := function(n, p, k) local G; G := UnitriangularPcpGroup(n, p); Reset(GlobalMersenneTwister, 1); return List([1..k], i -> Random(G)); end;; gap> dt := function(n, k) local G; G := UnitriangularPcpGroup(n, 0); AddHallPolynomials(Collector(G)); Reset(GlobalMersenneTwister, 1); return List([1..k], i -> Random(G)); end;; gap> sift := function(gens) local t; t := Runtime(); AddToIgs([], gens); return Runtime() - t; end;; Before: gap> sift(ut(20, 5, 30)); 19245 gap> sift(ut(20, 5, 1000)); 18643 gap> sift(ut(25, 5, 1000)); 163026 gap> sift(ut(6, 0, 1000)); 200 gap> sift(dt(10, 1000)); 13987 After: gap> sift(ut(20, 5, 30)); 4424 gap> sift(ut(20, 5, 1000)); 7095 gap> sift(ut(25, 5, 1000)); 42009 gap> sift(ut(6, 0, 1000)); 7 gap> sift(dt(10, 1000)); 410 Assisted-by: Generative AI --- CHANGES.md | 2 + gap/basic/collect.gd | 2 + gap/basic/collect.gi | 1 + gap/basic/pcppara.gi | 121 ++++++++++++++++++++++----------- gap/basic/pcppcps.gi | 157 +++++++++++++++++++++++++++++++++++-------- tst/AddToIgs.tst | 28 ++++++++ 6 files changed, 245 insertions(+), 66 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index be458a3..9553c0e 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -48,6 +48,8 @@ This file describes changes in the GAP package 'polycyclic'. - Allow a finite group fallback for `RepresentativeAction` and `IsConjugate` for subgroups - Optimise `NormalIntersection` and `EfaSeries` +- Speed up `AddToIgs` and `AddToIgsParallel`, in particular for many + generators and for large groups - Avoid needlessly large powers when norming and sifting elements at a depth of finite relative order; this also changes which power `NormedPcpElement` returns diff --git a/gap/basic/collect.gd b/gap/basic/collect.gd index 70a0bd2..d9c77de 100644 --- a/gap/basic/collect.gd +++ b/gap/basic/collect.gd @@ -128,3 +128,5 @@ BindGlobal( "PC_INVERSECOMMUTATORSINVERSE", 27 ); BindGlobal( "PC_NILPOTENT_COMMUTE", 28 ); BindGlobal( "PC_WEIGHTS", 29 ); BindGlobal( "PC_ABELIAN_START", 30 ); + +BindGlobal( "PC_COMMUTATOR_DEPTH_BOUNDS", 31 ); diff --git a/gap/basic/collect.gi b/gap/basic/collect.gi index 5753cfd..706d410 100644 --- a/gap/basic/collect.gi +++ b/gap/basic/collect.gi @@ -940,6 +940,7 @@ function( coll ) fi; FromTheLeftCollector_SetCommute( coll ); + Unbind( coll![PC_COMMUTATOR_DEPTH_BOUNDS] ); ## We have to declare the collector up to date now because the following ## functions need to collect and are careful enough. diff --git a/gap/basic/pcppara.gi b/gap/basic/pcppara.gi index 4503573..bc55b5d 100644 --- a/gap/basic/pcppara.gi +++ b/gap/basic/pcppara.gi @@ -113,18 +113,40 @@ function ( ind, gen, indd, pgen, rel, pows, ppows ) return [ci, cg]; end ); +############################################################################# +## +#F ReduceExpoElmPara( ind, g, indd, gg, rel, pows, ppows ) +## +## Parallel version of ReduceExpoElm. +## +BindGlobal( "ReduceExpoElmPara", +function ( ind, g, indd, gg, rel, pows, ppows ) + local i, q; + for i in [ Depth( g ) .. Length( ind ) ] do + if not IsBool( ind[i] ) and rel[i] = 0 then + q := QuoInt( Exponents( g )[i], LeadingExponent( ind[i] ) ); + if q <> 0 then + g := g * IgsNegativePower(ind, pows, i, q); + gg := gg * IgsNegativePower(indd, ppows, i, q); + fi; + fi; + od; + return [ g, gg ]; +end ); + ############################################################################# ## #F AddToIgsParallel( , , , ) ## ## This function adds the elements in to the induced pcs . ## It acts simultaneously on and as well as and . +## See AddToIgs for the algorithm. ## InstallGlobalFunction( AddToIgsParallel, function( pcs, gens, ppcs, pgens ) - local coll, rels, n, inf, todo, tododo, ind, indd, pows, ppows, g, gg, - d, h, hh, a, b, q, k, c, i, r, sub, keep, val, j, f, nrmd, pairs, - added, oldc, chg, t; + local coll, rels, n, inf, bnd, todo, tododo, ind, indd, pows, ppows, g, + gg, d, h, hh, a, b, q, k, c, i, r, sub, keep, idx, queue, pqueue, + qpos, val, j, f, nrmd, pairs, added, oldc, chg, t; if Length( gens ) = 0 then return [pcs, ppcs]; fi; @@ -133,7 +155,7 @@ function( pcs, gens, ppcs, pgens ) rels := RelativeOrders( coll ); n := NumberOfGenerators( coll ); inf := 0 in rels; - c := n+1; + bnd := fail; # set up; pows[d] and ppows[d] cache negative powers of ind[d], indd[d] ind := ListWithIdenticalEntries(n, false); @@ -145,29 +167,51 @@ function( pcs, gens, ppcs, pgens ) od; pows := []; ppows := []; + c := TailLimit(ind, n+1); - # initial gens still need to be reduced - added := true; - - # do a reduction step - c := TailLimit(ind, c); + # gens wait in queue, keeping the first of equal elements and sorted by + # IGSValFun; todo holds powers and commutators sub := Filtered([1..Length(gens)], i -> Depth(gens[i]) < c); - # Keep the first occurrence of each generator, in input order. - todo := gens{sub}; - StableSortParallel(todo, sub); - keep := Filtered([1..Length(sub)], i -> i = 1 or todo[i] <> todo[i-1]); + queue := gens{sub}; + StableSortParallel(queue, sub); + keep := Filtered([1..Length(sub)], i -> i = 1 or queue[i] <> queue[i-1]); sub := Set(sub{keep}); - todo := gens{sub}; - tododo := pgens{sub}; - val := List(todo, IGSValFun); - - # loop over to-do list until it is empty - while Length( todo ) > 0 and c > 1 do - j := PositionMinimum(val); - g := Remove(todo, j); - gg := Remove(tododo, j); - Remove(val, j); - d := Depth( g ); + queue := gens{sub}; + pqueue := pgens{sub}; + idx := [1..Length(sub)]; + StableSortParallel(List(queue, IGSValFun), idx); + queue := queue{idx}; + pqueue := pqueue{idx}; + qpos := 1; + todo := []; + tododo := []; + val := []; + added := false; + + while c > 1 do + + # take the next element, from queue only once todo is empty + if Length( todo ) > 0 then + j := PositionMinimum(val); + g := Remove(todo, j); + gg := Remove(tododo, j); + Remove(val, j); + else + while qpos <= Length(queue) and Depth(queue[qpos]) >= c do + qpos := qpos+1; + od; + if qpos > Length(queue) then break; fi; + g := queue[qpos]; + gg := pqueue[qpos]; + qpos := qpos+1; + # sifting g unreduced can blow up its exponents + if inf then + pairs := ReduceExpoElmPara(ind, g, indd, gg, rels, pows, ppows); + g := pairs[1]; + gg := pairs[2]; + fi; + fi; + d := Depth( g ); f := []; # shift g into ind @@ -254,27 +298,28 @@ function( pcs, gens, ppcs, pgens ) fi; fi; for j in [1..n] do - if j = d or Minimum( d, j ) >= c-1 then + # skip trivial commutators and those in the tail + if j = d or Minimum( d, j ) >= c-1 or IsBool(ind[j]) then continue; fi; - if not IsBool(ind[j]) then - k := Comm(g, ind[j]); + if bnd = fail then bnd := CommutatorDepthBounds(coll); fi; + if bnd[Minimum(d, j)][Maximum(d, j)] >= c then continue; fi; + k := Comm(g, ind[j]); + if Depth(k) < c then + Add(todo, k); + Add(tododo, Comm(gg, indd[j])); + Add(val, IGSValFun(k)); + added := true; + fi; + if rels[j] = 0 then + k := Comm(g, IgsNegativePower(ind, pows, j, 1)); if Depth(k) < c then Add(todo, k); - Add(tododo, Comm(gg, indd[j])); + Add(tododo, Comm(gg, + IgsNegativePower(indd, ppows, j, 1))); Add(val, IGSValFun(k)); added := true; fi; - if rels[j] = 0 then - k := Comm(g, IgsNegativePower(ind, pows, j, 1)); - if Depth(k) < c then - Add(todo, k); - Add(tododo, Comm(gg, - IgsNegativePower(indd, ppows, j, 1))); - Add(val, IGSValFun(k)); - added := true; - fi; - fi; fi; od; od; diff --git a/gap/basic/pcppcps.gi b/gap/basic/pcppcps.gi index 1623cc0..15b82fc 100644 --- a/gap/basic/pcppcps.gi +++ b/gap/basic/pcppcps.gi @@ -23,6 +23,72 @@ BindGlobal( "TailLimit", function( ind, c ) return i+1; end ); +############################################################################# +## +#F CommutatorDepthBounds( ) +## +## Returns a table B such that Comm(x, y) has depth at least B[d][e] whenever +## d < e, Depth(x) >= d and Depth(y) >= e. The table is cached in . +## +## Modulo the normal closure N of all [g_s, g_t] with s >= d and t >= e, the +## subgroups G_d and G_e commute, so Comm(x, y) lies in N. If these [g_s, g_t] +## lie in G_k, then N lies in the largest normal G_j with j <= k. +## +BindGlobal( "CommutatorDepthBounds", function( coll ) + local n, conj, low, norm, prev, row, bnd, w, k, s, t; + + if IsBound( coll![PC_COMMUTATOR_DEPTH_BOUNDS] ) then + return coll![PC_COMMUTATOR_DEPTH_BOUNDS]; + fi; + + n := NumberOfGenerators( coll ); + conj := coll![PC_CONJUGATES]; + + # low[t]: least depth of g_u^g_s over s < u, t <= u. Thus G_t is normal + # if and only if low[t] = t. + low := [1..n+1]; + for t in [2..n] do + for s in [1..t-1] do + if IsBound( conj[t][s] ) and conj[t][s][1] < low[t] then + low[t] := conj[t][s][1]; + fi; + od; + od; + for t in [n-1, n-2 .. 1] do low[t] := Minimum( low[t], low[t+1] ); od; + + # norm[k]: the largest j <= k with G_j normal + norm := []; + for k in [1..n+1] do + t := k; + while low[t] < t do t := low[t]; od; + norm[k] := t; + od; + + # row[t] for row s: least depth of [g_u, g_v] over s <= v < u, t <= u + bnd := []; + prev := ListWithIdenticalEntries( n+1, n+1 ); + for s in [n, n-1 .. 1] do + row := ListWithIdenticalEntries( n+1, n+1 ); + for t in [n, n-1 .. s+1] do + k := n+1; + if IsBound( conj[t][s] ) then + w := conj[t][s]; + if w[1] <> t or w[2] <> 1 then + k := w[1]; + elif Length( w ) > 2 then + k := w[3]; + fi; + fi; + row[t] := Minimum( k, prev[t], row[t+1] ); + od; + prev := row; + bnd[s] := norm{ row }; + od; + + coll![PC_COMMUTATOR_DEPTH_BOUNDS] := bnd; + return bnd; +end ); + ############################################################################# ## #F CheckIgs( [, ] ) @@ -202,13 +268,35 @@ BindGlobal( "ReduceExpo", function( ind, gen, rel, pows ) return [ci, cg]; end ); +############################################################################# +## +#F ReduceExpoElm( , , , ) . . . ReduceExpo for one element +## +BindGlobal( "ReduceExpoElm", function( ind, g, rel, pows ) + local i, q; + for i in [Depth(g)..Length(ind)] do + if not IsBool(ind[i]) and rel[i]=0 then + q := QuoInt(Exponents(g)[i], LeadingExponent(ind[i])); + if q <> 0 then g := g*IgsNegativePower(ind, pows, i, q); fi; + fi; + od; + return g; +end ); + ############################################################################# ## #F AddToIgs( , ) ## +## Sifts the elements of into one at a time. After each, the +## powers and commutators of the changed entries are sifted in until ind is +## closed under them again. The entries from c on have leading exponent 1 +## and so generate the tail G_c of the pc series: elements of depth c or +## more need no sifting, and the loop stops at c = 1, often long before +## all of is used. +## InstallGlobalFunction(AddToIgs, function(igs, gens) - local coll, rels, n, inf, c, ind, pows, g, d, todo, val, i, j, f, h, - a, b, q, k, pair, t, added, oldc, chg; + local coll, rels, n, inf, bnd, c, ind, pows, queue, qpos, g, d, todo, + val, i, j, f, h, a, b, q, k, pair, t, added, oldc, chg; if Length(gens) = 0 then return igs; fi; @@ -217,26 +305,39 @@ InstallGlobalFunction(AddToIgs, function(igs, gens) rels := RelativeOrders(coll); n := NumberOfGenerators(coll); inf := 0 in rels; - c := n+1; + bnd := fail; # set up; pows[d] caches negative powers of ind[d] ind := ListWithIdenticalEntries(n, false); for g in igs do ind[Depth(g)] := g; od; pows := []; - - # initial gens still need to be reduced - added := true; - - # do a reduction step - c := TailLimit(ind, c); - todo := Set(Filtered(gens, x -> Depth(x) < c)); - val := List(todo, IGSValFun); - - # loop over to-do list until it is empty - while Length(todo) > 0 and c > 1 do - j := PositionMinimum(val); - g := Remove(todo, j); - Remove(val, j); + c := TailLimit(ind, n+1); + + # gens wait in queue; todo holds powers and commutators + queue := Set(Filtered(gens, x -> Depth(x) < c)); + StableSortParallel(List(queue, IGSValFun), queue); + qpos := 1; + todo := []; + val := []; + added := false; + + while c > 1 do + + # take the next element, from queue only once todo is empty + if Length(todo) > 0 then + j := PositionMinimum(val); + g := Remove(todo, j); + Remove(val, j); + else + while qpos <= Length(queue) and Depth(queue[qpos]) >= c do + qpos := qpos+1; + od; + if qpos > Length(queue) then break; fi; + g := queue[qpos]; + qpos := qpos+1; + # sifting g unreduced can blow up its exponents + if inf then g := ReduceExpoElm(ind, g, rels, pows); fi; + fi; d := Depth(g); f := []; @@ -306,24 +407,24 @@ InstallGlobalFunction(AddToIgs, function(igs, gens) fi; for j in [1..n] do # skip trivial commutators and those in the tail - if j = d or Minimum( d, j ) >= c-1 then + if j = d or Minimum( d, j ) >= c-1 or IsBool(ind[j]) then continue; fi; - if not IsBool(ind[j]) then - k := Comm(g, ind[j]); + if bnd = fail then bnd := CommutatorDepthBounds(coll); fi; + if bnd[Minimum(d, j)][Maximum(d, j)] >= c then continue; fi; + k := Comm(g, ind[j]); + if Depth(k) < c then + Add(todo, k); + Add(val, IGSValFun(k)); + added := true; + fi; + if rels[j] = 0 then + k := Comm(g, IgsNegativePower(ind, pows, j, 1)); if Depth(k) < c then Add(todo, k); Add(val, IGSValFun(k)); added := true; fi; - if rels[j] = 0 then - k := Comm(g, IgsNegativePower(ind, pows, j, 1)); - if Depth(k) < c then - Add(todo, k); - Add(val, IGSValFun(k)); - added := true; - fi; - fi; fi; od; od; diff --git a/tst/AddToIgs.tst b/tst/AddToIgs.tst index fee0af8..bff6ee0 100644 --- a/tst/AddToIgs.tst +++ b/tst/AddToIgs.tst @@ -216,5 +216,33 @@ gap> v := S.2 * S.1^-2 * S.3 * S.2^2;; gap> Index(S, Subgroup(S, [u, v])); 1 +# +# Many generators, most of them redundant +# +gap> Reset(GlobalMersenneTwister, 1);; +gap> G := UnitriangularPcpGroup(10, 5);; +gap> gens := List([1..200], i -> Random(G));; +gap> igs := AddToIgs([], gens);; +gap> Cgs(igs) = Cgs(G); +true +gap> res := AddToIgsParallel([], gens, [], List(gens, g -> g^G.1));; +gap> Cgs(res[1]) = Cgs(G); +true +gap> res[2] = List(res[1], g -> g^G.1); +true +gap> G := UnitriangularPcpGroup(5, 0);; +gap> gens := List([1..50], i -> Random(G));; +gap> igs := AddToIgs([], gens);; +gap> Cgs(igs) = Cgs(G); +true +gap> CheckIgs(igs, gens); +true +gap> gens := [G.1^2*G.5, G.2^3*G.4^-1, G.3^2*G.1];; +gap> igs := AddToIgs([], gens);; +gap> CheckIgs(igs, gens); +true +gap> Index(G, SubgroupByIgs(G, igs)); +infinity + # gap> STOP_TEST( "AddToIgs.tst", 1);