Add Frank Lübeck SL2 recognition algorithm

Author: Till Müller

gap/projective/constructive_recognition/SL/BaseCase.gi

@@ -475,3 +475,359 @@ RECOG.FindStdGensUsingBSGS := function(g,stdgens,projective,large)
   od;
   return SLPOfElms(gens);
 end;
+
+
+###################################################################################################
+###################################################################################################
+######## Algorithm provided by Frank Lübeck #######################################################
+###################################################################################################
+###################################################################################################
+
+
+
+####################   (c) 2025  Frank Lübeck   ##########################
+##  
+##  G must be SL(2,q) generated by 2x2 matrices over GF(q).
+##  Returns [slps, mat] where slps is list of SLPs in given generators of G 
+##  to elements which are conjugate by mat to standard/nice generators: 
+##     t, u1, u2 = diag(1/Z(q),Z(q)) [[1,0],[1,1]], [[1,1],[0,1]]
+##  
+##  We construct elements with the strategy of Section 2. from the article
+##    Conder, Leedham-Green, Fast recognition of classical groups
+##    over large fields, Groups and Computation III, 2001.
+##  and then adjust them to the form described above.
+##  
+##  Some comments:
+##    - To find short SLP to the nice generators we avoid 'PseudoRandom(g)'
+##      and instead just start with the trivial element and multiply random
+##      generators to it.
+##      This should be good enough because in both cases (element xm and twice
+##      element cm in the code) almost half of the elements in G are suitable.
+##  
+##    - The list of SLPs in the result is probably what is needed for the
+##      SLPforNiceGens for the RecogNode.
+##  
+##    - To find for a given element A in G an SLP in the nice generators compute
+##       Astd := A^mat and use the function 'SLPinStandardSL2' from the file 
+##       "SL2.g".
+##    
+##    - If q = 2^m with odd m then the computation of the eigenvalues of xm
+##      needs the quadratic extension GF(2^(2m)).
+##  
+##  
+##  Example (takes about 15 seconds on my notebook):
+##      Size(SL(2,23^15));
+##      G:=Group(List([1..4],i->PseudoRandom(SL(2,23^15))));
+##      l:=RecogNaturalSL2(G,23^15);
+##      nicegens:=List(l[1],a->ResultOfStraightLineProgram(a,GeneratorsOfGroup(G)));
+##      List(last, a->a^l[2]);
+##  
+RecogNaturalSL2 := function(G, q)
+  local GM, one, zero, qm1fac, c, m, gens, xm, x, pol, v, z, exp, a, 
+        mat, tm, ym, y, ymat, tr, d, cm, r1, r2, r, log, i, trupm, 
+        smm, trlowm, F, a2, bas, e, l, emax, tmp;
+  GM := GroupWithMemory(G);
+  one := One(G.1[1,1]);
+  zero := Zero(one);
+  
+  # find element x of order q-1
+  # (a power of x will become the nice generator t later)
+  qm1fac := Factors(q-1);
+  c := Product(Set(qm1fac));
+  m := (q-1)/c;
+  gens := GeneratorsOfGroup(GM);
+  xm := One(GM);
+  repeat 
+    #xm := PseudoRandom(GM);
+    xm := xm*Random(gens);
+    x := StripMemory(xm);
+    pol := [one, -x[1,1]-x[2,2], one];
+    v := [zero, one];
+    v := PowerModCoeffs(v, m, pol);
+  until PowerModCoeffs(v, c, pol) = [one] and
+            ForAll(qm1fac, p-> PowerModCoeffs(v, c/p, pol) <> [one]);
+  # eigenvalues and eigenvectors of x (zeroes of pol)
+  # we use Cantor-Zassenhaus
+  z := Z(q);
+  if q mod 2 = 0 then
+    if q mod 3 = 1 then
+      exp := (q-1)/3;
+    else
+      exp := (q^2-1)/3;
+      z := Z(q^2);
+    fi;
+  else
+    exp := (q-1)/2;
+  fi;
+  repeat
+    v := [z^Random(0,q-1), one];
+    v := PowerModCoeffs(v, exp, pol);
+  until Length(v) = 2 and ValuePol(pol, (-v[1]+one)/v[2]) = zero;
+  a := (-v[1]+one)/v[2];
+  # colums of mat are eigenvectors for a and 1/a (x^mat is diagonal)
+  mat := [[-x[1,2], -x[1,2]], [x[1,1]-a, x[1,1]-1/a]];
+  if mat[1,1] = zero and mat[2,1] = zero then
+    mat[1,1] := x[2,2]-a; mat[2,1] := -x[2,1];
+  fi;
+  if mat[1,2] = zero and mat[2,2] = zero then
+    mat[1,2] := x[2,2]-1/a; mat[2,2] := -x[2,1];
+  fi;
+
+
+  # find conjugate of x with different eigenspaces
+  # (almost all conjugates will do)
+  tm := One(GM);
+  repeat
+    tm := tm * Random(gens);
+    ym := tm*xm*tm^-1;
+    y := StripMemory(ym);
+    ymat := y*mat;
+  until ymat[1,1]*mat[2,1]-ymat[2,1]*mat[1,1] <> zero and
+        ymat[1,2]*mat[2,2]-ymat[2,2]*mat[1,2] <> zero;
+  # now y^(tm * mat) = diag(a, a^-1) 
+  tr := tm*mat;
+
+  # a-eigenvector of x in new basis
+  d := tr^-1 * [mat[1,1],mat[2,1]];
+  # can be scaled to [1,d]
+  d := d[2]/d[1];
+  cm := One(GM);
+  repeat
+    # look for cm with non-trivial conditions (i <> 0, (q-1)/2)
+    repeat 
+      cm := cm*Random(gens);
+      c := StripMemory(cm)^tr;
+      r1 := c[2,1]+d*c[2,2];
+      r2 := d^2*c[1,2]+d*c[1,1];
+    until r2 <> zero and r1 <> zero and r1 <> r2 and r1 <> -r2;;
+    r := r1 / r2;
+    log := DLog(a, r, qm1fac);
+    i := false;
+    if log mod 2 = 0 then
+      i := log/2;
+    elif q mod 2 = 0 then
+      i := (q-1-log)/2;
+    fi;
+    if IsInt(i) then
+      # this will in most cases be a transvection normalized by x
+      trupm := Comm(xm, ym^i*cm);
+      smm := trupm^mat;
+      if smm[1,2] = zero then
+        i := false;
+      else
+        # rescale first column of mat such that trupm^mat = [[1,1],[0,1]]
+        mat[1,1] := mat[1,1]*smm[1,2];
+        mat[2,1] := mat[2,1]*smm[1,2];
+        tr := tm*mat;
+      fi;
+    fi;
+  until IsInt(i);
+
+  # same for the other eigenvector of x:
+  # 1/a-eigenvector of x in new basis
+  d := tr^-1 * [mat[1,2],mat[2,2]];
+  # can be scaled to [1,d]
+  d := d[2]/d[1];
+  cm := One(GM);
+  repeat
+    # look for cm with non-trivial conditions (i <> 0, (q-1)/2)
+    repeat 
+      cm := cm*Random(gens);
+      c := StripMemory(cm)^tr;
+      r1 := c[2,1]+d*c[2,2];
+      r2 := d^2*c[1,2]+d*c[1,1];
+    until r2 <> zero and r1 <> zero and r1 <> r2 and r1 <> -r2;;
+    r := r1 / r2;
+    log := DLog(a, r, qm1fac);
+    i := false;
+    if log mod 2 = 0 then
+      i := log/2;
+    elif q mod 2 = 0 then
+      i := (q-1-log)/2;
+    fi;
+    if IsInt(i) then
+      # in most cases a transvection which becomes conjugated by mat
+      # lower triangular (here it is more difficult to rescale such 
+      # that the conjugate matrix is [[1,0],[1,1]]).
+      trlowm := Comm(xm, ym^i*cm);
+      smm := trlowm^mat;
+      if smm[2,1] = zero then
+        i := false;
+      fi;
+    fi;
+  until IsInt(i);
+
+  # adjust lower left entry of trlowm^mat to one
+  # (we use F_p linear algebra in F_q to find the nice element
+  # of products of trlowm^(x^i) for some small i)
+  if smm[2,1] <> one then
+    F := GF(q);
+    a2 := a^2;
+    bas := [smm[2,1]];
+    e := DegreeOverPrimeField(F);
+    for i in [1..e-1] do
+      Add(bas, bas[i]*a2);
+    od;
+    bas := Basis(F, bas);
+    l := List(Coefficients(bas, one), IntFFE);
+    emax := e;
+    while l[emax] = 0 do
+      emax := emax-1;
+    od;
+    tmp := trlowm;
+    if l[1] = 0 then
+      trlowm := One(GM);
+    else
+      trlowm := tmp^l[1];
+    fi;
+    for i in [2..emax] do
+      tmp := tmp^xm;
+      if l[i] <> 0 then
+        trlowm := trlowm*tmp^l[i];
+      fi;
+    od;
+  fi;
+
+  # finally power x to change a to 1/Z(q)
+  if a <> 1/Z(q) then
+    log := DLog(a, 1/Z(q), qm1fac);
+    xm := xm^log;
+  fi;
+
+  # return SLPs of elements mapped by mat to 
+  # diag(1/Z(q),Z(q))[[1,0],[1,1]], [[1,1],[0,1]],
+  # and mat
+  return [List([xm, trlowm, trupm], SLPOfElm), mat];
+end;
+
+## The following code is provided by Till Eisenbrand.
+## It integrates the algorithm above to produce the same output format
+## as RECOG.RecogniseSL2NaturalEvenChar and
+## RECOG.RecogniseSL2NaturalOddCharUsingBSGS.
+## G must be SL(2,q) generated by 2x2 matrices over GF(q).
+## Note: does not work for q = 2, 3, 5.
+ConRecogNaturalSL2 := function(G, f)
+local q, res, nicegens, diag, u1, u2, umat, lmat, k, j, i, el, result, bas, true_diag, true_u1, true_u2, basi, p, basis, coeffs, m, c, l;
+q := Size(f);
+p := Characteristic(f);
+j := DegreeOverPrimeField(GF(q));
+## standard generators of SL(2,q) in the natural representation
+true_diag  := [[Z(q)^(-1),0*Z(q)],[0*Z(q),Z(q)]];
+true_u1 := [[Z(q)^0,0*Z(q)],[Z(q)^0,Z(q)^0]];
+true_u2 := [[Z(q)^0,Z(q)^0],[0*Z(q),Z(q)^0]];
+## retry until RecogNaturalSL2 returns generators matching the standard form
+for i in [1..100] do
+res := RecogNaturalSL2(G,q);
+nicegens :=List(res[1],a->ResultOfStraightLineProgram(a,GeneratorsOfGroup(G)));
+diag := nicegens[1];
+u1 := nicegens[2];
+u2 := nicegens[3];
+## check if we found the correct generators
+if diag^res[2] = true_diag and u1^res[2] = true_u1 and u2^res[2] = true_u2 then 
+break; 
+fi;
+od;
+Print("i = " , i ,"\n");
+if IsEvenInt(q) then
+## even characteristic: conjugation by diag generates all of GF(q)* directly
+lmat := [];
+for k in [0..j-1] do
+    i := (q-1-k)*Int(2)^-1 mod (q-1);
+    el := u1^(diag^i);
+Add(lmat,el);
+od;
+umat := [];
+for k in [0..j-1] do
+    i := k * Int(2)^-1 mod (q-1);
+    el := u2^(diag^i);
+Add(umat, el);
+od;
+else
+## odd characteristic: conjugation only yields squares, so express z^l
+## in the Fp-basis {z^0, z^2, ..., z^(2(j-1))} and multiply conjugates
+basis := List([0..j-1], i -> Z(q)^(2*i));
+lmat := [];
+for l in [0..j-1] do
+  coeffs := Coefficients(Basis(GF(q), basis), Z(q)^l);
+  m := u1^0;
+for i in [0..j-1] do
+    c := IntFFE(coeffs[i+1]);
+    m := m * (diag^i * u1 * diag^(-i))^c;
+od;
+Add(lmat, m);
+od;
+umat := [];
+for l in [0..j-1] do
+  coeffs := Coefficients(Basis(GF(q), basis), Z(q)^l);
+  m := u2^0;
+for i in [0..j-1] do
+    c := IntFFE(coeffs[i+1]);
+    m := m * (diag^(-i) * u2 * diag^i)^c;
+od;
+Add(umat, m);
+od;
+fi;
+basi := res[2];
+bas := basi^(-1);
+Print("umats: \n");
+for i in [1..Length(umat)] do
+Display(umat[i]^basi);
+od;
+Print("lmats: \n");
+for i in [1..Length(lmat)] do
+Display(lmat[i]^basi);
+od;
+result :=  rec( g := G, t := lmat, s := umat, bas := bas, basi := basi,
+              one := One(f), a := umat[1]*lmat[1]*umat[1], b := One(umat[1]),
+One := One(umat[1]), f := f, q := q, p := Characteristic(f), ext := j,
+              d := 2 );
+  result.all := Concatenation(result.s,result.t,[result.a],[result.b]);
+return result;
+end;
+
+
+## test function: compares ConRecogNaturalSL2 against the built-in recognition
+## note: does not work for q = 2, 3, 5
+## Input: either a list of prime powers or the empty list 
+## (then a preset list of prime powers is tested).
+test_ConRecogNaturalSL2 := function(input)
+local i, G, list, qlist, res_old, res, f, q, valid;
+if input = [] then
+    qlist := [2^3, 2^5, 3^4, 25, 17^3, 9967]; 
+elif Size(input)>0 then
+    qlist := [];
+for q in input do
+if IsPrimePowerInt(q) then
+Add(qlist, q);
+fi;
+od;
+fi;
+valid := true;
+for q in qlist do
+Print("testing q = ", q, "\n");
+    f := GF(q);
+    list := [];
+for i in [1..5] do
+Add(list, Random(SL(2,q)));
+od;
+    G := GroupWithGenerators(list);
+if IsEvenInt(q) then
+        res_old := RECOG.RecogniseSL2NaturalEvenChar(G,f,false);
+else
+        res_old := RECOG.RecogniseSL2NaturalOddCharUsingBSGS(G,f);
+fi;
+    res := ConRecogNaturalSL2(G,f);
+## compare all generators after change of basis
+for i in [1..Length(res.all)] do
+if res.all[i]^res.basi <> res_old.all[i]^res_old.basi then
+Print("Test failed for q = ", q, ", index i = ", i, "\n");
+            valid := false;
+fi;
+od;
+od;
+if valid then
+Print("All tests passed.\n");
+fi;
+end;
+
+