1//Immediate descendants of order p^7 of algebra 5.14
2//Modified to cut down complexity 20 August 2013
3
4readi p,"Input the prime p";
5//Get a primitive element
6w:=0;
7F:=FiniteField(p);
8for i in [2..p-1] do
9a:=F!i;
10S:={a};
11for j in [2..p-1] do
12  Include(~S,a^j);
13end for;
14if #S eq p-1 then
15  w:=i;
16  break;
17end if;
18end for;
19print "p equals",p;
20print "w equals",w;
21
22//Get the squares, and the least non-square lns
23SQ:={};
24for x in [1..((p-1) div 2)] do
25  y:=x^2 mod p;
26  Include(~SQ,y);
27end for;
28for i in [2..p-1] do
29  if i notin SQ then lns:=i; break; end if;
30end for;
31
32//Get cube roots of 1
33CU:={1};
34for x in [2..p-1] do
35  if x^3 mod p eq 1 then Include(~CU,x); end if;
36end for;
37
38//If p=1 mod 3 compute a transversal for the cube roots of 1
39transversal3:=[];
40if p mod 3 eq 1 then
41  transversal3:=[];
42  sofar:={};
43  for i in [1..p-1] do
44    if i^3 mod p in sofar then continue; end if;
45    Append(~transversal3,i);
46    Include(~sofar,i^3 mod p);
47  end for;
48end if;
49
50//If p=1 mod 4 then compute a transversal for the fourth roots of 1
51transversal4:=[];
52if p mod 4 eq 1 then
53  sofar:={};
54  for i in [1..(p-1) div 2] do
55    if i^4 mod p in sofar then continue; end if;
56    Append(~transversal4,i);
57    Include(~sofar,i^4 mod p);
58  end for;
59end if;
60
61//If p=1 mod 5 then compute a transversal for the fifth roots of 1
62transversal5:=[];
63if p mod 5 eq 1 then
64  sofar:={};
65  for i in [1..p-1] do
66    if i^5 mod p in sofar then continue; end if;
67    Append(~transversal5,i);
68    Include(~sofar,i^5 mod p);
69  end for;
70end if;
71
72//If p=1 mod 3 then compute a transversal for the sixth roots of 1
73transversal6:=[];
74if p mod 3 eq 1 then
75  sofar:={};
76  for i in [1..p-1] do
77    if i^6 mod p in sofar then continue; end if;
78    Append(~transversal6,i);
79    Include(~sofar,i^6 mod p);
80  end for;
81end if;
82
83Z:=Integers();
84V1:=VectorSpace(F,1);
85V2:=VectorSpace(F,2);
86V3:=VectorSpace(F,3);
87H12:=Hom(V1,V2);
88H22:=Hom(V2,V2);
89H32:=Hom(V3,V2);
90H33:=Hom(V3,V3);
91
92gtotal:=0;
93GR:=[];
94//These store the parameters in each of the 24 cases
95//Currently they are stored as integer sequences of length 6,
96//but an alternative (which takes less space) is to store them
97//as vectors or matrices over GF(p).  Two have size p^5+.
98params1:=[];
99params2:=[];
100params3:=[];
101params4:=[];
102params5:=[];
103params6:=[];
104params7:=[];
105params8:=[];
106params9:=[];
107params10:=[];
108params11:=[];
109params12:=[];
110params13:=[];
111params14:=[];
112params15:=[];
113params16:=[];
114params17:=[];
115params18:=[];
116params19:=[];
117params20:=[];
118params21:=[];
119params22:=[];
120params23:=[];
121params24:=[];
122
123tt:=Cputime();
124
125/*
126Descendants of 5.14.
127
128Case 1: cb=caa=cab=cac=0
129*/
130
131Append(~params1,[0,0,0,0,0,0]);
132Append(~params1,[1,0,0,0,0,0]);
133Append(~params1,[0,1,0,0,0,0]);
134Append(~params1,[0,w,0,0,0,0]);
135
136Append(~params1,[0,0,1,0,0,0]);
137Append(~params1,[0,0,w,0,0,0]);
138Append(~params1,[0,1,1,0,0,0]);
139Append(~params1,[0,1,w,0,0,0]);
140Append(~params1,[0,w,1,0,0,0]);
141Append(~params1,[0,w,w,0,0,0]);
142for y in [0..p-1] do
143  Append(~params1,[1,y,1,0,0,0]);
144  Append(~params1,[1,y,w,0,0,0]);
145end for;
146
147Append(~params1,[1,1,0,1,0,0]);
148Append(~params1,[1,w,0,1,0,0]);
149for x in [0..p-1] do
150  Append(~params1,[x,0,0,1,0,0]);
151end for;
152
153Append(~params1,[0,0,0,0,1,0]);
154Append(~params1,[0,1,0,0,1,0]);
155Append(~params1,[0,w,0,0,1,0]);
156Append(~params1,[0,0,0,1,1,0]);
157Append(~params1,[0,1,0,1,1,0]);
158Append(~params1,[0,w,0,1,1,0]);
159
160Append(~params1,[1,0,0,0,0,1]);
161Append(~params1,[0,0,0,0,0,1]);
162Append(~params1,[0,0,1,0,0,1]);
163Append(~params1,[0,0,w,0,0,1]);
164
165print #params1,3*p+22;
166gtotal:=gtotal+#params1;
167
168
169/*
170Descendants of 5.14.
171
172Case 2: caa=cab=cac=0, cb=baa
173*/
174
175Append(~params2,[0,0,0,0,0,0]);
176Append(~params2,[1,0,0,0,0,0]);
177Append(~params2,[0,1,0,0,0,0]);
178Append(~params2,[0,w,0,0,0,0]);
179
180Append(~params2,[0,0,1,0,0,0]);
181Append(~params2,[0,0,w,0,0,0]);
182Append(~params2,[0,1,1,0,0,0]);
183Append(~params2,[0,1,w,0,0,0]);
184Append(~params2,[0,w,1,0,0,0]);
185Append(~params2,[0,w,w,0,0,0]);
186for y in [0..p-1] do
187  Append(~params2,[1,y,1,0,0,0]);
188  Append(~params2,[1,y,w,0,0,0]);
189end for;
190
191Append(~params2,[1,1,0,1,0,0]);
192Append(~params2,[1,w,0,1,0,0]);
193for x in [0..p-1] do
194  Append(~params2,[x,0,0,1,0,0]);
195end for;
196
197for x in [0..p-1] do
198  Append(~params2,[0,x,0,0,1,0]);
199  Append(~params2,[0,x,0,1,1,0]);
200end for;
201
202Append(~params2,[0,0,0,0,0,1]);
203Append(~params2,[1,0,0,0,0,1]);
204if p mod 3 eq 1 then
205  Append(~params2,[w,0,0,0,0,1]);
206  Append(~params2,[w^2,0,0,0,0,1]);
207end if;
208Append(~params2,[0,0,1,0,0,1]);
209Append(~params2,[0,0,w,0,0,1]);
210if p mod 4 eq 1 then
211  Append(~params2,[0,0,w^2,0,0,1]);
212  Append(~params2,[0,0,w^3,0,0,1]);
213end if;
214
215print #params2,5*p+13+Gcd(p-1,3)+Gcd(p-1,4);
216gtotal:=gtotal+#params2;
217
218/*
219Descendants of 5.14.
220
221Case 3: cb=bab=bac=cac=0
222*/
223
224Append(~params3,[0,0,0,0,0,0]);
225Append(~params3,[1,0,0,0,0,0]);
226Append(~params3,[0,0,0,0,1,0]);
227Append(~params3,[0,1,0,0,1,0]);
228Append(~params3,[0,0,0,0,1,1]);
229Append(~params3,[1,0,0,0,1,1]);
230Append(~params3,[0,0,0,0,1,w]);
231Append(~params3,[1,0,0,0,1,w]);
232Append(~params3,[0,0,1,0,0,1]);
233Append(~params3,[0,0,w,0,0,w]);
234for i in [1..p-1] do
235  Append(~params3,[0,0,0,i,1,1]);
236  Append(~params3,[0,0,0,i,1,w]);
237end for;
238Append(~params3,[0,0,0,1,1,0]);
239Append(~params3,[0,0,0,w,1,0]);
240if p mod 4 eq 1 then
241  Append(~params3,[0,0,0,w^2 mod p,1,0]);
242  Append(~params3,[0,0,0,w^3 mod p,1,0]);
243end if;
244
245print #params3,2*p+8+Gcd(p-1,4);
246gtotal:=gtotal+#params3;
247
248/*
249Descendants of 5.14.
250
251Case 4: bab=bac=cac=0, cb=baa
252*/
253
254Append(~params4,[0,0,0,0,0,0]);
255Append(~params4,[1,0,0,0,0,0]);
256Append(~params4,[0,1,0,0,0,0]);
257if p mod 3 eq 1 then
258  Append(~params4,[0,w,0,0,0,0]);
259  Append(~params4,[0,w^2,0,0,0,0]);
260end if;
261
262Append(~params4,[0,0,0,0,1,0]);
263Append(~params4,[0,1,0,0,1,0]);
264if p mod 3 eq 1 then
265  Append(~params4,[0,w,0,0,1,0]);
266  Append(~params4,[0,w^2,0,0,1,0]);
267end if;
268
269Append(~params4,[0,0,0,0,0,1]);
270Append(~params4,[1,0,0,0,0,1]);
271Append(~params4,[0,0,0,0,0,w]);
272Append(~params4,[1,0,0,0,0,w]);
273
274for y in [0..p-1] do
275  Append(~params4,[0,0,1,0,0,y]);
276  Append(~params4,[0,0,w,0,0,y]);
277end for;
278Append(~params4,[0,0,1,0,1,1]);
279Append(~params4,[0,0,w,0,1,w]);
280for z in [1..(p-1) div 2] do
281  Append(~params4,[0,z,1,0,0,0]);
282  Append(~params4,[0,z,w,0,0,0]);
283end for;
284
285Append(~params4,[0,0,0,1,0,0]);
286Append(~params4,[1,0,0,1,0,0]);
287if p mod 5 eq 1 then
288  Append(~params4,[w,0,0,1,0,0]);
289  Append(~params4,[w^2,0,0,1,0,0]);
290  Append(~params4,[w^3,0,0,1,0,0]);
291  Append(~params4,[w^4,0,0,1,0,0]);
292end if;
293
294Append(~params4,[0,0,0,1,1,0]);
295Append(~params4,[0,0,0,1,w,0]);
296if p mod 4 eq 1 then
297  Append(~params4,[0,0,0,1,w^2,0]);
298  Append(~params4,[0,0,0,1,w^3,0]);
299end if;
300for x in [0..p-1] do
301  Append(~params4,[0,0,0,1,x,1]);
302  Append(~params4,[0,0,0,1,x,w]);
303end for;
304for z in [1..(p-1) div 2] do
305  Append(~params4,[z,0,0,1,0,1]);
306  Append(~params4,[z,0,0,1,0,w]);
307end for;
308
309print #params4,6*p+8+2*Gcd(p-1,3)+Gcd(p-1,4)+Gcd(p-1,5);
310gtotal:=gtotal+#params4;
311
312/*
313Descendants of 5.14.
314
315Case 5: cb=bac=cac=0, caa=bab
316*/
317
318Append(~params5,[0,0,0,0,0,0]);
319Append(~params5,[1,0,0,0,0,0]);
320Append(~params5,[0,1,0,0,0,0]);
321Append(~params5,[0,w,0,0,0,0]);
322
323Append(~params5,[0,0,1,0,0,0]);
324Append(~params5,[0,0,w,0,0,0]);
325Append(~params5,[0,1,1,0,0,0]);
326Append(~params5,[0,1,w,0,0,0]);
327Append(~params5,[0,w,1,0,0,0]);
328Append(~params5,[0,w,w,0,0,0]);
329
330for x in [0..p-1] do
331  Append(~params5,[x,0,0,1,0,0]);
332end for;
333
334Append(~params5,[0,0,0,0,1,0]);
335Append(~params5,[0,1,0,0,1,0]);
336Append(~params5,[0,w,0,0,1,0]);
337if p mod 3 eq 1 then
338  Append(~params5,[0,w^2,0,0,1,0]);
339  Append(~params5,[0,w^3,0,0,1,0]);
340  Append(~params5,[0,w^4,0,0,1,0]);
341  Append(~params5,[0,w^5,0,0,1,0]);
342end if;
343
344if p mod 4 eq 1 then
345for x in [0..(p-1) div 2] do
346  Append(~params5,[0,x,0,1,1,0]);
347  Append(~params5,[0,x,0,w,1,0]);
348  Append(~params5,[0,x,0,w^2,1,0]);
349  Append(~params5,[0,x,0,w^3,1,0]);
350end for;
351end if;
352
353if p mod 4 eq 3 then
354for x in [0..p-1] do
355  Append(~params5,[0,x,0,1,1,0]);
356  Append(~params5,[0,x,0,w,1,0]);
357end for;
358end if;
359
360for x in [0..p-1] do
361Append(~params5,[0,0,x,0,0,1]);
362Append(~params5,[0,0,x,0,0,w]);
363end for;
364Append(~params5,[1,0,0,0,0,1]);
365Append(~params5,[1,0,0,0,0,w]);
366Append(~params5,[1,0,1,0,0,1]);
367Append(~params5,[1,0,w,0,0,w]);
368
369print #params5,5*p+13+2*Gcd(p-1,3)+Gcd(p-1,4);
370gtotal:=gtotal+#params5;
371
372/*
373Descendants of 5.14.
374
375Case 6: bac=cac=0, caa=bab, cb=baa
376*/
377
378Append(~params6,[0,0,0,0,0,0]);
379Append(~params6,[1,0,0,0,0,0]);
380if p mod 5 eq 1 then
381  Append(~params6,[w,0,0,0,0,0]);
382  Append(~params6,[w^2,0,0,0,0,0]);
383  Append(~params6,[w^3,0,0,0,0,0]);
384  Append(~params6,[w^4,0,0,0,0,0]);
385end if;
386Append(~params6,[0,1,0,0,0,0]);
387Append(~params6,[0,w,0,0,0,0]);
388if p mod 3 eq 1 then
389  Append(~params6,[0,w^2,0,0,0,0]);
390  Append(~params6,[0,w^3,0,0,0,0]);
391  Append(~params6,[0,w^4,0,0,0,0]);
392  Append(~params6,[0,w^5,0,0,0,0]);
393end if;
394if p mod 4 eq 3 then
395  for x in [0..p-1] do
396    Append(~params6,[0,x,1,0,0,0]);
397    Append(~params6,[0,x,w,0,0,0]);
398  end for;
399end if;
400if p mod 4 eq 1 then
401  for x in [0..(p-1) div 2] do
402    Append(~params6,[0,x,1,0,0,0]);
403    Append(~params6,[0,x,w,0,0,0]);
404    Append(~params6,[0,x,w^2,0,0,0]);
405    Append(~params6,[0,x,w^3,0,0,0]);
406  end for;
407end if;
408for x in [0..p-1] do
409  Append(~params6,[x,0,0,1,0,0]);
410  if p mod 5 eq 1 then
411    Append(~params6,[x,0,0,w,0,0]);
412    Append(~params6,[x,0,0,w^2,0,0]);
413    Append(~params6,[x,0,0,w^3,0,0]);
414    Append(~params6,[x,0,0,w^4,0,0]);
415  end if;
416end for;
417if p mod 3 eq 2 then
418  for x in [0..p-1] do
419  for y in [0..p-1] do
420    Append(~params6,[0,y,0,x,1,0]);
421  end for;
422  end for;
423end if;
424if p mod 3 eq 1 then
425  xrange:=transversal3;
426  Append(~xrange,0);
427  for x in xrange do
428  for y in [0..p-1] do
429    Append(~params6,[0,y,0,x,1,0]);
430    Append(~params6,[0,y,0,x,w,0]);
431    Append(~params6,[0,y,0,x,w^2,0]);
432  end for;
433  end for;
434end if;
435if p mod 4 eq 3 then
436  for y in [0..(p-1) div 2] do
437    Append(~params6,[y,0,0,0,0,1]);
438    Append(~params6,[y,0,0,0,0,w]);
439  end for;
440  for x in [1..p-1] do
441    Append(~params6,[0,0,x,0,0,1]);
442    Append(~params6,[0,0,x,0,0,w]);
443  end for;
444  for y in [1..(p-1) div 2] do
445    Append(~params6,[y,0,1,0,0,1]);
446    Append(~params6,[y,0,w,0,0,w]);
447  end for;
448end if;
449if p mod 4 eq 1 then
450  yrange:=transversal4;
451  Append(~yrange,0);
452  for y in yrange do
453    Append(~params6,[y,0,0,0,0,1]);
454    Append(~params6,[y,0,0,0,0,w]);
455    Append(~params6,[y,0,0,0,0,w^2]);
456    Append(~params6,[y,0,0,0,0,w^3]);
457  end for;
458  for x in [1..p-1] do
459    Append(~params6,[0,0,x,0,0,1]);
460    Append(~params6,[0,0,x,0,0,w]);
461    Append(~params6,[0,0,x,0,0,w^2]);
462    Append(~params6,[0,0,x,0,0,w^3]);
463  end for;
464  for y in transversal4 do
465    Append(~params6,[y,0,1,0,0,1]);
466    Append(~params6,[y,0,w,0,0,w]);
467    Append(~params6,[y,0,w^2,0,0,w^2]);
468    Append(~params6,[y,0,w^3,0,0,w^3]);
469  end for;
470end if;
471
472
473print #params6,
474p^2+3*p-3+(p+2)*Gcd(p-1,3)+(p+1)*Gcd(p-1,4)+(p+1)*Gcd(p-1,5);
475gtotal:=gtotal+#params6;
476
477
478/*
479Descendants of 5.14.
480
481Case 7: cb=baa=bac=cac=0
482*/
483
484for u in [0,1,w] do
485for t in [0,1] do
486for x in [0,1] do
487  Append(~params7,[u,0,t,x,0,0]);
488end for;
489end for;
490end for;
491
492for u in [0,1,w] do
493for x in [0,1] do
494  Append(~params7,[u,1,0,x,0,0]);
495end for;
496end for;
497
498for x in [0,1,w] do
499  Append(~params7,[0,1,1,x,0,0]);
500end for;
501
502for u in [1,w] do
503for x in [0..p-1] do
504  Append(~params7,[u,1,1,x,0,0]);
505end for;
506end for;
507
508for u in [0,1,w] do
509for x in [0,1] do
510for z in [1,w] do
511  Append(~params7,[u,0,0,x,0,z]);
512end for;
513end for;
514end for;
515
516for u in [0..p-1] do
517for x in [0,1] do
518for z in [1,w] do
519  Append(~params7,[u,0,1,x,0,z]);
520end for;
521end for;
522end for;
523
524for v in [0,1,w] do
525for z in [0,1,w] do
526  Append(~params7,[0,v,0,0,1,z]);
527end for;
528end for;
529
530vrange:=[0,1,w];
531if p mod 4 eq 1 then vrange:=[0,1,w,w^2,w^3]; end if;
532for v in vrange do
533  Append(~params7,[0,v,0,1,1,0]);
534end for;
535
536for v in [0..p-1] do
537for z in [1,w] do
538  Append(~params7,[0,v,0,1,1,z]);
539end for;
540end for;
541
542for v in [0..p-1] do
543for z in [0,1,w] do
544  Append(~params7,[0,v,1,0,1,z]);
545end for;
546end for;
547
548for v in [0..p-1] do
549for x in [1,w] do
550for z in [0..p-1] do
551  Append(~params7,[0,v,1,x,1,z]);
552end for;
553end for;
554end for;
555
556print #params7,2*p^2+11*p+43+Gcd(p-1,4);
557gtotal:=gtotal+#params7;
558
559/*
560Descendants of 5.14.
561
562Case 8: cb=caa, baa=bac=cac=0
563*/
564
565
566urange:=[0,1,w];
567if p mod 4 eq 1 then urange:=[0,1,w,w^2,w^3]; end if;
568for u in urange do
569  Append(~params8,[u,0,0,0,0,0]);
570  Append(~params8,[u,0,0,1,0,0]);
571  Append(~params8,[u,1,0,0,0,0]);
572end for;
573
574if p mod 3 eq 1 then
575  urange:=transversal3;
576  Append(~urange,0);
577else;
578  urange:=[0..p-1];
579end if;
580for u in urange do
581  Append(~params8,[u,0,1,0,0,0]);
582  Append(~params8,[u,0,1,1,0,0]);
583  Append(~params8,[u,1,1,0,0,0]);
584  if p mod 3 eq 1 then
585    Append(~params8,[u,0,w,0,0,0]);
586    Append(~params8,[u,0,w,1,0,0]);
587    Append(~params8,[u,1,w,0,0,0]);
588    Append(~params8,[u,0,w^2,0,0,0]);
589    Append(~params8,[u,0,w^2,1,0,0]);
590    Append(~params8,[u,1,w^2,0,0,0]);
591  end if;
592end for;
593
594for u in [0..p-1] do
595for t in [0..p-1] do
596  Append(~params8,[u,1,t,1,0,0]);
597end for;
598end for;
599
600for u in [0..p-1] do
601for t in [0..(p-1) div 2] do
602for x in [0,1] do
603for z in [1,w] do
604  Append(~params8,[u,0,t,x,0,z]);
605end for;
606end for;
607end for;
608end for;
609
610vrange:=[0,1,w];
611if p mod 3 eq 1 then vrange:=[0,1,w,w^2,w^3,w^4,w^5]; end if;
612for v in vrange do
613  Append(~params8,[0,v,0,0,1,0]);
614end for;
615
616if p mod 5 ne 1 then
617  for v in [0..p-1] do
618    Append(~params8,[0,v,0,1,1,0]);
619  end for;
620end if;
621
622if p mod 5 eq 1 then
623  vrange:=transversal5;
624  Append(~vrange,0);
625  for v in vrange do
626  for x in [1,w,w^2,w^3,w^4] do
627    Append(~params8,[0,v,0,x,1,0]);
628  end for;
629  end for;
630end if;
631
632if p mod 3 eq 2 then
633  for v in [0..p-1] do
634  for x in [0..p-1] do
635    Append(~params8,[0,v,1,x,1,0]);
636  end for;
637  end for;
638end if;
639
640if p mod 3 eq 1 then
641  xrange:=transversal3;
642  Append(~xrange,0);
643  for v in [0..p-1] do
644  for x in xrange do
645    Append(~params8,[0,v,1,x,1,0]);
646    Append(~params8,[0,v,w,x,1,0]);
647    Append(~params8,[0,v,w^2,x,1,0]);
648  end for;
649  end for;
650end if;
651
652for v in [0..p-1] do
653for x in [0..(p-1) div 2] do
654for z in [1,w] do
655  Append(~params8,[0,v,0,x,1,z]);
656end for;
657end for;
658end for;
659
660for v in [0..p-1] do
661for t in [1..(p-1) div 2] do
662for x in [0..p-1] do
663for z in [1,w] do
664  Append(~params8,[0,v,t,x,1,z]);
665end for;
666end for;
667end for;
668end for;
669
670print #params8,
671p^3+4*p^2+6*p+(p+5)*Gcd(p-1,3)+3*Gcd(p-1,4)+Gcd(p-1,5);
672gtotal:=gtotal+#params8;
673
674/*
675Descendants of 5.14.
676
677Cases 9 & 10: cb=bac=caa=0, cac=kbab (k=1,w)
678*/
679
680Append(~params9,[0,0,0,0,0,0]);
681Append(~params9,[1,0,0,0,0,0]);
682Append(~params9,[0,1,0,0,0,0]);
683Append(~params9,[0,w,0,0,0,0]);
684
685for y in [0,1,w] do
686  Append(~params9,[0,y,0,0,1,0]);
687  if p mod 4 eq 1 then
688    Append(~params9,[0,y,0,0,w,0]);
689   end if;
690end for;
691yrange:=[0..p-1];
692if p mod 4 eq 1 then yrange:=[0..(p-1) div 2]; end if;
693for y in yrange do
694  Append(~params9,[1,y,0,0,1,0]);
695  if p mod 4 eq 1 then
696    Append(~params9,[1,y,0,0,w,0]);
697  end if;
698end for;
699
700for x in [0..(p-1) div 2] do
701  Append(~params9,[x,0,0,0,0,1]);
702end for;
703Append(~params9,[0,1,0,0,0,1]);
704Append(~params9,[0,w,0,0,0,1]);
705
706for y in [0,1,w] do
707for z in [0..(p-1) div 2] do
708  Append(~params9,[0,y,1,0,z,0]);
709  Append(~params9,[0,y,w,0,z,0]);
710end for;
711end for;
712
713for y in [0..p-1] do
714for z in [0..(p-1) div 2] do
715  Append(~params9,[0,y,1,0,z,1]);
716  Append(~params9,[0,y,w,0,z,1]);
717end for;
718end for;
719
720for y in [0..p-1] do
721for t in [0..(p-1) div 2] do
722  Append(~params9,[1,y,1,0,0,t]);
723  Append(~params9,[1,y,w,0,0,t]);
724end for;
725end for;
726
727for y in [0..p-1] do
728for z in [1..(p-1) div 2] do
729for t in [0..p-1] do
730  Append(~params9,[1,y,1,0,z,t]);
731  Append(~params9,[1,y,w,0,z,t]);
732end for;
733end for;
734end for;
735
736for x in [0..(p-1) div 2] do
737for y in [0..p-1] do
738  Append(~params9,[y,0,0,1,0,x]);
739end for;
740end for;
741for x in [0..(p-1) div 2] do
742  Append(~params9,[1,1,0,1,0,x]);
743  Append(~params9,[1,w,0,1,0,x]);
744end for;
745
746if p mod 4 eq 1 then
747  for x in [0..p-1] do
748  for y in [0..(p-1) div 2] do
749    Append(~params9,[x,y,0,1,1,0]);
750    Append(~params9,[x,y,0,1,w,0]);
751  end for;
752  end for;
753end if;
754
755if p mod 4 eq 3 then
756  for x in [0..p-1] do
757  for y in [0..p-1] do
758    Append(~params9,[x,y,0,1,1,0]);
759  end for;
760  end for;
761end if;
762
763print #params9,
764p^3+5*p^2/2+7*p+19/2+(p+4)*Gcd(p-1,4)/2;
765gtotal:=gtotal+#params9;
766
767params10:=params9;
768gtotal:=gtotal+#params9;
769
770/*
771Descendants of 5.14.
772
773Cases 11 & 12: bac=caa=0, cb=baa, cac=kbab (k=1,w)
774*/
775
776Append(~params11,[0,0,0,0,0,0]);
777Append(~params11,[1,0,0,0,0,0]);
778if p mod 3 eq 1 then
779  Append(~params11,[w,0,0,0,0,0]);
780  Append(~params11,[w^2,0,0,0,0,0]);
781end if;
782Append(~params11,[0,1,0,0,0,0]);
783Append(~params11,[0,w,0,0,0,0]);
784if p mod 4 eq 1 then
785  Append(~params11,[0,w^2,0,0,0,0]);
786  Append(~params11,[0,w^3,0,0,0,0]);
787end if;
788
789if p mod 4 eq 1 then
790  xrange:=transversal4;
791  Append(~xrange,0);
792  for x in xrange do
793  for y in [0..p-1] do
794    Append(~params11,[x,y,0,0,1,0]);
795    Append(~params11,[x,y,0,0,w,0]);
796  end for;
797  end for;
798end if;
799if p mod 4 eq 3 then
800  for x in [0..(p-1) div 2] do
801  for y in [0..p-1] do
802    Append(~params11,[x,y,0,0,1,0]);
803  end for;
804  end for;
805end if;
806
807if p mod 3 eq 1 then
808  for x in [0..(p-1) div 2] do
809    Append(~params11,[x,0,0,0,0,1]);
810    Append(~params11,[x,0,0,0,0,w]);
811    Append(~params11,[x,0,0,0,0,w^2]);
812  end for;
813  for x in transversal3 do
814    Append(~params11,[0,x,0,0,0,1]);
815    Append(~params11,[0,x,0,0,0,w]);
816    Append(~params11,[0,x,0,0,0,w^2]);
817  end for;
818end if;
819if p mod 3 eq 2 then
820  for x in [0..(p-1) div 2] do
821    Append(~params11,[x,0,0,0,0,1]);
822  end for;
823  for x in [1..p-1] do
824    Append(~params11,[0,x,0,0,0,1]);
825  end for;
826end if;
827
828for y in [0..p-1] do
829for z in [0..(p-1) div 2] do
830for t in [0..(p-1) div 2] do
831  Append(~params11,[0,y,1,0,z,t]);
832  Append(~params11,[0,y,w,0,z,t]);
833end for;
834end for;
835end for;
836for y in [0..p-1] do
837for x in [1..(p-1) div 2] do
838for z in [0..(p-1) div 2] do
839  Append(~params11,[x,y,1,0,z,0]);
840  Append(~params11,[x,y,w,0,z,0]);
841end for;
842end for;
843end for;
844for y in [0..p-1] do
845for x in [1..(p-1) div 2] do
846for t in [1..(p-1) div 2] do
847for z in [0..p-1] do
848  Append(~params11,[x,y,1,0,z,t]);
849  Append(~params11,[x,y,w,0,z,t]);
850end for;
851end for;
852end for;
853end for;
854
855
856if p mod 3 eq 1 then
857  for x in [0..p-1] do
858    Append(~params11,[x,0,0,1,0,0]);
859    Append(~params11,[x,0,0,w,0,0]);
860    Append(~params11,[x,0,0,w^2,0,0]);
861  end for;
862  for y in transversal3 do
863    Append(~params11,[1,y,0,1,0,0]);
864    Append(~params11,[w,y,0,w,0,0]);
865    Append(~params11,[w^2,y,0,w^2,0,0]);
866  end for;
867  for x in transversal6 do
868  for y in [0..p-1] do
869  for z in [0..p-1] do
870    Append(~params11,[y,z,0,1,x,0]);
871    Append(~params11,[y,z,0,w,x,0]);
872    Append(~params11,[y,z,0,w^2,x,0]);
873  end for;
874  end for;
875  end for;
876  for x in [0..p-1] do
877  for z in [1..(p-1) div 2] do
878    Append(~params11,[x,0,0,1,0,z]);
879    Append(~params11,[x,0,0,w,0,z]);
880    Append(~params11,[x,0,0,w^2,0,z]);
881  end for;
882  end for;
883  for y in transversal3 do
884  for z in [1..(p-1) div 2] do
885    Append(~params11,[1,y,0,1,0,z]);
886    Append(~params11,[w,y,0,w,0,z]);
887    Append(~params11,[w^2,y,0,w^2,0,z]);
888  end for;
889  end for;
890end if;
891
892if p mod 3 eq 2 then
893  for x in [0..p-1] do
894    Append(~params11,[x,0,0,1,0,0]);
895  end for;
896  for y in [1..p-1] do
897    Append(~params11,[1,y,0,1,0,0]);
898  end for;
899  for x in [1..(p-1) div 2] do
900  for y in [0..p-1] do
901  for z in [0..p-1] do
902    Append(~params11,[y,z,0,1,x,0]);
903  end for;
904  end for;
905  end for;
906  for x in [0..p-1] do
907  for z in [1..(p-1) div 2] do
908    Append(~params11,[x,0,0,1,0,z]);
909  end for;
910  end for;
911  for y in [1..p-1] do
912  for z in [1..(p-1) div 2] do
913    Append(~params11,[1,y,0,1,0,z]);
914  end for;
915  end for;
916end if;
917
918
919print #params11,
920(p^4+p^3+4*p^2+p-1+(p^2+2*p+3)*Gcd(p-1,3)+(p+2)*Gcd(p-1,4))/2;
921gtotal:=gtotal+#params11;
922
923
924params12:=params11;
925gtotal:=gtotal+#params11;
926
927/*
928Descendants of 5.14.
929
930Case 13, cb=bac=0, caa=baa, cac=-bab
931*/
932
933//pb=pc=0
934Append(~params13,[0,0,0,0,0,0]);
935Append(~params13,[1,0,0,0,0,0]);
936Append(~params13,[0,1,0,0,0,0]);
937
938//pb=0, pc=baa or wbaa
939Append(~params13,[0,0,0,0,1,0]);
940Append(~params13,[0,0,0,0,w,0]);
941Append(~params13,[0,1,0,0,1,0]);
942Append(~params13,[0,1,0,0,w,0]);
943Append(~params13,[0,w,0,0,1,0]);
944Append(~params13,[0,w,0,0,w,0]);
945for y in [0..p-1] do
946  Append(~params13,[1,y,0,0,1,0]);
947  Append(~params13,[1,y,0,0,w,0]);
948end for;
949
950//pb=0, pc=bab
951for x in [0..p-1] do
952  Append(~params13,[x,0,0,0,0,1]);
953end for;
954Append(~params13,[-1,1,0,0,0,1]);
955Append(~params13,[-1,w,0,0,0,1]);
956
957//pb=pc=bab
958Append(~params13,[0,0,0,1,0,1]);
959Append(~params13,[1,0,0,1,0,1]);
960Append(~params13,[0,1,0,1,0,1]);
961
962//pb=bab, pc=-bab
963for x in [0..p-1] do
964  Append(~params13,[x,0,0,1,0,-1]);
965end for;
966Append(~params13,[2,1,0,1,0,-1]);
967
968//pb=bab, pc=baa or wbaa
969for x in [0..p-1] do
970for y in [0..p-1] do
971  Append(~params13,[x,y,0,1,1,0]);
972  Append(~params13,[x,y,0,1,w,0]);
973end for;
974end for;
975
976//pb=pc=baa
977Append(~params13,[0,0,1,0,1,0]);
978Append(~params13,[1,0,1,0,1,0]);
979Append(~params13,[0,1,1,0,1,0]);
980Append(~params13,[1,1,1,0,1,0]);
981
982//pb=baa, pc=baa+bab
983for x in [0..p-2] do //Note that -1 is excluded
984  Append(~params13,[x,0,1,0,1,1]);
985  y:=F!((x+1)/2); y:=Z!y;
986  Append(~params13,[x,y,1,0,1,1]);
987end for;
988Append(~params13,[-1,0,1,0,1,1]);
989Append(~params13,[-1,1,1,0,1,1]);
990
991//pb=baa, pc=-baa
992Append(~params13,[0,0,1,0,-1,0]);
993Append(~params13,[0,1,1,0,-1,0]);
994Append(~params13,[0,w,1,0,-1,0]);
995for y in [0..p-1] do
996  Append(~params13,[1,y,1,0,-1,0]);
997end for;
998
999//pb=baa, pc=-baa+bab
1000half:=(p+1) div 2;
1001Append(~params13,[1,-half,1,0,-1,1]);
1002Append(~params13,[1,half,1,0,-1,1]);
1003Append(~params13,[1,w-half,1,0,-1,1]);
1004if p mod 4 eq 1 then
1005  Append(~params13,[1,w^2-half,1,0,-1,1]);
1006  Append(~params13,[1,w^3-half,1,0,-1,1]);
1007end if;
1008for y in [0..p-1] do
1009  Append(~params13,[0,y,1,0,-1,1]);
1010  Append(~params13,[1-w,y,1,0,-1,1]);
1011end for;
1012
1013//pb=pc=wbaaa
1014Append(~params13,[0,0,w,0,w,0]);
1015Append(~params13,[1,0,w,0,w,0]);
1016Append(~params13,[0,1,w,0,w,0]);
1017Append(~params13,[1,1,w,0,w,0]);
1018
1019//pb=wbaa, pc=wbaa+bab
1020for x in [0..p-2] do //Note that -1 is excluded
1021  Append(~params13,[x,0,w,0,w,1]);
1022  y:=F!((x+1)/(2*w)); y:=Z!y;
1023  Append(~params13,[x,y,w,0,w,1]);
1024end for;
1025Append(~params13,[-1,0,w,0,w,1]);
1026Append(~params13,[-1,1,w,0,w,1]);
1027
1028print #params13,2*p^2+11*p+27+Gcd(p-1,4);
1029gtotal:=gtotal+#params13;
1030
1031/*
1032Descendants of 5.14.
1033
1034Case 14, bac=0, cb=caa=baa, cac=-bab
1035*/
1036
1037//1. v ne +/- y
1038for y in [1,w] do
1039for u in [0..p-1] do
1040for x in [0..(p-1) div 2] do
1041  zrange:=[0..p-1];
1042  if x eq 0 then zrange:=[0..(p-1) div 2]; end if;
1043  for z in zrange do
1044    Append(~params14,[0,u,0,x,y,z]);
1045  end for;
1046end for;
1047end for;
1048end for;
1049
1050//1.5 t=0, v=1, y=-1, u ne 0
1051for u in [1,w] do
1052for x in [0..(p-1) div 2] do
1053  zrange:=[0..p-1];
1054  if x eq 0 then zrange:=[0..(p-1) div 2]; end if;
1055  for z in zrange do
1056    Append(~params14,[0,u,1,x,-1,z]);
1057  end for;
1058end for;
1059end for;
1060
1061//1.5 t=0, v=1, y=-1, u=0
1062for x in [0..p-1] do
1063for z in [0..p-1] do
1064  new:=true;
1065  for a in [1..p-1] do
1066    x1:=F!((x+z+x*a^2-z*a^2)/(2*a^3)); x1:=Z!x1;
1067    z1:=F!((x+z-x*a^2+z*a^2)/(2*a^3)); z1:=Z!z1;
1068    if [x1,z1] lt [x,z] then new:=false; break; end if;
1069  end for;
1070  if new then Append(~params14,[0,0,1,x,-1,z]); end if;
1071end for;
1072end for;
1073
1074
1075
1076//2. t ne 0, v=y, x ne z+1, v ne 0
1077for v in [1,w] do
1078for x in [0..p-1] do
1079for z in [0..p-1] do
1080  if x eq (z+1) mod p then continue; end if;
1081  x1:=(-x) mod p; z1:=(-z) mod p;
1082  if [z1,x1] lt [x,z] then continue; end if;
1083  Append(~params14,[1,0,v,x,v,z]);
1084end for;
1085end for;
1086end for;
1087
1088//2. t ne 0, v=y=0, x ne z+1
1089for x in [0..p-1] do
1090for z in [0..p-1] do
1091  if x eq (z+1) mod p then continue; end if;
1092  new:=true;
1093  for a in [1..p-1] do
1094    x1:=F!((x+z+x*a^3-z*a^3)/(2*a^3)); x1:=Z!x1;
1095    z1:=F!((x+z-x*a^3+z*a^3)/(2*a^3)); z1:=Z!z1;
1096    if [x1,z1] lt [x,z] then new:=false; break; end if;
1097  end for;
1098  if new then Append(~params14,[1,0,0,x,0,z]); end if;
1099end for;
1100end for;
1101
1102//3. t ne 0, v=y, x eq z+1, v ne 0
1103for v in [1,w] do
1104for u in [0..(p-1) div 2] do
1105  Append(~params14,[1,u,v,0,v,-1]);
1106end for;
1107end for;
1108
1109//3. t ne 0, v=y, x eq z+1, v eq 0, u ne 0
1110for z in [0..p-1] do
1111  Append(~params14,[1,1,0,z+1,0,z]);
1112end for;
1113
1114//3. t ne 0, v=y, x eq z+1, v eq 0, u eq 0
1115for z in [0..p-1] do
1116  new:=true;
1117  for a in [1..p-1] do
1118    z1:=F!((-a^3+2*z+1)/(2*a^3)); z1:=Z!z1;
1119    if z1 lt z then new:=false; break; end if;
1120  end for;
1121  if new then Append(~params14,[1,0,0,z+1,0,z]); end if;
1122end for;
1123
1124//4. t=0, v=y, x ne +/- z, v ne 0
1125for v in [1,w] do
1126for z in [1..(p-1) div 2] do
1127  Append(~params14,[0,0,v,0,v,z]);
1128end for;
1129end for;
1130
1131//4. t=0, v=y=0, x ne +/- z
1132Append(~params14,[0,0,0,0,0,1]);
1133if p mod 3 eq 1 then
1134  Append(~params14,[0,0,0,0,0,w]);
1135  Append(~params14,[0,0,0,0,0,w^2]);
1136end if;
1137
1138//4.5 t=0, v=y, x=-z ne 0
1139Append(~params14,[0,0,0,1,0,-1]);
1140Append(~params14,[0,0,w,1,w,-1]);
1141Append(~params14,[0,0,w^2,1,w^2,-1]);
1142
1143//5. t=0, v=y ne 0, x=z
1144for v in [1,w] do
1145for u in [0,1] do
1146  Append(~params14,[0,u,v,0,v,0]);
1147end for;
1148end for;
1149
1150
1151//6. t=0, v=y=0, x=z
1152xrange:=[0,1];
1153if p mod 3 eq 1 then xrange:=[0,1,w,w^2]; end if;
1154for u in [0,1] do
1155for x in xrange do
1156  Append(~params14,[0,u,0,x,0,x]);
1157end for;
1158end for;
1159
1160print #params14,p^3+2*p^2+6*p+10+(p+4)*Gcd(p-1,3);
1161gtotal:=gtotal+#params14;
1162
1163/*
1164Descendants of 5.14.
1165
1166Case 15, cb=baa=bac=caa=0
1167*/
1168
1169//x=y=z=v=0
1170Append(~params15,[0,0,0,0,0,0]);
1171Append(~params15,[0,0,0,0,1,0]);
1172Append(~params15,[0,0,0,1,1,0]);
1173if p mod 3 eq 1 then
1174  Append(~params15,[0,0,0,1,w,0]);
1175end if;
1176
1177//x=z=v=0, y=1,w
1178Append(~params15,[0,1,0,0,0,0]);
1179Append(~params15,[0,w,0,0,0,0]);
1180Append(~params15,[0,1,0,0,1,0]);
1181Append(~params15,[0,w,0,0,1,0]);
1182Append(~params15,[0,1,0,1,0,0]);
1183Append(~params15,[0,w,0,1,0,0]);
1184Append(~params15,[0,1,0,1,1,0]);
1185Append(~params15,[0,w,0,1,1,0]);
1186if  p mod 3 eq 1 then
1187  Append(~params15,[0,1,0,1,w,0]);
1188  Append(~params15,[0,w,0,1,w,0]);
1189  Append(~params15,[0,1,0,1,w^2,0]);
1190  Append(~params15,[0,w,0,1,w^2,0]);
1191end if;
1192
1193//z=v=0, x=y=1,w
1194Append(~params15,[1,1,0,0,0,0]);
1195Append(~params15,[w,w,0,0,0,0]);
1196Append(~params15,[1,1,0,0,1,0]);
1197Append(~params15,[w,w,0,0,1,0]);
1198for u in [1..(p-1) div 2] do
1199  u1:=F!u^-1; u1:=Z!u1;
1200  if u1 lt u or (p-u1) lt u then continue; end if;
1201  Append(~params15,[1,1,0,1,u,0]);
1202  Append(~params15,[w,w,0,1,u,0]);
1203end for;
1204
1205//z=v=0, x=1, y=w
1206Append(~params15,[1,w,0,0,0,0]);
1207Append(~params15,[1,w,0,0,1,0]);
1208Append(~params15,[1,w,0,1,0,0]);
1209for u in [1..(p-1) div 2] do
1210  Append(~params15,[1,w,0,1,u,0]);
1211end for;
1212
1213//z=0, v=1, t=0
1214Append(~params15,[0,0,0,0,0,1]);
1215Append(~params15,[0,0,0,0,1,1]);
1216Append(~params15,[0,0,0,0,w,1]);
1217Append(~params15,[0,1,0,0,0,1]);
1218Append(~params15,[0,1,0,0,1,1]);
1219Append(~params15,[0,1,0,0,w,1]);
1220Append(~params15,[0,w,0,0,0,1]);
1221Append(~params15,[0,w,0,0,1,1]);
1222Append(~params15,[0,w,0,0,w,1]);
1223Append(~params15,[1,0,0,0,0,1]);
1224Append(~params15,[1,0,0,0,1,1]);
1225Append(~params15,[1,0,0,0,w,1]);
1226Append(~params15,[w,0,0,0,0,1]);
1227Append(~params15,[w,0,0,0,1,1]);
1228Append(~params15,[w,0,0,0,w,1]);
1229for x in [1,w] do
1230for y in [1,w] do
1231for u in [0..p-1] do
1232  Append(~params15,[x,y,0,0,u,1]);
1233end for;
1234end for;
1235end for;
1236
1237//z=0, v=1, t=1
1238for x in [1,w] do
1239for y in [0..p-1] do
1240for u in [0..p-1] do
1241  Append(~params15,[x,y,0,1,u,1]);
1242end for;
1243end for;
1244end for;
1245for y in [0,1,w] do
1246for u in [0..p-1] do
1247  Append(~params15,[0,y,0,1,u,1]);
1248end for;
1249end for;
1250
1251//z=v=1
1252for t in [0..p-2] do
1253for u in [t+1..p-1] do
1254  Append(~params15,[0,0,1,t,u,1]);
1255  Append(~params15,[0,1,1,t,u,1]);
1256  Append(~params15,[0,w,1,t,u,1]);
1257  for y in [0..p-1] do
1258    Append(~params15,[1,y,1,t,u,1]);
1259    Append(~params15,[w,y,1,t,u,1]);
1260  end for;
1261end for;
1262end for;
1263
1264//Get set of non-zero squares with k^2 ~ k^-2
1265range:={};
1266for k in [1..(p-1) div 2] do
1267  k2:=k^2 mod p;
1268  km2:=F!(k2^-1); km2:=Z!km2;
1269  if k2 le km2 then Include(~range,k2); end if;
1270end for;
1271for t in [0..p-1] do
1272   Append(~params15,[0,0,1,t,t,1]);
1273   Append(~params15,[0,1,1,t,t,1]);
1274   Append(~params15,[0,w,1,t,t,1]);
1275   for k2 in range do
1276      Append(~params15,[1,k2,1,t,t,1]);
1277      Append(~params15,[w,w*k2,1,t,t,1]);
1278   end for;
1279   for k in [1..(p-1) div 2] do
1280      Append(~params15,[1,w*k^2,1,t,t,1]);
1281   end for;
1282end for;
1283
1284print #params15,
1285p^3+(7*p^2+17*p+59+5*Gcd(p-1,3)+(p+1)*Gcd(p-1,4))/2;
1286gtotal:=gtotal+#params15;
1287
1288/*
1289Descendants of 5.14.
1290
1291Case 16, cb=bac=caa=0, baa=cac
1292*/
1293
1294//u=v=z=0
1295for y in [0,1,w] do
1296  Append(~params16,[0,y,0,0,0,0]);
1297  Append(~params16,[0,y,0,1,0,0]);
1298  Append(~params16,[0,y,0,w,0,0]);
1299  for t in [0..p-1] do
1300    Append(~params16,[1,y,0,t,0,0]);
1301    Append(~params16,[w,y,0,t,0,0]);
1302  end for;
1303end for;
1304
1305//u=v=0, z=1,w
1306Append(~params16,[0,0,1,0,0,0]);
1307Append(~params16,[0,0,w,0,0,0]);
1308Append(~params16,[0,1,1,0,0,0]);
1309Append(~params16,[0,1,w,0,0,0]);
1310Append(~params16,[0,w,1,0,0,0]);
1311Append(~params16,[0,w,w,0,0,0]);
1312if p mod 4 eq 1 then
1313  Append(~params16,[0,w^2,1,0,0,0]);
1314  Append(~params16,[0,w^2,w,0,0,0]);
1315  Append(~params16,[0,w^3,1,0,0,0]);
1316  Append(~params16,[0,w^3,w,0,0,0]);
1317end if;
1318for y in [0..p-1] do
1319  Append(~params16,[0,y,1,1,0,0]);
1320  Append(~params16,[0,y,w,1,0,0]);
1321  Append(~params16,[0,y,1,w,0,0]);
1322  Append(~params16,[0,y,w,w,0,0]);
1323end for;
1324for y in [0..p-1] do
1325for t in [0..p-1] do
1326  Append(~params16,[1,y,1,t,0,0]);
1327  Append(~params16,[1,y,w,t,0,0]);
1328  Append(~params16,[w,y,1,t,0,0]);
1329  Append(~params16,[w,y,w,t,0,0]);
1330end for;
1331end for;
1332
1333//u=0, v=1
1334Append(~params16,[0,0,0,0,0,1]);
1335Append(~params16,[0,0,1,0,0,1]);
1336Append(~params16,[0,0,w,0,0,1]);
1337if p mod 3 eq 1 then
1338  Append(~params16,[0,0,w^2,0,0,1]);
1339  Append(~params16,[0,0,w^3,0,0,1]);
1340  Append(~params16,[0,0,w^4,0,0,1]);
1341  Append(~params16,[0,0,w^5,0,0,1]);
1342end if;
1343for z in [0..p-1] do
1344  Append(~params16,[0,0,z,1,0,1]);
1345  Append(~params16,[0,0,z,w,0,1]);
1346end for;
1347for z in [0..p-1] do
1348for t in [0..p-1] do
1349  Append(~params16,[0,1,z,t,0,1]);
1350  Append(~params16,[0,w,z,t,0,1]);
1351end for;
1352end for;
1353for y in [0..p-1] do
1354for z in [0..p-1] do
1355for t in [0..p-1] do
1356  Append(~params16,[1,y,z,t,0,1]);
1357  Append(~params16,[w,y,z,t,0,1]);
1358end for;
1359end for;
1360end for;
1361
1362//u=1
1363Append(~params16,[0,0,0,0,1,0]);
1364Append(~params16,[0,1,0,0,1,0]);
1365Append(~params16,[0,w,0,0,1,0]);
1366if p mod 3 eq 1 then
1367  Append(~params16,[0,w^2,0,0,1,0]);
1368  Append(~params16,[0,w^3,0,0,1,0]);
1369  Append(~params16,[0,w^4,0,0,1,0]);
1370  Append(~params16,[0,w^5,0,0,1,0]);
1371end if;
1372if p mod 4 eq 1 then
1373  for v in [1,w,w^2,w^3] do
1374  for y in [0..(p-1) div 2] do
1375    Append(~params16,[0,y,0,0,1,v]);
1376  end for;
1377  end for;
1378end if;
1379if p mod 4 eq 3 then
1380  for v in [1,w] do
1381  for y in [0..p-1] do
1382    Append(~params16,[0,y,0,0,1,v]);
1383  end for;
1384  end for;
1385end if;
1386for v in [0..p-1] do
1387for y in [0..p-1] do
1388  Append(~params16,[0,y,0,1,1,v]);
1389  Append(~params16,[0,y,0,w,1,v]);
1390end for;
1391end for;
1392for t in [0..p-1] do
1393for v in [0..p-1] do
1394for y in [0..p-1] do
1395  Append(~params16,[0,y,1,t,1,v]);
1396  Append(~params16,[0,y,w,t,1,v]);
1397end for;
1398end for;
1399end for;
1400for z in [0..p-1] do
1401for t in [0..p-1] do
1402for v in [0..p-1] do
1403for y in [0..p-1] do
1404  Append(~params16,[1,y,z,t,1,v]);
1405  Append(~params16,[w,y,z,t,1,v]);
1406end for;
1407end for;
1408end for;
1409end for;
1410
1411print #params16,
14122*p^4+4*p^3+8*p^2+14*p+11+4*Gcd(p-1,3)+3*Gcd(p-1,4);
1413gtotal:=gtotal+#params16;
1414
1415/*
1416Descendants of 5.14.
1417
1418Case 17, cb=bac=0, cac=baa, caa=bab
1419*/
1420
1421mats:=[];
1422//get pc,pb
1423
1424for y5 in [0,1,lns] do
1425for y6 in [0..p-1] do
1426for y3 in [0..p-1] do
1427for y4 in [0..p-1] do
1428
1429A:=H22![y3,y4,y5,y6];
1430if A eq 0 then
1431  Append(~mats,A);
1432  continue;
1433end if;
1434
1435new:=1;
1436index:=p^3*y5+p^2*y6+p*y3+y4;
1437
1438for a in [1..p-1] do
1439for x in CU do
1440c:=a*x;
1441
1442B:=H22![a^-1*c^2,0,0,c];
1443C:=H22![a*c^2,0,0,a^2*c];
1444
1445D:=B*A*C^-1;
1446
1447z3:=Z!(D[1][1]);
1448z4:=Z!(D[1][2]);
1449z5:=Z!(D[2][1]);
1450z6:=Z!(D[2][2]);
1451
1452
1453ind1:=p^3*z5+p^2*z6+p*z3+z4;
1454
1455if ind1 lt index then new:=0; end if;
1456
1457B:=H22![0,a^-1*c^2,c,0];
1458C:=H22![0,a*c^2,a^2*c,0];
1459
1460D:=B*A*C^-1;
1461
1462z3:=Z!(D[1][1]);
1463z4:=Z!(D[1][2]);
1464z5:=Z!(D[2][1]);
1465z6:=Z!(D[2][2]);
1466
1467
1468ind1:=p^3*z5+p^2*z6+p*z3+z4;
1469
1470if ind1 lt index then new:=0; end if;
1471
1472
1473if new eq 0 then break; end if;
1474end for;
1475if new eq 0 then break; end if;
1476end for;
1477
1478if new eq 1 then
1479  Append(~mats,A);
1480  //print y3,y4,y5,y6;
1481end if;
1482
1483end for;
1484end for;
1485end for;
1486end for;
1487
1488
1489for AS in mats do
1490//Get stabilizer of AS
1491stab1:=[];
1492stab2:=[];
1493for a in [1..p-1] do
1494for x in CU do
1495c:=a*x;
1496
1497B:=H22![a^-1*c^2,0,0,c];
1498C:=H22![a*c^2,0,0,a^2*c];
1499
1500D:=B*AS*C^-1;
1501
1502if D eq AS then
1503  Append(~stab1,H33![a,0,0,0,a^-1*c^2,0,0,0,c]);
1504  Append(~stab2,C);
1505end if;
1506
1507B:=H22![0,a^-1*c^2,c,0];
1508C:=H22![0,a*c^2,a^2*c,0];
1509
1510D:=B*AS*C^-1;
1511
1512if D eq AS then
1513  Append(~stab1,H33![a,0,0,0,0,a^-1*c^2,0,c,0]);
1514  Append(~stab2,C);
1515end if;
1516
1517end for;
1518end for;
1519
1520y3:=Z!(AS[1][1]); y4:=Z!(AS[1][2]); y5:=Z!(AS[2][1]); y6:=Z!(AS[2][2]);
1521
1522for y1 in [0..p-1] do
1523for y2 in [0..p-1] do
1524
1525A:=H32![y1,y2,y3,y4,y5,y6];
1526if A eq 0 then
1527  Append(~params17,[0,0,0,0,0,0]);
1528  continue;
1529end if;
1530
1531new:=1;
1532index:=p*y1+y2;
1533
1534for ii in [1..#stab1] do
1535
1536B:=stab1[ii];
1537C:=stab2[ii];
1538
1539D:=B*A*C^-1;
1540
1541z1:=Z!(D[1][1]);
1542z2:=Z!(D[1][2]);
1543
1544ind1:=p*z1+z2;
1545
1546if ind1 lt index then new:=0; end if;
1547
1548if new eq 0 then break; end if;
1549end for;
1550
1551
1552if new eq 1 then
1553  Append(~params17,[y1,y2,y3,y4,y5,y6]);
1554end if;
1555
1556end for;
1557end for;
1558
1559end for;
1560
1561print #params17,
1562(p^4+2*p^3+3*p^2+4*p+2)*(p-1)/Gcd(p-1,3)+3*p+4+(p^2+p+1)*Gcd(p-1,4)/2;
1563gtotal:=gtotal+#params17;
1564
1565/*
1566Descendants of 5.14.
1567
1568Case 18, cb=bac=0, cac=baa, caa=wbab (p=1mod3 only)
1569*/
1570if p mod 3 eq 1 then
1571
1572mats:=[];
1573//get pc,pb
1574
1575for y5 in [0,1,lns] do
1576for y6 in [0..p-1] do
1577for y3 in [0..p-1] do
1578for y4 in [0..p-1] do
1579
1580A:=H22![y3,y4,y5,y6];
1581if A eq 0 then
1582  Append(~mats,A);
1583  continue;
1584end if;
1585
1586new:=1;
1587index:=p^3*y5+p^2*y6+p*y3+y4;
1588
1589for a in [1..p-1] do
1590for x in CU do
1591c:=a*x;
1592
1593B:=H22![a^-1*c^2,0,0,c];
1594C:=H22![a*c^2,0,0,a^2*c];
1595
1596D:=B*A*C^-1;
1597
1598z3:=Z!(D[1][1]);
1599z4:=Z!(D[1][2]);
1600z5:=Z!(D[2][1]);
1601z6:=Z!(D[2][2]);
1602
1603
1604ind1:=p^3*z5+p^2*z6+p*z3+z4;
1605
1606if ind1 lt index then new:=0; end if;
1607
1608
1609if new eq 0 then break; end if;
1610end for;
1611if new eq 0 then break; end if;
1612end for;
1613
1614if new eq 1 then
1615  Append(~mats,A);
1616  //print y3,y4,y5,y6;
1617end if;
1618
1619end for;
1620end for;
1621end for;
1622end for;
1623
1624for AS in mats do
1625//Get stabilizer of AS
1626stab1:=[];
1627stab2:=[];
1628for a in [1..p-1] do
1629for x in CU do
1630c:=a*x;
1631
1632B:=H22![a^-1*c^2,0,0,c];
1633C:=H22![a*c^2,0,0,a^2*c];
1634
1635D:=B*AS*C^-1;
1636
1637if D eq AS then
1638  Append(~stab1,H33![a,0,0,0,a^-1*c^2,0,0,0,c]);
1639  Append(~stab2,C);
1640end if;
1641
1642end for;
1643end for;
1644
1645y3:=Z!(AS[1][1]); y4:=Z!(AS[1][2]); y5:=Z!(AS[2][1]); y6:=Z!(AS[2][2]);
1646
1647for y1 in [0..p-1] do
1648for y2 in [0..p-1] do
1649
1650A:=H32![y1,y2,y3,y4,y5,y6];
1651if A eq 0 then
1652  Append(~params18,[0,0,0,0,0,0]);
1653  continue;
1654end if;
1655
1656new:=1;
1657index:=p*y1+y2;
1658
1659for ii in [1..#stab1] do
1660
1661B:=stab1[ii];
1662C:=stab2[ii];
1663
1664D:=B*A*C^-1;
1665
1666z1:=Z!(D[1][1]);
1667z2:=Z!(D[1][2]);
1668
1669ind1:=p*z1+z2;
1670
1671if ind1 lt index then new:=0; end if;
1672
1673if new eq 0 then break; end if;
1674end for;
1675
1676
1677if new eq 1 then
1678  Append(~params18,[y1,y2,y3,y4,y5,y6]);
1679end if;
1680
1681end for;
1682end for;
1683
1684end for;
1685
1686print #params18,
1687(2*p^5+2*p^4+2*p^3+2*p^2+14*p+17)/3;
1688gtotal:=gtotal+#params18;
1689
1690end if;
1691
1692/*
1693Descendants of 5.14.
1694
1695Case 19, cb=baa=caa=cac=0
1696*/
1697
1698Append(~params19,[0,0,0,0,0,0]);
1699Append(~params19,[1,0,0,0,0,0]);
1700Append(~params19,[w,0,0,0,0,0]);
1701Append(~params19,[0,1,0,0,0,0]);
1702
1703for x in [0,1,w] do
1704for y in [0,1] do
1705Append(~params19,[x,y,1,0,0,0]);
1706end for;
1707end for;
1708
1709Append(~params19,[0,0,0,1,0,0]);
1710Append(~params19,[1,0,0,1,0,0]);
1711Append(~params19,[w,0,0,1,0,0]);
1712Append(~params19,[0,1,0,1,0,0]);
1713
1714for t in [0..p-1] do
1715  Append(~params19,[0,0,1,t,1,0]);
1716  Append(~params19,[0,1,1,t,1,0]);
1717  Append(~params19,[0,w,1,t,1,0]);
1718  for y in [0..p-1] do
1719    Append(~params19,[1,y,1,t,1,0]);
1720    Append(~params19,[w,y,1,t,1,0]);
1721  end for;
1722end for;
1723
1724for x in [0,1,w] do
1725for t in [0,1,w] do
1726  Append(~params19,[x,0,0,t,1,0]);
1727end for;
1728end for;
1729
1730for x in [1,w] do
1731for t in [0..p-1] do
1732  Append(~params19,[x,1,0,t,1,0]);
1733end for;
1734end for;
1735
1736Append(~params19,[0,1,0,0,1,0]);
1737Append(~params19,[0,1,0,1,1,0]);
1738Append(~params19,[0,1,0,w,1,0]);
1739if p mod 4 eq 1 then
1740  Append(~params19,[0,1,0,w^2,1,0]);
1741  Append(~params19,[0,1,0,w^3,1,0]);
1742end if;
1743
1744half:=(p+1) div 2;
1745for x in [0..p-1] do
1746  if x eq half then continue; end if;
1747  for z in [0,1,w] do
1748  for t in [0,1] do
1749    Append(~params19,[z,t,x,0,0,1]);
1750  end for;
1751  end for;
1752end for;
1753
1754Append(~params19,[1,0,half,0,0,1]);
1755Append(~params19,[w,0,half,0,0,1]);
1756Append(~params19,[1,0,half,1,0,1]);
1757Append(~params19,[w,0,half,1,0,1]);
1758
1759Append(~params19,[0,0,half,0,0,1]);
1760Append(~params19,[0,0,half,1,0,1]);
1761Append(~params19,[0,1,half,0,0,1]);
1762Append(~params19,[0,1,half,1,0,1]);
1763Append(~params19,[0,w,half,1,0,1]);
1764
1765
1766print #params19,2*p^2+11*p+27+Gcd(p-1,4);
1767gtotal:=gtotal+#params19;
1768
1769/*
1770Descendants of 5.14.
1771
1772Case 20, cb=baa=cac, caa=bab
1773*/
1774
1775//z=v=y=0
1776for u in [0,1,w] do
1777  Append(~params20,[0,0,0,0,u,0]);
1778  Append(~params20,[0,0,0,1,u,0]);
1779  Append(~params20,[0,0,0,w,u,0]);
1780  for t in [0..p-1] do
1781    Append(~params20,[1,0,0,t,u,0]);
1782    Append(~params20,[w,0,0,t,u,0]);
1783  end for;
1784end for;
1785
1786//z=v=0, y=1, x=t=0
1787Append(~params20,[0,1,0,0,0,0]);
1788Append(~params20,[0,1,0,0,1,0]);
1789Append(~params20,[0,1,0,0,w,0]);
1790if p mod 3 eq 1 then
1791  Append(~params20,[0,1,0,0,w^2,0]);
1792  Append(~params20,[0,1,0,0,w^3,0]);
1793  Append(~params20,[0,1,0,0,w^4,0]);
1794  Append(~params20,[0,1,0,0,w^5,0]);
1795end if;
1796
1797//z=v=0, y=1, x=0, t=1,w
1798for u in [0..p-1] do
1799  Append(~params20,[0,1,0,1,u,0]);
1800  Append(~params20,[0,1,0,w,u,0]);
1801end for;
1802
1803//z=v=0, y=1, x=1,w
1804for t in [0..p-1] do
1805for u in [0..p-1] do
1806  Append(~params20,[1,1,0,t,u,0]);
1807  Append(~params20,[w,1,0,t,u,0]);
1808end for;
1809end for;
1810
1811//z=0, v=1, x=t=u=0
1812Append(~params20,[0,0,0,0,0,1]);
1813Append(~params20,[0,1,0,0,0,1]);
1814Append(~params20,[0,w,0,0,0,1]);
1815if p mod 4 eq 1 then
1816  Append(~params20,[0,w^2,0,0,0,1]);
1817  Append(~params20,[0,w^3,0,0,0,1]);
1818end if;
1819
1820//z=0, v=1, x=t=0, u=1,w
1821for y in [0..p-1] do
1822  Append(~params20,[0,y,0,0,1,1]);
1823  Append(~params20,[0,y,0,0,w,1]);
1824end for;
1825
1826//z=0, v=1, x=0, t=1,w
1827for y in [0..p-1] do
1828for u in [0..p-1] do
1829  Append(~params20,[0,y,0,1,u,1]);
1830  Append(~params20,[0,y,0,w,u,1]);
1831end for;
1832end for;
1833
1834//z=0, v=1, x=1,w
1835for y in [0..p-1] do
1836for t in [0..p-1] do
1837for u in [0..p-1] do
1838  Append(~params20,[1,y,0,t,u,1]);
1839  Append(~params20,[w,y,0,t,u,1]);
1840end for;
1841end for;
1842end for;
1843
1844//z=1
1845for v in [0..p-1] do
1846
1847Append(~params20,[0,0,1,0,0,v]);
1848Append(~params20,[0,1,1,0,0,v]);
1849Append(~params20,[0,w,1,0,0,v]);
1850if p mod 4 eq 1 then
1851  Append(~params20,[0,w^2,1,0,0,v]);
1852  Append(~params20,[0,w^3,1,0,0,v]);
1853end if;
1854
1855for y in [0..p-1] do
1856  Append(~params20,[0,y,1,0,1,v]);
1857  Append(~params20,[0,y,1,0,w,v]);
1858end for;
1859
1860for y in [0..p-1] do
1861for u in [0..p-1] do
1862  Append(~params20,[0,y,1,1,u,v]);
1863  Append(~params20,[0,y,1,w,u,v]);
1864end for;
1865end for;
1866
1867for y in [0..p-1] do
1868for t in [0..p-1] do
1869for u in [0..p-1] do
1870  Append(~params20,[1,y,1,t,u,v]);
1871  Append(~params20,[w,y,1,t,u,v]);
1872end for;
1873end for;
1874end for;
1875
1876end for;
1877
1878
1879print #params20,
18802*p^4+4*p^3+6*p^2+11*p+11+2*Gcd(p-1,3)+(p+1)*Gcd(p-1,4);
1881gtotal:=gtotal+#params20;
1882
1883/*
1884Descendants of 5.14.
1885
1886Case 21, cb=caa=cac=0, bab=baa
1887*/
1888
1889Append(~params21,[0,0,0,0,0,0]);
1890Append(~params21,[1,0,0,0,0,0]);
1891Append(~params21,[w,0,0,0,0,0]);
1892Append(~params21,[0,1,0,0,0,0]);
1893
1894for x in [0..p-1] do
1895  Append(~params21,[x,0,1,0,0,0]);
1896  Append(~params21,[x,0,w,0,0,0]);
1897  Append(~params21,[x,1,1,0,0,0]);
1898  Append(~params21,[x,1,w,0,0,0]);
1899end for;
1900
1901Append(~params21,[1,0,0,1,0,0]);
1902Append(~params21,[w,0,0,1,0,0]);
1903for x in [0..p-1] do
1904  Append(~params21,[0,x,0,1,0,0]);
1905end for;
1906
1907//pc=baa
1908for y in [0..p-1] do
1909for z in [0..p-1] do
1910for t in [0..p-1] do
1911  Append(~params21,[1,y,z,t,1,0]);
1912  Append(~params21,[w,y,z,t,1,0]);
1913end for;
1914end for;
1915end for;
1916
1917for y in [0..p-1] do
1918for t in [0..p-1] do
1919  Append(~params21,[0,y,1,t,1,0]);
1920  Append(~params21,[0,y,w,t,1,0]);
1921end for;
1922end for;
1923
1924for t in [0..p-1] do
1925  Append(~params21,[0,1,0,t,1,0]);
1926  Append(~params21,[0,w,0,t,1,0]);
1927  if p mod 4 eq 1 then
1928    Append(~params21,[0,w^2,0,t,1,0]);
1929    Append(~params21,[0,w^3,0,t,1,0]);
1930  end if;
1931end for;
1932
1933Append(~params21,[0,0,0,0,1,0]);
1934Append(~params21,[0,0,0,1,1,0]);
1935Append(~params21,[0,0,0,w,1,0]);
1936if p mod 4 eq 1 then
1937  Append(~params21,[0,0,0,w^2,1,0]);
1938  Append(~params21,[0,0,0,w^3,1,0]);
1939end if;
1940
1941for x in [0..p-1] do
1942if x eq 1 then continue; end if;
1943for z in [0..p-1] do
1944for t in [0,1] do
1945Append(~params21,[x,0,z,t,0,1]);
1946end for;
1947end for;
1948end for;
1949
1950for x in [0..p-1] do
1951if x eq w then continue; end if;
1952for z in [0..p-1] do
1953for t in [0,1] do
1954Append(~params21,[x,0,z,t,0,w]);
1955end for;
1956end for;
1957end for;
1958
1959half:=(p+1) div 2;
1960wover2:=F!(w/2); wover2:=Z!wover2;
1961for z in [0..p-1] do
1962  if z eq half then continue; end if;
1963  Append(~params21,[1,0,z,0,0,1]);
1964  Append(~params21,[1,1,z,0,0,1]);
1965end for;
1966for z in [0..p-1] do
1967  if z eq wover2 then continue; end if;
1968  Append(~params21,[w,0,z,0,0,w]);
1969  Append(~params21,[w,1,z,0,0,w]);
1970end for;
1971
1972Append(~params21,[1,0,half,0,0,1]);
1973Append(~params21,[w,0,wover2,0,0,w]);
1974Append(~params21,[1,0,half,1,0,1]);
1975Append(~params21,[w,0,wover2,1,0,w]);
1976for t in [0..p-1] do
1977  Append(~params21,[1,1,half,t,0,1]);
1978  Append(~params21,[w,1,wover2,t,0,w]);
1979end for;
1980
1981print #params21,
19822*p^3+6*p^2+7*p+7+(p+1)*Gcd(p-1,4);
1983gtotal:=gtotal+#params21;
1984
1985/*
1986Descendants of 5.14.
1987
1988Case 22, cb=baa=caa=0, cac=wbab
1989*/
1990
1991mats:=[];
1992//get pb,pc
1993
1994bcrange:=[[0,1]];
1995for i in [0..p-1] do
1996  Append(~bcrange,[1,i]);
1997end for;
1998
1999for y3 in [0,1] do
2000for y4 in [0..p-1] do
2001for y5 in [0..p-1] do
2002for y6 in [0..p-1] do
2003
2004A:=H22![y3,y4,y5,y6];
2005if A eq 0 then
2006  Append(~mats,A);
2007  continue;
2008end if;
2009
2010new:=1;
2011index:=p^3*y3+p^2*y4+p*y5+y6;
2012
2013for bc in bcrange do
2014b:=bc[1]; c:=bc[2];
2015for m in [-1,1] do
2016
2017B:=H22![w*b,m*c,w*c,w*m*b];
2018C:=H22![w*(w*b^2+c^2),m*2*w*b*c,2*w^2*b*c,m*w*(w*b^2+c^2)];
2019
2020D:=B*A*C^-1;
2021if D[1][1] ne 0 then D:=D[1][1]^-1*D; end if;
2022if D[1][1] eq 0 and D[1][2] ne 0 then D:=D[1][2]^-1*D; end if;
2023if D[1][1] eq 0 and D[1][2] eq 0 and D[2][1] ne 0 then D:=D[2][1]^-1*D; end if;
2024if D[1][1] eq 0 and D[1][2] eq 0 and D[2][1] eq 0 and D[2][2] ne 0 then D:=D[2][2]^-1*D; end if;
2025
2026z3:=Z!(D[1][1]);
2027z4:=Z!(D[1][2]);
2028z5:=Z!(D[2][1]);
2029z6:=Z!(D[2][2]);
2030
2031
2032ind1:=p^3*z3+p^2*z4+p*z5+z6;
2033
2034if ind1 lt index then new:=0; end if;
2035
2036
2037if new eq 0 then break; end if;
2038end for;
2039if new eq 0 then break; end if;
2040end for;
2041
2042if new eq 1 then
2043  Append(~mats,A);
2044  //print y3,y4,y5,y6;
2045end if;
2046
2047end for;
2048end for;
2049end for;
2050end for;
2051
2052for AS in mats do
2053//Get stabilizer of AS
2054stab1:=[];
2055stab2:=[];
2056for b in [0..p-1] do
2057for c in [0..p-1] do
2058for m in [-1,1] do
2059if b+c eq 0 then continue; end if;
2060
2061B:=H22![w*b,m*c,w*c,w*m*b];
2062C:=H22![w*(w*b^2+c^2),m*2*w*b*c,2*w^2*b*c,m*w*(w*b^2+c^2)];
2063
2064D:=B*AS*C^-1;
2065if D[1][1] ne 0 then D:=D[1][1]^-1*D; end if;
2066if D[1][1] eq 0 and D[1][2] ne 0 then D:=D[1][2]^-1*D; end if;
2067if D[1][1] eq 0 and D[1][2] eq 0 and D[2][1] ne 0 then D:=D[2][1]^-1*D; end if;
2068if D[1][1] eq 0 and D[1][2] eq 0 and D[2][1] eq 0 and D[2][2] ne 0 then D:=D[2][2]^-1*D; end if;
2069
2070if D eq AS then
2071  Append(~stab1,H33![1,0,0,0,w*b,m*c,0,w*c,w*m*b]);
2072  Append(~stab2,C);
2073end if;
2074
2075end for;
2076end for;
2077end for;
2078
2079y3:=Z!(AS[1][1]); y4:=Z!(AS[1][2]); y5:=Z!(AS[2][1]); y6:=Z!(AS[2][2]);
2080
2081for y1 in [0..p-1] do
2082for y2 in [0..p-1] do
2083
2084A:=H32![y1,y2,y3,y4,y5,y6];
2085if A eq 0 then
2086  Append(~params22,[0,0,0,0,0,0]);
2087  continue;
2088end if;
2089
2090new:=1;
2091index:=p*y1+y2;
2092
2093for ii in [1..#stab1] do
2094
2095B:=stab1[ii];
2096C:=stab2[ii];
2097
2098D:=B*A*C^-1;
2099
2100z1:=Z!(D[1][1]);
2101z2:=Z!(D[1][2]);
2102
2103ind1:=p*z1+z2;
2104
2105if ind1 lt index then new:=0; end if;
2106
2107if new eq 0 then break; end if;
2108end for;
2109
2110
2111if new eq 1 then
2112  Append(~params22,[y1,y2,y3,y4,y5,y6]);
2113end if;
2114
2115end for;
2116end for;
2117
2118end for;
2119
2120print #params22,
2121(2*p^3+3*p^2+3*p+13-Gcd(p-1,3)+(p+1)*Gcd(p-1,4))/2;
2122gtotal:=gtotal+#params22;
2123
2124/*
2125Descendants of 5.14.
2126
2127Case 23, cb=baa=0, caa=bac, cac=wbab
2128*/
2129sol:=0;
2130if p mod 3 eq 2 then
2131  //look for solution of 12wx^2=-1
2132  for x in [1..p-1] do
2133    if F!(12*w*x^2+1) eq 0 then sol:=x; break; end if;
2134  end for;
2135end if;
2136
2137mats:=[];
2138//get pc,pb
2139
2140for y5 in [0,1,lns] do
2141for y6 in [0..p-1] do
2142for y3 in [0..p-1] do
2143for y4 in [0..p-1] do
2144
2145A:=H22![y3,y4,y5,y6];
2146if A eq 0 then
2147  Append(~mats,A);
2148  continue;
2149end if;
2150
2151new:=1;
2152index:=p^3*y5+p^2*y6+p*y3+y4;
2153
2154for a in [1..p-1] do
2155for m in [-1,1] do
2156
2157B:=H22![a,0,0,m*a];
2158C:=H22![a^3,0,0,m*a^3];
2159
2160D:=B*A*C^-1;
2161
2162z3:=Z!(D[1][1]);
2163z4:=Z!(D[1][2]);
2164z5:=Z!(D[2][1]);
2165z6:=Z!(D[2][2]);
2166
2167
2168ind1:=p^3*z5+p^2*z6+p*z3+z4;
2169
2170if ind1 lt index then new:=0; end if;
2171
2172if p mod 3 eq 2 and new eq 1 then
2173for b in [sol,-sol] do
2174
2175B:=H22![-2*w*a*b,a,m*w*a,-2*m*w*a*b];
2176C:=H22![8/3*w^2*a^3*b,4/3*w*a^3,4/3*m*w^2*a^3,8/3*m*w^2*a^3*b];
2177
2178D:=B*A*C^-1;
2179
2180z3:=Z!(D[1][1]);
2181z4:=Z!(D[1][2]);
2182z5:=Z!(D[2][1]);
2183z6:=Z!(D[2][2]);
2184
2185ind1:=p^3*z5+p^2*z6+p*z3+z4;
2186
2187if ind1 lt index then new:=0; end if;
2188
2189if new eq 0 then break; end if;
2190end for;
2191
2192end if;
2193
2194if new eq 0 then break; end if;
2195end for;
2196if new eq 0 then break; end if;
2197end for;
2198
2199if new eq 1 then
2200  Append(~mats,A);
2201  //print y3,y4,y5,y6;
2202end if;
2203
2204end for;
2205end for;
2206end for;
2207end for;
2208
2209for AS in mats do
2210//Get stabilizer of AS
2211stab1:=[];
2212stab2:=[];
2213for a in [1..p-1] do
2214for m in [-1,1] do
2215
2216B:=H22![a,0,0,m*a];
2217C:=H22![a^3,0,0,m*a^3];
2218
2219D:=B*AS*C^-1;
2220
2221if D eq AS then
2222  Append(~stab1,H33![a,0,0,0,a,0,0,0,m*a]);
2223  Append(~stab2,C);
2224end if;
2225
2226if p mod 3 eq 2 then
2227for b in [sol,-sol] do
2228
2229B:=H22![-2*w*a*b,a,m*w*a,-2*m*w*a*b];
2230C:=H22![8/3*w^2*a^3*b,4/3*w*a^3,4/3*m*w^2*a^3,8/3*m*w^2*a^3*b];
2231
2232D:=B*AS*C^-1;
2233
2234if D eq AS then
2235  Append(~stab1,H33![4*w*a*b,-3*w*a*b,a/2,0,-2*w*a*b,a,0,m*w*a,-2*m*w*a*b]);
2236  Append(~stab2,C);
2237end if;
2238
2239end for;
2240end if;
2241
2242end for;
2243end for;
2244
2245y3:=Z!(AS[1][1]); y4:=Z!(AS[1][2]); y5:=Z!(AS[2][1]); y6:=Z!(AS[2][2]);
2246
2247for y1 in [0..p-1] do
2248for y2 in [0..p-1] do
2249
2250A:=H32![y1,y2,y3,y4,y5,y6];
2251if A eq 0 then
2252  Append(~params23,[0,0,0,0,0,0]);
2253  continue;
2254end if;
2255
2256new:=1;
2257index:=p*y1+y2;
2258
2259for ii in [1..#stab1] do
2260
2261B:=stab1[ii];
2262C:=stab2[ii];
2263
2264D:=B*A*C^-1;
2265
2266z1:=Z!(D[1][1]);
2267z2:=Z!(D[1][2]);
2268
2269ind1:=p*z1+z2;
2270
2271if ind1 lt index then new:=0; end if;
2272
2273if new eq 0 then break; end if;
2274end for;
2275
2276
2277if new eq 1 then
2278  Append(~params23,[y1,y2,y3,y4,y5,y6]);
2279end if;
2280
2281end for;
2282end for;
2283
2284end for;
2285
2286if p mod 3 eq 1 then
2287  expect:=p^5+p^4+p^3+p^2+p+2+(p^2+p+1)*Gcd(p-1,4)/2;
2288else;
2289  expect:=p^5/3+p^4/3+p^3/3+p^2/3+p+2+(p^2+p+1)*Gcd(p-1,4)/2;
2290end if;
2291print #params23,expect;
2292gtotal:=gtotal+#params23;
2293
2294/*
2295Descendants of 5.14.
2296
2297Case 24, cb=baa=0, caa=xbab+bac, cac=wbab
2298where x is not a value of y(y^2+3w)/(3y^2+w)
2299
2300p=2 mod 3 only
2301*/
2302
2303if p mod 3 eq 2 then
2304
2305val:={};
2306for x in [0..p-1] do
2307  if F!(3*x^2+w) ne 0 then
2308    a:=F!(x*(x^2+3*w))*(F!(3*x^2+w))^-1;
2309    Include(~val,Z!a);
2310  end if;
2311end for;
2312for a in [1..p-1] do
2313  if a notin val then s:=a; break; end if;
2314end for;
2315//This is the value of s we need for the presentation
2316
2317//look for solution of wb^2=-3
2318for y in [1..p-1] do
2319  if F!(w*y^2+3) eq 0 then b:=y; break; end if;
2320end for;
2321
2322B1:=H22![2,2*b,2*w*b,2];
2323C1:=H22![32,-32*b,-32*w*b,32];
2324B2:=H22![2,-2*b,-2*w*b,2];
2325C2:=H22![32,32*b,32*w*b,32];
2326BB1:=H33![-4,s*b+3,3*s*w^-1+b,0,2,2*b,0,2*w*b,2];
2327BB2:=H33![-4,-s*b+3,3*s*w^-1-b,0,2,-2*b,0,-2*w*b,2];
2328
2329mats:=[];
2330//get pc,pb
2331
2332for y5 in [0,1,lns] do
2333for y6 in [0..p-1] do
2334for y3 in [0..p-1] do
2335for y4 in [0..p-1] do
2336
2337A:=H22![y3,y4,y5,y6];
2338if A eq 0 then
2339  Append(~mats,A);
2340  continue;
2341end if;
2342
2343new:=true;
2344index:=p^3*y5+p^2*y6+p*y3+y4;
2345
2346A1:=B1*A*C1^-1;
2347A2:=B2*A*C2^-1;
2348
2349for a in [1..(p-1) div 2] do
2350
2351D:=(F!a)^-2*A;
2352
2353z3:=Z!(D[1][1]);
2354z4:=Z!(D[1][2]);
2355z5:=Z!(D[2][1]);
2356z6:=Z!(D[2][2]);
2357
2358
2359ind1:=p^3*z5+p^2*z6+p*z3+z4;
2360
2361if ind1 lt index then new:=false; break; end if;
2362
2363D:=(F!a)^-2*A1;
2364
2365z3:=Z!(D[1][1]);
2366z4:=Z!(D[1][2]);
2367z5:=Z!(D[2][1]);
2368z6:=Z!(D[2][2]);
2369
2370ind1:=p^3*z5+p^2*z6+p*z3+z4;
2371
2372if ind1 lt index then new:=false; break; end if;
2373
2374D:=(F!a)^-2*A2;
2375
2376z3:=Z!(D[1][1]);
2377z4:=Z!(D[1][2]);
2378z5:=Z!(D[2][1]);
2379z6:=Z!(D[2][2]);
2380
2381ind1:=p^3*z5+p^2*z6+p*z3+z4;
2382
2383if ind1 lt index then new:=false; break; end if;
2384
2385end for;
2386
2387if new then
2388  Append(~mats,A);
2389  //print y3,y4,y5,y6;
2390end if;
2391
2392end for;
2393end for;
2394end for;
2395end for;
2396
2397
2398for AS in mats do
2399//Get stabilizer of AS
2400stab1:=[];
2401stab2:=[];
2402
2403AS1:=B1*AS*C1^-1;
2404AS2:=B2*AS*C2^-1;
2405
2406for a in [1..(p-1) div 2] do
2407
2408D:=(F!a)^-2*AS;
2409
2410if D eq AS then
2411  Append(~stab1,H33![1,0,0,0,1,0,0,0,1]);
2412  Append(~stab2,H22![a^2,0,0,a^2]);
2413end if;
2414
2415D:=(F!a)^-2*AS1;
2416
2417if D eq AS then
2418  Append(~stab1,BB1);
2419  Append(~stab2,a^2*C1);
2420end if;
2421
2422D:=(F!a)^-2*AS2;
2423
2424if D eq AS then
2425  Append(~stab1,H33!BB2);
2426  Append(~stab2,a^2*C2);
2427end if;
2428
2429end for;
2430
2431y3:=Z!(AS[1][1]); y4:=Z!(AS[1][2]); y5:=Z!(AS[2][1]); y6:=Z!(AS[2][2]);
2432
2433for y1 in [0..p-1] do
2434for y2 in [0..p-1] do
2435
2436A:=H32![y1,y2,y3,y4,y5,y6];
2437if A eq 0 then
2438  Append(~params24,[0,0,0,0,0,0]);
2439  continue;
2440end if;
2441
2442new:=1;
2443index:=p*y1+y2;
2444
2445for ii in [1..#stab1] do
2446
2447B:=stab1[ii];
2448C:=stab2[ii];
2449
2450D:=B*A*C^-1;
2451
2452z1:=Z!(D[1][1]);
2453z2:=Z!(D[1][2]);
2454
2455ind1:=p*z1+z2;
2456
2457if ind1 lt index then new:=0; end if;
2458
2459if new eq 0 then break; end if;
2460end for;
2461
2462
2463if new eq 1 then
2464  Append(~params24,[y1,y2,y3,y4,y5,y6]);
2465end if;
2466
2467end for;
2468end for;
2469
2470end for;
2471
2472print #params24,2*(p^5+p^4+p^3+p^2)/3+2*p+3;
2473gtotal:=gtotal+#params24;
2474
2475end if;
2476
2477print "Algebra 5.14 has",gtotal,"descendants of order p^7 and class 3";
2478print "2p^5+7p^4+19p^3+49p^2+128p+256+(p^2+7p+29)gcd(p-1,3)+
2479(p^2+7p+24)gcd(p-1,4)+(p+3)gcd(p-1,5) =",
24802*p^5+7*p^4+19*p^3+49*p^2+128*p+256+(p^2+7*p+29)*Gcd(p-1,3)+
2481(p^2+7*p+24)*Gcd(p-1,4)+(p+3)*Gcd(p-1,5);
2482
2483
2484print Cputime(tt);
2485
2486