|
#############################################################################
##
## 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 <i> in position <i>
## 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 <i>th standard basis vector
# (equivalently, of the orbit of `mat[<i>]',
# 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>[<i>]' divides '<d>[<i>+1]' and <mat> is equivalent
## to a diagonal matrix with the elements '<d>[<i>]' 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,
--> --------------------
--> maximum size reached
--> --------------------
[ Dauer der Verarbeitung: 0.42 Sekunden
(vorverarbeitet)
]
|
