|
|
|
|
Impressum matrix.gi
Interaktion und
Portierbarkeitunbekannt
|
|
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, M artin 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,
"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 <p>
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 <U>
## factorzero: zero of complement
## subspace: triangulized basis of <mat>
## factorspace: base of a complement of <U> 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)<i or mat[i,j]=z) then
# Add vector from bas with j-th component not zero (if any exists)
v:=PositionProperty(bas,k->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( <B>, <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( <B>, <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>, <b>, <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 < i do
val := mat * val;
i := i - 1;
for j in [1..NrRows(mat)] do
val[j,j] := val[j,j]+c[i];
od;
od;
if 0 <> 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);
[Verzeichnis aufwärts0.83unsichere VerbindungÜbersetzung europäischer Sprachen durch Browser2026-05-06]
|
2026-05-26
|
|
|
|
|
