
(* Quasifields.m *)

(* Vedran Krcadinac, April, 1996. *)


(* Initial setup *)

<<FiniteFieldsV2.m;
Unprotect[NonCommutativeMultiply,v,var];
ClearAttributes[NonCommutativeMultiply,Flat];
Clear[var,v];
Protect[var,v];

(********************************************

           Hall quasifields                    
 
 ********************************************)

(* Irreducible quadratic polynomials over finite fields *)

ReducibleQ[el_,i_,j_] :=
  Module[ {m, reducible,k},
          m = Length[el];
          If[ i > m || j > m,True,
              reducible = False;
              k = 1;
              While[ k<= m && Not[reducible] ,
 reducible = SameQ[el[[k]]^2+el[[i]]*el[[k]]+el[[j]], el[[1]] ];
 ++k 
                   ];
              reducible
            ]  
        ] 

IrreducibleQuadratic[{p_,n_}] := 
  IrreducibleQuadratic[{p,n}] =
  Module[ {el,m,i,j},
          el = AllElements[v,{p,n}];
          m = Length[el];
          i = 1;
          j = 1;
          While[ i<=m && ReducibleQ[el,i,j] ,
                 j = 1;
                 While[j<=m && ReducibleQ[el,i,j], ++j];
                 ++i;
               ];
          --i;
 If[ReducibleQ[el,i,j],
    Message[IrreducibleQuadratic::nonexsistent,p^n]; $Failed,
    var^2 + el[[i]]*var+el[[j]] ]
   ]               
  
IrreducibleQuadratic::nonexsistent =
 "Irreducible quadratic over the field with `1` elements not found."

IrreducibleQuadratic[x_,{p_,n_}] :=
  IrreducibleQuadratic[{p,n}] /. var -> x

AllIrreducibleQuadratics[{p_,n_}] :=
  AllIrreducibleQuadratics[{p,n}] = 
  Module[ {el,m,i,j,t},
          el = AllElements[v,{p,n}];
          m = Length[el];
          t = Table[ If[ReducibleQ[el,i,j],Null,var^2+el[[i]]*var+el[[j]] ],
                       {i,1,m},{j,1,m} ]  // Flatten;
          DeleteCases[t,Null]             
        ]  

(* Definition of addition and multiplication *)

Hall /: Hall[a_List] + Hall[b_List] := Hall[ a + b ]

a[{p_,n_}] :=
 -Coefficient[IrreducibleQuadratic[{p,n}],var,1]

b[{p_,n_}] :=
 -Coefficient[IrreducibleQuadratic[{p,n}],var,0]

Hall /: Hall[{x_FiniteField , FiniteField[0,v_,{p_,n_}] }] ** 
        Hall[{z_FiniteField , t_FiniteField }]  :=  Hall[{x*z,x*t}]

Hall /: Hall[{x_FiniteField , FiniteField[y_,v_,{p_,n_}] }] **
        Hall[{z_FiniteField , t_FiniteField}]   :=
        Hall[{x*z-(FiniteField[y,v,{p,n}]^(-1))*t*
	IrreducibleQuadratic[x,{p,n}] ,
        FiniteField[y,v,{p,n}]*z - x*t + a[{p,n}]*t }]

(* Find all elements in a Hall quasifield *)

HallAll[{p_,n_}] := HallAll[{p,n}] =
 Map[Hall,Flatten[Table[{AllElements[v,{p,n}] [[j]],
     AllElements[v,{p,n}] [[i]]},{i,1,p^n},{j,1,p^n}],1]]

(* Define addition and multiplication in a Hall quasifiled *)

HallOperations[{p_,n_}] :=
  Module[ {ff,i,j},
	   Clear[FFPlus,FFTimes];
	   Attributes[FFPlus]={Flat,Orderless,Listable};
	   Attributes[FFTimes]={Listable};
	   FFOrder = p^(2*n);
           ff = HallAll[{p,n}];
           For[i=1, i<=p^(2*n), ++i,
	       For[j=1, j<=p^(2*n), ++j,
	           FFPlus[i-1,j-1] = Position[ff,ff[[i]]+ff[[j]]] [[1,1]]-1;
		   FFTimes[i-1,j-1] = Position[ff,ff[[i]]**ff[[j]]] [[1,1]]-1
                  ]
              ]
	]

(* Save a Hall quasifileld to a file *)

SaveHall[{p_,n_},name_] :=
 ( HallOperations[{p,n}]; Put[Definition[FFOrder,FFPlus,FFTimes],name] )

(*********************************************

          Andre quasifields

 *********************************************)


(* Definition of Andre objects. They should behave as FiniteField
   objects with respect to the operations * and +                    *)

Andre[f_, {p_, s_, d_}] :=
          Andre[PolynomialMod[PolynomialRemainder[PolynomialMod[f, p],
                   Generator[v, p, s*d], v], p], {p, s, d}] /;
     Exponent[f, v] >= s*d || Select[CoefficientList[f, v], # >= p&] != {}

Andre /: Andre[f_,  {p_, s_, d_}] + Andre[g_,{p_, s_, d_}] :=
   Andre[PolynomialMod[f + g, p], {p, s, d }]

Andre /: Andre[f_,{p_, s_, d_}] + m_Integer :=
   Andre[PolynomialMod[f + m,p],{p,s,d}]

Andre /: Andre[f_, {p_,s_,d_}] * Andre[g_,{p_, s_,d_}] :=
   Andre[PolynomialMod[f * g, p] , {p, s,d}]

Andre /: Andre[f_, {p_, s_,d_}] * m_Integer :=
   Andre[PolynomialMod[f * m, p], {p, s,d}]

Andre /: Andre[0, {p_, s_,d_}] * x_ :=
   Andre[0, {p, s,d}]


Andre /: Power[Andre[f_, {p_, s_,d_}], -1] :=
   Andre[Last[FFExtendedGCD[Generator[v,p,s*d],f,v,p]],{p, s,d}]  /;
   PolynomialGCD[Generator[v, p, s*d], f, Modulus -> p] == 1

Andre /: Andre[f_, {p_, s_,d_}]^k_ :=
   Andre[PolynomialMod[f^k, p], {p, s,d}]    /;  k >= 0


(* All Andre objects for given p,s,d   *)

AndreAll[{p_,s_,d_}] := AndreAll[{p,s,d}] =
   Array[ Andre[Expand[Poly[v,##]],{p,s,d}] &,
	  Table[p,{s*d}] ]  // Flatten


(* The primitive root in the Andre finite field *)

PrimitiveQ[Andre[f_,{p_,s_,d_}]] :=
 Length[Union[Table[Andre[f,{p,s,d}]^i,{i,0,p^(s*d)-2}]]] == p^(s*d)-1

PrimitiveRoot[{p_,s_,d_}] := Andre[v,{p,s,d}] 

CyclicOrder[{p_,s_,d_}] := CyclicOrder[{p,s,d}] =
 Table[Andre[v,{p,s,d}]^i,{i,0,p^(s*d)-2}]


(* The norm mapping *)

Norm[ Andre[f_,{p_,s_,d_}] ] :=
 Mod[ Position[CyclicOrder[{p,s,d}] , Andre[f,{p,s,d}]] [[1,1]] -1,
      p^s-1]

Norm2[ Andre[f_,{p_,s_,d_}] ] :=
  Andre[f,{p,s,d}] ^ Sum[p^(i*s),{i,0,d-1}]

Norm3[  Andre[f_,{p_,s_,d_}] ] :=
  Position[CyclicOrder[{p,s,d}],Norm2[Andre[f,{p,s,d}]]][[1,1]]-1


(* Find all possible functions f from {0,1,..,p^s-2} to {0,1,..,d-1}
   such that f[0] = 0. There are altogether d^(p^s-2) of them.
   One should choose one and assign it to Phi[{p,s,d}].              *)

MakeFunction = Function[i, #1[[i + 1]]] & 

AllPhis1[{2,1,d_}] = {{0}}

AllPhis1[{p_, s_, d_}] :=
  Flatten[Array[{1,##1} - 1 & , Table[d, {p^s - 2}]],p^s-3]

AllPhis[{p_,s_,d_}] :=
 Map[MakeFunction,AllPhis1[{p,s,d}]]

(* If it is not already defined, randomly choose one Phi[{p,s,d}] *)

Phi[{p_,s_,d_}] := Phi[{p,s,d}] =
 AllPhis[{p,s,d}] [[ Random[Integer,{1,Length[AllPhis[{p,s,d}]]}] ]]

(* An isomorphism between the group {0,1,..,d-1} with addition (mod d)
   and the authomorphism group of the field with p^(s*d) elements
   that fix it's subfield of order p^s.                            *)

Aut[i_Integer, {p_,s_,d_}] :=
  Function[x,x^(p^(s*i))]


(* The mapping Alpha, from the multiplicative group of the field
   of order p^(s*d) to the above described automorphism group.
   The Andre quasifield will be associative if and only if this
   mapping is a homomorphism, i.e. if the chosen Phi is. Norm and  
   Aut are always homomorphisms.                               *)

Alpha[ Andre[f_,{p_,s_,d_}] ] := 
  Aut[ Phi[{p,s,d}] [ Norm[Andre[f,{p,s,d}]] ], {p,s,d}]

(* Finally, the definition of multiplication in an Andre quasifield. *)

Andre /: Andre[0,{p_,s_,d_}] ** Andre[y_,{p_,s_,d_}] :=
 Andre[0,{p,s,d}]

Andre /: Andre[x_,{p_,s_,d_}] ** Andre[y_,{p_,s_,d_}] :=
 Andre[x,{p,s,d}] * Alpha[Andre[x,{p,s,d}]][Andre[y,{p,s,d}]]


(* Define operations in an Andre quasifield *)

AndreOperations[{p_,s_,d_}] :=
  Module[ {ff,i,j},
	   Clear[FFPlus,FFTimes];
	   Attributes[FFPlus]={Flat,Orderless,Listable};
	   Attributes[FFTimes]={Listable};
	   FFOrder = p^(s*d);
           ff = AndreAll[{p,s,d}];
           For[i=1, i<=p^(s*d), ++i,
	       For[j=1, j<=p^(s*d), ++j,
	           FFPlus[i-1,j-1] = Position[ff,ff[[i]]+ff[[j]]] [[1,1]]-1;
		   FFTimes[i-1,j-1] = Position[ff,ff[[i]]**ff[[j]]] [[1,1]]-1
                  ]
              ]
	]


(* Save an Andre quasifield to a file *)

SaveAndre[{p_,s_,d_},name_] :=
 ( AndreOperations[{p,s,d}]; Put[Definition[FFOrder,FFPlus,FFTimes],name] )



(*********************************************

          Hughes nearfields

 *********************************************)

(*
    Hughes nearfields are in fact special Andre quasifields.
    They are nearfields of order q^2, where q=p^s is an odd prime power.
    Phi is chosen to be the homomorphism
 
        Phi[i] = i (mod 2) , i=0,1,..,p^s-2

    The kernel F of a Hughes nearfield N is a finite field of order q.
    Elements of the kernel commute with all elements in N.

    We need Hughes nearfields for the construction of Hughes planes.
    What we also need is a 3*3 matrix A over F such that A^(q^2+q+1)
    is scalar and no smaller power of A has this property.
    
 *)


(* The subfield of order q=p^s in the field of order q^d = p^(s*d).
   This will be the kernel F.                                         *)

Subfield[{p_,1,d_}] := Subfield[{p,1,d}] = Take[AndreAll[{p,1,d}],{1,p}]

Subfield[{p_,s_,d_}] := Subfield[{p,s,d}] =
 Prepend[ Table[ Evaluate[Andre[v,{p,s,d}]^(j*Sum[p^(s*i),{i,0,d-01}])],
		{j,0,p^s-2}],
         Andre[0,{p,s,d}] ]


(* Ordering of the Andre finite field in which its subfield comes first *)

SubfieldFirst[{p_,1,d_}] := SubfieldFirst[{p,1,d}] = AndreAll[{p,1,d}]

SubfieldFirst[{p_,s_,d_}] := SubfiledFirst[{p,s,d}] =
 Join[ Subfield[{p,s,d}],
       DeleteCases[AndreAll[{p,s,d}] , 
		   x_Andre /; MemberQ[ Subfield[{p,s,d}],x] 
		  ]
     ]

 
(* To obtain the matrix A we need irreducible cubic polynomials over
    the kernel F                                                        *)

ReducibleQ[el_,i_,j_,k_] :=
  Module[ {m, reducible,l},
          m = Length[el];
          If[ i > m || j > m || k >m,True,
              reducible = False;
              l = 1;
              While[ l<= m && Not[reducible] ,
 reducible = SameQ[el[[l]]^3+el[[i]]*(el[[l]]^2)+el[[j]]*el[[l]]+el[[k]],
                   el[[1]] ];
                     ++l 
                   ];
              reducible
            ]  
        ] 

AllIrreducibleCubics[{p_,s_,d_}] :=
  AllIrreducibleCubics[{p,s,d}] = 
  Module[ {el,m,i,j,k,t},
          el = Subfield[{p,s,d}];
          m = Length[el];
          t = Table[ If[ReducibleQ[el,i,j,k],Null,
                        var^3+el[[i]]*var^2+el[[j]]*var+el[[k]]
                       ],
                     {i,1,m},{j,1,m},{k,1,m} 
                   ]  // Flatten;
          DeleteCases[t,Null]             
        ]  

 
(* A function that checks whether a 3*3 matrix is scalar *)

ScalarQ[{{s_,a_,b_},{c_,s_,d_},{e_,f_,s_}}] := True /; 
  Apply[ And , Map[ # == 0 || MatchQ[#,Andre[0,{_,_,_}]] & , 
                    {a,b,c,d,e,f} ] ]

ScalarQ[x_] := False


(* This function returns the least positive integer i such that M^i is scalar *)

MatrixOrder[M_] :=
  Module[ {P,i},
          i=1;  P=M;
          While[Not[ScalarQ[P]], ++i; P = P.M];
          i
        ]


(* Definition of the matrix A *)

MatrixA[pol_] := 
  Transpose[{Reverse[Take[CoefficientList[-pol,var], {1, 3}]],
             {1, 0, 0}, {0, 1, 0}}]


MatrixA[{p_,s_,d_}] :=  MatrixA[{p,s,d}] =
  Module[{i,A},
         i = 1;
         A = MatrixA[AllIrreducibleCubics[{p,s,d}][[1]]];
         While[ MatrixOrder[A] != (p^(2*s)+p^s+1) , ++i;
                A =  MatrixA[AllIrreducibleCubics[{p,s,d}][[i]]]
              ];
         A
        ]


(* Define operatons in a Hughes nearfield and save it to a file,
   together with A and all its relevant powers. In order to construct
   a Hughes plane we will have to take left multiplicative inverses,
   so they are computed and saved, too.                                      *)

DefineA[{p_,s_,d_}] :=
 ( Clear[A]; Do[A[i] = MatrixPower[MatrixA[{p,s,d}],i] /.
 x_ /; MemberQ[Subfield[{p,s,d}],x] :> Position[Subfield[{p,s,d}],x][[1,1]]-1,
                { i, 0, p^(2*s)+p^s}] )

HughesOperations[{p_,s_}] :=
  Module[ {ff,i,j},
           Phi[{p,s,2}] = Function[i,Mod[i,2]];
	   Clear[FFPlus,FFTimes];
	   Attributes[FFPlus]={Flat,Orderless,Listable};
	   Attributes[FFTimes]={Listable};
	   FFOrder = p^(s*2);
           ff = SubfieldFirst[{p,s,2}];
           For[i=1, i<=p^(s*2), ++i,
	       For[j=1, j<=p^(s*2), ++j,
	           FFPlus[i-1,j-1] = Position[ff,ff[[i]]+ff[[j]]] [[1,1]]-1;
		   FFTimes[i-1,j-1] = Position[ff,ff[[i]]**ff[[j]]] [[1,1]]-1
                  ]
              ]
	]

DefineInverses[{p_,s_}] :=
  Module[ {el,i},
          Clear[FFInverse];
          el = Table[i,{i,1,p^(2*s)-1}];
          Do[ FFInverse[i] = Position[Map[FFTimes[#,i]&,el],1][[1,1]],
              {i,1,p^(2*s)-1} ]
        ]

SaveHughes[{p_,s_},name_] :=
 ( DefineA[{p,s,2}]; HughesOperations[{p,s}]; DefineInverses[{p,s}]; 
   Put[Definition[FFOrder,FFPlus,FFTimes,FFInverse,A],name] )

