Quelle matrix.gi Sprache: unbekannt

Spracherkennung für: .gi vermutete Sprache: Unknown {[0] [0] [0]} [Methode: Schwerpunktbildung, einfache Gewichte, sechs Dimensionen]

#############################################################################
##
## This file is part of GAP, a system for computational discrete algebra.
## This file's authors include Thomas Breuer, Frank Celler, Alexander Hulpke, Heiko Theißen, Martin Schönert.
##
## Copyright of GAP belongs to its developers, whose names are too numerous
## to list here. Please refer to the COPYRIGHT file for details.
##
## SPDX-License-Identifier: GPL-2.0-or-later
##
## This file contains methods for matrices.
##

#
# Kernel method for computing
#

InstallMethod(Zero,
 [IsRectangularTable and IsAdditiveElementWithZeroCollColl and IsInternalRep],
 ZERO_ATTR_MAT);

# Fallback methods for matrix entry access.
InstallOtherMethod( \[\], [ IsMatrix, IsPosInt, IsPosInt ],
 1-RankFilter(IsMatrix),
 function (m,i,j) return m[i][j]; end );
InstallOtherMethod( \[\]\:\=, [ IsMatrix and IsMutable, IsPosInt, IsPosInt, IsObject ],
 1-RankFilter(IsMatrix),
 function (m,i,j,o) m[i][j]:=o; end );

# Note that installing 'ShallowCopy' is not allowed here, for several reasons.
# Firstly, the result has to be independent of the input, and fully mutable.
# Secondly, a shallow copy of an 8bit matrix would be again an 8bit matrix.
# Also installing 'M -> List( M, ShallowCopy )' is not allowed
# because the entries of 'M' can be proper vector objects.
InstallOtherMethod( Unpack,
 [ IsMatrix ],
 M -> List( M, Unpack ) );

# generic unpack method for row list matrices
InstallMethod( Unpack,
 [ IsRowListMatrix ],
 M -> List( M, Unpack ) );

#############################################################################
##
#F PrintArray( <array> ) . . . . . . . . . . . . . . . . pretty print matrix
##
InstallGlobalFunction(PrintArray,function( array )
 local arr, max, l, k,maxp,compact,bl;

 compact:=ValueOption("compact")=true;
 if compact then
 maxp:=0;
 bl:="";
 else
 maxp:=1;
 bl:=" ";
 fi;
 if not IsDenseList( array ) then
 Error( "<array> must be a dense list" );
 elif Length( array ) = 0 then
 Print( "[ ]\n" );
 elif array = [[]] then
 Print( "[ [ ] ]\n" );
 elif not ForAll( array, IsList ) then
 arr := List( array, String );
 max := Maximum( List( arr, Length ) );
 Print( "[ ", String( arr[ 1 ], max + maxp ) );
 for l in [ 2 .. Length( arr ) ] do
 Print( ", ", String( arr[ l ], max + 1 ) );
 od;
 Print( " ]\n" );
 else
 arr := List( array, x -> List( x, String ) );
 if compact then
 max:=List([1..Length(arr[1])],
 x->Maximum(List([1..Length(arr)],y->Length(arr[y][x]))));
 else
 max := Maximum( List( arr,
 function(x)
 if Length(x) = 0 then
 return 1;
 else
 return Maximum( List(x,Length) );
 fi;
 end) );
 fi;

 Print( "[",bl );
 for l in [ 1 .. Length( arr ) ] do
 if l > 1 then
 Print(bl," ");
 fi;
 Print( "[",bl );
 if Length(arr[ l ]) = 0 then
 Print(bl,bl,"]" );
 else
 for k in [ 1 .. Length( arr[ l ] ) ] do
 if compact then
 Print( String( arr[ l ][ k ], max[k] + maxp ) );
 else
 Print( String( arr[ l ][ k ], max + maxp ) );
 fi;
 if k = Length( arr[ l ] ) then
 Print( bl,"]" );
 else
 Print( ", " );
 fi;
 od;
 fi;
 if l = Length( arr ) then
 Print( bl,"]\n" );
 else
 Print( ",\n" );
 fi;
 od;
 fi;
end);

##########################################################################
##
#M Display( <mat> )
##
InstallMethod( Display,
 "for a matrix",
 [IsMatrix ],
 PrintArray );

#############################################################################
##
#M IsGeneralizedCartanMatrix( <A> )
##
InstallMethod( IsGeneralizedCartanMatrix,
 "for a matrix",
 [ IsMatrixOrMatrixObj ],
 function( A )

 local n, i, j;

 if NrRows( A ) <> NrCols( A ) then
 Error( "<A> must be a square matrix" );
 fi;

 n:= NrRows( A );
 for i in [ 1 .. n ] do
 if A[i,i] <> 2 then
 return false;
 fi;
 od;
 for i in [ 1 .. n ] do
 for j in [ i+1 .. n ] do
 if not IsInt( A[i,j] ) or not IsInt( A[j,i] )
 or 0 < A[i,j] or 0 < A[j,i] then
 return false;
 elif ( A[i,j] = 0 and A[j,i] <> 0 )
 or ( A[j,i] = 0 and A[i,j] <> 0 ) then
 return false;
 fi;
 od;
 od;
 return true;
 end );

#############################################################################
##
#M IsDiagonalMatrix(<mat>)
##
InstallMethod( IsDiagonalMatrix,
 "for a matrix",
 [ IsMatrixOrMatrixObj ],
 function( mat )
 local i, j, z;
 z:=ZeroOfBaseDomain(mat);
 for i in [ 1 .. NrRows( mat ) ] do
 for j in [ 1 .. NrCols( mat ) ] do
 if mat[i,j] <> z and i <> j then
 return false;
 fi;
 od;
 od;
 return true;
 end);

InstallTrueMethod( IsDiagonalMatrix, IsMatrixOrMatrixObj and IsEmptyMatrix );

#############################################################################
##
## Note that an empty list is not a matrix,
## but 'IsDiagonalMat(rix)' used to return 'true' for it.
##
## Note also that there was no such hack for other operations such as
## 'IsUpperTriangularMat(rix)' or 'IsLowerTriangularMat(rix)'.
##
## And note that installing an implication from 'IsList and IsEmpty' to
## 'IsDiagonalMatrix' does not work because the type of an empty plain list
## cannot store the information.
## Furthermore, we cannot simply install 'ReturnTrue' as a method because
## then the method installation would complain that one should just install
## an implication.
##
InstallOtherMethod( IsDiagonalMatrix, [ IsList and IsEmpty ], x -> true );

#############################################################################
##
#M IsUpperTriangularMatrix(<mat>)
##
InstallMethod( IsUpperTriangularMatrix,
 "for a matrix",
 [ IsMatrixOrMatrixObj ],
 function( mat )
 local i, j, z;
 z:=ZeroOfBaseDomain(mat);
 for j in [ 1 .. NrCols( mat ) ] do
 for i in [ j+1 .. NrRows( mat ) ] do
 if mat[i,j] <> z then
 return false;
 fi;
 od;
 od;
 return true;
 end);

#############################################################################
##
#M IsLowerTriangularMatrix(<mat>)
##
InstallMethod( IsLowerTriangularMatrix,
 "for a matrix",
 [ IsMatrixOrMatrixObj ],
 function( mat )
 local i, j, z;
 z:=ZeroOfBaseDomain(mat);
 for i in [ 1 .. NrRows( mat ) ] do
 for j in [ i+1 .. NrCols( mat ) ] do
 if mat[i,j] <> z then
 return false;
 fi;
 od;
 od;
 return true;
 end);

#############################################################################
##
#M DiagonalOfMatrix(<mat>) . . . . . . . . . . . . . . . diagonal of matrix
##
InstallGlobalFunction( DiagonalOfMatrix, function ( mat )
 local diag, i;

 diag := [];
 i := 1;
 while i <= NrRows(mat) and i <= NrCols(mat) do
 diag[i] := mat[i,i];
 i := i + 1;
 od;
 while 1 <= NrRows(mat) and i <= NrCols(mat) do
 diag[i] := mat[1,1] - mat[1,1];
 i := i + 1;
 od;
 return diag;
end );

#############################################################################
##
#R IsNullMapMatrix . . . . . . . . . . . . . . . . . . . null map as matrix
##
DeclareRepresentation( "IsNullMapMatrix", IsMatrix and IsPositionalObjectRep, [ ] );

BindGlobal( "NullMapMatrix",
 Objectify( NewType( ListsFamily, IsNullMapMatrix ), [ ] ) );

InstallOtherMethod( Length,
 "for null map matrix",
 [ IsNullMapMatrix ],
 null -> 0 );

InstallMethod( IsZero,
 "for null map matrix",
 [ IsNullMapMatrix ],
 x -> true);

InstallMethod( ZeroSameMutability,
 "for null map matrix",
 [ IsNullMapMatrix ],
 null -> null );

InstallMethod( \+,
 "for two null map matrices",
 [ IsNullMapMatrix, IsNullMapMatrix ],
 ReturnFirst );

InstallMethod( AdditiveInverseSameMutability,
 "for a null map matrix",
 [ IsNullMapMatrix ],
 null -> null );

InstallMethod( AdditiveInverseOp,
 "for a null map matrix",
 [ IsNullMapMatrix ],
 null -> null );

InstallMethod( \*,
 "for two null map matrices",
 [ IsNullMapMatrix, IsNullMapMatrix ],
 ReturnFirst );

InstallMethod( \*,
 "for a scalar and a null map matrix",
 [ IsScalar, IsNullMapMatrix ],
 function(s,null)
 return null;
end );

InstallMethod( \*,
 "for a null map matrix and a scalar",
 [ IsNullMapMatrix, IsScalar ],
 ReturnFirst );

InstallMethod( \*,
 "for vector and null map matrix",
 [ IsVector, IsNullMapMatrix ],
 function( v, null )
 return [ ];
end );

InstallOtherMethod( \*,
 "for empty list and null map matrix",
 [ IsList and IsEmpty, IsNullMapMatrix ],
 function( v, null )
 return [ ];
end );

InstallMethod( \*,
 "for matrix and null map matrix",
 [ IsMatrix, IsNullMapMatrix ],
 function( A, null )
 return List( A, row -> [ ] );
end );

InstallOtherMethod( \*,
 "for list and null map matrix",
 [ IsList, IsNullMapMatrix ],
 function( A, null )
 return List( A, row -> [ ] );
end );

InstallMethod( ViewObj,
 "for null map matrix",
 [ IsNullMapMatrix ],
 function( null )
 Print( "<null map>" );
end );

InstallMethod( PrintObj,
 "for null map matrix",
 [ IsNullMapMatrix ],
 function( null )
 Print( "NullMapMatrix" );
end );

#############################################################################
##
#F Matrix_OrderPolynomialInner( <fld>, <mat>, <vec>, <spannedspace> )
##
## Returns the coefficients of the order polynomial of <mat> at <vec>
## modulo <spannedspace>. No conversions are attempted on <mat> or
## <vec>, which should usually be immutable and compressed for best
## results. <spannedspace> should be a semi-echelonized basis, stored
## as a list with holes with the vector with pivot in position 
## Vectors are added to <spannedspace> so that it also spans all the images
## of <vec> under the algebra generated by <mat>
##
## The result, and any vectors added to <spannedspace> are compressed
## and immutable
##
#N In characteristic zero, or for structured sparse matrices, the naive
#N Gaussian elimination here may not be optimal
##
#N Shift to using ClearRow once we have kernel methods that give a
#N performance benefit
##
BindGlobal( "Matrix_OrderPolynomialInner", function( fld, mat, vec, vecs)
 local d, w, p, one, zero, zeroes, piv, pols, x;
 Info(InfoMatrix,2,"Order Polynomial Inner on ",NrRows(mat),
 " x ",NrCols(mat)," matrix over ",fld," with ",
 Number(vecs)," basis vectors already given");
 d := Length(vec);
 pols := [];
 one := One(fld);
 zero := Zero(fld);
 zeroes := [];

 # this loop runs images of <vec> under powers of <mat>
 # trying to reduce them with smaller powers (and tracking the polynomial)
 # or with vectors from <spannedspace> as passed in
 # when we succeed, we know the order polynomial

 repeat
 w := ShallowCopy(vec);
 p := ShallowCopy(zeroes);
 Add(p,one);
 ConvertToVectorRepNC(p,fld);
 piv := PositionNonZero(w,0);

 #
 # Strip as far as we can
 #

 while piv <= d and IsBound(vecs[piv]) do
 x := -w[piv];
 if IsBound(pols[piv]) then
 AddCoeffs(p, pols[piv], x);
 fi;
 AddRowVector(w, vecs[piv], x, piv, d);
 piv := PositionNonZero(w,piv);
 od;

 #
 # if something is left then we don't have the order poly yet
 # update tables etc.
 #

 if piv <=d then
 x := Inverse(w[piv]);
 MultVector(p, x);
 MakeImmutable(p);
 pols[piv] := p;
 MultVector(w, x );
 MakeImmutable(w);
 vecs[piv] := w;
 vec := vec*mat;
 Add(zeroes,zero);
 fi;
 until piv > d;
 MakeImmutable(p);
 Info(InfoMatrix,2,"Order Polynomial returns ",p);
 return p;
end );

#############################################################################
##
#F Matrix_OrderPolynomialSameField( <fld>, <mat>, <vec>, <ind> )
##
## Compute the order polynomial, all the work is done in the
## routine above
##
BindGlobal( "Matrix_OrderPolynomialSameField", function( fld, mat, vec, ind )
 local imat, ivec, coeffs;
 imat:=ImmutableMatrix(fld,mat);
 ivec := Immutable(vec);
 ConvertToVectorRepNC(ivec, fld);
 coeffs := Matrix_OrderPolynomialInner( fld, imat, ivec, []);
 return UnivariatePolynomialByCoefficients(ElementsFamily(FamilyObj(fld)), coeffs, ind );
end );

#############################################################################
##
#F Matrix_CharacteristicPolynomialSameField( <fld>, <mat>, <ind> )
##
BindGlobal( "Matrix_CharacteristicPolynomialSameField",
 function( fld, mat, ind)
 local i, n, base, imat, vec, one,cp,op,zero,fam;
 Info(InfoMatrix,1,"Characteristic Polynomial called on ",
 NrRows(mat)," x ",NrCols(mat)," matrix over ",fld);
 imat := ImmutableMatrix(fld,mat);
 n := NrRows(mat);
 base := [];
 vec := ZeroOp(mat[1]);
 one := One(fld);
 zero := Zero(fld);
 fam := ElementsFamily(FamilyObj(fld));
 cp:=[one];
 if Is8BitMatrixRep(mat) and NrRows(mat)>0 then
 # stay in the same field as matrix
 ConvertToVectorRepNC(cp,Q_VEC8BIT(mat[1]));
 fi;
 cp := UnivariatePolynomialByCoefficients(fam,cp,ind);
 for i in [1..n] do
 if not IsBound(base[i]) then
 vec[i] := one;
 op := Matrix_OrderPolynomialInner( fld, imat, vec, base);
 cp := cp * UnivariatePolynomialByCoefficients( fam,op,ind);
 vec[i] := zero;
 fi;
 od;
 Assert(2, Length(CoefficientsOfUnivariatePolynomial(cp)) = n+1);
 if AssertionLevel()>=3 then
 # cannot use Value(cp,imat), as this uses characteristic polynomial
 n:=Zero(imat);
 one:=One(imat);
 for i in Reversed(CoefficientsOfUnivariatePolynomial(cp)) do
 n:=n*imat+(i*one);
 od;
 Assert(3,IsZero(n));
 fi;
 Info(InfoMatrix,1,"Characteristic Polynomial returns ", cp);
 return cp;
end );

##########################################################################
##
#F Matrix_MinimalPolynomialSameField( <fld>, <mat>, <ind> )
##
BindGlobal( "Matrix_MinimalPolynomialSameField", function( fld, mat, ind )
 local i, n, base, imat, vec, one, zero, fam,
 mp, dim, span,op,w, piv,j;

 Info(InfoMatrix,1,"Minimal Polynomial called on ",
 NrRows(mat)," x ",NrCols(mat)," matrix over ",fld);
 imat := ImmutableMatrix(fld,mat);
 n := NrRows(imat);
 base := [];
 dim := 0; # should be number of bound positions in base
 one := One(fld);
 zero := Zero(fld);
 fam := ElementsFamily(FamilyObj(fld));
 mp:=[one];
 if Is8BitMatrixRep(mat) and NrRows(mat)>0 then
 # stay in the same field as matrix
 ConvertToVectorRepNC(mp,Q_VEC8BIT(mat[1]));
 fi;
 #keep coeffs
 #mp := UnivariatePolynomialByCoefficients( fam, mp,ind);
 while dim < n do
 vec := ShallowCopy(mat[1]);
 for i in [1..n] do
 #Add(vec,Random([one,zero]));
 vec[i]:=Random([one,zero]);
 od;
 vec[Random(1,n)] := one; # make sure it's not zero
 #ConvertToVectorRepNC(vec,fld);
 MakeImmutable(vec);
 span := [];
 op := Matrix_OrderPolynomialInner( fld, imat, vec, span);
 #op := UnivariatePolynomialByCoefficients(fam, op, ind);
 #mp := Lcm(mp, op);
 # this command takes much time since a polynomial ring is created.
 # Instead use the quick gcd-based method (avoiding the dispatcher):
 #mp := (mp*op)/GcdOp(mp, op);
 #mp:=mp/LeadingCoefficient(mp);
 mp:=QUOTREM_LAURPOLS_LISTS(ProductCoeffs(mp,op),GcdCoeffs(mp,op))[1];
 mp:=mp/Last(mp);

 for j in [1..Length(span)] do
 if IsBound(span[j]) then
 if dim < n then
 if not IsBound(base[j]) then
 base[j] := span[j];
 dim := dim+1;
 else
 w := ShallowCopy(span[j]);
 piv := j;
 repeat
 AddRowVector(w,base[piv],-w[piv], piv, n);
 piv := PositionNonZero(w, piv);
 until piv > n or not IsBound(base[piv]);
 if piv <= n then
 MultVector(w,Inverse(w[piv]));
 MakeImmutable(w);
 base[piv] := w;
 dim := dim+1;
 fi;
 fi;
 fi;
 fi;
 od;
 od;
 mp := UnivariatePolynomialByCoefficients( fam, mp,ind);
 Assert(3, IsZero(Value(mp,imat)));
 Info(InfoMatrix,1,"Minimal Polynomial returns ", mp);
 return mp;
end );

##########################################################################
##
#M Display( <ffe-mat> )
##
InstallMethod( Display,
 "for matrix of FFEs",
 [ IsFFECollColl and IsMatrix ],
function( m )
 local deg, chr, zero, w, t, x, v, f, z, y;

 if NrCols(m) = 0 then
 TryNextMethod();
 fi;
 if IsZmodnZObj(m[1,1]) then
 # this case mostly applies to large characteristic, in
 # which the regular code for FFE elements does not work
 # (e.g. it tries to create the range [2..chr], which
 # means chr may be at most 2^28 resp. 2^60).
 Print("ZmodnZ matrix:\n");
 t:=List(m,i->List(i,i->i![1]));
 Display(t);
 Print("modulo ",Characteristic(m[1,1]),"\n");
 else
 # get the degree and characteristic
 deg := Lcm( List( m, DegreeFFE ) );
 chr := Characteristic(m[1,1]);
 zero := ZeroOfBaseDomain(m);

 # if it is a finite prime field, use integers for display
 if deg = 1 then

 # compute maximal width
 w := LogInt( chr, 10 ) + 2;

 # create strings
 t := [];
 for x in [ 2 .. chr ] do
 t[x] := String( x-1, w );
 od;
#T useful only for (very) small characteristic, or?
 t[1] := String( ".", w );

 # print matrix
 for v in m do
 for x in List( v, IntFFE ) do
#T !
 Print( t[x+1] );
 od;
 Print( "\n" );
 od;

 # if it a finite, use mixed integers/z notation
#T ...
 else
 Print( "z = Z(", chr^deg, ")\n" );

 # compute maximal width
 w := LogInt( chr^deg-1, 10 ) + 4;

 # create strings
 t := [];
 f := GF(chr^deg);
 z := Z(chr^deg);
 for x in [ 0 .. Size(f)-2 ] do
 y := z^x;
 if DegreeFFE(y) = 1 then
 t[x+2] := String( IntFFE(y), w );
#T !
 else
 t[x+2] := String(Concatenation("z^",String(x)),w);
 fi;
 od;
 t[1] := String( ".", w );

 # print matrix
 for v in m do
 for x in v do
 if x = zero then
 Print( t[1] );
 else
 Print( t[LogFFE(x,z)+2] );
 fi;
 od;
 Print( "\n" );
 od;

 fi;
 fi;
end );

##########################################################################
##
#M Display( <ZmodnZ-mat> )
##
InstallMethod( Display,
 "for matrix over Integers mod n",
 [ IsZmodnZObjNonprimeCollColl and IsMatrix ],
 function( m )
 Print( "matrix over Integers mod ", Characteristic( m[1,1] ),
 ":\n" );
 Display( List( m, i -> List( i, i -> i![1] ) ) );
 end );

#############################################################################
##
#M CharacteristicPolynomial( <mat> )
##
InstallMethod( CharacteristicPolynomial,
 "supply field and indeterminate 1",
 [ IsMatrix ],
 mat -> CharacteristicPolynomialMatrixNC(
 DefaultFieldOfMatrix( mat ), mat, 1 ) );

#############################################################################
##
#M CharacteristicPolynomial( <F>, <E>, <mat> )
##
InstallMethod( CharacteristicPolynomial,
 "supply indeterminate 1",
 function (famF, famE, fammat)
 local fam;
 if HasElementsFamily (fammat) then
 fam := ElementsFamily (fammat);
 return IsIdenticalObj (famF, fam) and IsIdenticalObj (famE, fam);
 fi;
 return false;
 end,
 [ IsField, IsField, IsMatrix ],
 function( F, E, mat )
 return CharacteristicPolynomial( F, E, mat, 1);
 end );

#############################################################################
##
#M CharacteristicPolynomial( <mat>, <indnum> )
##
InstallMethod( CharacteristicPolynomial,
 "supply field",
 [ IsMatrix, IsPosInt ],
 function( mat, indnum )
 local F;
 F := DefaultFieldOfMatrix( mat );
 return CharacteristicPolynomial( F, F, mat, indnum );
 end );

#############################################################################
##
#M CharacteristicPolynomial( <subfield>, <field>, <matrix>, <indnum> )
##
InstallMethod( CharacteristicPolynomial, "spinning over field",
 function (famF, famE, fammat, famid)
 local fam;
 if HasElementsFamily (fammat) then
 fam := ElementsFamily (fammat);
 return IsIdenticalObj (famF, fam) and IsIdenticalObj (famE, fam);
 fi;
 return false;
 end,
 [ IsField, IsField, IsOrdinaryMatrix, IsPosInt ],
 function( F, E, mat, inum )
 local B;

 if not IsSubset (E, F) then
 Error ("<F> must be a subfield of <E>.");
 elif F <> E then
 # Replace the matrix by a matrix with the same char. polynomial
 # but with entries in `F'.
 B:= Basis( AsVectorSpace( F, E ) );
 mat:= BlownUpMat( B, mat );
 fi;

 return CharacteristicPolynomialMatrixNC( F, mat, inum);
 end );

InstallMethod( CharacteristicPolynomialMatrixNC, "spinning over field",
 IsElmsCollsX,
 [ IsField, IsOrdinaryMatrix, IsPosInt ],
 Matrix_CharacteristicPolynomialSameField);

#############################################################################
##
#M MinimalPolynomial( <field>, <matrix>, <indnum> )
##
InstallMethod( MinimalPolynomial,
 "spinning over field",
 IsElmsCollsX,
 [ IsField, IsOrdinaryMatrix, IsPosInt ],
function( F, mat,inum )
 local fld, B;

 fld:= DefaultFieldOfMatrix( mat );

 if fld <> fail and not IsSubset( F, fld ) then

 # Replace the matrix by a matrix with the same minimal polynomial
 # but with entries in `F'.
 if not IsSubset( fld, F ) then
 fld:= ClosureField( fld, F );
 fi;
 B:= Basis( AsField( F, fld ) );
 mat:= BlownUpMat( B, mat );

 fi;

 return MinimalPolynomialMatrixNC( F, mat,inum);
end );

InstallOtherMethod( MinimalPolynomial,
 "supply field",
 [ IsMatrix,IsPosInt ],
function(m,n)
 return MinimalPolynomial( DefaultFieldOfMatrix( m ), m, n );
end);

InstallOtherMethod( MinimalPolynomial,
 "supply field and indeterminate 1",
 [ IsMatrix ],
function(m)
 return MinimalPolynomial( DefaultFieldOfMatrix( m ), m, 1 );
end);

InstallMethod( MinimalPolynomialMatrixNC, "spinning over field",
 IsElmsCollsX,
 [ IsField, IsOrdinaryMatrix, IsPosInt ],
 Matrix_MinimalPolynomialSameField);

#############################################################################
##
#M Order( <mat> ) . . . . . . . . . . . . . . . . . . . . order of a matrix
##
OrderMatLimit := 1000;

InstallOtherMethod( Order,
 "generic method for ordinary matrices",
 [ IsOrdinaryMatrix ],
function ( mat )
 local m, rank;

 # check that the argument is an invertible square matrix
 m := NrRows(mat);
 if m <> NrCols(mat) then
 Error( "Order: <mat> must be a square matrix" );
 fi;
 rank:= RankMat( mat );
 if rank = fail then
 if not IsUnit( DeterminantMat( mat ) ) then
 Error( "Order: <mat> must be invertible" );
 fi;
 elif rank <> m then
 Error( "Order: <mat> must be invertible" );
#T also test here that the determinant is in fact a unit in the ring
#T that is generated by the matrix entries?
#T (Do we need `IsPossiblyInvertibleMat' and `IsSurelyInvertibleMat',
#T the first meaning that the inverse in some ring exists,
#T the second meaning that the inverse exists in the ring generated by the
#T matrix entries?
#T For `Order', it is `IsSurelyInvertibleMat' that one wants to check;
#T then one can return `infinity' if the determinant is not a unit in the
#T ring generated by the matrix entries.)
 fi;

 # loop over the standard basis vectors
 return OrderMatTrial(mat,infinity);
end );

#############################################################################
##
#F OrderMatTrial( <mat>,<lim> )
##
InstallGlobalFunction(OrderMatTrial,function(mat,lim)
local ord,i,vec,v,o;

 # loop over the standard basis vectors
 ord := 1;
 for i in [1..NrRows(mat)] do

 # compute the length of the orbit of the th standard basis vector
 # (equivalently, of the orbit of `mat[]',
 # the image of the standard basis vector under `mat')
 vec := mat[i];
 v := vec * mat;
 o := 1;
 while v <> vec do
 v := v * mat;
 o := o + 1;
 if o>lim then
 return fail;
 elif OrderMatLimit = o and Characteristic(v[1])=0 then
 Info( InfoWarning, 1,
 "Order: warning, order of <mat> might be infinite" );
 fi;
 od;

 # raise the matrix to this length (new mat will fix basis vector)
 if o>1 then
 mat := mat ^ o;
 ord := ord * o;
 fi;
 od;
 if IsOne(mat) then return ord; else return fail; fi;
end);

# #############################################################################
# ##
# #M Order( <cycmat> ) . . . . . . . . . . . order of a matrix of cyclotomics
# ##
# ## The idea is to compute the minimal polynomial of the matrix,
# ## and to decompose this polynomial into cyclotomic polynomials.
# ## This is due to R. Beals, who used it in his `grim' package for {\GAP}~3.
# ##
# InstallMethod( Order,
# "ordinary matrix of cyclotomics",
# [ IsOrdinaryMatrix and IsCyclotomicCollColl ],
# function( cycmat )
# local m, # dimension of the matrix
# trace, # trace of the matrix
# minpol, # minimal polynomial of the matrix
# n, # degree of `minpol'
# p, # loop over small primes
# t, # product of the primes `p'
# l, # product of the values `p-1'
# ord, # currently known factor of the order
# d, # loop over the indices of cyclotomic polynomials
# phi, # `Phi( d )'
# c, # `d'-th cyclotomic polynomial
# q; # quotient and remainder
#
# # Before we start with expensive calculations,
# # we check whether the matrix has a *small* order.
# ord:= OrderMatTrial( cycmat, OrderMatLimit - 1 );
# if ord <> fail then
# return ord;
# fi;
#
# # Check that the argument is an invertible square matrix.
# m:= NrRows( cycmat );
# if m <> NrCols( cycmat ) then
# Error( "Order: <cycmat> must be a square matrix" );
# elif RankMat( cycmat ) <> m then
# Error( "Order: <cycmat> must be invertible" );
# fi;
# #T Here I could compute the inverse;
# #T its trace could be checked, too.
# #T Additionally, if `mat' consists of (algebraic) integers
# #T and the inverse does not then the order of `mat' is infinite.
#
# # If the order is finite then the trace must be an algebraic integer.
# trace:= TraceMat( cycmat );
# if not IsIntegralCyclotomic( trace ) then
# return infinity;
# fi;
#
# # If the order is finite then the absolute value of the trace
# # is bounded by the dimension of the matrix.
# #T compute this (approximately) for arbitrary cyclotomics
# #T (by the way: why isn't this function called `AbsRat'?)
# if IsInt( trace ) and NrRows( cycmat ) < AbsInt( trace ) then
# return infinity;
# fi;
#
# # Compute the minimal polynomial of the matrix.
# minpol:= MinimalPolynomial( Rationals, cycmat );
# n:= DegreeOfLaurentPolynomial( minpol );
#
# # The order is finite if and only if the minimal polynomial
# # is a product of cyclotomic polynomials.
# # (Note that cyclotomic polynomials over the rationals are irreducible.)
#
# # A necessary condition is that the constant term of the polynomial
# # is $\pm 1$, since this holds for every cyclotomic polynomial.
# if AbsInt( Value( minpol, 0 ) ) <> 1 then
# return infinity;
# fi;
#
# # Another necessary condition is that no irreducible factor
# # may occur more than once.
# # (Note that the minimal polynomial divides $X^{ord} - 1$.)
# if not IsOne( Gcd( minpol, Derivative( minpol ) ) ) then
# return infinity;
# fi;
#
# # Compute an upper bound `t' for the numbers $i$ with the property
# # that $\varphi(i) \leq n$ holds.
# # (Let $p_k$ denote the $k$-th prime divisor of $i$,
# # and $q_k$ the $k$-th prime; then clearly $q_k \leq p_k$ holds.
# # Now let $h$ be the smallest *positive* integer --be careful that the
# # products considered below are not empty-- such that
# # $\prod_{k=1}^h ( q_k - 1 ) \geq n$, and set $t = \prod_{k=1}^h q_k$.
# # If $i$ has the property $\varphi(i) \leq n$ then
# # $i \leq n \frac{i}{\varphi(i)} = n \prod_{k} \frac{p_k}{p_k-1}$.
# # Replacing $p_k$ by $q_k$ means to replace the factor
# # $\frac{p_k}{p_k-1}$ by a larger factor,
# # and if $i$ has less than $h$ prime divisors then
# # running over the first $h$ primes increases the value of the product
# # again, so we get $i \leq n \prod_{k=1}^h \frac{q_k}{q_k-1} \leq t$.)
# p:= 2;
# t:= 2;
# l:= 1;
# while l < n do
# p:= NextPrimeInt( p );
# t:= t * p;
# l:= l * ( p - 1 );
# od;
#
# # Divide by each possible cyclotomic polynomial.
# ord:= 1;
# for d in [ 1 .. t ] do
#
# phi:= Phi( d );
# if phi <= n then
# c:= CyclotomicPolynomial( Rationals, d );
# q:= QuotientRemainder( minpol, c );
# if IsZero( q[2] ) then
# minpol:= q[1];
# n:= n - phi;
# ord:= Lcm( ord, d );
# if n = 0 then
#
# # The minimal polynomial is a product of cyclotomic polynomials.
# return ord;
#
# fi;
# fi;
# fi;
#
# od;
#
# # The matrix has infinite order.
# return infinity;
# end );

#############################################################################
##
#M Order( <mat> ) . . . . . . . . . . . . order of a matrix of cyclotomics
##
InstallMethod( Order,
 "for a matrix of cyclotomics, with Minkowski kernel",
 [ IsOrdinaryMatrix and IsCyclotomicCollColl ],

 function ( mat )

 local dim, F, tracemat, lat, red, det, trace, order, orddet, powdet,
 ordpowdet;

 # Check that the argument is an invertible square matrix.
 dim:= NrRows( mat );
 if dim <> NrCols( mat ) then
 Error( "Order: <mat> must be a square matrix" );
 fi;

 # Before we start with expensive calculations,
 # we check whether the matrix has a *very small* order.
 order:= OrderMatTrial( mat, 6 );
 if order <> fail then return order; fi;

 # We compute the determinant <det>, issue an error message in case <mat>
 # is not invertible, compute the order <orddet> of <det> and check
 # whether <mat>^<orddet> has small order.
 det := DeterminantMat( mat );
 if det = 0 then Error( "Order: <mat> must be invertible" ); fi;
 orddet := Order(det);
 if orddet = infinity then return infinity; fi;
 powdet := mat^orddet;
 ordpowdet := OrderMatTrial( powdet, 12 );
 if ordpowdet <> fail then return orddet * ordpowdet; fi;

 # If the order is finite then the trace must be an algebraic integer.
 trace := TraceMat( mat );
 if not IsIntegralCyclotomic( trace ) then return infinity; fi;

 # If the order is finite then the absolute value of the trace
 # is bounded by the dimension of the matrix.
 if IsInt( trace ) and NrRows( mat ) < AbsInt( trace ) then
 return infinity;
 fi;

 F:= DefaultFieldOfMatrix( mat );

 # Convert to a rational matrix if necessary.
 if 1 < Conductor( F ) then

 # Check whether the trace is larger than the dimension.
 tracemat := BlownUpMat( Basis(F), [[ trace ]] );
 if AbsInt(Trace(tracemat)) > NrRows(mat) * NrRows(tracemat)
 then return infinity; fi;

 mat:= BlownUpMat( Basis( F ), mat );
 dim:= NrRows( mat );
 fi;

 # Convert to an integer matrix if necessary.
 if not ForAll( mat, row -> ForAll( row, IsInt ) ) then

 # The following checks trace and determinant.
 lat:= InvariantLattice( GroupWithGenerators( [ mat ] ) );
 if lat = fail then
 return infinity;
 fi;
 mat:= lat * mat * Inverse( lat );

 fi;

 # Compute the order of the reduction modulo $2$.
 red:= mat * Z(2);
 ConvertToMatrixRep(red,2);
 order:= Order( red );
#T if OrderMatTrial was used above then call `ProjectiveOrder' directly?

 # Now use the theorem (see Morris Newman, Integral Matrices)
 # that `mat' has infinite order if the `2 * order'-th
 # power is not equal to the identity matrix.
 mat:= mat ^ order;
 if IsOne(mat) then
 return order;
 elif IsOne(mat ^ 2) then
 return 2 * order;
 else
 return infinity;
 fi;
 end );

#############################################################################
##
#M Order( <ffe-mat> ) . . . . . order of a matrix of finite field elements
##
InstallMethod( Order, "ordinary matrix of finite field elements", true,
 [ IsOrdinaryMatrix and IsFFECollColl ], 0,
 function( mat )
 local o;
 # catch the (unlikely in GL but likely in group theory...) case that mat
 # has a small order

 # the following limit is very crude but seems to work OK. It picks small
 # orders but still does not cost too much if the order gets larger.
 if NrRows(mat) <> NrCols(mat) then
 Error("Order of non-square matrix is not defined");
 fi;
 o:=Characteristic(mat[1,1])^DegreeFFE(mat[1,1]); # size of field of
 # first entry
 o:=QuoInt(NrRows(mat),o)*5;

 o:=OrderMatTrial(mat,o);
 if o<>fail then
 return o;
 fi;

 o := ProjectiveOrder(mat);
 return o[1] * Order( o[2] );
end );

#############################################################################
##
#M IsZero( <mat> )
##
InstallMethod( IsZero,
 "method for a matrix",
 [ IsMatrix ],
 function( mat )
 local ncols, # number of columns
 row; # loop over rows in 'obj'

 ncols:= NrCols( mat );
 for row in mat do
 if PositionNonZero( row ) <= ncols then
 return false;
 fi;
 od;
 return true;
 end );

#############################################################################
##
#M IsOne( <mat> )
##
InstallMethod( IsOne,
 "method for a matrix",
 [ IsMatrix ],
 function( mat )
 local ncols, # number of columns
 i,
 row; # loop over rows in 'obj'

 ncols:= NrCols( mat );
 for i in [1 .. NrRows( mat )] do
 row := mat[i];
 if PositionNonZero( row ) <> i or not IsOne( row[i] ) then
 return false;
 fi;
 if PositionNonZero( row, i ) <= ncols then
 return false;
 fi;
 od;
 return true;
 end );

#############################################################################
##
#M BaseMat( <mat> ) . . . . . . . . . . base for the row space of a matrix
##
InstallMethod( BaseMatDestructive,
 "generic method for matrices",
 [ IsMatrix ],
 mat -> SemiEchelonMatDestructive( mat ).vectors );

InstallMethod( BaseMat,
 "generic method for matrices",
 [ IsMatrix ],
 mat -> BaseMatDestructive( MutableCopyMatrix( mat ) ) );

#############################################################################
##
#M DefaultFieldOfMatrix( <mat> )
##
InstallMethod( DefaultFieldOfMatrix,
 "default method for a matrix (return `fail')",
 [ IsMatrix ],
 ReturnFail );

#############################################################################
##
#M DefaultFieldOfMatrix( <ffe-mat> )
##
InstallMethod( DefaultFieldOfMatrix,
 "method for a matrix over a finite field",
 [ IsMatrix and IsFFECollColl ],
function( mat )
 local deg, j;

 deg := 1;
 for j in mat do
 deg := LcmInt( deg, DegreeFFE(j) );
 od;
 return GF( Characteristic(mat[1]), deg );
end );

#############################################################################
##
#M DefaultFieldOfMatrix( <cyc-mat> )
##
InstallMethod( DefaultFieldOfMatrix,
 "method for a matrix over the cyclotomics",
 [ IsMatrix and IsCyclotomicCollColl ],
function( mat )
 local deg, j;

 deg := 1;
 for j in mat do
 deg := LcmInt( deg, Conductor(j) );
 od;
 return CF( deg );
end );

#############################################################################
##
#M BaseDomain( <v> )
#M OneOfBaseDomain( <v> )
#M ZeroOfBaseDomain( <v> )
#M ConstructingFilter( <v> )
#M BaseDomain( <mat> )
#M OneOfBaseDomain( <mat> )
#M ZeroOfBaseDomain( <mat> )
#M RowsOfMatrix( <mat> )
#M NumberRows( <mat> )
#M NumberColumns( <mat> )
#M ConstructingFilter( <mat> )
##
## Note that a row vector or a matrix that is a plain list
## is necessarily nonempty and dense,
## hence it is safe to access the first entry.
## Since the rows of a matrix that is a plain list are homogeneous and
## nonempty, one can also access the first entry in the first row.
##
InstallOtherMethod( BaseDomain,
 "generic method for a row vector",
 [ IsRowVector ],
 -SUM_FLAGS,
 DefaultRing );

InstallOtherMethod( OneOfBaseDomain,
 "generic method for a row vector",
 [ IsRowVector ],
 -SUM_FLAGS,
 v -> One( v[1] ) );

InstallOtherMethod( ZeroOfBaseDomain,
 "generic method for a row vector",
 [ IsRowVector ],
 -SUM_FLAGS,
 v -> Zero( v[1] ) );

InstallOtherMethod( ConstructingFilter,
 "generic method for a row vector",
 [ IsRowVector ],
 -SUM_FLAGS,
 v -> IsPlistRep );

InstallOtherMethod( BaseDomain,
 "generic method for a plain list matrix",
 [ IsMatrix and IsPlistRep ],
 mat -> BaseDomain( Concatenation( mat ) ) );

InstallOtherMethod( BaseDomain,
 "generic method for a plain list matrix over a finite field",
 [ IsMatrix and IsPlistRep and IsFFECollColl ],
 DefaultFieldOfMatrix );

InstallOtherMethod( BaseDomain,
 "generic method for a plain list over the cyclotomics",
 [ IsMatrix and IsPlistRep and IsCyclotomicCollColl ],
 DefaultFieldOfMatrix );

InstallOtherMethod( OneOfBaseDomain,
 "generic method for a matrix that is a plain list",
 [ IsMatrix and IsPlistRep ],
 mat -> One( mat[1,1] ) );

InstallOtherMethod( ZeroOfBaseDomain,
 "generic method for a matrix that is a plain list",
 [ IsMatrix and IsPlistRep ],
 mat -> Zero( mat[1,1] ) );

InstallMethod( RowsOfMatrix,
 "generic method for a matrix that is a plain list",
 [ IsMatrix and IsPlistRep ],
 Immutable );

InstallMethod( NumberRows,
 "generic method for a (perhaps empty) matrix",
 [ IsMatrix ],
 -SUM_FLAGS,
 Length );

InstallMethod( NumberColumns,
 "generic method for a (perhaps empty) matrix",
 [ IsMatrix ],
 -SUM_FLAGS,
 function( mat )
 if Length( mat ) = 0 then
 return 0;
 fi;
 return Length( mat[1] );
 end );

InstallOtherMethod( ConstructingFilter,
 "generic method for a matrix that is a plain list",
 [ IsMatrix and IsPlistRep ],
 M -> IsPlistRep );

#############################################################################
##
## The following "auxiliary methods" are needed when 'DiagonalMatrix' etc.
## is called with filter 'IsPlistRep'.
## (Strictly speaking, 'NewMatrix' and 'NewZeroMatrix' are not
## defined for this situation, since they promise to return matrix objects.)
##
InstallTagBasedMethod( NewMatrix,
 IsPlistRep,
 function( filter, basedomain, ncols, list )
 local copied, nd;

 # If applicable then replace a flat list 'list' by a nested list
 # of lists of length 'ncols'.
 copied:= false;
 if Length( list ) > 0 and not IsVectorObj( list[1] ) then
 nd:= NestingDepthA( list );
 if nd < 2 or nd mod 2 = 1 then
 if Length( list ) mod ncols <> 0 then
 Error( "NewMatrix: length of <list> is not a multiple of <ncols>" );
 fi;
 list:= List( [ 0, ncols .. Length( list )-ncols ],
 i -> list{ [ i+1 .. i+ncols ] } );
 copied:= true;
 fi;
 fi;

 if ValueOption( "check" ) <> false then
 if ForAny( list, l -> Length( l ) <> ncols ) then
 Error( "NewMatrix: all entries in <list> must have length <ncols>" );
 elif ForAny( list, l -> not IsSubset( basedomain, l ) ) then
 Error( "NewMatrix: entries of all in <list> must lie in <basedomain>" );
 fi;
 fi;

 if not copied then
 list:= List( list, ShallowCopy );
 fi;

 return list;
 end );

# This is faster than the default method.
InstallTagBasedMethod( NewZeroMatrix,
 IsPlistRep,
 function( filter, basedomain, nrows, ncols )
 return NullMat( nrows, ncols, basedomain );
 end );

#############################################################################
##
#M DepthOfUpperTriangularMatrix( <mat> )
##
InstallMethod( DepthOfUpperTriangularMatrix,
 [ IsMatrix ],
function( mat )
 local dim, zero, i, j;

 # find the correct layer of <m>
 dim := NrRows(mat);
 zero := ZeroOfBaseDomain(mat);
 for i in [ 1 .. dim-1 ] do
 for j in [ 1 .. dim-i ] do
 if mat[j,i+j] <> zero then
 return i;
 fi;
 od;
 od;
 return dim;

end);

InstallOtherMethod( SumIntersectionMat,
 [ IsEmpty, IsMatrix ],
function(a,b)
 b:=MutableCopyMatrix(b);
 TriangulizeMat(b);
 b:=Filtered(b,i->not IsZero(i));
 return [b,a];
end);

InstallOtherMethod( SumIntersectionMat,
 [ IsMatrix, IsEmpty ],
function(a,b)
 a:=MutableCopyMatrix(a);
 TriangulizeMat(a);
 a:=Filtered(a,i->not IsZero(i));
 return [a,b];
end);

InstallOtherMethod( SumIntersectionMat,
 IsIdenticalObj,
 [ IsEmpty, IsEmpty ],
function(a,b)
 return [a,b];
end);

#############################################################################
##
#M DeterminantMat( <mat> )
##
## Fractions free method, will never introduce denominators
##
## This method is better for cyclotomics, but pivoting is really needed
##
InstallMethod( DeterminantMatDestructive,
 "fraction-free method",
 [ IsOrdinaryMatrix and IsMutable],
 function ( mat )
 local det, sgn, row, zero, m, i, j, k, mult, row2, piv;

 # check that the argument is a square matrix and get the size
 m := NrRows(mat);
 zero := ZeroOfBaseDomain(mat);
 if m <> NrCols(mat) then
 Error("DeterminantMat: <mat> must be a square matrix");
 fi;

 # run through all columns of the matrix
 i := 0; det := 1; sgn := 1;
 for k in [1..m] do

 # find a nonzero entry in this column
 #N 26-Oct-91 martin if <mat> is a rational matrix look for a pivot
 j := i + 1;
 while j <= m and mat[j,k] = zero do j := j+1; od;

 # if there is a nonzero entry
 if j <= m then

 # increment the rank
 i := i + 1;

 # make its row the current row
 row := mat[j];
 if i <> j then
 SwapMatrixRows(mat, j, i);
 sgn := -sgn;
 fi;
 piv := row[k];

 # clear all entries in this column
 # Then divide through by det, this, amazingly, works, due
 # to a theorem about 3x3 determinants
 for j in [i+1..m] do
 row2 := mat[j];
 mult := -row2[k];
 if mult <> zero then
 MultVector( row2, piv );
 AddRowVector( row2, row, mult, k, m );
 MultVector( row2, Inverse(det) );
 else
 MultVector( row2, piv/det);
 fi;
 od;

 det := piv;
 else
 return zero;
 fi;

 od;

 # return the determinant
 return sgn * det;
end);

#############################################################################
##
#M DeterminantMat( <mat> )
##
## direct Gaussian elimination, not avoiding denominators
#T This method at the moment is better for finite fields
## another method is installed for cyclotomics. Anything else falls
## through here also.
##
InstallMethod( DeterminantMatDestructive,"non fraction free",
 [ IsOrdinaryMatrix and IsFFECollColl and IsMutable],
function( mat )
 local m, zero, det, sgn, k, j, row, l, row2, x;

 Info( InfoMatrix, 1, "DeterminantMat called" );

 # check that the argument is a square matrix, and get the size
 m := NrRows(mat);
 if m = 0 or m <> NrCols(mat) then
 Error( "<mat> must be a square matrix at least 1x1" );
 fi;
 zero := ZeroOfBaseDomain(mat);

 det := One(zero);
 sgn := det;

 # run through all columns of the matrix
 for k in [ 1 .. m ] do

 # look for a nonzero entry in this column
 j := k;
 while j <= m and mat[j,k] = zero do
 j := j+1;
 od;

 # if there is a nonzero entry
 if j <= m then

 # increment the rank, ...
 Info( InfoMatrix, 2, " nonzero columns: ", k );

 # ... make its row the current row, ...
 row := mat[j];
 if k <> j then
 SwapMatrixRows(mat, j, k);
 sgn := -sgn;
 fi;

 # ... and normalize the row.
 x := row[k];
 det := det * x;
 MultVector( mat[k], Inverse(x) );

 # clear all entries in this column, adjust only columns > k
 # (Note that we need not clear the rows from 'k+1' to 'j'.)
 for l in [ j+1 .. m ] do
 row2 := mat[l];
 x := row2[k];
 if x <> zero then
 AddRowVector( row2, row, -x, k+1, m );
 fi;
 od;

 # the determinant is zero
 else
 Info( InfoMatrix, 1, "DeterminantMat returns ", zero );
 return zero;
 fi;
 od;
 det := sgn * det;
 Info( InfoMatrix, 1, "DeterminantMat returns ", det );

 # return the determinant
 return det;

end );

InstallMethod( DeterminantMat,
 "for matrices",
 [ IsMatrix ],
 function( mat )
 return DeterminantMatDestructive( MutableCopyMatrix( mat ) );
 end );

InstallMethod( DeterminantMatDestructive,"nonprime residue rings",
 [ IsOrdinaryMatrix and
 CategoryCollections(CategoryCollections(IsZmodnZObjNonprime)) and IsMutable],
 DeterminantMatDivFree);

#############################################################################
##
#M DeterminantMatDivFree( <M> )
##
## Division free method. This is an alternative to the fraction free method
## when division of matrix entries is expensive or not possible.
##
## This method implements a division free algorithm found in
## Mahajan and Vinay \cite{MV97}.
##
## The run time is $O(n^4)$
## Auxiliary storage size $n^2+n + C$
##
## Our implementation has two runtime optimizations (both noted
## by Mahajan and Vinay)
## 1. Partial monomial sums, subtractions, and products are done at
## each level.
## 2. Prefix property is maintained allowing for a pruning of many
## vertices at each level
##
## and two auxiliary storage size optimizations
## 1. only the upper triangular and diagonal portion of the
## auxiliary storage is used.
## 2. Level information storage is reused (2 levels).
##
## This code was implemented by:
## Timothy DeBaillie
## Robert Morse
## Marcus Wassmer
##
InstallMethod( DeterminantMatDivFree,
 "Division-free method",
 [ IsMatrix ],
 function ( M )
 local u,v,w,i, ## indices
 a,b,c,x,y, ## temp indices
 temp, ## temp variable
 nlevel, ## next level
 clevel, ## current level
 pmone, ## plus or minus one
 zero, ## zero of the ring
 n, ## size of the matrix
 Vs, ## final sum
 V; ## graph

 # check that the argument is a square matrix and set the size
 n := NumberRows( M );
 if not n = NumberColumns( M ) or not IsRectangularTable(M) then
 Error("DeterminantMatDivFree: <mat> must be a square matrix");
 fi;

 ## initialize the final sum, the vertex set, initial parity
 ## and level indexes
 ##
 zero := ZeroOfBaseDomain( M );
 Vs := zero;
 V := [];
 pmone := (-OneOfBaseDomain( M ))^((n mod 2)+1);
#T better if/else
 clevel := 1; nlevel := 2;

 ## Vertices are indexed [u,v,i] holding the (partial) monomials
 ## whose sums will form the determinant
 ## where i = depth in the tree (current and next reusing
 ## level storage)
 ## u,v indices in the matrix
 ##
 ## Only the upper triangular portion of the storage space is
 ## needed. It is easier to create lower triangular data type
 ## which we do here and index via index arithmetic.
 ##
 for u in [1..n] do
 Add(V,[]);
 for v in [1..u] do
 Add(V[u],[zero,zero]);
 od;
 ## Initialize the level 0 nodes with +/- one, depending on
 ## the initial parity determined by the size of the matrix
 ##
 V[u][u][clevel] := pmone;
 od;

 ## Here are the $O(n^4)$ edges labeled by the elements of
 ## the matrix $M$. We build up products of the labels which form
 ## the monomials which make up the determinant.
 ##
 ## 1. Parity of monomials are maintained implicitly.
 ## 2. Partial sums for some vertices are not part of the final
 ## answer and can be pruned.
 ##
 for i in [0..n-2] do
 for u in [1..i+2] do ## <---- pruning of vertices
 for v in [u..n] do ## (maintains the prefix property)
 for w in [u+1..n] do
#T for b in [ n-u, n-u-1 .. 1 ] do

 ## translate indices to lower triangluar coordinates
 ##
 a := n-u+1; b := n-w+1; c := n-v+1;
#T move a to for u ...
#T move c to for v ...
 V[a][b][nlevel]:= V[a][b][nlevel]+
 V[a][c][clevel] * M[ v, w ];
 V[b][b][nlevel]:= V[b][b][nlevel]-
 V[a][c][clevel] * M[ v, u ];
 od;
 od;
 od;

 ## set the new current and next level. The new next level
 ## is initialized to zero
 ##
 temp := nlevel; nlevel := clevel; clevel := temp;
 for x in [1..n] do
 for y in [1..x] do
 V[x][y][nlevel] := zero;
 od;
 od;
 od;

 ## with the final level, we form the last monomial product and then
 ## sum these monomials (parity has been accounted for)
 ## to find the determinant.
 ##
 for u in [1..n] do
 for v in [u..n] do
 Vs := Vs + V[n-u+1][n-v+1][clevel] * M[ v, u ];
 od;
 od;

 ## Return the final sum
 ##
 return Vs;

 end);

#############################################################################
##
#M DimensionsMat( <mat> )
##
InstallOtherMethod( DimensionsMat,
 [ IsMatrix ],
 function( A )
 if IsRectangularTable(A) then
 return [ NrRows(A), NrCols(A) ];
 else
 return fail;
 fi;
 end );

BindGlobal("DoDiagonalizeMat",function(arg)
local R,M,transform,divide,swaprow, swapcol, addcol, addrow, multcol, multrow, l, n, start, d,
 ed, posi,posj, a, b, qr, c, i,j,left,right,cleanout, alldivide,basmat,origtran;

 R:=arg[1];
 M:=arg[2];
 transform:=arg[3];
 divide:=arg[4];

 l:=NrRows(M);
 n:=NrCols(M);

 basmat:=fail;
 if transform then
 left:=IdentityMat(l,R);
 right:=IdentityMat(n,R);
 if Length(arg)>4 then
 origtran:=arg[5];
 basmat:=IdentityMat(l,R); # for RCF -- transpose of P' in D&F, sec. 12.2
 fi;
 fi;

 swaprow:=function(a,b)
 SwapMatrixRows(M, a, b);
 if transform then
 SwapMatrixRows(left, a, b);
 if basmat<>fail then
 SwapMatrixRows(basmat, a, b);
 fi;
 fi;
 end;

 swapcol:=function(a,b)
 SwapMatrixColumns(M, a, b);
 if transform then
 SwapMatrixColumns(right, a, b);
 fi;
 end;

 addcol:=function(a,b,m)
 local i;
 for i in [1..l] do
 M[i,a]:=M[i,a]+m*M[i,b];
 od;
 if transform then
 for i in [1..n] do
 right[i,a]:=right[i,a]+m*right[i,b];
 od;
 fi;
 end;

 addrow:=function(a,b,m)
 AddCoeffs(M[a],M[b],m);
 if transform then
 AddCoeffs(left[a],left[b],m);
 if basmat<>fail then
 basmat[b]:=basmat[b]-basmat[a]*Value(m,origtran);
 fi;
 fi;
 end;

 multcol:=function(a,m)
 local i;
 for i in [1..l] do
 M[i,a]:=M[i,a]*m;
 od;
 if transform then
 for i in [1..n] do
 right[i,a]:=right[i,a]*m;
 od;
 fi;
 end;

 multrow:=function(a,m)
 MultVector(M[a],m);
 if transform then
 MultVector(left[a],m);
 if basmat<>fail then
 MultVector(basmat[a],1/m);
 fi;
 fi;
 end;

 # clean out row and column
 cleanout:=function()
 local a,i,b,c,qr;
 repeat
 # now do the GCD calculations only in row/column
 for i in [start+1..n] do
 a:=i;
 b:=start;
 if not IsZero(M[start,b]) then
 repeat
 qr:=QuotientRemainder(R,M[start,a],M[start,b]);
 addcol(a,b,-qr[1]);
 c:=a;a:=b;b:=c;
 until IsZero(qr[2]);
 if b=start then
 swapcol(start,i);
 fi;
 fi;

 # normalize
 qr:=StandardAssociateUnit(R,M[start,start]);
 multcol(start,qr);

 od;

 for i in [start+1..l] do
 a:=i;
 b:=start;
 if not IsZero(M[b,start]) then
 repeat
 qr:=QuotientRemainder(R,M[a,start],M[b,start]);
 addrow(a,b,-qr[1]);
 c:=a;a:=b;b:=c;
 until IsZero(qr[2]);
 if b=start then
 swaprow(start,i);
 fi;
 fi;

 qr:=StandardAssociateUnit(R,M[start,start]);
 multrow(start,qr);

 od;
 until ForAll([start+1..n],i->IsZero(M[start,i]));
 end;

 start:=1;
 while start<=NrRows(M) and start<=n do

 # find element of lowest degree and move it into pivot
 # hope is this will reduce the total number of iterations by making
 # it small in the first place
 d:=infinity;

 for i in [start..l] do
 for j in [start..n] do
 if not IsZero(M[i,j]) then
 ed:=EuclideanDegree(R,M[i,j]);
 if ed<d then
 d:=ed;
 posi:=i;
 posj:=j;
 fi;
 fi;
 od;
 od;

 if d<>infinity then # there is at least one nonzero entry

 if posi<>start then
 swaprow(start,posi);
 fi;
 if posj<>start then
 swapcol(start,posj);
 fi;
 cleanout();

 if divide then
 repeat
 alldivide:=true;
 #make sure the pivot also divides the rest
 for i in [start+1..l] do
 for j in [start+1..n] do
 if Quotient(M[i,j],M[start,start])=fail then
 alldivide:=false;
 # do gcd
 addrow(start,i,One(R));
 cleanout();
 fi;
 od;
 od;
 until alldivide;

 fi;

 # normalize
 qr:=StandardAssociateUnit(R,M[start,start]);
 multcol(start,qr);

 fi;
 start:=start+1;
 od;

 if transform then
 M:=rec(rowtrans:=left,coltrans:=right,normal:=M);
 if basmat<>fail then
 M.basmat:=basmat;
 fi;
 return M;
 else
 return M;
 fi;
end);

#############################################################################
##
#M DiagonalizeMat(<euclring>,<mat>)
##
# this is a very naive implementation but it should work for any euclidean
# ring.
InstallMethod( DiagonalizeMat,
 "method for general Euclidean Ring",
 true, [ IsEuclideanRing,IsMatrix and IsMutable], 0,function(R,M)
 return DoDiagonalizeMat(R,M,false,false);
end);

#############################################################################
##
#M ElementaryDivisorsMat(<mat>) . . . . . . elementary divisors of a matrix
##
## 'ElementaryDivisors' returns a list of the elementary divisors, i.e., the
## unique <d> with '<d>[]' divides '<d>[+1]' and <mat> is equivalent
## to a diagonal matrix with the elements '<d>[]' on the diagonal.
##

InstallGlobalFunction(ElementaryDivisorsMatDestructive,function(ring,mat)
 # diagonalize the matrix
 DoDiagonalizeMat(ring, mat,false,true );

 # get the diagonal elements
 return DiagonalOfMat(mat);
end );

InstallMethod( ElementaryDivisorsMat,
 "generic method for euclidean rings",
 [ IsEuclideanRing,IsMatrix ],
function ( ring,mat )
 # make a copy to avoid changing the original argument
 mat := MutableCopyMatrix( mat );
 if IsIdenticalObj(ring,Integers) then
 DiagonalizeMat(Integers,mat);
 return DiagonalOfMat(mat);
 fi;
 return ElementaryDivisorsMatDestructive(ring,mat);
end);

InstallOtherMethod( ElementaryDivisorsMat,
 "compatibility method -- supply ring",
 [ IsMatrix ],
function(mat)
local ring;
 if ForAll(mat,row->ForAll(row,IsInt)) then
 return ElementaryDivisorsMat(Integers,mat);
 fi;
 ring:=DefaultRing(Flat(mat));
 return ElementaryDivisorsMat(ring,mat);
end);

#############################################################################
##
#M ElementaryDivisorsTransformationsMat(<mat>) elem. divisors of a matrix
##
## 'ElementaryDivisorsTransformationsMat' does not only compute the
## elementary divisors, but also transforming matrices.

InstallGlobalFunction(ElementaryDivisorsTransformationsMatDestructive,
function(ring,mat)

 # diagonalize the matrix
 return DoDiagonalizeMat(ring, mat,true,true );

end );

InstallMethod( ElementaryDivisorsTransformationsMat,
 "generic method for euclidean rings",
 [ IsEuclideanRing,IsMatrix ],
function ( ring,mat )
 # make a copy to avoid changing the original argument
 mat := MutableCopyMatrix( mat );
 return ElementaryDivisorsTransformationsMatDestructive(ring,mat);
end);

InstallOtherMethod( ElementaryDivisorsTransformationsMat,
 "compatibility method -- supply ring",
 [ IsMatrix ],
function(mat)
local ring;
 if ForAll(mat,row->ForAll(row,IsInt)) then
 return ElementaryDivisorsTransformationsMat(Integers,mat);
 fi;
 ring:=DefaultRing(Flat(mat));
 return ElementaryDivisorsTransformationsMat(ring,mat);
end);

#############################################################################
##
#M MutableCopyMatrix( <mat> )
##
InstallOtherMethod( MutableCopyMatrix, "generic method", [ IsMatrix ],
 -SUM_FLAGS,
 mat -> List( mat, ShallowCopy ) );

InstallOtherMethod( MutableCopyMatrix, "for empty lists", [ IsList and IsEmpty ],
 mat -> [ ] );

#############################################################################
##
#M MutableTransposedMat( <mat> ) . . . . . . . . . . transposed of a matrix
##
InstallOtherMethod( MutableTransposedMat,
 "generic method",
 [ IsRectangularTable and IsMatrix ],
 function( mat )
 local trn, n, m, j;

 m:= NrRows( mat );
 if m = 0 then return []; fi;

 # initialize the transposed
 m:= [ 1 .. m ];
 n:= [ 1 .. NrCols( mat ) ];
 trn:= [];

 # copy the entries
 for j in n do
 trn[j]:= mat{ m }[j];
# ConvertToVectorRepNC( trn[j] );
 od;

 # return the transposed
 return trn;
end );

#############################################################################
##
#M MutableTransposedMat( <mat> ) . . . . . . . . . . transposed of a matrix
##
InstallOtherMethod( MutableTransposedMat,
 "for arbitrary lists of lists",
 [ IsList ],
 function( t )
 local res, m, i, j, row;
 res := [];
 if Length(t)>0 and IsDenseList(t) and ForAll(t, IsDenseList) then
 # special case with dense list of dense lists
 m := Maximum(List(t, Length));
 for i in [m,m-1..1] do
 res[i] := [];
 od;
 for i in [1..Length(t)] do
 res{[1..Length(t[i])]}[i] := t[i];
 od;
 else
 # general case, non dense lists allowed
 for i in [1..Length(t)] do
 if IsBound(t[i]) then
 row := t[i];
 if IsList(row) then
 for j in [1..Length(row)] do
 if IsBound(row[j]) then
 if not IsBound(res[j]) then
 res[j] := [];
 fi;
 res[j,i] := row[j];
 fi;
 od;
 else
 Error("bound entries must be lists");
 fi;
 fi;
 od;
 fi;
 return res;
end);

#############################################################################
##
#M MutableTransposedMatDestructive( <mat> ) . . . . . transposed of a matrix
## may destroy `mat'.
##
InstallMethod( MutableTransposedMatDestructive,
 "generic method",
 [IsMatrix and IsMutable],
 function( mat )

 local m, n, min, i, j, store;

 m:= NrRows( mat );
 if m = 0 then return []; fi;

 n:= NrCols( mat );
 min:= Minimum( m, n );

 # swap the entries in the "square part" of the matrix.
 for i in [1..min] do
 for j in [i+1..min] do
 store:= mat[i,j];
 mat[i,j]:= mat[j,i];
 mat[j,i]:= store;
 od;
 od;

 # if the matrix is not square, then we have to adjust some rows or
 # columns.
 if m < n then
 for i in [1..n-m] do
 store:= [ ];
 for j in [1..m] do
 store[j]:= mat[j,m+i];
 Unbind( mat[j][m+i] );
 od;
 Add( mat, store );
 od;
 for i in [1..m] do
 mat[i]:= Compacted( mat[i] );
 od;
 fi;

 if m > n then
 for i in [n+1..m] do
 for j in [1..n] do
 mat[j,i]:= mat[i,j];
 od;
 Unbind( mat[i] );
 od;
 mat:= Compacted( mat );
 fi;

 # return the transposed
 return mat;
end );

#############################################################################
##
#M NullspaceMat( <mat> ) . . . . . . basis of solutions of <vec> * <mat> = 0
##
InstallMethod( NullspaceMat,
 "generic method for ordinary matrices",
 [ IsOrdinaryMatrix ],
 mat -> SemiEchelonMatTransformation(mat).relations );

InstallOtherMethod(NullspaceMat,"matrix objects",[IsMatrixObj],
function(mat)
local r,nr,nc,n,i,j;
 r:=SemiEchelonMatTransformation(mat).relations;
 if Length(r)=0 then return r;fi;
 nr:=Length(r);
 nc:=Length(r[1]);
 n:=ZeroMatrix(BaseDomain(mat),nr,nc);
 for i in [1..nr] do
 for j in [1..nc] do
 n[i,j]:=r[i][j];
 od;
 od;
 return n;
end);

InstallMethod( NullspaceMatDestructive,
 "generic method for ordinary matrices",
 [ IsOrdinaryMatrix and IsMutable],
 mat -> SemiEchelonMatTransformationDestructive(mat).relations );

InstallOtherMethod( TriangulizedNullspaceMat,
 "generic method for matrices",
 [ IsMatrixOrMatrixObj ],
 mat -> TriangulizedNullspaceMatDestructive( MutableCopyMatrix( mat ) ) );

InstallOtherMethod( TriangulizedNullspaceMatDestructive,
 "generic method for matrices",
 [ IsMatrixOrMatrixObj and IsMutable],
 function( mat )
 local ns;
 ns := SemiEchelonMatTransformationDestructive(mat).relations;
 if IsPlistRep(ns) and Length(ns)>0
 # filter for vector objects, not compressed FF vectors
 and ForAll(ns,x->IsVectorObj(x) and not IsDataObjectRep(x)) then
 ns:=Matrix(BaseDomain(mat),ns);
 fi;
 TriangulizeMat(ns);
 return ns;
end );

InstallMethod( TriangulizedNullspaceMatNT,
 "generic method",
 [ IsOrdinaryMatrix ],
 function( mat )
 local nullspace, n, empty, i, k, row, zero, one;#

 TriangulizeMat( mat );
 n := NrCols(mat);
 zero := ZeroOfBaseDomain( mat );
 one := OneOfBaseDomain( mat );

 # insert empty rows to bring the leading term of each row on the diagonal
 empty := 0*mat[1];
 i := 1;
 while i <= NrRows(mat) do
 if i < n and mat[i,i] = zero then
 for k in Reversed([i..Minimum(NrRows(mat),n-1)]) do
 mat[k+1] := mat[k];
 od;
 mat[i] := empty;
 fi;
 i := i+1;
 od;
 for i in [ NrRows(mat)+1 .. n ] do
 mat[i] := empty;
 od;

 # 'mat' now looks like [ [1,2,0,2], [0,0,0,0], [0,0,1,3], [0,0,0,0] ],
 # and the solutions can be read in those columns with a 0 on the diagonal
 # by replacing this 0 by a -1, in this example [2,-1,0,0], [2,0,3,-1].
 nullspace := [];
 for k in Reversed([1..n]) do
 if mat[k,k] = zero then
 row := [];
 for i in [1..k-1] do row[n-i+1] := -mat[i,k]; od;
 row[n-k+1] := one;
 for i in [k+1..n] do row[n-i+1] := zero; od;
 ConvertToVectorRepNC( row );
 Add( nullspace, row );
 fi;
 od;

 return nullspace;
end );

#InstallMethod(TriangulizedNullspaceMat,"generic method",
# [IsOrdinaryMatrix],
# function ( mat )
# # triangulize the transposed of the matrix
# return TriangulizedNullspaceMatNT(
# MutableTransposedMat( Reversed( mat ) ) );
#end );

#InstallMethod(TriangulizedNullspaceMatDestructive,"generic method",
# [IsOrdinaryMatrix],
# function ( mat )
# # triangulize the transposed of the matrix
# return TriangulizedNullspaceMatNT(
# MutableTransposedMatDestructive( Reversed( mat ) ) );
#end );

#############################################################################
##
#M GeneralisedEigenvalues( <F>, <A> )
##
InstallMethod( GeneralisedEigenvalues,
 "for a matrix or matrix obj",
 [ IsField, IsMatrixOrMatrixObj ],
 function( F, A )
 return Set( Factors( UnivariatePolynomialRing(F), MinimalPolynomial(F, A,1) ) );
 end );

#############################################################################
##
#M GeneralisedEigenspaces( <F>, <A> )
##
InstallMethod( GeneralisedEigenspaces,
 "for a matrix or matrix obj",
 [ IsField, IsMatrixOrMatrixObj ],
 function( F, A )
 local M,S,eval;
 S:=[];
 for eval in GeneralisedEigenvalues(F,A) do
 M:=TriangulizedNullspaceMat( Value( eval, A ) );
 if IsMatrixObj(M) and not IsMatrix(M) then
 M:=RowsOfMatrix(M); fi;
 Add(S,VectorSpace(F,M));
 od;
 return S;
 end );

#############################################################################
##
#M Eigenvalues( <F>, <A> )
##
InstallMethod( Eigenvalues,
 "for a matrix",
 [ IsField, IsMatrixOrMatrixObj ],
 function( F, A )
 return List( Filtered( GeneralisedEigenvalues(F,A),
 eval -> DegreeOfLaurentPolynomial(eval) = 1 ),
 eval -> -1 * Value(eval,0) );
 end );

#############################################################################
##
#M Eigenspaces( <F>, <A> )
##
InstallMethod( Eigenspaces,
 "for a matrix",
 [ IsField, IsMatrix ],
 function( F, A )
 return List( Eigenvalues(F,A), eval ->
 VectorSpace( F, TriangulizedNullspaceMat(A - eval*One(A)) ) );
 end );

#############################################################################
##
#M Eigenvectors( <F>, <A> )
##
InstallMethod( Eigenvectors,
 "for a matrix",
 [ IsField, IsMatrix ],
 function( F, A )
 return Concatenation( List( Eigenspaces(F,A),
 esp -> AsList(Basis(esp)) ) );
 end );

#############################################################################
##
#M ProjectiveOrder( <mat> ) . . . . . . . . . . . . . . . order mod scalars
##
InstallMethod( ProjectiveOrder,
 "ordinary matrix of finite field elements",
 [ IsOrdinaryMatrix and IsFFECollColl ],
function( mat )
 local p, c;

 # construct the minimal polynomial of <A>
 p := MinimalPolynomialMatrixNC( DefaultFieldOfMatrix(mat), mat,1 );

 # check if <A> is invertible
 c := CoefficientsOfUnivariatePolynomial(p);
 if IsZero(c[1]) then
 Error( "matrix <mat> must be invertible" );
 fi;

 # compute the order of 
 return ProjectiveOrder(p);
end );

#############################################################################
##
#M RankMat( <mat> ) . . . . . . . . . . . . . . . . . . . rank of a matrix
##
InstallOtherMethod( RankMatDestructive,
 "generic method for mutable matrices",
 [ IsMatrix and IsMutable ],
 function( mat )
 mat:= SemiEchelonMatDestructive( mat );
 if mat <> fail then
 mat:= Length( mat.vectors );
 fi;
 return mat;
 end );

InstallOtherMethod( RankMat,
 "generic method for matrices",
 [ IsMatrix ],
 mat -> RankMatDestructive( MutableCopyMatrix( mat ) ) );

#############################################################################
##
#M SemiEchelonMat( <mat> )
##
InstallOtherMethod( SemiEchelonMatDestructive,
 "generic method for matrices",
 [ IsMatrixOrMatrixObj and IsMutable ],
 function( mat )
 local zero, # zero of the field of <mat>
 nrows, # number of rows in <mat>
 ncols, # number of columns in <mat>
 vectors, # list of basis vectors
 heads, # list of pivot positions in `vectors'
 i, # loop over rows
 j, # loop over columns
 x, # a current element
 nzheads, # list of non-zero heads
 row, # the row of current interest
 inv; # inverse of a matrix entry

 nrows:= NrRows( mat );
 ncols:= NrCols( mat );

 zero:= ZeroOfBaseDomain( mat );

 heads:= ListWithIdenticalEntries( ncols, 0 );
 nzheads := [];
 vectors := [];

 for i in [ 1 .. nrows ] do

 row := mat[i];
 # Reduce the row with the known basis vectors.
 for j in [ 1 .. Length(nzheads) ] do
 x := row[nzheads[j]];
 if x <> zero then
 AddRowVector( row, vectors[ j ], - x );
 fi;
 od;

 j := PositionNonZero( row );
 if j <= ncols then

 # We found a new basis vector.
 inv:= Inverse( row[j] );
 if inv = fail then
 return fail;
 fi;
 MultVector( row, inv );
 Add( vectors, row );
 Add( nzheads, j );
 heads[j]:= Length( vectors );

 fi;

 od;

 return rec( heads := heads,
 vectors := vectors );
 end );

InstallMethod( SemiEchelonMat,
 "generic method for matrices",
 [ IsMatrix ],
 function( mat )
 local copymat, v, vc, f;
 copymat := [];
 f := DefaultFieldOfMatrix(mat);
 for v in mat do
 vc := ShallowCopy(v);
 ConvertToVectorRepNC(vc,f);
 Add(copymat, vc);
 od;
 return SemiEchelonMatDestructive( copymat );
end );

InstallOtherMethod( SemiEchelonMat,
 "generic method for list of vector objects",
 [ IsList ],
 function( mat )
 local copymat, v, vc, f;
 # filter for vector objects, not compressed FF vectors
 if not (ForAll(mat,x->IsVectorObj(x) and not IsDataObjectRep(x))
 and Length(mat)>0) then
 TryNextMethod();
 fi;
 copymat := [];
 if Length(mat)>0 then
 f := BaseDomain(mat[1]);
 for v in mat do
 vc := ShallowCopy(v);
 Add(copymat, vc);
 od;
 fi;
 return SemiEchelonMatDestructive(Matrix(f, copymat ));
end );

#############################################################################
##
#M SemiEchelonMatTransformation( <mat> )
##
InstallMethod( SemiEchelonMatTransformation,
 "generic method for matrices",
 [ IsMatrix ],
 function( mat )
 local copymat, v, vc, f;
 copymat := [];
 f := DefaultFieldOfMatrix(mat);
 for v in mat do
 vc := ShallowCopy(v);
 ConvertToVectorRepNC(vc,f);
 Add(copymat, vc);
 od;
 return SemiEchelonMatTransformationDestructive( copymat );
end);

InstallOtherMethod(SemiEchelonMatTransformation,"matrix objects",[IsMatrixObj],
function( mat )
 local copymat;
 copymat:=MutableCopyMatrix(mat);
 return SemiEchelonMatTransformationDestructive( copymat );
end);

BindGlobal("DoSemiEchelonMatTransformationDestructive", function(f, mat )
 local zero, # zero of the field of <mat>
 nrows, # number of rows in <mat>
 ncols, # number of columns in <mat>
 vectors, # list of basis vectors
 heads, # list of pivot positions in 'vectors'
 i, # loop over rows
 j, # loop over columns
 T, # transformation matrix
 coeffs, # list of coefficient vectors for 'vectors'
 relations, # basis vectors of the null space of 'mat'
 row, head, x, row2;

 nrows := NrRows( mat );
 ncols := NrCols( mat );

 if f = fail then
 f := mat[1,1];
 fi;
 zero := Zero(f);

 heads := ListWithIdenticalEntries( ncols, 0 );
 vectors := [];

 T := IdentityMat( nrows, f );
 if IsMatrixObj(mat) then
 T:=Matrix(f,T);
 fi;
 coeffs := [];
 relations := [];

 for i in [ 1 .. nrows ] do

 row := mat[i];
 row2 := T[i];

 # Reduce the row with the known basis vectors.
 for j in [ 1 .. ncols ] do
 head := heads[j];
 if head <> 0 then
 x := - row[j];
 if x <> zero then
 AddRowVector( row2, coeffs[ head ], x );
 AddRowVector( row, vectors[ head ], x );
 fi;
 fi;
 od;

 j:= PositionNonZero( row );
 if j <= ncols then

 # We found a new basis vector.
 x:= Inverse( row[j] );
 if x = fail then
 TryNextMethod();
 fi;
 Add( coeffs, row2 * x );
 Add( vectors, row * x );
 heads[j]:= Length( vectors );

 else
 Add( relations, row2 );
 fi;

 od;

 return rec( heads := heads,
 vectors := vectors,
 coeffs := coeffs,
 relations := relations );
end );

InstallOtherMethod( SemiEchelonMatTransformationDestructive,
 "generic method for matrices",
 [ IsMatrix and IsMutable],
 mat->DoSemiEchelonMatTransformationDestructive(
 DefaultFieldOfMatrix(mat),mat));

InstallOtherMethod( SemiEchelonMatTransformationDestructive,
 "generic method for matrix objects over fields",
 [ IsMatrixObj and IsRowListMatrix and IsMutable],function(mat)
local f;
 f:=BaseDomain(mat);
 if not IsField(f) then TryNextMethod();fi;
 return DoSemiEchelonMatTransformationDestructive(f,mat);
end);

#############################################################################
##
#M SemiEchelonMats( <mats> )
##
InstallGlobalFunction( SemiEchelonMatsNoCo, function( mats )
 local zero, # zero coefficient
 m, # number of rows
 n, # number of columns
 v,
 vectors, # list of matrices in the echelonized basis
 heads, # list with info about leading entries
 mat, # loop over generators of 'V'
 i, j, # loop over rows and columns of the matrix
 k, l,
 mij,
 scalar,
 x;

 zero:= ZeroOfBaseDomain( mats[1] );
 m:= NrRows( mats[1] );
 n:= NrCols( mats[1] );

 # Compute an echelonized basis.
 vectors := [];
 heads := ListWithIdenticalEntries( n, 0 );
 heads := List( [ 1 .. m ], x -> ShallowCopy( heads ) );

 for mat in mats do

 # Reduce the matrix modulo 'ech'.
 for i in [ 1 .. m ] do
 for j in [ 1 .. n ] do
 if heads[i][j] <> 0 and mat[i,j] <> zero then

 # Compute 'mat:= mat - mat[i,j] * vectors[ heads[i][j] ];'
 scalar:= - mat[i,j];
 v:= vectors[ heads[i][j] ];
 for k in [ 1 .. m ] do
 AddRowVector( mat[k], v[k], scalar );
 od;

 fi;
 od;
 od;

 # Get the first nonzero column.
 i:= 1;
 j:= PositionNonZero( mat[1] );
 while n < j and i < m do
 i:= i + 1;
 j:= PositionNonZero( mat[i] );
 od;

 if j <= n then

 # We found a new basis vector.
 mij:= mat[i,j];
 for k in [ 1 .. m ] do
 for l in [ 1 .. n ] do
 x:= Inverse( mij );
 if x = fail then
 TryNextMethod();
 fi;
 mat[k,l]:= mat[k,l] * x;
 od;
 od;

 Add( vectors, mat );
 heads[i][j]:= Length( vectors );

 fi;

 od;

 # Return the result.
 return rec(
 vectors := vectors,
 heads := heads
 );
end );

InstallMethod( SemiEchelonMats,
 "for list of matrices",
 [ IsList ],
 mats -> SemiEchelonMatsNoCo( List( mats, MutableCopyMatrix ) ) );

InstallMethod( SemiEchelonMatsDestructive,
 "for list of matrices",
 [ IsList ],
 SemiEchelonMatsNoCo );

#############################################################################
##
#M TransposedMat( <mat> ) . . . . . . . . . . . . . transposed of a matrix
##
InstallOtherMethod( TransposedMat,
 "generic method for matrices and lists",
 [ IsList ],
 MutableTransposedMat );

#############################################################################
##
#M TransposedMatDestructive( <mat> ) . . . . . . . . transposed of a matrix
##
InstallMethod( TransposedMatDestructive,
 "generic method for matrices",
 [ IsMatrix ],
 MutableTransposedMatDestructive );

InstallOtherMethod(TransposedMatDestructive,
 "method for empty matrices",[IsList],
function(mat)
 if mat<>[] and mat<>[[]] then
 TryNextMethod();
 fi;
 return mat;
end);

############################################################################
##

#M IsMonomialMatrix( <mat> )
##
InstallMethod( IsMonomialMatrix,
 "generic method for matrices",
 [ IsMatrix ],
 function( M )
 local len, # length of rows
 found, # store positions of nonzero elements
 row, # loop over rows
 j; # position of first non-zero element

 len:= NrCols( M );
 if NrRows( M ) <> len then
 return false;
 fi;
 found:= ListWithIdenticalEntries( len, false );
 for row in M do
 j := PositionNonZero( row );
 if len < j or found[j] then
 return false;
 fi;
 if PositionNonZero( row, j ) <= len then
 return false;
 fi;
 found[j] := true;
 od;
 return true;
end );

##########################################################################
##
#M InverseMatMod( <cyc-mat>, <integer> )
##
InstallMethod( InverseMatMod,
 "generic method for matrix and integer",
 IsCollCollsElms,
 [ IsMatrix, IsInt ],
function( mat, m )
 local n, MM, inv, perm, i, pj, elem, liste, l;

 if NrRows(mat) <> NrCols(mat) then
 Error( "<mat> must be a square matrix" );
 fi;

 MM := List( mat, x -> List( x, y -> y mod m ) );
 n := NrRows(MM);

 # construct the identity matrix
 inv := IdentityMat( n, Cyclotomics );
 perm := [];

 # loop over the rows
 for i in [ 1 .. n ] do

 pj := 1;
 while MM[i,pj] = 0 do
 pj := pj + 1;
 if pj > n then
 # <mat> is not invertible mod <m>
 return fail;
 fi;
 od;
 perm[pj] := i;
 elem := MM[i,pj];
 MM[i] := List( MM[i], x -> (x/elem) mod m );
 inv[i] := List( inv[i], x -> (x/elem) mod m );

 liste := [ 1 .. i-1 ];
 Append( liste, [i+1..n] );
 for l in liste do
 elem := MM[l,pj];
 MM[l] := MM[l] - MM[i] * elem;
 MM[l] := List( MM[l], x -> x mod m );
 inv[l] := inv[l] - inv[i] * elem;
 inv[l] := List( inv[l], x -> x mod m );
 od;
 od;
 return List( perm, i->inv[i] );
end );

#############################################################################
##
#M KroneckerProduct( <mat1>, <mat2> )
##
InstallOtherMethod( KroneckerProduct,
 "generic method for matrices",
 IsIdenticalObj,
 [ IsMatrix, IsMatrix ],
function ( mat1, mat2 )
 local i, row1, row2, row, kroneckerproduct;
 kroneckerproduct := [];
 for row1 in mat1 do
 for row2 in mat2 do
 row := [];
 for i in row1 do
 Append( row, i * row2 );
 od;
#T application of the new 'AddRowVector' function?
 ConvertToVectorRepNC( row );
 Add( kroneckerproduct, row );
 od;
 od;
 return kroneckerproduct;
end );

## symmetric and alternating parts of the kronecker product
## code is due to Jean Michel

#############################################################################
##
#M ExteriorPower( <mat>, <m> )
##
InstallOtherMethod(ExteriorPower,
 "for matrices", true,[IsMatrix,IsPosInt],
function ( A, m )
local basis;
 basis := Combinations( [ 1 .. NrRows( A ) ], m );
 return List( basis, i->List( basis, j->DeterminantMat( A{i}{j} )));
end);

#############################################################################
##
#M SymmetricPower( <mat>, <m> )
##
InstallOtherMethod(SymmetricPower,
 "for matrices", true,[IsMatrix,IsPosInt],
function ( A, m )
local basis, f;
 f := j->Product( List( Collected( j ), x->x[2]), Factorial );
 basis := UnorderedTuples( [ 1 .. NrRows( A ) ], m );
 return List( basis, i-> List( basis, j->Permanent( A{i}{j}) / f( i )));
#T This code makes sense only in characteristic zero.
end);

#############################################################################
##
#M SolutionMat( <mat>, <vec> ) . . . . . . . . . . one solution of equation
##
## One solution <x> of <x> * <mat> = <vec> or `fail'.
##
InstallOtherMethod( SolutionMatDestructive,
 "generic method",
 IsCollsElms,
 [ IsMatrixOrMatrixObj and IsMutable,
 IsListOrCollection and IsMutable],
 function( mat, vec )
 local i,ncols,sem, vno, z,x, sol;
 ncols := Length(vec);
 z := ZeroOfBaseDomain(mat);
 sol := ListWithIdenticalEntries(NrRows(mat),z);
 ConvertToVectorRepNC(sol);
 if ncols <> NrCols(mat) then
 Error("SolutionMat: matrix and vector incompatible");
 fi;
 sem := SemiEchelonMatTransformationDestructive(mat);
 for i in [1..ncols] do
 vno := sem.heads[i];
 if vno <> 0 then
 x := vec[i];
 if x <> z then
 AddRowVector(vec, sem.vectors[vno], -x);
 AddRowVector(sol, sem.coeffs[vno], x);
 fi;
 fi;
 od;
 if IsZero(vec) then
 return sol;
 else
 return fail;
 fi;
end);

#InstallMethod( SolutionMatNoCo,
# "generic method for ordinary matrix and vector",
# IsCollsElms,
# [ IsOrdinaryMatrix,
# IsRowVector ],
#function ( mat, vec )
# local h, v, tmp, i, l, r, s, c, zero;
#
# # solve <mat> * x = <vec>.
# vec := ShallowCopy( vec );
# l := NrRows( mat );
# r := 0;
# zero := ZeroOfBaseDomain( mat );
# Info( InfoMatrix, 1, "SolutionMat called" );
#
# # Run through all columns of the matrix.
# c := 1;
# while c <= NrCols( mat ) and r < l do
#
# # Find a nonzero entry in this column.
# s := r + 1;
# while s <= l and mat[ s , c ] = zero do s := s + 1; od;
#
# # If there is a nonzero entry,
# if s <= l then
#
# # increment the rank.
# Info( InfoMatrix, 2, " nonzero columns ", c );
# r := r + 1;
#
# # Make its row the current row and normalize it.
# tmp := mat[ s , c ] ^ -1;
# v := mat[ s ]; mat[ s ] := mat[ r ]; mat[ r ] := tmp * v;
# v := vec[ s ]; vec[ s ] := vec[ r ]; vec[ r ] := tmp * v;
#
# # Clear all entries in this column.
# for s in [ 1 .. NrRows( mat ) ] do
# if s <> r and mat[ s , c ] <> zero then
# tmp := mat[ s , c ];
# mat[ s ] := mat[ s ] - tmp * mat[ r ];
# vec[ s ] := vec[ s ] - tmp * vec[ r ];
# fi;
# od;
# fi;
# c := c + 1;
# od;
#
# # Find a solution.
# for i in [ r + 1 .. l ] do
# if vec[ i ] <> zero then return fail; fi;
# od;
# h := [];
# s := NrCols( mat );
# v := Zero( mat[ 1 , 1 ] );
# r := 1;
# c := 1;
# while c <= s and r <= l do
# while c <= s and mat[ r , c ] = zero do
# c := c + 1;
# Add( h, v );
# od;
# if c <= s then
# Add( h, vec[ r ] );
# r := r + 1;
# c := c + 1;
# fi;
# od;
# while c <= s do Add( h, v ); c := c + 1; od;
#
# Info( InfoMatrix, 1, "SolutionMat returns" );
# return h;
#end );
#

InstallOtherMethod( SolutionMat,
 "generic method for ordinary matrix and vector",
 IsCollsElms,
 [ IsMatrixOrMatrixObj,IsListOrCollection ],
 function ( mat, vec )
 return SolutionMatDestructive( MutableCopyMatrix( mat ),
 ShallowCopy( vec ) );
end );

#InstallMethod( SolutionMatDestructive,
# "generic method for ordinary matrix and vector",
# IsCollsElms,
# [ IsOrdinaryMatrix,
# IsRowVector ],
# function ( mat, vec )
# return SolutionMatNoCo( MutableTransposedMatDestructive( mat ),
# vec );
#end );

############################################################################
##
#M SumIntersectionMat( <M1>, <M2> ) . . sum and intersection of two spaces
##
## performs Zassenhaus' algorithm to compute bases for the sum and the
## intersection of spaces generated by the rows of the matrices <M1>, <M2>.
##
## returns a list of length 2, at first position a base of the sum, at
## second position a base of the intersection. Both bases are in
## semi-echelon form.
##
InstallMethod( SumIntersectionMat,"ordinary matrices",
 IsIdenticalObj,
 [ IsMatrix, IsMatrix ],
function( M1, M2 )
 local n, # number of columns
 mat, # matrix for Zassenhaus algorithm
 zero, # zero vector
 v, # loop over 'M1' and 'M2'
 heads, # list of leading positions
 sum, # base of the sum
 i, # loop over rows of 'mat'
 int; # base of the intersection

 if Length( M1 ) = 0 then
 return [ M2, M1 ];
 elif Length( M2 ) = 0 then
 return [ M1, M2 ];
 elif NrCols( M1 ) <> NrCols( M2 ) then
 Error( "dimensions of matrices are not compatible" );
 elif ZeroOfBaseDomain( M1 ) <> ZeroOfBaseDomain( M2 ) then
 Error( "fields of matrices are not compatible" );
 fi;

 n:= Length( M1[1] );
 mat:= [];
 zero:= Zero( M1[1] );

 # Set up the matrix for Zassenhaus' algorithm.
 mat:= [];
 for v in M1 do
 v:= ShallowCopy( v );
 Append( v, v );
 ConvertToVectorRepNC( v );
 Add( mat, v );
 od;
 for v in M2 do
 v:= ShallowCopy( v );
 Append( v, zero );
 ConvertToVectorRepNC( v );
 Add( mat, v );
 od;

 # Transform `mat' into semi-echelon form.
 mat := SemiEchelonMatDestructive( mat );
 heads := mat.heads;
 mat := mat.vectors;

 # Extract the bases for the sum \ldots
 sum:= [];
 for i in [ 1 .. n ] do
 if heads[i] <> 0 then
 Add( sum, mat[ heads[i] ]{ [ 1 .. n ] } );
 fi;
 od;

 # \ldots and the intersection.
 int:= [];
 for i in [ n+1 .. Length( heads ) ] do
 if heads[i] <> 0 then
 Add( int, mat[ heads[i] ]{ [ n+1 .. 2*n ] } );
 fi;
 od;

 # return the result
 return [ sum, int ];
end );

InstallOtherMethod( SumIntersectionMat,"MatrixObject", IsIdenticalObj,
 [ IsMatrixObj, IsMatrixObj ],
function( M1, M2 )
 local n, # number of columns
 mat, # matrix for Zassenhaus algorithm
 v, # loop over 'M1' and 'M2'
 heads, # list of leading positions
 sum, # base of the sum
 i, # loop over rows of 'mat'
 int; # base of the intersection

 if NrRows( M1 ) = 0 then
 return [ M2, M1 ];
 elif NrRows( M2 ) = 0 then
 return [ M1, M2 ];
 elif NrCols( M1 ) <> NrCols( M2 ) then
 Error( "dimensions of matrices are not compatible" );
 elif ZeroOfBaseDomain( M1 ) <> ZeroOfBaseDomain( M2 ) then
 Error( "fields of matrices are not compatible" );
 fi;

 n:= NrCols( M1 );
 mat:= [];

 # Set up the matrix for Zassenhaus' algorithm.
 mat:= ZeroMatrix(BaseDomain(M1),NrRows(M1)+NrRows(M2),2*n);

 i:=1;
 for v in List(M1) do
 CopySubVector(v,mat[i],[1..n],[1..n]);
 CopySubVector(v,mat[i],[1..n],[n+1..2*n]);
 #v:= ShallowCopy( v );
 #Append( v, v );
 i:=i+1;
 od;
 for v in List(M2) do
 CopySubVector(v,mat[i],[1..n],[1..n]);
 #mat[i]{[n+1..2*n]}:=zero; # not needed, as initially 0
 #v:= ShallowCopy( v );
 #Append( v, zero );
 i:=i+1;
 od;

 # Transform `mat' into semi-echelon form.
 mat := SemiEchelonMatDestructive( mat );
 heads := mat.heads;
 mat := mat.vectors;

 # Extract the bases for the sum \ldots
 sum:= [];
 for i in [ 1 .. n ] do
 if heads[i] <> 0 then
 Add( sum, mat[ heads[i] ]{ [ 1 .. n ] } );
 fi;
 od;

 # \ldots and the intersection.
 int:= [];
 for i in [ n+1 .. Length( heads ) ] do
 if heads[i] <> 0 then
 Add( int, mat[ heads[i] ]{ [ n+1 .. 2*n ] } );
 fi;
 od;

 # return the result
 return [ sum, int ];
end );

#############################################################################
#end
#M TriangulizeMat( <mat> ) . . . . . bring a matrix in upper triangular form
##

BindGlobal( "TRIANGULIZE_MAT_GENERIC", function ( mat, m, n, zero )
local i, j, k, row, x, row2;

 Info( InfoMatrix, 1, "TriangulizeMat called" );

 if m>0 then

 # make sure that the rows are mutable
 for i in [ 1 .. m ] do
 if not IsMutable( mat[i] ) then
 mat[i]:= ShallowCopy( mat[i] );
 fi;
 od;

 # run through all columns of the matrix
 i := 0;
 for k in [1..n] do

 # find a nonzero entry in this column
 j := i + 1;
 while j <= m and mat[j,k] = zero do j := j + 1; od;

 # if there is a nonzero entry
 if j <= m then

 # increment the rank
 Info( InfoMatrix, 2, " nonzero columns: ", k );
 i := i + 1;

 # make its row the current row and normalize it
 row := mat[j];
 mat[j] := mat[i];
 x:= Inverse( row[k] );
 if x = fail then
 TryNextMethod();
 fi;
 MultVector( row, x );
 mat[i] := row;

 # clear all entries in this column
 for j in [1..i-1] do
 row2 := mat[j];
 x := row2[k];
 if x <> zero then
 AddRowVector( row2, row, - x );
 fi;
 od;
 for j in [i+1..m] do
 row2 := mat[j];
 x := row2[k];
 if x <> zero then
 AddRowVector( row2, row, - x );
 fi;
 od;

 fi;

 od;

 fi;

 Info( InfoMatrix, 1, "TriangulizeMat returns" );
end );

InstallMethod( TriangulizeMat,
 "generic method for mutable matrices",
 [ IsMatrix and IsMutable ],
function(mat)
 TRIANGULIZE_MAT_GENERIC(mat,NrRows(mat),NrCols(mat),
 ZeroOfBaseDomain(mat));
end);

InstallOtherMethod( TriangulizeMat,
 "generic method for mutable matrix obj",
 [ IsMatrixObj and IsMutable ],
function(mat)
 TRIANGULIZE_MAT_GENERIC(mat,NrRows(mat),NrCols(mat),
 ZeroOfBaseDomain(mat));
end);

InstallOtherMethod( TriangulizeMat,
 "list of mutable vectors",
 [ IsList and IsMutable ],
function(mat)
 if Length(mat)=0 or not ForAll(mat,x->IsRowVector(x) or IsVectorObj(x))
 then TryNextMethod();fi;
 TRIANGULIZE_MAT_GENERIC(mat,Length(mat),Length(mat[1]),
 ZeroOfBaseDomain(mat[1]));
end);

InstallOtherMethod( TriangulizeMat,
 "for an empty list",
 [ IsList and IsEmpty],
 function( m ) return; end );

InstallMethod( TriangulizedMat, "generic method for matrices", [ IsMatrix ],
function ( mat )
local m;
 m:=List(mat,ShallowCopy);
 TriangulizeMat(m);
 return m;
end);

InstallOtherMethod( TriangulizedMat, "matrix objects", [IsMatrixObj],
function ( mat )
local m;
 m:=MutableCopyMatrix(mat);
 TriangulizeMat(m);
 return m;
end);

InstallOtherMethod( TriangulizedMat,"list of vectors", [IsList],
function ( mat )
local m;
 m:=MutableCopyMatrix(mat);
 TriangulizeMat(m);
 return m;
end);

InstallOtherMethod( TriangulizedMat, "for an empty list", [ IsList and IsEmpty ],
mat -> []);

#############################################################################
##
#M UpperSubdiagonal( <mat>, <pos> )
##
InstallMethod( UpperSubdiagonal,
 [ IsMatrixOrMatrixObj, IsPosInt ],
function( mat, l )
 return List( [ 1 .. Minimum( NrRows( mat ), NrCols( mat ) - l ) ],
 i -> mat[ i, l+i ] );
end );

#############################################################################
##
#F BaseFixedSpace( <mats> ) . . . . . . . . . . . . calculate fixed points
##
## 'BaseFixedSpace' returns a base of the vector space $V$ such that
## $M v = v$ for all $v$ in $V$ and all matrices $M$ in the list <mats>.
##
InstallGlobalFunction( BaseFixedSpace, function( matrices )
 local I, # identity matrix
 size, # dimension of vector space
 E, # linear system
 M, # one matrix of 'matrices'
 N, # M - I
 j;

 I := matrices[1]^0;
 size := NrRows(I);
 E := List( [ 1 .. size ], x -> [] );
 for M in matrices do
 N := M - I;
 for j in [ 1..size ] do
 Append( E[ j ], N[ j ] );
 od;
 od;
 return NullspaceMatDestructive( E );
end );

##########################################################################
##
#F BaseSteinitzVectors( <bas>, <mat> )
##
## find vectors extending mat to a basis spanning the span of <bas>.
## 'BaseSteinitz' returns a
## record describing a base for the factorspace and ways to decompose
## vectors:
##
## zero: zero of <V> and 
## factorzero: zero of complement
## subspace: triangulized basis of <mat>
## factorspace: base of a complement of in <V>
## heads: a list of integers i_j, such that if i_j>0 then a vector
## with head j is at position i_j in factorspace. If i_j<0
## then the vector is in subspace.
##
InstallGlobalFunction( BaseSteinitzVectors, function(bas,mat)
local z,l,b,i,j,k,stop,v,dim,h,zv;

 # catch trivial case
 if Length(bas)=0 then
 return rec(subspace:=[],factorspace:=[]);
 fi;

 z:=Zero(bas[1][1]);
 zv:=Zero(bas[1]);
 if Length(mat)>0 then
 mat:=MutableCopyMatrix(mat);
 TriangulizeMat(mat);
 fi;
 bas:=MutableCopyMatrix(bas);
 dim:=Length(bas[1]);
 l:=Length(bas)-Length(mat); # missing dimension
 b:=[];
 h:=[];
 i:=1;
 j:=1;
 while Length(b)<l do
 stop:=false;
 repeat
 if j<=dim and (Length(mat)k[j]<>z);
 if v<>fail then
 # add the vector
 v:=bas[v];
 v:=1/v[j]*v; # normed
 Add(b,v);
 h[j]:=Length(b);
 # if fail, then this dimension is only dependent (and not needed)
 fi;
 else
 stop:=true;
 # check whether we are running to fake zero columns
 if i<=Length(mat) then
 # has a step, clean with basis vector
 v:=mat[i];
 v:=1/v[j]*v; # normed
 h[j]:=-i;
 else
 v:=fail;
 fi;
 fi;
 if v<>fail then
 # clean j-th component from bas with v
 for k in [1..Length(bas)] do
 if not IsZero(bas[k][j]) then
 bas[k]:=bas[k]-bas[k][j]/v[j]*v;
 fi;
 od;
 v:=Zero(v);
 bas:=Filtered(bas,k->k<>v);
 fi;
 j:=j+1;
 until stop;
 i:=i+1;
 od;
 # add subspace indices
 while i<=Length(mat) do
 if mat[i,j]<>z then
 h[j]:=-i;
 i:=i+1;
 fi;
 j:=j+1;
 od;
 return rec(factorspace:=b,
 factorzero:=zv,
 subspace:=mat,
 heads:=h);
end );

#############################################################################
##
#F BlownUpMat( , <mat> )
##
InstallGlobalFunction( BlownUpMat, function ( B, mat )
 local result, # blown up matrix, result
 vectors, # basis vectors of 'B'
 row, # loop over rows of 'mat'
 b, # loop over 'vectors'
 resrow, # one row of 'result'
 entry; # loop over 'row'

 vectors:= BasisVectors( B );
 result:= [];
 for row in mat do
 for b in vectors do
 resrow:= [];
 for entry in row do
 entry := Coefficients( B, entry * b );
 if entry = fail then return fail; fi;
 Append( resrow, entry );
 od;
 ConvertToVectorRepNC( resrow );
 Add( result, resrow );
 od;
 od;

 # Return the result.
 return result;
end );

#############################################################################
##
#F BlownUpVector( , <vector> )
##
InstallGlobalFunction( BlownUpVector, function ( B, vector )
 local result, # blown up vector, result
 entry; # loop over 'vector'

 result:= [];
 for entry in vector do
 Append( result, Coefficients( B, entry ) );
 od;
 ConvertToVectorRepNC( result );

 # Return the result.
 return result;
end );

#############################################################################
##
#F IdentityMat( <m>[, <F>] ) . . . . . . . . identity matrix of a given size
##
InstallGlobalFunction( IdentityMat, function ( arg )
 local id, m, zero, one, row, i, f;

 # check the arguments and get dimension, zero and one
 if Length(arg) = 1 then
 m := arg[1];
 zero := 0;
 one := 1;
 f := Rationals;
 elif Length(arg) = 2 and IsRing(arg[2]) then
 m := arg[1];
 zero := Zero( arg[2] );
 one := One( arg[2] );
 f := arg[2];
 if one = fail then
 Error( "ring must be a ring-with-one" );
 fi;
 elif Length(arg) = 2 then
 m := arg[1];
 zero := Zero( arg[2] );
 one := One( arg[2] );
 f := Ring( one, arg[2] );
 else
 Error("usage: IdentityMat( <m>[, <R>] )");
 fi;

 # special treatment for 0-dimensional spaces
 if m=0 then
 return NullMapMatrix;
 fi;

 # make an empty row
 row := ListWithIdenticalEntries(m,zero);
 ConvertToVectorRepNC(row,f);

 # make the identity matrix
 id := [];
 for i in [1..m] do
 id[i] := ShallowCopy( row );
 id[i,i] := one;
 od;

 # We do *not* call ConvertToMatrixRep here, as that can cause
 # unexpected problems for the user (e.g. if a matrix over GF(2) is
 # created, and the user then tries to change an entry to Z(4),

 return id;
end );

#############################################################################
##
#F NullMat( <m>, <n> [, <F>] ) . . . . . . . . . null matrix of a given size
##
InstallGlobalFunction( NullMat, function ( arg )
 local null, m, n, zero, row, i, f;

 if Length(arg) = 2 then
 m := arg[1];
 n := arg[2];
 f := Rationals;
 elif Length(arg) = 3 and IsRing(arg[3]) then
 m := arg[1];
 n := arg[2];
 f := arg[3];
 elif Length(arg) = 3 then
 m := arg[1];
 n := arg[2];
 f := Ring(One(arg[3]), arg[3]);
 else
 Error("usage: NullMat( <m>, <n> [, <R>] )");
 fi;
 zero := Zero(f);

# # special treatment for 0-dimensional spaces
# if m = 0 or n = 0 then
# return NullMapMatrix;
# fi;
#T Adding this would break a test.

 # make an empty row
 row := ListWithIdenticalEntries(n,zero);
 ConvertToVectorRepNC( row, f );

 # make the null matrix
 null := [];
 for i in [1..m] do
 null[i] := ShallowCopy( row );
 od;

 # We do *not* call ConvertToMatrixRep here, as that can cause
 # unexpected problems for the user (e.g. if a matrix over GF(2) is
 # created, and the user then tries to change an entry to Z(4),

 return null;
end );

#############################################################################
##
#F NullspaceModN( <M>, <n> ) . . . . . . . . . . . nullspace of <M> mod <n>
##
InstallGlobalFunction( NullspaceModN, function( M, n )
 local B, coeffs;

 B := BasisNullspaceModN(M, n);
 coeffs := Cartesian(ListWithIdenticalEntries(Length(B), [0..n-1]));
 return Set(coeffs, c -> c * B mod n);
end );

#############################################################################
##
#F BasisNullspaceModN( <M>, <n> ) . . basis of the nullspace of <M> mod <n>
##
InstallGlobalFunction( BasisNullspaceModN, function( M, n )
 local snf, null, nullM, i;

 # if n is a prime, Gaussian elimination is fastest
 if IsPrimeInt (n) then
 return List (NullspaceMat (M*One(GF(n))),
 v -> List (v, IntFFE));
 fi;

 # compute the Smith normal form S for M, i.e., S = R M C
 snf := NormalFormIntMat (M, 1+4);

 # compute the nullspace of S mod n
 null := IdentityMat (Length (M));

 for i in [1..snf.rank] do
 null[i,i] := n/GcdInt (n, snf.normal[i,i]);
 od;

 # nullM = null*R is the nullspace of M C mod n
 # since solutions do not change under elementary matrices
 # nullM is also the nullspace for M

 nullM := null*snf.rowtrans mod n;
 Assert (1, ForAll (nullM, v -> v*M mod n =0*M[1]));
 return nullM;
end);

#############################################################################
##
#F PermutationMat( <perm>, <dim> [, <F> ] ) . . . . . . permutation matrix
##
InstallGlobalFunction( PermutationMat, function( arg )
 local i, # loop variable
 perm, # permutation
 dim, # desired dimension of the permutation matrix
 F, # field of the matrix entries (defauled to 'Rationals')
 mat; # matrix corresponding to 'perm', result

 if not ( ( Length( arg ) = 2 or Length( arg ) = 3 )
 and IsPerm( arg[1] ) and IsInt( arg[2] ) ) then
 Error( "usage: PermutationMat( <perm>, <dim> [, <F> ] )" );
 fi;

 perm:= arg[1];
 dim:= arg[2];
 if Length( arg ) = 2 then
 F:= Rationals;
 else
 F:= arg[3];
 fi;

 mat:= NullMat( dim, dim, F );

 for i in [ 1 .. dim ] do
 mat[i, i^perm ]:= One( F );
 od;

 return mat;
end );

#############################################################################
##
#M DiagonalMatrix( [<filt>, ]<R>, <vector> )
#M DiagonalMatrix( <vector>[, <M>] )
##
InstallTagBasedMethod( DiagonalMatrix,
 function( filt, R, vector )
 local n, mat, i;

 n:= Length( vector );
 mat:= NewZeroMatrix( filt, R, n, n );
 for i in [ 1 .. n ] do
 mat[i, i]:= vector[i];
 od;

 return mat;
 end );

InstallMethod( DiagonalMatrix,
 [ "IsSemiring", "IsRowVectorOrVectorObj" ],
 function( R, vector )
 return DiagonalMatrix( DefaultMatrixRepForBaseDomain( R ), R, vector );
 end );

InstallMethod( DiagonalMatrix,
 [ "IsRowVectorOrVectorObj", "IsMatrixOrMatrixObj" ],
 function( vector, M )
 return DiagonalMatrix( ConstructingFilter( M ), BaseDomain( M ), vector );
 end );

InstallMethod( DiagonalMatrix,
 [ "IsRowVectorOrVectorObj" ],
 function( vector )
 local R;

 if Length( vector ) = 0 then
 Error( "do not know over which ring the matrix shall be defined" );
 fi;
 R:= DefaultScalarDomainOfMatrixList( [ [ vector ] ] );
 return DiagonalMatrix( DefaultMatrixRepForBaseDomain( R ), R, vector );
 end );
#T Alternatively, we could use the code of 'DiagonalMat' in this method,
#T turn 'DiagonalMat' into an obsolescent synonym of 'DiagonalMatrix',
#T and turn the documentation of 'DiagonalMat' into a link to the
#T documentation of 'DiagonalMatrix'.
#T The behaviour is slightly different for example if the argument is a list
#T of integers or rationals.

#############################################################################
##
#F DiagonalMat( <vector> )
##
InstallGlobalFunction( DiagonalMat, function( vector )
 local zerovec,
 M,
 i;

 M:= [];
 zerovec:= Zero( vector[1] );
 zerovec:= List( vector, x -> zerovec );

 for i in [ 1 .. Length( vector ) ] do
 M[i]:= ShallowCopy( zerovec );
 M[i,i]:= vector[i];
 ConvertToVectorRepNC(M[i]);
 od;

 # We do *not* call ConvertToMatrixRep here, as that can cause
 # unexpected problems for the user (e.g. if a matrix over GF(2) is
 # created, and the user then tries to change an entry to Z(4),

 return M;
end );

#############################################################################
##
#F ReflectionMat( <coeffs> )
#F ReflectionMat( <coeffs>, <root> )
#F ReflectionMat( <coeffs>, <conj> )
#F ReflectionMat( <coeffs>, <conj>, <root> )
##
InstallGlobalFunction( ReflectionMat, function( arg )
 local coeffs, # coefficients vector, first argument
 w, # root of unity, second argument (optional)
 conj, # conjugation function, third argument (optional)
 M, # matrix of the reflection, result
 n, # length of 'coeffs'
 one, # identity of the ring over that 'coeffs' is written
 c, # coefficient of 'M'
 i, # loop over rows of 'M'
 j, # loop over columns of 'M'
 row; # one row of 'M'

 # Get and check the arguments.
 if Length( arg ) < 1 or 3 < Length( arg )
 or not IsList( arg[1] ) then
 Error( "usage: ReflectionMat( <coeffs> [, <conj> ] [, <k> ] )" );
 fi;
 coeffs:= arg[1];
 if Length( arg ) = 1 then
 w:= -1;
 conj:= List( coeffs, ComplexConjugate );
 elif Length( arg ) = 2 then
 if not IsFunction( arg[2] ) then
 w:= arg[2];
 conj:= List( coeffs, ComplexConjugate );
 else
 w:= -1;
 conj:= arg[2];
 if not IsFunction( conj ) then
 Error( "<conj> must be a function" );
 fi;
 conj:= List( coeffs, conj );
 fi;
 elif Length( arg ) = 3 then
 conj:= arg[2];
 if not IsFunction( conj ) then
 Error( "<conj> must be a function" );
 fi;
 conj:= List( coeffs, conj );
 w:= arg[3];
 fi;

 # Construct the matrix.
 M:= [];
 one:= coeffs[1] ^ 0;
 w:= w * one;
 n:= Length( coeffs );
 c:= ( w - one ) / ( coeffs * conj );
 for i in [ 1 .. n ] do
 row:= [];
 for j in [ 1 .. n ] do
 row[j]:= conj[i] * c * coeffs[j];
 od;
 row[i]:= row[i] + one;
 ConvertToVectorRepNC( row );
 M[i]:= row;
 od;

 # We do *not* call ConvertToMatrixRep here, as that can cause
 # unexpected problems for the user (e.g. if a matrix over GF(2) is
 # created, and the user then tries to change an entry to Z(4),

 return M;
end );

#########################################################################
##
#F RandomInvertibleMat( [rs ,] <m> [, <R>] ) . . . make a random invertible matrix
##
## 'RandomInvertibleMat' returns a invertible random matrix with <m> rows
## and columns with elements taken from the ring <R>, which defaults to
## 'Integers'.
##
InstallGlobalFunction( RandomInvertibleMat, function ( arg )
 local rs, mat, m, R, i, row, k;

 if Length(arg) >= 1 and IsRandomSource(arg[1]) then
 rs := Remove(arg, 1);
 else
 rs := GlobalMersenneTwister;
 fi;

 # check the arguments and get the list of elements
 if Length(arg) = 1 then
 m := arg[1];
 R := Integers;
 elif Length(arg) = 2 then
 m := arg[1];
 R := arg[2];
 else
 Error("usage: RandomInvertibleMat( [rs ,] <m> [, <R>] )");
 fi;

 # now construct the random matrix
 mat := [];
 repeat
 for i in [1..m] do
 row := [];
 for k in [1..m] do
 row[k] := Random( rs, R );
 od;
 ConvertToVectorRepNC( row, R );
 mat[i] := row;
 od;
 until RankMat( mat ) = m;

 # We do *not* call ConvertToMatrixRep here, as that can cause
 # unexpected problems for the user (e.g. if a matrix over GF(2) is
 # created, and the user then tries to change an entry to Z(4),

 return mat;
end );

#############################################################################
##
#F RandomMat( [rs ,] <m>, <n> [, <R>] ) . . . . . . . . make a random matrix
##
## 'RandomMat' returns a random matrix with <m> rows and <n> columns with
## elements taken from the ring <R>, which defaults to 'Integers'.
##
InstallGlobalFunction( RandomMat, function ( arg )
 local rs, mat, m, n, R, i, row, k;

 if Length(arg) >= 1 and IsRandomSource(arg[1]) then
 rs := Remove(arg, 1);
 else
 rs := GlobalMersenneTwister;
 fi;

 # check the arguments and get the list of elements
 if Length(arg) = 2 then
 m := arg[1];
 n := arg[2];
 R := Integers;
 elif Length(arg) = 3 then
 m := arg[1];
 n := arg[2];
 R := arg[3];
 else
 Error("usage: RandomMat( [rs ,] <m>, <n> [, <F>] )");
 fi;

 # now construct the random matrix
 mat := [];
 for i in [1..m] do
 row := [];
 for k in [1..n] do
 row[k] := Random( rs, R );
 od;
 ConvertToVectorRepNC( row, R );
 mat[i] := row;
 od;

 # We do *not* call ConvertToMatrixRep here, as that can cause
 # unexpected problems for the user (e.g. if a matrix over GF(2) is
 # created, and the user then tries to change an entry to Z(4),

 return mat;
end );

#############################################################################
##
#F RandomUnimodularMat( [rs ,] <m> ) . . . . . . . random unimodular matrix
##
InstallGlobalFunction( RandomUnimodularMat, function ( arg )
local rs, m, mat, c, i, j, k, l, a, b, v, w, gcd,dom;

 dom:=ValueOption("domain");
 if dom=fail or dom=false then
 dom:=Integers;
 fi;

 if Length(arg) >= 1 and IsRandomSource(arg[1]) then
 rs := Remove(arg, 1);
 else
 rs := GlobalMersenneTwister;
 fi;

 if Length(arg) = 1 then
 m := arg[1];
 else
 Error("usage: RandomUnimodularMat( [rs ,] <m> )");
 fi;

 if not IsPosInt( m ) then
 Error("<m> must be a positive integer");
 fi;

 # start with the identity matrix
 mat := IdentityMat( m );

 if m = 1 then
 return Random( rs, [ 1, -1 ] ) * mat;
 fi;

 for c in [1..m] do

 # multiply two random rows with a random? unimodular 2x2 matrix
 i := Random(rs, 1, m);
 repeat
 j := Random(rs, 1, m);
 until j <> i;
 repeat
 a := Random( rs, dom );
 b := Random( rs, dom );
 gcd := Gcdex( a, b );
 until gcd.gcd = 1;
 v := mat[i]; w := mat[j];
 mat[i] := ShallowCopy(gcd.coeff1 * v + gcd.coeff2 * w);
 mat[j] := ShallowCopy(gcd.coeff3 * v + gcd.coeff4 * w);

 # multiply two random cols with a random? unimodular 2x2 matrix
 k := Random(rs, 1, m);
 repeat
 l := Random(rs, 1, m);
 until l <> k;
 repeat
 a := Random( rs, dom );
 b := Random( rs, dom );
 gcd := Gcdex( a, b );
 until gcd.gcd = 1;
 for i in [1..m] do
 v := mat[i,k]; w := mat[i,l];
 mat[i,k] := gcd.coeff1 * v + gcd.coeff2 * w;
 mat[i,l] := gcd.coeff3 * v + gcd.coeff4 * w;
 ConvertToVectorRepNC( mat[i] );
 od;

 od;

 return mat;
end );

#############################################################################
##
#F SimultaneousEigenvalues( <matlist>, <expo> ) . . . . . . . . .eigenvalues
##
## The matgroup generated by <matlist> must be an abelian p-group of
## exponent <expo>. The matrices in matlist must be matrices over GF(<q>)
## for some prime <q>. Then the eigenvalues of <mat> in the splitting field
## GF(<q>^r) for some r are powers of an element ksi in the splitting field,
## which is of order <expo>.
##
## 'SimultaneousEigenspaces' returns a matrix of integers mod <expo>, say
## (a_{i,j}), such that the power ksi^a_{i,j} is an eigenvalue of the i-th
## matrix in <matlist> and the eigenspaces of the different matrices to the
## eigenvalues ksi^a_{i,j} for fixed j are equal.
##
InstallGlobalFunction( SimultaneousEigenvalues,
 function( arg )
 local matlist, expo,
 q, # characteristic of field of matrices
 r, # such that <q>^r is splitting field
 field, # GF(<q>^r)
 ksi, # <expo>-root of field
 eival, # exponents of eigenvalues of the matrices
 eispa, # eigenspaces of the matrices
 eigen, # exponents of simultaneous eigenvalues
 I, # identity matrix
 w, # ksi^w is candidate for a eigenvalue
 null, # basis of nullspace
 i, Split;

 Split := function( space, i )
 local int, # intersection of two row spaces
 j;

 for j in [1..Length(eival[i])] do
 if 0 < Length( eispa[i][j] ) then
 int := SumIntersectionMat( space, eispa[i][j] )[2];
 if 0 < Length( int ) then
 Append( eigen[i],
 List( int, x -> eival[i][j] ) );
 if i < Length( matlist ) then
 Split( int, i+1 );
 fi;
 fi;
 fi;
 od;
 end;

 matlist := arg[1];
 expo := arg[2];

 # compute ksi
 if Length( arg ) = 2 then
 q := Characteristic( matlist[1][1,1] );

 # get splitting field
 r := 1;
 while EuclideanRemainder( q^r - 1, expo ) <> 0 do
 r := r+1;
 od;
 field := GF(q^r);
 ksi := GeneratorsOfField(field)[1]^((q^r - 1) / expo);
 else
 ksi := arg[3];
 fi;

 # set up eigenvalues and spaces and Idmat
 eival := List( matlist, x -> [] );
 eispa := List( matlist, x -> [] );
 I := matlist[1]^0;

 # calculate eigenvalues and spaces for each matrix
 for i in [1..Length(matlist)] do
 for w in [0..expo-1] do
 null := NullspaceMat( matlist[i] - (ksi^w * I) );
 if 0 < Length(null) then
 Add( eival[i], w );
 Add( eispa[i], null );
 fi;
 od;
 od;

 # now make the eigenvalues simultaneous
 eigen := List( matlist, x -> [] );
 for i in [1..Length(eival[1])] do
 Append( eigen[1], List( eispa[1][i], x -> eival[1][i] ) );
 if Length( matlist ) > 1 then
 Split( eispa[1][i], 2 );
 fi;
 od;

 # return
 return eigen;
end );

#############################################################################
##
#F FlatBlockMat( <blockmat> ) . . . . . . . . . . . convert block mat to mat
##
InstallGlobalFunction( FlatBlockMat, function( block )
 local d, l, mat, i, j, h, k, a, b;

 d := Length( block );
 l := Length( block[1][1] );
 mat := List( [1..d*l], x -> List( [1..d*l], y -> false ) );
 for i in [1..d] do
 for j in [1..d] do
 for h in [1..l] do
 for k in [1..l] do
 a := (i-1)*l + h;
 b := (j-1)*l + k;
 mat[a][b] := block[i][j][h][k];
 od;
 od;
 od;
 od;
 return mat;
end );

#############################################################################
##
#F DirectSumMat( <matlist> ) . . . . . . . . . . . create block diagonal mat
#F DirectSumMat( <mat1>, ..., <matn> ) . . . . . . create block diagonal mat
##
InstallGlobalFunction( DirectSumMat, function (arg)
 local m, F, res, B, r, c;

 if Length( arg ) = 1 and not IsMatrixOrMatrixObj( arg[1] ) then
 # We expect that 'arg[1]' is a list of matrices or matrix objects.
 arg:= arg[1];
 fi;

 # Distinguish:
 # - If we are given a list of 'IsMatrix' that are *not* 'IsMatrixObj'
 # or of empty plain lists, we want to return an 'IsMatrix',
 # and we do not care about base domains.
 # - If we are given a list of 'IsMatrixObj',
 # we want to return an 'IsMatrixObj',
 # and eventually we want to require that all given inputs have the
 # same base domain and the same internal representation.
 # For backwards compatibility reasons,
 # we are currently on the one hand less restrictive,
 # because we admit different base domains,
 # and on the other hand we require the base domains to be fields.
 # (This situation is not at all satisfactory.)
 arg:= Filtered( arg, x -> x <> [] );
 if Length( arg ) = 0 then
 return [];
 elif ForAll( arg, m -> IsMatrix( m ) and not IsMatrixObj( m ) ) then
 m:= arg[1];
 F:= DefaultField( m[1,1] );
 res:= NullMat( Sum( arg, NrRows ), Sum( arg, NrCols ), F );
 elif ForAny( arg, IsMatrixObj ) then
 # Find a common 'BaseDomain'.
 F:= fail;
 for m in arg do
 if HasBaseDomain( m ) then
 B:= BaseDomain( m );
 else
 B:= DefaultFieldOfMatrix( m );
 fi;
 if F = fail then
 F:= B;
 elif not IsSubset( F, B ) then
 if IsSubset( B, F ) then
 F:= B;
 else
 F:= Field( Concatenation( GeneratorsOfField( F ),
 GeneratorsOfField( B ) ) );
 fi;
 fi;
 od;
 res:= ZeroMatrix( F, Sum( arg, NrRows ), Sum( arg, NrCols ) );
 else
 Error( "<arg> must be a list of IsMatrix or of IsMatrixObj" );
 fi;

 r:= 0;
 c:= 0;
 for m in arg do
 CopySubMatrix( m, res, [ 1 .. NrRows( m ) ], r + [ 1 .. NrRows( m ) ],
 [ 1 .. NrCols( m ) ], c + [ 1 .. NrCols( m ) ] );
 r:= r + NrRows( m );
 c:= c + NrCols( m );
 od;

 return res;
end );

#############################################################################
##
#M Trace( <mat> ) . . . . . . . . . . . . . . . . . . . . . . for a matrix
##
InstallOtherMethod( Trace,
 "generic method for matrices",
 [ IsMatrix ],
 TraceMat );

#############################################################################
##
#M JordanDecomposition( <mat> )
##
InstallMethod( JordanDecomposition,
 "method for matrices",
 [IsMatrix],
function( mat )

 local F,p,B,f,g,fac,ff,h,w;

# The algorithm is due to R. Beals

 F:= DefaultFieldOfMatrix( mat );
 if F = fail then
 TryNextMethod();
 fi;
 p:= Characteristic( F );

# First we determine a squarefree polynomial 'g' such that 'g^d(mat)=0'.

 f:= CharacteristicPolynomial( F, F, mat );
 if p = 0 or p > Length( mat ) then
 g:= f/Gcd( f, Derivative( f ) );
 else
 fac:= Factors(f);
 g:= One( F );
 ff:= [ ];
 for h in fac do
 if not h in ff then
 g:= g*h;
 Add( ff, h );
 fi;
 od;
 fi;

 if f=g then return [ mat, 0*mat ]; fi;

# Now 'B' will be the semisimple part of the matrix 'mat'.

 w:= GcdRepresentation( g, Derivative( g ) )[2];
 w:= w*g;
 B:= ShallowCopy( mat );
 while Value( g, B ) <> 0*B do
 B:= B - Value( w, B );
 od;

 return [ B, mat-B ];

end );

#############################################################################
##
#F OnSubspacesByCanonicalBasis(<bas>,<mat>)
##
InstallGlobalFunction(OnSubspacesByCanonicalBasis,function( mat, obj )
 local row;
 mat:=mat*obj;
 if not IsMutable(mat) then
 mat := MutableCopyMatrix(mat);
 else
 for row in [1..Length(mat)] do
 if not IsMutable(mat[row]) then
 mat[row] := ShallowCopy(mat[row]);
 fi;
 od;
 fi;
 TriangulizeMat(mat);
 return mat;
end);

#############################################################################
##
#F OnSubspacesByCanonicalBasisConcatenations(<basvec>,<mat>)
##
InstallGlobalFunction(OnSubspacesByCanonicalBasisConcatenations,
function( bvec, obj )
 local n,a,mat,r;
 n:=Length(obj); # acting dimension
 mat:=[];
 a:=1;
 while a<Length(bvec) do
 r:=bvec{[a..a+n-1]}*obj;
 if not IsMutable(r) then r:=ShallowCopy(r);fi;
 Add(mat,r);
 a:=a+n;
 od;
 TriangulizeMat(mat);
 return Concatenation(mat);
end);

#############################################################################
##
#M FieldOfMatrixList
##
InstallMethod(FieldOfMatrixList,"generic: form field",
 [IsListOrCollection],
function(l)
local i,j,k,fg,f;
 # try to find out the field
 if Length(l)=0 or ForAny(l,i->not IsMatrix(i)) then
 Error("<l> must be a list of matrices");
 fi;
 fg:=[l[1][1,1]];
 f:=Field(fg);
 for i in l do
 for j in i do
 for k in j do
 if not k in f then
 Add(fg,k);
 f:=Field(fg);
 fi;
 od;
 od;
 od;
 return f;
end);

#############################################################################
##
#M DefaultScalarDomainOfMatrixList
##
InstallMethod( DefaultScalarDomainOfMatrixList,
 "generic: form ring",
 [ IsList ],
 function(l)
 local B, i,j,k,fg,f,char;
 # try to find out the field
 if Length( l ) = 0 or ForAny( l, i -> not IsMatrixOrMatrixObj( i ) ) then
 Error( "<l> must be a nonempty list of matrices or matrix objects" );
 elif ForAll( l, HasBaseDomain ) then
 B:= BaseDomain( l[1] );
 if ForAll( l, x -> B = BaseDomain( x ) ) then
 return B;
 fi;
 elif ForAll( l, IsMatrix ) then
 fg:=[l[1][1,1]];
 # This fixes Issue #5134 on GitHub.
 # See https://github.com/gap-system/gap/issues/5134
 # The problem is that Characteristic([1.0]) = fail
 # which is catched by the new else case
 # Note that DefaultRing should be defined for each element with
 # a characteristic
 char := Characteristic(fg);
 if char=0 then
 f:=DefaultField(fg);
 elif char = fail then
 Error("DefaultScalarDomainOfMatrixList: Cannot find a base domain for the list entries.");
 else
 f:=DefaultRing(fg);
 fi;
 for i in l do
 for j in i do
 for k in j do
 if not k in f then
 Add(fg,k);
 f:=DefaultRing(fg);
 fi;
 od;
 od;
 od;
 return f;
 fi;
 Error( "all entries in <l> must either be lists of lists ",
 "or have the same BaseDomain" );
end);

#############################################################################
##
#F BaseOrthogonalSpaceMat( <mat> )
##
InstallMethod( BaseOrthogonalSpaceMat,
 "for a matrix",
 [ IsMatrix ],
 mat -> NullspaceMat( TransposedMat( mat ) ) );

# simplex method, code by Ken Monks, AH

#in matrix M, row reduce to get 1s
#in exactly the columns given by
#L a list of indices
BindGlobal("TriangulizeMatPivotColumns",function(M,L)
local idx,i;

 if L=[1..Length(L)] then
 TriangulizeMat(M);
 else
 idx:=Concatenation(L,Filtered([1..NrCols(M)],x->not x in L));
 for i in [1..Length(M)] do M[i]:=M[i]{idx}; od;
 TriangulizeMat(M);
 idx:=ListPerm(PermList(idx)^-1,NrCols(M));
 for i in [1..Length(M)] do M[i]:=M[i]{idx}; od;
 fi;

end);

#inputs a linear form c and maximizes it subject to
#the constraints Ax <= b where all entries of b are nonnegative.
InstallGlobalFunction(SimplexMethod,function(A,b,c)
local M, n, p, vars, slackVars, i, id, bestMove,
 newNonzero, len, ratios, newZero, positiveRatios, point, value,Val;

 Val:=function(M,vars,slackVars,len,x)
 if x in vars then
 return 0;
 else return M[Position(slackVars,x)+1,len];
 fi;
 end;

 #check the size of the data is legit

 n:=Size(c);
 p:=Size(b);
 if not (IsMatrix(A) and Size(A)=p and ForAny(A,R->Size(R)=n)) then
 Error( "usage: SimplexMethod( <A>, , <c>)");
 fi;

 id:=IdentityMat(p,Rationals);

 #build the augmented matrix

 #first row
 M:=[Concatenation([1],-c,List([1..p+1],x->0))];
 #the rest of the rows
 for i in [1..p] do
 Add(M,Concatenation([0],A[i],id[i],[b[i]]));
 od;

 len:=Size(M[1]);

 #initialize the feasible starting vertex
 if ForAll(b,x->not x<0) then
 vars:=[2..2+n-1];
 slackVars:=[2+n..n+p+1];
 else
 return "Invalid data: not all constraints nonnegative!";
 fi;

 #Print("slackVars are ",slackVars ,"\n");
 #Print("vars are ",vars,"\n");
 #Display(M);

 TriangulizeMatPivotColumns(M,Concatenation([1],slackVars));
 #Display(M);
 #bestMove is the coeff var that will become nonzero
 bestMove:=Minimum(List(vars,i->M[1,i]));
 newNonzero:=vars[Position(List(vars,i->M[1,i]),bestMove)];

 #Print(newNonzero, " is the new nonzero guy \n");

 while bestMove<0 do

 #see if figure is unbounded
 #Print("about to do some ratios \n");
 #Print(List([1..p],x-> M[x+1,newNonzero]), " is what we're going to divide by \n");
 ratios:=List([1..p], function(x) if M[x+1,newNonzero]=0 then return infinity; else return M[x+1,len]/M[x+1,newNonzero]; fi; end);
 #Print("done doing some ratios");
 positiveRatios:=Filtered(ratios,x -> x>0);
 if Size(positiveRatios)=0 then return "Feasible region unbounded!"; fi;

 #Print("Feasible region still looks bounded. \n");

 #figure out who will become zero
 newZero:=slackVars[Position(ratios,Minimum(positiveRatios))];

 #Print(newZero, " is the new zero guy \n");

 Remove(slackVars,Position(slackVars,newZero));
 Remove(vars,Position(vars,newNonzero));
 Add(vars,newZero);
 Add(slackVars,newNonzero);

 slackVars:=Set(slackVars);
 vars:=Set(vars);

 #Print("slackVars are ",slackVars,"\n");
 #Print("vars are ",vars,"\n");

 TriangulizeMatPivotColumns(M,Concatenation([1],slackVars));
 #Display(M);
 bestMove:=Minimum(List(vars,i->M[1,i]));

 newNonzero:=vars[Position(List(vars,i->M[1,i]),bestMove)];
 #Print(newNonzero," is the new nonzero guy");

 od;

 #calculate the original point and the max value there

 point:=List([2..2+n-1],x -> Val(M,vars,slackVars,len,x));
 value:=point*c;

 return [point,value];
end);

# can do better for matrices and large exponents, preliminary improvement,
# will be further improved (FL)
## InstallMethod( \^,
## "for matrices, use char. poly. for large exponents",
## [ IsMatrix, IsPosInt ],
## function(mat, n)
## local pol, indet;
## # generic method for small n, break even point probably a bit lower,
## # needs rethinking and some experiments.
## if n < 2^Length(mat) then
## return POW_OBJ_INT(mat, n);
## fi;
## pol := CharacteristicPolynomial(mat);
## indet := IndeterminateOfUnivariateRationalFunction(pol);
## # cost of this needs to be investigated
## pol := PowerMod(indet, n, pol);
## # now we are sure that we need at most Length(mat) matrix multiplications
## return Value(pol, mat);
## end);
##
# next iteration, conjugate matrix such that it is often very sparse
# (a companion matrix), could still be improved, maybe with kernel functions
# for compact matrices (FL)
BindGlobal("POW_MAT_INT", function(mat, n)
 local d, addb, trafo, value, t, ti, mm, pol, ind;
 d := NrRows(mat);
 # finding a better break even point probably also depends on q
 if n < 2^QuoInt(3*d,4) or not IsField(DefaultFieldOfMatrix(mat)) then
 return POW_OBJ_INT(mat, n);
 fi;
 # helper function to build up a semi-echelon basis
 addb := function(seb, v)
 local rows, pivots, len, vv, c, pos, i;
 rows := seb.vectors;
 pivots := seb.pivots;
 len := Length(rows);
 vv := ShallowCopy(v);
 for i in [1..len] do
 c := vv[pivots[i]];
 if not IsZero(c) then
 AddRowVector(vv, rows[i], -c);
 fi;
 od;
 pos := PositionNonZero(vv);
 if pos <= Length(vv) then
 if not IsOne(vv[pos]) then
 vv := vv/vv[pos];
 fi;
 Add(rows, vv);
 Add(pivots, pos);
 seb.heads[pos] := len + 1;
 return true;
 else
 return false;
 fi;
 end;
 # this returns a base change matrix such that t*m*t^-1 is block triangular
 # with companion matrices along the diagonal
 # (could/should? be improved to return t^-1, t*m*t^-1 and the
 # characteristic polynomial of m at the same time)
 trafo := function(m)
 local id, b, t, r, a;
 id := m^0;
 b := rec(vectors := [], pivots := [], heads := []);
 t := [];
 # maybe better start with a random vector?
 for a in id do
 r := addb(b,a);
 if r = true then
 repeat
 Add(t, a);
 a := a*m;
 r := addb(b,a);
 until r <> true;
 fi;
 od;
 t := Matrix(t, m);
 return t;
 end;
 # compared to standard method, we avoid some zero or identity matrices
 # and we multiply with mat from left to take advantage of sparseness of mat
 value := function(pol, mat)
 local f, c, i, val, j;
 f := CoefficientsOfLaurentPolynomial(pol);
 c := f[1];
 i := Length(c);
 if i = 0 then
 return 0*mat;
 fi;
 if i = 1 then
 val := POW_OBJ_INT(mat, f[2]);
 return c[1] * val;
 fi;
 val := c[i] * mat;
 if not IsMutable(val[1]) then
 val := MutableCopyMatrix(val);
 fi;
 i := i-1;
 for j in [1..NrRows(mat)] do
 val[j,j] := val[j,j]+c[i];
 od;
 while 1 f[2] then
 val := val * POW_OBJ_INT(mat, f[2]);
 fi;
 return val;
 end;
 t := trafo(mat);
 ti := t^-1;
 mm := t * mat * ti;
 pol := CharacteristicPolynomial(mm);
 ind := IndeterminateOfUnivariateRationalFunction(pol);
 pol := PowerMod(ind, n, pol);
 mm := value(pol, mm);
 return ti * mm * t;
end);

InstallMethod( \^,
 "for matrices, use char. poly. for large exponents",
 [ IsMatrixOrMatrixObj, IsPosInt ], POW_MAT_INT );

InstallGlobalFunction(RationalCanonicalFormTransform,function(mat)
local cr,R,x,com,nf,matt,p,i,j,di,d,v;
 matt:=TransposedMat(mat);
 cr:=DefaultFieldOfMatrix(mat);
 R:=PolynomialRing(cr,1);
 x:=IndeterminatesOfPolynomialRing(R)[1];
 com:=x*mat^0-mat;
 com:=List(com,ShallowCopy);
 nf:=DoDiagonalizeMat(R,com,true,true,matt);
 di:=DiagonalOfMat(nf.normal);
 p:=[];
 for i in [1..Length(di)] do
 d:=DegreeOfUnivariateLaurentPolynomial(di[i]);
 if d>0 then
 v:=List(nf.basmat[i],x->Value(x,Zero(cr))); # move in base ring
 Add(p,v);
 for j in [1..d-1] do
 v:=v*matt;
 Add(p,v);
 od;
 fi;
 od;
 return TransposedMat(p);
end);

Quelle matrix.gi Sprache: unbekannt

[zur Elbe Produktseite wechseln0.87QuellennavigatorsAnalyse erneut starten2026-05-06]