Speed up IsPowerfulPGroup

The speed up is achieved by carefully caching some more
intermediate results. Also avoid an Agemo call, and drop a size
check that I thought was beneficial, but in benchmarking harder
examples turned out to be not helpful.

Some timeing on my M1 MacBook Pro. Before this patch:

    gap> IsRegularPGroup(DirectProduct(SmallGroup(3^5,22),SmallGroup(3^5,22)));; time;
    27609
    gap> IsRegularPGroup(DirectProduct(SmallGroup(3^5,22),SmallGroup(3^5,39)));; time;
    77479

After this patch:

    gap> IsRegularPGroup(DirectProduct(SmallGroup(3^5,22),SmallGroup(3^5,22)));; time;
    1488
    gap> IsRegularPGroup(DirectProduct(SmallGroup(3^5,22),SmallGroup(3^5,39)));; time;
    3957

Also fix some references to Huppert's book.

lib/grp.gi

@@ -433,43 +433,44 @@ InstallMethod( IsPowerfulPGroup,
 InstallMethod( IsRegularPGroup,
     [ IsGroup ],
 function( G )
-local p, hom, reps, a, b, ap_bp, ab_p, H;
+local p, hom, reps, as, bs, a, b, ap, bp, ab, ap_bp, ab_p, g, h, H;
 
   if not IsPGroup(G) then
     return false;
   fi;
 
   p:=PrimePGroup(G);
   if p = 2 then
-    # see [Hup67, Satz 10.3 a)]
+    # see [Hup67, Satz III.10.3 a)]
     return IsAbelian(G);
   elif p = 3 and DerivedLength(G) > 2 then
-    # see [Hup67, Satz 10.3 b)]
+    # see [Hup67, Satz III.10.3 b)]
     return false;
   elif Size(G) <= p^p then
-    # see [Hal34, Corollary 14.14], [Hall, p. 183], [Hup67, Satz 10.2 b)]
+    # see [Hal34, Corollary 14.14], [Hall, p. 183], [Hup67, Satz III.10.2 b)]
     return true;
   elif NilpotencyClassOfGroup(G) < p then
-    # see [Hal34, Corollary 14.13], [Hall, p. 183], [Hup67, Satz 10.2 a)]
+    # see [Hal34, Corollary 14.13], [Hall, p. 183], [Hup67, Satz III.10.2 a)]
     return true;
   elif IsCyclic(DerivedSubgroup(G)) then
-    # see [Hup67, Satz 10.2 c)]
+    # see [Hup67, Satz III.10.2 c)]
     return true;
   elif Exponent(G) = p then
-    # see [Hup67, Satz 10.2 d)]
+    # see [Hup67, Satz III.10.2 d)]
     return true;
   elif p = 3 and RankPGroup(G) = 2 then
     # see [Hup67, Satz 10.3 b)]: at this point we know that the derived
     # subgroup is not cyclic, hence G is not regular
     return false;
   elif Size(G) < p^p * Size(Agemo(G,p)) then
-    # see [Hal36, Theorem 2.3], [Hup67, Satz 10.13]
+    # see [Hal36, Theorem 2.3], [Hup67, Satz III.10.13]
     return true;
   elif Index(DerivedSubgroup(G),Agemo(DerivedSubgroup(G),p)) < p^(p-1) then
-    # see [Hal36, Theorem 2.3], [Hup67, Satz 10.13]
+    # see [Hal36, Theorem 2.3], [Hup67, Satz III.10.13]
     return true;
   fi;
 
+
   # Fallback to actually check the defining criterion, i.e.:
   # for all a,b in G, we must have that a^p*b^p/(a*b)^p in (<a,b>')^p
 
@@ -483,27 +484,35 @@ local p, hom, reps, a, b, ap_bp, ab_p, H;
   reps := Filtered(reps, g -> not IsOne(g));
   reps := List(reps, g -> PreImagesRepresentative(hom, g));
 
-  for b in Image(hom) do
-    b := PreImagesRepresentative(hom, b);
-    for a in reps do
+  as := List(reps, a -> [a,a^p]);;
+  bs := List(Image(hom), function(b)
+      b := PreImagesRepresentative(hom, b);
+      return [b,b^p];
+      end);;
+
+  for b in bs do
+    bp := b[2]; b := b[1];
+    for a in as do
+      ap := a[2]; a := a[1];
       # if a and b commute the regularity condition automatically holds
-      if a*b = b*a then continue; fi;
+      ab := a*b;
+      if ab = b*a then continue; fi;
 
       # regularity is also automatic if a^p * b^p = (a*b)^p
-      ap_bp := a^p * b^p;
-      ab_p := (a*b)^p;
+      ap_bp := ap * bp;
+      ab_p := ab^p;
       if ap_bp = ab_p then continue; fi;
 
       # if the subgroup generated H by a and b is itself regular, we are also
-      # done. However we don't use recursion, here, as H may be equal to G;
-      # and also we have to be careful to not use too expensive code here.
-      # But a quick size check is certainly fine.
+      # done. However we don't use recursion here, as it is too expensive.
+      # we just check the direct definition, with a minor twist to avoid
+      # calling Agemo
+      g := ap_bp / ab_p;
+      h := Comm(a,b)^p;
+      if g = h or IsOne(g*h) then continue; fi;
       H := Subgroup(G, [a,b]);
-      if Size(H) <= p^p then continue; fi;
-
-      # finally the full check
-      H := DerivedSubgroup(H);
-      if not (ap_bp / ab_p) in Agemo(H, p) then
+      H := NormalClosure(H, [h]); # faster than Agemo(DerivedSubgroup(H), p);
+      if not g in H then
         return false;
       fi;
     od;