1readi p,"Input prime";
2
3w:=0;
4F:=FiniteField(p);
5for i in [2..p-1] do
6a:=F!i;
7S:={a};
8for j in [2..p-1] do
9  Include(~S,a^j);
10end for;
11if #S eq p-1 then
12  w:=i;
13  break;
14end if;
15end for;
16
17SQ:={};
18for x in [0..p-1] do
19  y:=x^2 mod p;
20  Include(~SQ,y);
21end for;
22
23//Get least non-square
24ns:=0;
25for x in [2..p-1] do
26  if x in SQ then continue; end if;
27  ns:=x;
28  break;
29end for;
30
31w2:=w^2 mod p;
32
33Z:=Integers();
34V2:=VectorSpace(F,2);
35H22:=Hom(V2,V2);
36
37//ca=bab, cb=baa
38mats1:=[];
39range1:=[0,1,ns];
40range2:=[0,1];
41y1range:=range1;
42if p mod 4 ne 1 then y1range:=range2; end if;
43
44for y1 in y1range do
45y2range:=[0..p-1];
46if y1 eq 0 then y2range:=range1; end if;
47for y2 in y2range do
48for y3 in [0..p-1] do
49for y4 in [0..p-1] do
50
51A:=H22![y1,y2,y3,y4];
52if A eq 0 then Append(~mats1,[0,0,0,0]); continue; end if;
53
54new:=1;
55index:=p^3*y1+p^2*y2+p*y3+y4;
56
57for a in [0..p-1] do
58for b in [0..p-1] do
59
60e:=F!(a^2-b^2);
61if e eq 0 then continue; end if;
62
63P:=H22![a,b,b,a];
64
65D:=e^-1*P*A*P^-1;
66
67z1:=Z!(D[1,1]);
68z2:=Z!(D[1,2]);
69z3:=Z!(D[2,1]);
70z4:=Z!(D[2,2]);
71
72ind1:=p^3*z1+p^2*z2+p*z3+z4;
73
74if ind1 lt index then
75  new:=0;
76end if;
77
78P:=H22![a,b,-b,-a];
79
80D:=-e^-1*P*A*P^-1;
81
82z1:=Z!(D[1,1]);
83z2:=Z!(D[1,2]);
84z3:=Z!(D[2,1]);
85z4:=Z!(D[2,2]);
86
87ind1:=p^3*z1+p^2*z2+p*z3+z4;
88
89if ind1 lt index then
90  new:=0;
91end if;
92
93if new eq 0 then break; end if;
94end for;
95if new eq 0 then break; end if;
96end for;
97
98if new eq 1 then
99  Append(~mats1,[y1,y2,y3,y4]);
100end if;
101
102end for;
103end for;
104end for;
105end for;
106
107print #mats1,p^2+(7*p+15)/2;
108
109//ca=wbab, cb=baa
110mats2:=[];
111
112expect:=p^2+(3*p+5)/2;
113
114for y1 in [0,1] do
115y2range:=[0..p-1];
116if y1 eq 0 then y2range:=[0,1,ns]; end if;
117for y2 in y2range do
118for y3 in [0..p-1] do
119for y4 in [0..p-1] do
120
121A:=H22![y1,y2,y3,y4];
122if A eq 0 then Append(~mats2,[0,0,0,0]); continue; end if;
123
124new:=1;
125index:=p^3*y1+p^2*y2+p*y3+y4;
126
127for a in [0..p-1] do
128for b in [0..p-1] do
129
130e:=F!(a^2-w*b^2);
131if e eq 0 then continue; end if;
132
133P:=H22![a,w*b,b,a];
134
135D:=e^-1*P*A*P^-1;
136
137z1:=Z!(D[1,1]);
138z2:=Z!(D[1,2]);
139z3:=Z!(D[2,1]);
140z4:=Z!(D[2,2]);
141
142ind1:=p^3*z1+p^2*z2+p*z3+z4;
143
144if ind1 lt index then
145  new:=0;
146end if;
147
148P:=H22![a,w*b,-b,-a];
149
150D:=-e^-1*P*A*P^-1;
151
152z1:=Z!(D[1,1]);
153z2:=Z!(D[1,2]);
154z3:=Z!(D[2,1]);
155z4:=Z!(D[2,2]);
156
157ind1:=p^3*z1+p^2*z2+p*z3+z4;
158
159if ind1 lt index then
160  new:=0;
161end if;
162
163if new eq 0 then break; end if;
164end for;
165if new eq 0 then break; end if;
166end for;
167
168if new eq 1 then
169  Append(~mats2,[y1,y2,y3,y4]);
170end if;
171
172end for;
173end for;
174end for;
175end for;
176
177print #mats2,p^2+(3*p+5)/2;
178
179