Quelle  group5.38a   Sprache: unbekannt

readi p,"Input the prime p";
w:=0;
F:=FiniteField(p);
for i in [2..p-1] do
a:=F!i;
S:={a};
for j in [2..p-1] do
  Include(~S,a^j);
end for;
if #S eq p-1 then
  w:=i;
  break;
end if;
end for;
print "p equals",p;
print "w equals",w;



//Take advantage of the fact that you can multiply any matrix by k^3
//for any non-zero k.  This means you can take the leading entry
//in a non-zero matrix to lie in a transversal for the cubes in Z_p\{0}.

tt:=Cputime();

range:=[0,1];

if (p mod 3) eq 1 then
  //Get a transversal for the cubes
  cubes:={i^3 mod p: i in [1..p-1]};
  for i in [2..p-1] do
    if i in cubes then continue; end if;
    u:=i;
    break;
  end for;
  cubesxu:={(x*u) mod p:x in cubes};
  for i in [2..p-1] do
    if i in cubes join cubesxu then continue; end if;
    v:=i;
    break;
  end for;
  range:=[0,1,u,v];
end if;


Z:=Integers();
V2:=VectorSpace(F,2);
H22:=Hom(V2,V2);

//Case babb=baaa
mats1:={[0,0,0,0]};

for y1 in range do
for y2 in [0..p-1] do
if y1 eq 0 and y2 notin range then continue; end if;
for y3 in [0..p-1] do
for y4 in [0..p-1] do

A:=H22![y1,y2,y3,y4];
if A eq 0 then continue; end if;

new:=true;
index:=p^3*y1+p^2*y2+p*y3+y4;

for a in [0..p-1] do
for b in [0..p-1] do

if a eq b or a eq p-b then continue; end if;

B:=H22![a,b,b,a];
C:=H22![a^4-b^4,2*a*b*(a^2-b^2),2*a*b*(a^2-b^2),a^4-b^4];
D:=B*A*C^-1;

z1:=Z!(D[1][1]);
z2:=Z!(D[1][2]);
z3:=Z!(D[2][1]);
z4:=Z!(D[2][2]);

ind1:=p^3*z1+p^2*z2+p*z3+z4;

if ind1 lt index then new:=false; break; end if;

B[2]:=-B[2];
C[1]:=-C[1];
D:=B*A*C^-1;

z1:=Z!(D[1][1]);
z2:=Z!(D[1][2]);
z3:=Z!(D[2][1]);
z4:=Z!(D[2][2]);

ind1:=p^3*z1+p^2*z2+p*z3+z4;

if ind1 lt index then new:=false; break; end if;


end for;
if not new then break; end if;
end for;

if new then
  Include(~mats1,[y1,y2,y3,y4]);
end if;

//These checks whether we have reached the expected size really help!
if #mats1 eq (Gcd(p-1,3)*(p^2+3*p+11)+1)/2 then break; end if;
end for;
if #mats1 eq (Gcd(p-1,3)*(p^2+3*p+11)+1)/2 then break; end if;
end for;
if #mats1 eq (Gcd(p-1,3)*(p^2+3*p+11)+1)/2 then break; end if;
end for;
if #mats1 eq (Gcd(p-1,3)*(p^2+3*p+11)+1)/2 then break; end if;
end for;

//Case babb=wbaaa
mats2:={[0,0,0,0]};

for y1 in range do
for y2 in [0..p-1] do
if y1 eq 0 and y2 notin range then continue; end if;
for y3 in [0..p-1] do
for y4 in [0..p-1] do

A:=H22![y1,y2,y3,y4];
if A eq 0 then continue; end if;

new:=true;
index:=p^3*y1+p^2*y2+p*y3+y4;

for a in [0..p-1] do
for b in [0..p-1] do

if a+b eq 0 then continue; end if;

B:=H22![a,b,w*b,a];
C:=H22![a^4-w^2*b^4,2*a*b*(a^2-w*b^2),2*w*a*b*(a^2-w*b^2),a^4-w^2*b^4];
D:=B*A*C^-1;

z1:=Z!(D[1][1]);
z2:=Z!(D[1][2]);
z3:=Z!(D[2][1]);
z4:=Z!(D[2][2]);

ind1:=p^3*z1+p^2*z2+p*z3+z4;

if ind1 lt index then new:=false; break; end if;


B[2]:=-B[2];
C[1]:=-C[1];
D:=B*A*C^-1;

z1:=Z!(D[1][1]);
z2:=Z!(D[1][2]);
z3:=Z!(D[2][1]);
z4:=Z!(D[2][2]);

ind1:=p^3*z1+p^2*z2+p*z3+z4;

if ind1 lt index then new:=false; break; end if;


end for;
if not new then break; end if;
end for;

if new then
  Include(~mats2,[y1,y2,y3,y4]);
end if;

if #mats2 eq (Gcd(p-1,3)*(p^2+p+1)+5)/2 then break; end if;
end for;
if #mats2 eq (Gcd(p-1,3)*(p^2+p+1)+5)/2 then break; end if;
end for;
if #mats2 eq (Gcd(p-1,3)*(p^2+p+1)+5)/2 then break; end if;
end for;
if #mats2 eq (Gcd(p-1,3)*(p^2+p+1)+5)/2 then break; end if;
end for;


print #mats1,(Gcd(p-1,3)*(p^2+3*p+11)+1)/2;
print #mats2,(Gcd(p-1,3)*(p^2+p+1)+5)/2;

print Cputime(tt);