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 5c6c48c..bc55b5d 100644 --- a/gap/basic/pcppara.gi +++ b/gap/basic/pcppara.gi @@ -73,13 +73,15 @@ 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 . +## Returns the positions changed in and in . ## -BindGlobal( "ReduceExpoPara", function ( ind, gen, indd, pgen, rel ) - local i, j, a, b, q, cg; +BindGlobal( "ReduceExpoPara", +function ( ind, gen, indd, pgen, rel, pows, ppows ) + 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 @@ -89,8 +91,10 @@ 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; od; @@ -98,14 +102,36 @@ 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; fi; od; - return cg; + 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 ); ############################################################################# @@ -114,11 +140,13 @@ end ); ## ## 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, 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, 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; @@ -127,9 +155,9 @@ function( pcs, gens, ppcs, pgens ) rels := RelativeOrders( coll ); n := NumberOfGenerators( coll ); inf := 0 in rels; - c := n+1; + bnd := fail; - # 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,29 +165,53 @@ function( pcs, gens, ppcs, pgens ) ind[d] := pcs[i]; indd[d] := ppcs[i]; 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 @@ -173,6 +225,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 +240,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 +252,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,12 +269,20 @@ 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, pows, ppows); + 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; # add powers and commutators + # no need to clean up pows/ppows here added := false; for d in f do g := ind[d]; @@ -234,26 +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, ind[j]^-1); - if Depth(k) < c then - Add(todo, k); - Add(tododo, Comm(gg, indd[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 b283a2f..15b82fc 100644 --- a/gap/basic/pcppcps.gi +++ b/gap/basic/pcppcps.gi @@ -25,37 +25,68 @@ end ); ############################################################################# ## -#F ReduceExpo +#F CommutatorDepthBounds( ) ## -## Returns the positions changed in . +## 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 . ## -BindGlobal( "ReduceExpo", function( ind, gen, rel ) - local i, j, a, b, q, cg; +## 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; - 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; - 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); + 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; - od; - fi; + fi; + row[t] := Minimum( k, prev[t], row[t+1] ); + od; + prev := row; + bnd[s] := norm{ row }; od; - return cg; + + coll![PC_COMMUTATOR_DEPTH_BOUNDS] := bnd; + return bnd; end ); ############################################################################# @@ -161,13 +192,111 @@ 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 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 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, g, d, todo, val, j, f, h, a, b, k, - pair, t, added, oldc, chg, i; + 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; @@ -176,25 +305,39 @@ InstallGlobalFunction(AddToIgs, function(igs, gens) rels := RelativeOrders(coll); n := NumberOfGenerators(coll); inf := 0 in rels; - c := n+1; + bnd := fail; - # 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; - - # 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); + pows := []; + 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 := []; @@ -204,6 +347,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 +359,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 +383,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, pows); + 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; @@ -257,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, ind[j]^-1); - if Depth(k) < c then - Add(todo, k); - Add(val, IGSValFun(k)); - added := true; - fi; - fi; fi; od; od; @@ -412,7 +562,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) 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);