Linear Algebra over Finite Fields and Generic Domains in Maple Simon Lo, Michael Monagan Department of Mathematics, Sim

1. Linear Algebra over Finite Fields and Generic Domains in Maple Simon Lo, Michael Monagan Department of Mathematics, Simon Fraser University Examples in Maple Introduction Most linear algebra algorithms are general algorithms that work for matrices with coefficients over any field. Therefore, it is useful to implement the algorithms in Maple such that they are independent of any particular fields and field representations. For some algorithms, it is possible to generalize them further to work for matrices over any arbitrary domains such as Euclidean domain, integral domains etc. that are more general. Here, we restrict our attention to finite fields because we want to avoid dealing with problems such as representation and numerical stability. We represent a domain in Maple as a table of procedures and constants. Special procedures like ConvertIn and ConvertOut are used to convert from Maple representation to our representation and vice versa. Rand is a procedure that generates a random element and is very useful for generating test cases. The following is a list of procedures/constants that must be additionally defined in the table for each type of domains. Integral domain (ID): 0, 1, +, , , UnitPart, NormalPart, Divide (exact division) Unique factorization domain (UFD): Gcd Euclidean domain (ED): EuclideanNorm, Quo, Rem, Gcdex Finite field (GF): Inv,  (division) Three types of base domains have been implemented. They are Z, Zp, and GF(pn) where p is a prime. We have also included procedures for constructing univariate polynomial domains over a general coefficient domain that has already been constructed, and procedures for constructing square matrix domains over a general coefficient domain that has already been constructed. Multiplication in a given square matrix domain in general is not commutative, so we have defined operations for a ring with identity (RI), and a commutative ring with identity (CRI). We will define operations for a general ring in the future. If the coefficient domain is a GF, then the polynomial is a ED. If the coefficient domain is a ED, then the polynomial is a UFD. If the coefficient domain is any other domains, say type D, then the polynomial domain D[x] is also of type D. A matrix domain is in general is of type RI, given that the coefficients are also of type RI. Operations for polynomial domains (Poly): CoefficientDomain, ConvertInD, Variable, Coeff, Degree, Diff, Eval, Lcoeff, Ldegree, Resultant, Sqrfree, Tcoeff. Currently, Factors, Splits, and Roots are only implemented for polynomial domains over finite fields. Operations for matrix domains (Matrix): CoefficientDomain, ConvertInD, Trace, Adjoint, Inverse, Charpoly, Ffgausselim, Det, Rank, Rowspan, Colspan, Hermite, Smith, Gausselim, Gaussjord, Rowspace, Colspace, Nullspace, Backsub, Forwardsub, Hessenberg, Linsolve, LUdecomp, Eigenvalues, Eigenvectors. > alias(a=RootOf(x^3+x^2+1)): > F := `makeGF/Fq`(a,2): > M := makeMatrix(F,4): > A := M['ConvertIn']([[a^2+a, a+1, 0, 0], [a^2+a, a^2+a, a+1, 1], [a^2+a, 0, a^2+1, a+1], [a+1, 0, a, a^2+1]]); > M['Charpoly'](A); > M['Ffgausselim'](A); > M['Gausselim'](A); > F := `makeGF/Zp`(2): > M := makeMatrix(F,4): > A := M['ConvertIn']([[1,1,1,0],[1,1,1,1], [1,1,1,1],[0,0,1,1]]); > M['Eigenvalues']('F2',A); > F2['Minpoly']; > M['Eigenvectors']('F2',A); > F := `makeGF/Zp`(2):> P := makePoly(F,x): > M := makeMatrix(P,4): > A := M['Rand'](); > M['Adjoint'](A); > M['Det'](A); > M['Ffgausselim'](A); > M['Hermite'](A); > M['Smith'](A); Implementation of Generic Algorithms In order to define generic algorithms for linear algebra suitable for matrices with arbitrary coefficient domains, it is necessary for the generic algorithm to take the table of operations/constants for the domain as an input parameter and do all operations using the appropriate operation from the table. This approach makes sure that the generic code will work for any domain with the appropriate operations defined. Generally, computation using this approach will be slower than the equivalent Maple code because of the increased number of table access and procedure calls. The procedure calls will generally execute slower than the equivalent builtin Maple functions. Efficient data structures such as modp1 and modp2 have been used to speed up computations for finite fields and polynomials over finite fields. Consider the example below of the simplified code for fraction-free Gaussian elimination (D is the coefficient domain, B is the input matrix). Ffgausselim := proc(D,B) A := Matrix(B); m,n := LinearAlgebra:-Dimensions(A); divisor := D[1]; r := 1; for c to n while r <= m do # find pivot for i from r to m while A[i,c] = D[0] do od; if i > n then next fi; if i <> r then # interchange row i with row r to make A[r,c] the pivot for j from c to m do t := A[i,j]; A[i,j] := A[r,j]; A[r,j] := t od fi; for i from r+1 to m do for j from c+1 to n do A[i,j] := D[`-`]( D[`*`](A[r,c],A[i,j]), D[`*`](A[i,c],A[r,j]) ); if divisor<>D[1] then D['Divide'](A[i,j],divisor,'q'); A[i,j] := q; fi; od; A[i,c] := D[0]; od; divisor := A[r,c]; r := r+1 # go to next row od; # go to next column Aend: The code for fraction-free Gaussian elimination requires only the operations 0, 1, , , Divide from the coefficient domain, so the same code would work for any coefficient domain with these operations properly defined. We have parameterized the code by D, so that the code works for any domain D, not just the integers we might have had in mind when we wrote the code. This main idea of passing as a parameter a collection of functions as a single unit came from the AXIOM system, which was formerly called Scratchpad II. The same idea is also used in the Domains package in Maple. Compare the code for the fraction-free Gaussian elimination with the code for Gaussian elimination below which generally requires the coefficient domain to be a field (since we need to compute inverses). Gausselim := proc(F,B) A := Matrix(B); m,n := LinearAlgebra:-Dimensions(A); for i to m do for j to n while A[i,j] = F[0] do od; if j > n or A[i,j] = F[1] then next fi; s := F['Inv'](A[i,j]); for k from j+1 to n do A[i,k] := F[`*`](s,A[i,k]) od; A[i,j] := F[1]; od; r := 1; for c to n while r <= m do for i from r to m while A[i,c] = F[0] do od; if i > m then next fi; if i <> r then # interchange row i with row r for j from c to n do t := A[i,j]; A[i,j] := A[r,j]; A[r,j] := t od fi; for i from r+1 to m do if A[i,c] = F[0] then next fi; A[i,c] := F[0]; for j from c+1 to n do A[i,j] := F[`-`](A[i,j],A[r,j]) od; for j from c+1 to n while A[i,j] = F[0] do od; if j > n then next fi; s := F['Inv'](A[i,j]); for k from j+1 to n do A[i,k] := F[`*`](s,A[i,k]) od; A[i,j] := F[1]; od; r := r + 1 # go to next row od; # go to next column A end: