Implement `Intersection(U,V)` for the case were neither subgroup normalizes the other

[fingolfin]

A good starting point might be to look at the "Handbook of computational group theory" section 8.8.1 which deals with finite polycyclic groups. For the general case and also a bit more details, look at section 8.4 of Bettina Eick's Habilitation thesis.

For the infinite case, one work along an efa series / "linear series". Then, besides elementary abelian factors, we also get free abelian factors. To compute the required stabilizers there, one could use StabilizerIntegralAction.

And of course when adding this feature, let's add many, many test cases for it :-)

[stertooy]

I've made an attempt at implementing this over the holiday. The code is still horribly unoptimised(*) and contains extra checks and a few Print commands for testing purposes, but it mostly(**) works, I think.

IntersectionNew := function( U, H )
    local G, L, n, A, p, UA_A, HA_A, UA_HA_A, gens, imgs, i, U_HA, K, HA, p1, A_H, p2, q, H_A_H, h_a_h, ha_a,
        secisom, secisominv, l, k, delta, x, px, py, y, z, A_A_H, T, r, ZZ, rl, rq,  mats, e, oper, prei, orbstabinfo, S;

    ###
    # Eliminate trivial cases
    ###

    # Remark 8.6: if either group normalises the other, use NormalIntersection instead
    if IsNormal( U, H ) then
        Print("H normalizes U\n");
        return NormalIntersection( H, U );
    elif IsNormal( H, U ) then
        Print("U normalizes H\n");
        return NormalIntersection( U, H );
    fi;

    ###
    # 8.4 Intersection of subgroups
    ###

    # Create parent group and find suitable abelian normal subgroup
    G := ClosureGroup( U, H );
    L := EfaSeries( G );
    n := Length( L );
    A := L[n-1];
    Print("Efa series has length ",n,"\n");

    # Calculate UA/A and HA/A
    p := NaturalHomomorphismByNormalSubgroup( G, A );
    UA_A := ImagesSet( p, U ); # UA/A
    HA_A := ImagesSet( p, H ); # HA/A

    # Induction step: calculate UA/A ∩ HA/A = (UA ∩ HA)/A (subgroups of GA)
    Print("Starting induction step\n");
    UA_HA_A := IntersectionNew( UA_A, HA_A );
    Print("Ending induction step\n");

    ###
    # 8.4.1 Determining U ∩ HA from UA ∩ HA
    ###

    # Calculate U ∩ HA as preimage of (UA ∩ HA)/A
    gens := SmallGeneratingSet( U );
    imgs := List( gens, u -> ImagesRepresentative( p, u ) );
    i := GroupHomomorphismByImages( U, UA_A, gens, imgs );
    U_HA := PreImagesSet( i, UA_HA_A ); # U ∩ HA

    # Correctness test
    HA := ClosureGroup( H, A );
    for x in GeneratorsOfGroup( U_HA ) do
        if ( not x in U ) or ( not x in HA ) then
            Error("Incorrect intersection U ∩ HA \n");
        fi;
    od;

    ###
    # 8.4.2 Determining U ∩ H from U ∩ HA
    ###

    K := U_HA;
    k := Pcp( K );

    # Preliminary work to find h in H, a in A s.t. k = h*a for given k in K
    p1 := RestrictedMapping( p, HA ); # HA -> HA/A
    A_H := NormalIntersection( A, H ); # A is normal so OK
    p2 := NaturalHomomorphismByNormalSubgroup( H, A_H ); # H -> H/(A ∩ H)
    q := NaturalHomomorphismByNormalSubgroup( A, A_H ); # A -> A/(A ∩ H)
    H_A_H := ImagesSource( p2 ); # H/(A ∩ H )
    h_a_h := SmallGeneratingSet( H_A_H );
    ha_a := List( h_a_h, x -> ImagesRepresentative( p1, PreImagesRepresentative( p2, x ) ) );
    secisom := GroupHomomorphismByImages( H_A_H, HA_A, h_a_h, ha_a ); # 2nd isom. thm. H/(A ∩ H) -> HA/A
    secisominv := InverseGeneralMapping( secisom ); # 2nd isom. thm. HA/A -> H/(A ∩ H)

    # Create derivation
    delta := function( x )
        local px, py, y, z;
        # Decompose x in K as x = y * z with y in H, z in A
        px := ImagesRepresentative( p1, x );
        py := ImagesRepresentative( secisominv, px );
        y := PreImagesRepresentative( p2, py );
        z := y^-1*x;
        # z(A ∩ H) is image of x under derivation delta
        return ImagesRepresentative( q, z );
    end;

    # Work inside A/(A ∩ H)
    A_A_H := ImagesSource( q ); # A/(A ∩ H)
    T := TorsionSubgroup( A_A_H );

    # If A/(A ∩ H) is infinite, more work is required
    if not IsFinite( A_A_H ) then
        ###
        # 7.5 Mod out torsion subgroup to get free abelian group (7.5)
        ###
        r := NaturalHomomorphismByNormalSubgroup( A_A_H, T );
        ZZ := ImagesSource( r );

        # Images of generators of z under derivation and quotient
        z := IndependentGeneratorsOfAbelianGroup( ZZ ); # Don't use PCP, may not be independent
        l := List( GeneratorsOfPcp( k ), x -> delta( x ) );
        rl := List( l, x -> ImagesRepresentative( r, x ) );
        rq := q*r;

        ###
        # Lemma 7.1
        ###

        # Create matrices corresponding to affine action
        mats := List( [1..Length(k)], i ->
            Concatenation(
                List( [1..Length(z)], j -> 
                    Concatenation( IndependentGeneratorExponents( 
                        ZZ, 
                        ImagesRepresentative( rq, PreImagesRepresentative( rq, z[j] )^k[i] )
                    ), [0] )
                ),
                [ Concatenation( IndependentGeneratorExponents( ZZ, rl[i] ), [1] ) ]
            )
        );

        # Calculate stabiliser of affine action, replace K
        Print("Calculating stabiliser on free abelian quotient\n");
        e := Concatenation( ListWithIdenticalEntries( Length( z ), 0 ), [1] );
        K := StabilizerIntegralAction( K, mats, e );
        k := Pcp( K );
        Print("Calculated stabiliser on free abelian quotient\n");
    fi;

    if IsTrivial( T ) then
        # Nothing to calculate
        S := K;
    else
        # According to 7.5, orbit of 0+(A ∩ H) under action is now finite
        Print("Calculating stabiliser with finite orbit\n");
        oper := function( pnt, g )
            local prei;
            prei := PreImagesRepresentative( q, pnt );
            return ImagesRepresentative( q, prei^g )*delta( g );
        end;
        orbstabinfo := PcpOrbitStabilizer( One( A_A_H ), k, k, oper );
        Print("Calculated stabiliser with finite orbit\n");
        S := Subgroup( G, orbstabinfo!.stab );
    fi;

    # Test if calculated stabiliser is, at the very least, a subgroup of the actual intersection
    for x in GeneratorsOfGroup( S ) do
        if not (x in U and x in H) then
            Error("Intersection wrong\n");
        fi;
    od;
    return S;

end;

(*) For optimisation, I can think at least the following things would improve runtime:

• Remove coding related to testing
• Turn functions into their NC variant
• Do the final stabiliser calculation without invoking delta so much, it's probably doable with only the images of the generators of K
• Catch more "trivial" cases, e.g. A contained in either U or H
• Maybe switch the role of U and H when this would lead to an improvement

(*) When I say the code mostly* works, I mean that I do encounter the following:

• StabilizerIntegralAction runs forever, because of #76
• StabilizerIntegralAction or PcpOrbitStabilizer can take a very long time to finish
• Incorrect IGS calculations (Error, igs is incorrect [...]), which is unrelated to this code but does pop up eventually

In case it's of any use, I've mostly been testing with the following code (for different values of n):

gap> CHECK_CENT@Polycyclic := true;;
gap> CHECK_IGS@Polycyclic := true;;
gap> CHECK_INTNORM@Polycyclic := true;;
gap> CHECK_INTSTAB@Polycyclic := true;;
gap> CHECK_NORM@Polycyclic := true;;
gap> SetInfoLevel(InfoIntStab,4);;

gap> n := 8;;
gap> G := ExamplesOfSomePcpGroups( n );;
gap> U := Subgroup(G,List([1..2],i->Random(G)));;
gap> H := Subgroup(G,List([1..2],i->Random(G)));;
gap> IntersectionNew( U, H ) = IntersectionNew( H, U );

[fingolfin]

By coincidence, a student of mine also implemented this just a few days before you added this comment (which happened while I was on vacation, hence the late reply).

Ideally I would like to review both his and your code. And then finish my new AddToIgs code, review your other PRs, and make a new release.... In reality, I don't know when I'll get to that, given that I also need to do a ton of other things :-(.