`PreservedSesquilinearForms` fails for forms preserved up to scalars

[chseeger]

This might be something that has been known for some time, but I could not find an issue documenting it.

As stated in the documentation of PreservedSesquilinearForms (see Section 4.8-2):

The function uses random methods to find all of the bilinear or unitary forms preserved by group [...] up to a scalar. Since the procedure relies on a pseudo-random generator, the user may need to execute the operation more than once to find all invariant sesquilinear forms.

In practice, however, I observe that forms preserved up to a scalar are sometimes not detected at all, even after repeated runs. This seems not just due to randomness, but caused by the scalar normalization step in ClassicalForms_GeneratorsWithBetterScalarsSesquilinear.

The following lines look problematic:

forms/lib/recognition_new.gi

Lines 126 to 134 in be8afc2

	root := NthRoot(field,a1[2],a1[1]);
	if a1[1] = 1 then # the matrix m1 has scalar a1[2]
	#if a1[2] has a square root, we can replace m1 with m1*sqrt{a1[2]};
	if LogFFE(a1[2],PrimitiveRoot(field)) mod 2 = 0 then
	return [m1/NthRoot(field,a1[2],2),[One(field)]];
	fi;
	return [m1,[a1[2]]]; #originally, the three lines above this return were not there. Those three lines make sure scalar becomes 1 if possible (basicaly if there is a sqrt).
	elif LogFFE(root,PrimitiveRoot(field)) mod (q+1) = 0 then
	return [m1/NthRoot(field,root,q+1),[One(field)]]; #either frob = id, then q+1 = 2, or frob is not trivial, then we take q+1-st root. In both cases, modify m1 to a matrix that has scalar one.

Two observations:

• In the case a1[1] = 1, normalization uses exponent 2. This is only correct in the bilinear case. I think a minimal fix would be to replace 2 by q+1:

      if a1[1] = 1 then
          if LogFFE(a1[2],PrimitiveRoot(field)) mod (q+1) = 0 then
              return [m1/NthRoot(field,a1[2],q+1),[One(field)]];
          fi;
          return [m1,[a1[2]]];

  Here is a reproducer (the unitary form is never found):

  Group([ [ [ Z(3^2)^3, Z(3^2)^6, Z(3^2)^7 ], [ Z(3^2)^7, Z(3^2)^2, Z(3)^0 ], [ Z(3^2)^7, Z(3^2)^7, 0*Z(3) ] ], [ [ 0*Z(3), 0*Z(3), Z(3) ], [ Z(3^2)^7, Z(3)^0, Z(3^2)^7 ], [ Z(3), Z(3^2), Z(3) ] ], [ [ Z(3^2), 0*Z(3), 0*Z(3) ], [ 0*Z(3), Z(3^2), 0*Z(3) ], [ 0*Z(3), 0*Z(3), Z(3^2) ] ] ]);

• It seems that the elif condition is too weak (but note that I haven't read [CLN+08]): It checks whether one (arbitrarily chosen) solution of $\lambda^{a1[1]}=a1[2]$ is a $(q+1)$-st power.

  The naive solution would be to simply remove the elif branch; in fact, that significantly improves the detection rate (at least in my use cases).

  Here is a reproducer (the bilinear form is found only occasionally):

  Group(Z(5)^0*[
  [ [ 0, 1, 0, 0, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 1, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 1 ],
    [ 2, 3, 4, 3, 1, 3, 2, 4 ],
    [ 1, 4, 0, 0, 1, 0, 0, 0 ],
    [ 0, 2, 2, 2, 4, 4, 2, 1 ],
    [ 2, 3, 4, 2, 1, 3, 3, 4 ],
    [ 2, 2, 3, 0, 1, 0, 0, 3 ] ],
  [ [ 0, 0, 1, 0, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 1, 0, 0 ],
    [ 2, 4, 0, 3, 2, 1, 2, 2 ],
    [ 2, 4, 0, 2, 4, 4, 4, 4 ],
    [ 4, 4, 1, 0, 1, 0, 3, 2 ],
    [ 1, 1, 2, 3, 2, 4, 1, 3 ],
    [ 3, 2, 0, 3, 1, 4, 1, 4 ],
    [ 3, 3, 1, 4, 0, 3, 4, 3 ] ],
  [ [ 0, 0, 0, 1, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 1, 0 ],
    [ 2, 2, 3, 1, 3, 3, 2, 4 ],
    [ 3, 0, 0, 0, 0, 0, 0, 0 ],
    [ 1, 0, 2, 1, 0, 3, 2, 1 ],
    [ 1, 3, 2, 4, 3, 1, 2, 2 ],
    [ 0, 3, 0, 0, 0, 0, 0, 0 ],
    [ 0, 1, 3, 0, 3, 1, 0, 1 ] ]
  ]);
  

If useful, I would be happy to open a small PR with the minimal changes described above – however, I'm aware of @GalaxyCo's ongoing work in PR #79, which may eventually replace the existing routine.

[fingolfin]

I was hoping that perhaps @jdebeule might chime in, or @GalaxyCo but so far no such luck.

In any case, please open a PR with your proposed changes. Then I can email all involved parties to figure out how we can move on!

[GalaxyCo]

Hello, sorry i was looking into this. My function to compute the non degenerate forms PreservedFormsWithScalars() is also using this routine, however i think i found a fix. The issue seems to lie in the improvegenerator() function. I have changed it to the following (in my branch):

improvegenerator := function(m1)
        local a1, scalars;
        a1 := ClassicalForms_PossibleScalarsSesquilinear(field,m1,frob);
        #Recall that the scalars satisfy the equation lambda^a[1] = a[2]
        if a1 = false then
            return false; #The group does not preserve a form modulo scalars
        fi;
                
        if IsList(a1) then
            if a1[1] = 1 then # the matrix m1 has scalar a1[2]
                return [a1[2]];
            else
                scalars := AsList(Group(NthRoot(field,a1[2],a1[1]))); # add all possible scalars for m1
                return scalars;
            fi;
        fi;
    end;

As you can see, I have taken out much of the logic that tries to make the scalar one, i am not sure why this is needed other than for aesthetic reasons. Furthermore, I have also taken out the code that tries to randomy multiply the group generator and calls the function recursively, this might result in a function that returns more possible scalars, which could result in a need to compute the preserved forms multiple times with respect to the different scalars, however it seems that in practice there is not mucht performance overhead added from removing this code. But it seems to have fixed the bug (at least for the two examples by @chseeger) Anyways, the performance of the functions should be tested carefully with this change. It might very well be the case that there is an example of a group where computing the preserved forms with this code in place becomes to slow.

Furthermore, longer term i am also planning on implementing a function to compute all scalars, even of degenerate forms because this function is not capable of that. However, i have not started with this yet.