Speed up computation of the intersection of group cosets in the generic cases and even more so for permutation groups

lib/csetgrp.gi

@@ -789,6 +789,48 @@ local a,r;
 end);
 
 
+InstallMethod(Intersection2, "general cosets", IsIdenticalObj,
+              [IsRightCoset,IsRightCoset],
+function(cos1,cos2)
+    local swap, H1, H2, x1, x2, sigma, U, rho;
+    if Size(cos1) < 10 then
+        TryNextMethod();
+    elif Size(cos2) < 10 then
+        return Intersection2(cos2, cos1);
+    fi;
+    if Size(cos1) > Size(cos2) then
+        swap := cos1;
+        cos1 := cos2;
+        cos2 := swap;
+    fi;
+    H1:=ActingDomain(cos1);
+    H2:=ActingDomain(cos2);
+    x1:=Representative(cos1);
+    x2:=Representative(cos2);
+    sigma := x1 / x2;
+    if Size(H1) = Size(H2) and H1 = H2 then
+        if sigma in H1 then
+            return cos1;
+        else
+            return [];
+        fi;
+    fi;
+    # We want to compute the intersection of cos1 = H1*x1 with cos2 = H2*x2.
+    # This is equivalent to intersecting H1 with H2*x2/x1, which is either empty
+    # or equal to a coset U*rho, where U is the intersection of H1 and H2.
+    # In the non-empty case, the overall result then is U*rho*x1.
+    #
+    # To find U*rho, we iterate over all cosets of U in H1 and for each test
+    # if it is contained in H2*x2/x1, which is the case if and only if rho is
+    # in H2*x2/x1, if and only if rho/(x2/x1) = rho*x1/x2 is in H2
+    U:=Intersection(H1, H2);
+    for rho in RightTransversal(H1, U) do
+        if rho * sigma in H2 then
+            return RightCoset(U, rho * x1);
+        fi;
+    od;
+    return [];
+end);
 
 # disabled because of comparison incompatibilities
 #InstallMethod(\<,"RightCosets",IsIdenticalObj,[IsRightCoset,IsRightCoset],0,

lib/csetperm.gi

@@ -548,3 +548,190 @@ function(a,b)
   return CanonicalRepresentativeOfExternalSet(a)
          <CanonicalRepresentativeOfExternalSet(b);
 end);
+
+
+
+InstallMethod(Intersection2, "perm cosets", IsIdenticalObj,
+  [IsRightCoset and IsPermCollection,IsRightCoset and IsPermCollection],0,
+function(cos1,cos2)
+    local H1, H2, x1, x2, shift, sigma, listMoved_H1, listMoved_H2,
+          listMoved_sigma, U, repr, H1_sigma, H2_sigma, H12, swap, rho, diff;
+    # We set cosInt = cos1 \cap cos2 = H1 x1 \cap H2 x2
+    H1:=ActingDomain(cos1);
+    H2:=ActingDomain(cos2);
+    x1:=Representative(cos1);
+    x2:=Representative(cos2);
+    if H1=H2 then
+        if cos1=cos2 then
+            return cos1;
+        else
+            return [];
+        fi;
+    fi;
+    # We are using that
+    #    H1*x1 \cap H2*x2 = (H1 \cap H2*x2/x1)*x1 = (H1 \cap H2*sigma)*shift,
+    # where shift and sigma are defined as below:
+    shift:=x1;
+    sigma:=x2 / x1;
+
+    # Reducing as much as possible in advance by using various relatively
+    # cheap to compute criteria
+    while true do
+        listMoved_H1:=MovedPoints(H1);
+        listMoved_H2:=MovedPoints(H2);
+        listMoved_sigma:=MovedPoints(sigma);
+
+        # If the coset intersection is non-empty, then there is h1 \in H1 and
+        # h2 \in H2 such that h1 = h2*sigma. Therefore sigma is contained in
+        # the group generated by H1 and H2. A necessary condition for this is
+        # that the points moved by sigma are a subset of the points moved by
+        # H1 and H2.
+        if not IsSubset(Union(listMoved_H1, listMoved_H2), listMoved_sigma) then
+            return [];
+        fi;
+
+        # Suppose x is an element of H1 \cap H2*sigma. Then for any positive
+        # integer n, we know that n^x is contained in n^H1 but also in
+        # (n^H2)^\sigma. Thus if the intersection of n^H1 and (n^H2)^\sigma is
+        # empty, then the intersection of the cosets is also empty. Clearly
+        # the orbit intersection contains n whenever n is fixed by sigma, so
+        # we only have to consider this for n moved by sigma.
+        if ForAny(listMoved_sigma, n -> IsEmpty(Intersection(Orbit(H1,n), OnTuples(Orbit(H2,n),sigma)))) then
+            return [];
+        fi;
+
+        # If there are points that are moved by sigma and by H2 but not by H1,
+        # then there must be an element in H2 which matches the action of sigma
+        # on these points, or else the intersection is empty
+        diff:=Difference(Intersection(listMoved_H2, listMoved_sigma), listMoved_H1);
+        if Length(diff) > 0 then
+            repr:=RepresentativeAction(H2, diff, OnTuples(diff, Inverse(sigma)), OnTuples);
+            if repr=fail then
+                return [];
+            fi;
+            # Since repr is in H2, we can replace sigma by repr*sigma
+            # without changing the coset H2*sigma. The new sigma then fixes
+            # all points in diff, as does H1. Hence replacing H2 by the
+            # stabilizer in H2 of diff does not change H1 \cap H2 sigma.
+            sigma:=repr * sigma;
+            H2:=Stabilizer(H2, diff, OnTuples);
+            continue;
+        fi;
+
+        # Mirror to the previous check:
+        # If there are points that are moved by sigma and by H1 but not by H2,
+        # then there must be an element in H1 which matches the action of sigma
+        # on these points, or else the intersection is empty
+        diff:=Difference(Intersection(listMoved_H1, listMoved_sigma), listMoved_H2);
+        if Length(diff) > 0 then
+            repr:=RepresentativeAction(H1, diff, OnTuples(diff, sigma), OnTuples);
+            if repr=fail then
+                return [];
+            fi;
+            # Again this is similar to the case before, except that we are adjusting
+            # H1 now and thus also need to take `shift` into account. The situation
+            # is as follows:
+            #
+            # cosInt = (H1 \cap H2 sigma) shift
+            #        = (H1 sigma^{-1} \cap H2) sigma shift
+            #        = (Stab_{H1}(diff) repr sigma^{-1} \cap H2) sigma shift
+            #        = (Stab_{H1}(diff) \cap H2 sigma repr^{-1} ) repr shift
+            H1:=Stabilizer(H1, diff, OnTuples);
+            sigma:=sigma / repr;
+            shift:=repr * shift;
+            continue;
+        fi;
+
+        # easy termination criterion: reduction to group intersection
+        if sigma in H2 then
+            U:=Intersection(H1, H2);
+            return RightCoset(U, shift);
+        fi;
+
+        # easy termination criterion: reduction to group intersection
+        if sigma in H1 then
+            # cosInt = (H1 \cap H2 sigma) shift
+            #        = (H1 sigma^{-1} \cap H2) sigma shift
+            U:=Intersection(H1, H2);
+            return RightCoset(U, sigma * shift);
+        fi;
+
+        # any element of H1 which moves points not moved by sigma or anything
+        # in H2 can not be in the intersection, so we may as well remove them
+        # by a stabilizer computation
+        diff:=Difference(listMoved_H1, Union(listMoved_H2, listMoved_sigma));
+        if Length(diff) > 0 then
+            H1:=Stabilizer(H1, diff, OnTuples);
+            continue;
+        fi;
+
+        # the same but with the roles of H1 and H2 reversed
+        diff:=Difference(listMoved_H2, Union(listMoved_H1, listMoved_sigma));
+        if Length(diff) > 0 then
+            H2:=Stabilizer(H2, diff, OnTuples);
+            continue;
+        fi;
+
+        # More general but more expensive than previous check
+        H1_sigma:=ClosureGroup(H1, sigma);
+        if not IsSubgroup(H1_sigma, H2) then
+            H2:=Intersection(H1_sigma, H2);
+            continue;
+        fi;
+        H2_sigma:=ClosureGroup(H2, sigma);
+        if not IsSubgroup(H2_sigma, H1) then
+            H1:=Intersection(H2_sigma, H1);
+            continue;
+        fi;
+        # No more reduction tricks available
+        break;
+    od;
+
+    # A final termination criterion
+    H12:=ClosureGroup(H1, H2);
+    if not sigma in H12 then
+        return [];
+    fi;
+
+    # We are now inspired by the algorithm from
+    # Lazlo Babai, Coset Intersection in Moderately Exponential Time
+    #
+    # We use the algorithm from Page 10 of coset analysis and we reformulate
+    # it here in order to avoid errors:
+    # --- The naive algorithm for computing H1 \cap H2 sigma is to iterate
+    # over elements of H1 and testing if one belongs to H2 sigma. If we find
+    # one such z then the result is the coset RightCoset(U, z). If not
+    # then it is empty.
+    # --- Since the result is independent of the cosets U, what we can
+    # do is iterate over the RightCosets(H1, U). The algorithm is the
+    # one of Proposition 3.2
+    # for r in RightCosets(H1, U) do
+    #   if r in H1*sigma then
+    #     return RightCoset(U, r * shift)
+    #   fi;
+    # od;
+    # --- (TODO for future work): The question is how to make it faster.
+    # One idea is to use an ascending chain between U and H1.
+    # Section 3.4 of above paper gives statement related to that but not a
+    # useful algorithm. The question deserves further exploration.
+    #
+    # We select the smallest group for that computation in order to have
+    # as few cosets as possible
+    if Order(H2) < Order(H1) then
+        # cosInt = (H1 \cap H2 sigma) shift
+        #        = (H1 sigma^{-1} \cap H2) sigma shift
+        swap:=H1;
+        H1:=H2;
+        H2:=swap;
+        shift:=sigma * shift;
+        sigma:=Inverse(sigma);
+    fi;
+    # So now Order(H1) <= Order(H2)
+    U:=Intersection(H1, H2);
+    for rho in RightTransversal(H1, U) do
+        if rho / sigma in H2 then
+            return RightCoset(U, rho * shift);
+        fi;
+    od;
+    return [];
+end);

tst/testextra/grpperm.tst

@@ -99,19 +99,6 @@ gap> IsNaturalAlternatingGroup(g);
 true
 gap> 2*Size(g) = Factorial(999);
 true
-gap> Intersection(SymmetricGroup([1..5]),SymmetricGroup([3..8]));
-Sym( [ 3 .. 5 ] )
-gap> Intersection(SymmetricGroup([1..5]),AlternatingGroup([3..8]));
-Alt( [ 3 .. 5 ] )
-gap> Intersection(AlternatingGroup([1..5]),AlternatingGroup([3..8]));
-Alt( [ 3 .. 5 ] )
-gap> Intersection(AlternatingGroup([1..5]),SymmetricGroup([3..8]));
-Alt( [ 3 .. 5 ] )
-gap> Intersection(Group( (1,2,3), (4,5,6), (11,12,13), (11,12,14) ),
-                  Group( (1,7,8), (4,9,11), (10,12), (13,14) ));
-Group(())
-gap> Intersection(Group((1,2), (3,4)), Group((3,4),(5,6)));
-Group( (3,4) )
 gap> s := SymmetricGroup(100);
 Sym( [ 1 .. 100 ] )
 gap> Stabilizer(s,3,OnPoints);