diff --git a/doc/ref/matobj.xml b/doc/ref/matobj.xml index fd66048984..1ebc4393de 100644 --- a/doc/ref/matobj.xml +++ b/doc/ref/matobj.xml @@ -4,8 +4,8 @@ - -Vector and matrix objects + +Vector and Matrix Objects This chapter is work in progress. It will eventually describe the new interface to vector and matrix objects. @@ -20,34 +20,46 @@ To eventually solve this problem, this chapter describes a new programming interface to vectors and matrices. -
- Fundamental ideas and rules +
+ Fundamental Ideas and Rules <#Include Label="MatObj_Overview">
-
- Categories of vectors and matrices +
+ Categories of Vectors and Matrices
-
- Constructing vector and matrix objects +
+ Constructing Vector and Matrix Objects
-
- Operations for row vector objects +
+ Operations for Vector Objects + + <#Include Label="MatObj_PositionNonZero"> + <#Include Label="MatObj_PositionLastNonZero"> + <#Include Label="MatObj_ListOp"> + <#Include Label="MatObj_UnpackVector"> + <#Include Label="MatObj_ConcatenationOfVectors"> + <#Include Label="MatObj_ExtractSubVector"> + <#Include Label="MatObj_ZeroVector"> + <#Include Label="MatObj_ConstructingFilter_Vector"> + <#Include Label="MatObj_Randomize_Vectors"> + <#Include Label="MatObj_WeightOfVector"> + <#Include Label="MatObj_DistanceOfVectors">
-
- Operations for row list matrix objects +
+ Operations for Row List Matrix Objects
-
- Operations for flat matrix objects +
+ Operations for Flat Matrix Objects
diff --git a/doc/ref/matrix.xml b/doc/ref/matrix.xml index 23ea02b98c..b852b70774 100644 --- a/doc/ref/matrix.xml +++ b/doc/ref/matrix.xml @@ -13,7 +13,7 @@ Matrices are represented in &GAP; by lists of row vectors (see ) (for future changes to this -policy see Chapter ). +policy see Chapter ). The vectors must all have the same length, and their elements must lie in a common ring. However, since checking rectangularness can be expensive functions and methods of operations for matrices often will not give an error diff --git a/hpcgap/lib/vecmat.gi b/hpcgap/lib/vecmat.gi index 255544d716..8c559c23f3 100644 --- a/hpcgap/lib/vecmat.gi +++ b/hpcgap/lib/vecmat.gi @@ -1844,18 +1844,6 @@ InstallMethod( Append, "for GF2 vectors", IsGF2VectorRep and IsList], 0, APPEND_GF2VEC); -############################################################################# -## -#M PositionCanonical( , ) . . for GF2 matrices -## -InstallMethod( PositionCanonical, - "for internally represented lists, fall back on `Position'", - true, # the list may be non-homogeneous. - [ IsList and IsGF2MatrixRep, IsObject ], 0, - function( list, obj ) - return Position( list, obj, 0 ); -end ); - ############################################################################# ## #M ShallowCopy( ) . . . for GF2 vectors @@ -2311,9 +2299,8 @@ InstallMethod(DeterminantMatDestructive, ## -InstallMethod(RankMatDestructive, +InstallOtherMethod(RankMatDestructive, "kernel method for plain list of GF2 vectors", - true, [IsMatrix and IsPlistRep and IsFFECollColl and IsMutable], GF2_AHEAD_OF_8BIT_RANK, RANK_LIST_GF2VECS); @@ -2422,27 +2409,15 @@ InstallMethod( Matrix, "for a list of vecs, an integer, and a gf2 mat", ConvertToMatrixRep(li,2); return li; end ); -BindGlobal( "PositionLastNonZeroFunc", + +InstallMethod( PositionLastNonZero, "for a row vector obj", + [IsVectorObj], function(l) local i; i := Length(l); while i >= 1 and IsZero(l[i]) do i := i - 1; od; return i; end ); -BindGlobal( "PositionLastNonZeroFunc2", - function(l,pos) - local i; - i := pos-1; - while i >= 1 and IsZero(l[i]) do i := i - 1; od; - return i; - end ); - -InstallMethod( PositionLastNonZero, "for a row vector obj", - [IsVectorObj], PositionLastNonZeroFunc ); -InstallMethod( PositionLastNonZero, "for a matrix obj", - [IsMatrixObj], PositionLastNonZeroFunc ); -InstallMethod( PositionLastNonZero, "for a matrix obj, and an index", - [IsMatrixObj, IsPosInt], PositionLastNonZeroFunc2 ); InstallMethod( ExtractSubMatrix, "for a gf2 matrix, and two lists", [IsGF2MatrixRep, IsList, IsList], diff --git a/lib/algebra.gi b/lib/algebra.gi index 79660edba5..c535061db0 100644 --- a/lib/algebra.gi +++ b/lib/algebra.gi @@ -930,7 +930,6 @@ InstallMethod( IsAnticommutative, function( A ) local n, # dimension of `A' T, # table of structure constants for `A' - zero, # zero coefficient i, j; # loop over rows and columns ot `T' if not IsFiniteDimensional( A ) then @@ -939,12 +938,11 @@ InstallMethod( IsAnticommutative, n:= Dimension( A ); T:= StructureConstantsTable( Basis( A ) ); - zero:= T[ n+2 ]; for i in [ 2 .. n ] do for j in [ 1 .. i-1 ] do if T[i][j][1] <> T[j][i][1] or ( not IsEmpty( T[i][j][1] ) - and PositionNot( T[i][j][2] + T[j][i][2], zero ) + and PositionNonZero( T[i][j][2] + T[j][i][2] ) <= Length( T[i][j][2] ) ) then return false; fi; @@ -2974,23 +2972,23 @@ InstallMethod( IsNilpotentElement, local B, # a basis of `L' A, # adjoint matrix of `x w.r. to `B' n, # dimension of `L' - i, # loop variable - zero; # zero coefficient + i; # loop variable B := Basis( L ); A := AdjointMatrix( B, x ); n := Dimension( L ); i := 1; - zero:= Zero( A[1][1] ); - if ForAll( A, x -> n < PositionNot( x, zero ) ) then + if ForAll( A, x -> n < PositionNonZero( x ) ) then +#T better ask IsZero? return true; fi; while i < n do i:= 2 * i; A:= A * A; - if ForAll( A, x -> n < PositionNot( x, zero ) ) then + if ForAll( A, x -> n < PositionNonZero( x ) ) then +#T better ask IsZero? return true; fi; od; diff --git a/lib/algsc.gi b/lib/algsc.gi index 7063f550a9..56b76c6e1a 100644 --- a/lib/algsc.gi +++ b/lib/algsc.gi @@ -162,8 +162,7 @@ InstallMethod( PrintObj, return; fi; - zero := Zero( elm[1] ); - depth := PositionNot( elm, zero ); + depth := PositionNonZero( elm ); if len < depth then @@ -174,6 +173,7 @@ InstallMethod( PrintObj, else one:= One( elm[1] ); + zero:= Zero( elm[1] ); if elm[ depth ] <> one then Print( "(", elm[ depth ], ")*" ); @@ -221,8 +221,7 @@ InstallMethod( String, return ""; fi; - zero := Zero( elm[1] ); - depth := PositionNot( elm, zero ); + depth := PositionNonZero( elm ); s:=""; if len < depth then @@ -235,6 +234,7 @@ InstallMethod( String, else one:= One( elm[1] ); + zero:= Zero( elm[1] ); if elm[ depth ] <> one then Add(s,'('); diff --git a/lib/ctblpope.gi b/lib/ctblpope.gi index b8ecac6b3d..42358bc34b 100644 --- a/lib/ctblpope.gi +++ b/lib/ctblpope.gi @@ -1384,7 +1384,8 @@ InstallGlobalFunction( PermCandidates, # delete zero rows: shrink:= []; for i in matrix do - if PositionNot( i, 0 ) <= Length( i ) then + if PositionNonZero( i ) <= Length( i ) then +#T better call IsZero? Add( shrink, i ); fi; od; @@ -2156,7 +2157,8 @@ InstallGlobalFunction( PermCandidatesFaithful, # delete zero rows: shrink:= []; for i in matrix do - if PositionNot( i, 0 ) <= Length( i ) then + if PositionNonZero( i ) <= Length( i ) then +#T better call IsZero? Add( shrink, i ); fi; od; diff --git a/lib/liefam.gi b/lib/liefam.gi index df75f0f06b..7706005a8a 100644 --- a/lib/liefam.gi +++ b/lib/liefam.gi @@ -368,28 +368,24 @@ InstallMethod( AdditiveInverseOp, #M IsBound\[\]( , ) #M Position( , ) ## -InstallMethod( \[\], +InstallOtherMethod( \[\], "for Lie matrix in default representation, and positive integer", - true, - [ IsLieMatrix and IsPackedElementDefaultRep, IsPosInt ], 0, + [ IsLieMatrix and IsPackedElementDefaultRep, IsPosInt ], function( mat, i ) return mat![1][i]; end ); -InstallMethod( Length, +InstallOtherMethod( Length, "for Lie matrix in default representation", - true, - [ IsLieMatrix and IsPackedElementDefaultRep ], 0, - mat -> Length( mat![1] ) ); + [ IsLieMatrix and IsPackedElementDefaultRep ], + mat -> NumberRows( mat![1] ) ); InstallMethod( IsBound\[\], "for Lie matrix in default representation, and integer", - true, - [ IsLieMatrix and IsPackedElementDefaultRep, IsPosInt ], 0, + [ IsLieMatrix and IsPackedElementDefaultRep, IsPosInt ], function( mat, i ) return IsBound( mat![1][i] ); end ); InstallMethod( Position, "for Lie matrix in default representation, row vector, and integer", - true, - [ IsLieMatrix and IsPackedElementDefaultRep, IsRowVector, IsInt ], 0, + [ IsLieMatrix and IsPackedElementDefaultRep, IsRowVector, IsInt ], function( mat, v, pos ) return Position( mat![1], v, pos ); end ); diff --git a/lib/list.gi b/lib/list.gi index b8e8bb43f7..f9485936e3 100644 --- a/lib/list.gi +++ b/lib/list.gi @@ -3807,7 +3807,7 @@ InstallMethod( PositionNot, "default method ", [IsList, IsObject, IsInt ], InstallOtherMethod( PositionNot, "default value of third argument ", [IsList, IsObject], function(l,x) - return PositionNot(l,x,0); + return POSITION_NOT(l,x,0); end ); @@ -3816,7 +3816,7 @@ InstallMethod( PositionNonZero, "default method", [IsHomogeneousList], if Length(l) = 0 then return 1; else - return PositionNot(l, Zero(l[1])); + return POSITION_NOT(l, Zero(l[1]), 0); fi; end); @@ -3825,7 +3825,7 @@ InstallMethod( PositionNonZero, "default method with start", [IsHomogeneousList, if Length(l) = 0 then return from+1; fi; - return PositionNot(l, Zero(l[1]), from); + return POSITION_NOT(l, Zero(l[1]), from); end); diff --git a/lib/mat8bit.gi b/lib/mat8bit.gi index 687e5bf221..40d7f52f53 100644 --- a/lib/mat8bit.gi +++ b/lib/mat8bit.gi @@ -1050,11 +1050,9 @@ InstallMethod(DeterminantMatDestructive, ## -InstallMethod(RankMatDestructive, +InstallOtherMethod(RankMatDestructive, "kernel method for plain list of GF2 vectors", - true, [IsMatrix and IsPlistRep and IsFFECollColl and IsMutable], - 0, RANK_LIST_VEC8BITS); InstallMethod(NestingDepthM, [Is8BitMatrixRep], m->2); diff --git a/lib/matblock.gi b/lib/matblock.gi index fac4901609..6f7504cbaa 100644 --- a/lib/matblock.gi +++ b/lib/matblock.gi @@ -132,7 +132,7 @@ end ); ## #M Length( ) . . . . . . . . . . . . . . . . for a block matrix ## -InstallMethod( Length, +InstallOtherMethod( Length, "for an ordinary block matrix", [ IsOrdinaryMatrix and IsBlockMatrixRep ], blockmat -> blockmat!.nrb * blockmat!.rpb ); @@ -142,7 +142,7 @@ InstallMethod( Length, ## #M \[\]( , ) . . . . . . . . . . . . . . . for a block matrix ## -InstallMethod( \[\], +InstallOtherMethod( \[\], "for an ordinary block matrix and a positive integer", [ IsOrdinaryMatrix and IsBlockMatrixRep, IsPosInt ], function( blockmat, n ) @@ -195,7 +195,7 @@ InstallOtherMethod( \[\], ## #M TransposedMat( ) . . . . . . . . . . . . . for a block matrix ## -InstallMethod( TransposedMat, +InstallOtherMethod( TransposedMat, "for an ordinary block matrix", [ IsOrdinaryMatrix and IsBlockMatrixRep ], m -> BlockMatrix( List( m!.blocks, i -> [ i[2], i[1], @@ -668,7 +668,7 @@ InstallMethod( PrintObj, ## #M DimensionsMat( ) . . . . . . . . . . . . . for a block matrix ## -InstallMethod( DimensionsMat, +InstallOtherMethod( DimensionsMat, "for an ordinary block matrix", [ IsOrdinaryMatrix and IsBlockMatrixRep ], m -> [ m!.nrb * m!.rpb, m!.ncb * m!.cpb ] ); diff --git a/lib/matobj.gi b/lib/matobj.gi index eff43509b4..b73d526375 100644 --- a/lib/matobj.gi +++ b/lib/matobj.gi @@ -12,6 +12,45 @@ ############################################################################ +# Install fallback methods for m[i,j] which delegate MatElm resp. SetMatElm, +# for old MatrixObj implementation which don't provide them. We lower the rank +# so that these are only used as a last resort. +InstallMethod( \[\], "for a matrix object and two positions", + [ IsMatrixObj, IsPosInt, IsPosInt ], + -RankFilter(IsMatrixObj), + function( m, row, col ) + return MatElm( m, row, col ); + end ); + +InstallMethod( \[\]\:\=, "for a matrix object, two positions, and an object", + [ IsMatrixObj and IsMutable, IsPosInt, IsPosInt, IsObject ], + -RankFilter(IsMatrixObj), + function( m, row, col, obj ) + SetMatElm( m, row, col, obj ); + end ); + + +InstallMethod( OneOfBaseDomain, + "generic method for IsVectorObj", + [ IsVectorObj ], + v -> One( BaseDomain( v ) ) ); + +InstallMethod( ZeroOfBaseDomain, + "generic method for IsVectorObj", + [ IsVectorObj ], + v -> Zero( BaseDomain( v ) ) ); + +InstallMethod( OneOfBaseDomain, + "generic method for IsMatrixObj", + [ IsMatrixObj ], + M -> One( BaseDomain( M ) ) ); + +InstallMethod( ZeroOfBaseDomain, + "generic method for IsMatrixObj", + [ IsMatrixObj ], + M -> Zero( BaseDomain( M ) ) ); + + InstallMethod( WeightOfVector, "generic method", [IsVectorObj], function(v) @@ -42,18 +81,226 @@ InstallMethod( DistanceOfVectors, "generic method", return n; end ); +# +# TODO: possibly rename the following +# +BindGlobal( "DefaultVectorRepForBaseDomain", +function( basedomain ) + if IsFinite(basedomain) and IsField(basedomain) and Size(basedomain) = 2 then + return IsGF2VectorRep; + elif IsFinite(basedomain) and IsField(basedomain) and Size(basedomain) <= 256 then + return Is8BitVectorRep; + fi; + return IsPlistVectorRep; +end); +BindGlobal( "DefaultMatrixRepForBaseDomain", +function( basedomain ) + if IsFinite(basedomain) and IsField(basedomain) and Size(basedomain) = 2 then + return IsGF2MatrixRep; + elif IsFinite(basedomain) and IsField(basedomain) and Size(basedomain) <= 256 then + return Is8BitMatrixRep; + fi; + return IsPlistMatrixRep; +end); + + + +# methods to create vectors + +InstallMethod( Vector, + [IsOperation, IsSemiring, IsList], +function(rep, base, l) + return NewVector(rep, base, l); +end); + +InstallMethod( Vector, + [IsOperation, IsSemiring, IsVectorObj], +function(rep, base, v) + return NewVector(rep, base, Unpack(v)); +end); + +InstallMethod( Vector, + [IsSemiring, IsList], +function(base, l) + return NewVector(DefaultVectorRepForBaseDomain(base), base, l); +end); + +InstallMethod( Vector, + [IsSemiring, IsVectorObj], +function(base, v) + return NewVector(DefaultVectorRepForBaseDomain(base), base, Unpack(v)); +end); + +InstallMethod( Vector, + [IsList, IsVectorObj], +function(l, example) + return NewVector(ConstructingFilter(example), BaseDomain(example), l); +end); + +InstallMethod( Vector, + [IsVectorObj, IsVectorObj], +function(v, example) + return NewVector(ConstructingFilter(example), BaseDomain(example), Unpack(v)); +end); + +InstallMethod( Vector, + [IsList], +function(l) + local dom; + dom := DefaultScalarDomainOfMatrixList([[l]]); + return NewVector(DefaultVectorRepForBaseDomain(dom), dom, l); +end); +# +# +# +InstallMethod( Matrix, + [IsOperation, IsSemiring, IsList, IsInt], + function( rep, basedomain, list, nrCols ) + # TODO: adjust NewMatrix to use same arg order as Matrix (or vice-versa) + return NewMatrix( rep, basedomain, nrCols, list ); + end ); + +InstallMethod( Matrix, + [IsOperation, IsSemiring, IsList], + function( rep, basedomain, list ) + if Length(list) = 0 then Error("list must be not empty; to create empty matrices, please specify nrCols"); fi; + return NewMatrix( rep, basedomain, Length(list[1]), list ); + end ); + +InstallMethod( Matrix, + [IsOperation, IsSemiring, IsMatrixObj], + function( rep, basedomain, mat ) + # TODO: can we do better? encourage MatrixObj implementors to overload this? + return NewMatrix( rep, basedomain, NrCols(mat), Unpack(mat) ); + end ); + +# +# +# +InstallMethod( Matrix, + [IsSemiring, IsList, IsInt], + function( basedomain, list, nrCols ) + local rep; + rep := DefaultMatrixRepForBaseDomain(basedomain); + return NewMatrix( rep, basedomain, nrCols, list ); + end ); + +InstallMethod( Matrix, + [IsSemiring, IsList], + function( basedomain, list ) + local rep; + if Length(list) = 0 then Error("list must be not empty"); fi; + rep := DefaultMatrixRepForBaseDomain(basedomain); + return NewMatrix( rep, basedomain, Length(list[1]), list ); + end ); + +InstallMethod( Matrix, + [IsSemiring, IsMatrixObj], + function( basedomain, mat ) + # TODO: can we do better? encourage MatrixObj implementors to overload this? + return NewMatrix( DefaultMatrixRepForBaseDomain(basedomain), basedomain, NrCols(mat), Unpack(mat) ); + end ); + +# +# +# +InstallMethod( Matrix, + [IsList, IsInt], + function( list, nrCols ) + local basedomain, rep; + if Length(list) = 0 then Error("list must be not empty"); fi; + if Length(list[1]) = 0 then Error("list[1] must be not empty, please specify base domain explicitly"); fi; + basedomain := DefaultScalarDomainOfMatrixList([list]); + rep := DefaultMatrixRepForBaseDomain(basedomain); + return NewMatrix( rep, basedomain, nrCols, list ); + end ); + +InstallMethod( Matrix, + [IsList], + function( list ) + local rep, basedomain; + if Length(list) = 0 then Error("list must be not empty"); fi; + if Length(list[1]) = 0 then Error("list[1] must be not empty, please specify base domain explicitly"); fi; + basedomain := DefaultScalarDomainOfMatrixList([list]); + rep := DefaultMatrixRepForBaseDomain(basedomain); + return NewMatrix( rep, basedomain, Length(list[1]), list ); + end ); + +# +# matrix constructors using example objects (as last argument) +# +InstallMethod( Matrix, + [IsList, IsInt, IsMatrixObj], + function( list, nrCols, example ) + return NewMatrix( ConstructingFilter(example), BaseDomain(example), nrCols, list ); + end ); + InstallMethod( Matrix, "generic convenience method with 2 args", - [IsList,IsMatrixObj], - function( l, m ) - if Length(l) = 0 then + [IsList, IsMatrixObj], + function( list, example ) + if Length(list) = 0 then ErrorNoReturn("Matrix: two-argument version not allowed with empty first arg"); fi; - if not (IsList(l[1]) or IsVectorObj(l[1])) then + if not (IsList(list[1]) or IsVectorObj(list[1])) then ErrorNoReturn("Matrix: flat data not supported in two-argument version"); fi; - return Matrix(l,Length(l[1]),m); + return NewMatrix( ConstructingFilter(example), BaseDomain(example), Length(list[1]), list ); + end ); + +InstallMethod( Matrix, + [IsMatrixObj, IsMatrixObj], + function( mat, example ) + # TODO: can we avoid using Unpack? resp. make this more efficient + # perhaps adjust NewMatrix to take an IsMatrixObj? + return NewMatrix( ConstructingFilter(example), BaseDomain(example), NrCols(mat), Unpack(mat) ); + end ); + +# +# +# +InstallMethod( ZeroMatrix, + [IsInt, IsInt, IsMatrixObj], + function( rows, cols, example ) + return ZeroMatrix( ConstructingFilter(example), BaseDomain(example), rows, cols ); + end ); + +InstallMethod( ZeroMatrix, + [IsSemiring, IsInt, IsInt], + function( basedomain, rows, cols ) + return ZeroMatrix( DefaultMatrixRepForBaseDomain(basedomain), basedomain, rows, cols ); + end ); + +InstallMethod( ZeroMatrix, + [IsOperation, IsSemiring, IsInt, IsInt], + function( rep, basedomain, rows, cols ) + # TODO: urge matrixobj implementors to overload this + return NewMatrix( rep, basedomain, cols, ListWithIdenticalEntries( rows * cols, Zero(basedomain) ) ); + end ); + +# +# +# +InstallMethod( IdentityMatrix, + [IsInt, IsMatrixObj], + function( dim, example ) + return IdentityMatrix( ConstructingFilter(example), BaseDomain(example), dim ); + end ); + +InstallMethod( IdentityMatrix, + [IsSemiring, IsInt], + function( basedomain, dim ) + return IdentityMatrix( DefaultMatrixRepForBaseDomain(basedomain), basedomain, dim ); end ); +InstallMethod( IdentityMatrix, + [IsOperation, IsSemiring, IsInt], + function( rep, basedomain, dim ) + # TODO: avoid using IdentityMat eventually + return NewMatrix( rep, basedomain, dim, IdentityMat( dim, basedomain ) ); + end ); + + + InstallMethod( Unfold, "for a matrix object, and a vector object", [ IsMatrixObj, IsVectorObj ], function( m, w ) @@ -86,9 +333,8 @@ InstallMethod( Fold, "for a vector, a positive int, and a matrix", InstallMethod( CompanionMatrix, "for a polynomial and a matrix", [ IsUnivariatePolynomial, IsMatrixObj ], function( po, m ) - local l, n, q, ll, i, bd, one; - bd := BaseDomain(m); - one := One(bd); + local l, n, q, ll, i, one; + one := OneOfBaseDomain( m ); l := CoefficientsOfUnivariatePolynomial(po); n := Length(l)-1; if not IsOne(l[n+1]) then @@ -182,35 +428,6 @@ InstallGlobalFunction( MakeVector, return NewVector(ty,bd,l); end ); -InstallGlobalFunction( MakeMatrix, - function( arg ) - local bd,l,len,rowlen,ty; - if Length(arg) = 1 then - l := arg[1]; - bd := DefaultFieldOfMatrix(l); - elif Length(arg) <> 2 then - Error("usage: MakeVector( [,] )"); - return fail; - else - l := arg[1]; - bd := arg[2]; - fi; - len := Length(l); - if len = 0 then - Error("does not work for matrices with zero rows"); - return fail; - fi; - rowlen := Length(l[1]); - if IsFinite(bd) and IsField(bd) and Size(bd) = 2 then - ty := IsGF2MatrixRep; - elif IsFinite(bd) and IsField(bd) and Size(bd) <= 256 then - ty := Is8BitMatrixRep; - else - ty := IsPlistMatrixRep; - fi; - return NewMatrix(ty,bd,rowlen,l); - end ); - InstallMethod( ExtractSubVector, "generic method", [ IsVectorObj, IsList ], function( v, l ) @@ -220,13 +437,12 @@ InstallMethod( ExtractSubVector, "generic method", InstallOtherMethod( ScalarProduct, "generic method", [ IsVectorObj, IsVectorObj ], function( v, w ) - local bd,i,s; - bd := BaseDomain(v); - s := Zero(bd); + local i,s; if Length(v) <> Length(w) then Error("vectors must have equal length"); return fail; fi; + s:= ZeroOfBaseDomain( v ); for i in [1..Length(v)] do s := s + v[i]*w[i]; od; @@ -236,19 +452,131 @@ InstallOtherMethod( ScalarProduct, "generic method", InstallMethod( TraceMat, "generic method", [ IsMatrixObj ], function( m ) - local bd,i,s; - bd := BaseDomain(m); - s := Zero(bd); + local i,s; if NumberRows(m) <> NumberColumns(m) then Error("matrix must be square"); return fail; fi; + s := ZeroOfBaseDomain( m ); for i in [1..NumberRows(m)] do - s := s + MatElm(m,i,i); + s := s + m[ i, i ]; od; return s; end ); +InstallMethod(PositionNonZero, + "generic method for a row vector", + [IsRowVector], + function(vec) + local i; + for i in [1..Length(vec)] do + if not IsZero(vec[i]) then return i; fi; + od; + return i+1; +end); +#T superfluous? + + +InstallMethod(PositionNonZero, + "generic method for a vector object", + [ IsVectorObj ], + function(vec) + local i; + for i in [1..Length(vec)] do + if not IsZero(vec[i]) then return i; fi; + od; + return i+1; +end); + +InstallMethod( ListOp, + "generic method for a vector object", + [ IsVectorObj ], + function(vec) + local result, i, len; + len := Length(vec); + result := []; + result[len] := vec[len]; + for i in [ 1 .. len - 1 ] do + result[i] := vec[i]; + od; + return result; +end ); + +InstallMethod( ListOp, + "generic method for a vector object and a function", + [ IsVectorObj, IsFunction ], + function(vec,func) + local result, i, len; + len := Length(vec); + result := []; + result[len] := func(vec[len]); + for i in [ 1 .. len - 1 ] do + result[i] := func(vec[i]); + od; + return result; +end ); + +InstallMethod( Unpack, + "generic method for a vector object", + [ IsVectorObj ], + ListOp ); ## Potentially slower than a direct implementation, + ## but avoids code duplication. + + +InstallMethod( \{\}, + "generic method for a vector object and a list", + [ IsVectorObj, IsList ], + function(vec,poss) + local vec_list; + vec_list := ListOp(vec); + vec_list := vec_list{poss}; + return Vector(vec_list,vec); +end ); + +InstallMethod( CopySubVector, + "generic method for vectors", + [ IsVectorObj and IsMutable, IsList, IsVectorObj, IsList ], + function(dst, dcols, src, scols) + local i; + if not Length( dcols ) = Length( scols ) then + Error( "source and destination index lists must be of equal length" ); + return; + fi; + for i in [ 1 .. Length( dcols ) ] do + dst[dcols[i]] := src[scols[i]]; + od; +end ); + +## Backwards compatible version +InstallMethod( CopySubVector, + "generic method for vectors", + [ IsVectorObj, IsVectorObj and IsMutable, IsList, IsList ], + function(src, dst, scols, dcols) + CopySubVector(dst,dcols,src,scols); +end ); + +InstallMethod( Randomize, + "generic method for a vector", + [ IsVectorObj and IsMutable ], + function(vec) + local basedomain, i; + basedomain := BaseDomain( vec ); + for i in [ 1 .. Length( vec ) ] do + vec[ i ] := Random( basedomain ); + od; +end ); + +InstallMethod( Randomize, + "generic method for a vector and a random source", + [ IsVectorObj and IsMutable, IsRandomSource ], + function(vec, rs) + local basedomain, i; + basedomain := BaseDomain( vec ); + for i in [ 1 .. Length( vec ) ] do + vec[ i ] := Random( rs, basedomain ); + od; +end ); + # # Compatibility code: Install MatrixObj methods for IsMatrix. # @@ -256,3 +584,9 @@ InstallOtherMethod( NumberRows, "for a plist matrix", [ IsMatrix ], Length); InstallOtherMethod( NumberColumns, "for a plist matrix", [ IsMatrix ], m -> Length(m[1]) ); + +# +# Compatibility code: Generic methods for IsMatrixObj +# +InstallOtherMethod( DimensionsMat, "for a matrix in IsMatrixObj", + [ IsMatrixObj ], m -> [ NumberRows( m ), NumberColumns( m ) ] ); diff --git a/lib/matobj1.gd b/lib/matobj1.gd index 7388fd0ffe..0c1af5caba 100644 --- a/lib/matobj1.gd +++ b/lib/matobj1.gd @@ -36,24 +36,40 @@ DeclareSynonym( "IsRowVectorObj", IsVectorObj ); # We should eventually remove this synonym. -# There are one main category for matrices and two disjoint sub-categories: +# There are one main category for matrices and one subcategory: DeclareCategory( "IsMatrixObj", IsVector and IsScalar and IsCopyable ); # All the arithmetical filters come from IsVector and IsScalar. # In particular, matrices are in "IsMultiplicativeElement" which defines # powering with a positive integer by the (kernel) method for POW_OBJ_INT. # Note that this is at least strange for non-associative base domains. +# The filter 'IsMatrixObj' for an object does *not* imply that the +# multiplication for this object is the usual matrix multiplication, +# one can specify this multiplication via the filter 'IsOrdinaryMatrix'. +# Also the associativity of the multiplication of matrices does not follow +# from 'IsMatrixObj' and the associativity of the base domain, +# one needs also 'IsOrdinaryMatrix' for this implication. +# Note that elements of matrix Lie algebras lie in 'IsMatrixObj' but not in +# 'IsOrdinaryMatrix'. # Matrices are no longer necessarily lists, since they do not promise all list # operations! Of course, in specific implementations the objects may # still be lists. # The family of an object in IsMatrixObj is the collections family of # the family of its base domain. -InstallTrueMethod(IsAssociativeElement, - IsMatrixObj and IsAssociativeElementCollColl); +InstallTrueMethod( IsMatrixObj, IsMatrix ); +# We want that those matrices that are plain lists of plain lists +# are also in 'IsMatrixObj', in order to be able to install a method for +# all representations of matrices with the same requirements. +# (With the current distribution of declarations to files, +# we cannot put the implication into the declaration of 'IsMatrix': +# Both 'IsVector' and 'IsMatrix' are declared in 'lib/arith.gd', +# and 'IsMatrixObj' --which shall be in the middle-- is declared in +# 'lib/matobj1.gd'.) DeclareCategory( "IsRowListMatrix", IsMatrixObj ); # The category of matrices behaving like lists of rows which are GAP objects. # Different matrices in this category can share rows and the same row can # occur more than once in a matrix. Row access just gives a reference # to the row object. + diff --git a/lib/matobj2.gd b/lib/matobj2.gd index 3b4b0075d2..62f3c0dcc3 100644 --- a/lib/matobj2.gd +++ b/lib/matobj2.gd @@ -13,6 +13,9 @@ # ############################################################################ +# TODO: make sure to document what exactly a new matrix obj implementation has to +# provide, and that we provide default implementations for everything else. + ############################################################################ # @@ -116,12 +119,54 @@ ############################################################################ -# The following are guaranteed to be always set or cheaply calculable: +############################################################################# +## +#A BaseDomain( ) +## +## <#GAPDoc Label="MatObj_BaseDomain_IsVectorObj"> +## +## +## +## +## The object vector is defined over. +## Note that not all entries of vector necessarily lie in +## BaseDomain( vector ) with respect to in. +## E.g. vector can be defined over a polynomial ring but also contain +## integers +## (see section ). +## +## +## <#/GAPDoc> +## DeclareAttribute( "BaseDomain", IsVectorObj ); # Typically, the base domain will be a ring, it need not be commutative # nor associative. For non-associative base domains, the behavior of # powering matrices is undefined. +DeclareAttribute( "OneOfBaseDomain", IsVectorObj ); +DeclareAttribute( "ZeroOfBaseDomain", IsVectorObj ); +# 'BaseDomain' shall work also for IsPlist objects (vectors and matrices), +# but it is regarded as expensive to call this function; +# fetching a one or zero is cheaper, and if one needs these values for a +# non-plist vector or matrix then calling 'OneOfBaseDomain' or +# 'ZeroOfBaseDomain' is not more expensive than calling 'BaseDomain' first. + + +############################################################################# +## +#A Length( ) +## +## <#GAPDoc Label="MatObj_Length_IsVectorObj"> +## +## +## +## +## returns the length of the vector object vector, which is defined to +## be the number of entries of vector. +## +## +## <#/GAPDoc> +## DeclareAttribute( "Length", IsVectorObj ); # can be zero # We have to declare this since a row vector is not necessarily # a list! Correspondingly we have to use InstallOtherMethod @@ -164,19 +209,68 @@ DeclareOperation( "{}", [IsVectorObj,IsList] ); # Of course the positions must all lie in [1..Length(VECTOR)]. # Returns a vector in the same representation! +## <#GAPDoc Label="MatObj_PositionNonZero"> +## +## +## An integer +## +## Returns the index of the first entry in the vector V which is not +## zero. If all entries are zero, the function +## returns Length(V) + 1. +## +## +## <#/GAPDoc> DeclareOperation( "PositionNonZero", [IsVectorObj] ); +## <#GAPDoc Label="MatObj_PositionLastNonZero"> +## +## +## An integer +## +## Returns the index of the last entry in the vector V which is not +## zero. If all entries are zero, the function +## returns 0. +## +## +## <#/GAPDoc> DeclareOperation( "PositionLastNonZero", [IsVectorObj] ); +## <#GAPDoc Label="MatObj_ListOp"> +## +## +## A plain list +## +## Applies func to each entry of the vector V and returns +## the results as a plain list. This allows for calling +## on vectors. +## If the argument func is not provided, applies +## to all entries. +## +## +## <#/GAPDoc> DeclareOperation( "ListOp", [IsVectorObj] ); DeclareOperation( "ListOp", [IsVectorObj,IsFunction] ); # This is an unpacking operation returning a mutable copy in form of a list. # It enables the "List" function to work. -# The following unwraps a vector to a list: -DeclareOperation( "Unpack", [IsVectorObj] ); +## <#GAPDoc Label="MatObj_UnpackVector"> +## +## +## A plain list +## +## Returns a new plain list containing the entries of V. +## Guarantees to return a new list which can be manipulated without +## changing V. The entries itself are not copied. +## +## +## <#/GAPDoc> +DeclareOperation( "Unpack", [IsVectorObj] ); # It guarantees to copy, that is changing the returned object does # not change the original object. +# TODO: replace by AsList ? +# TODO: this is already used by the fining package + # "PositionNot" is intentionally left out here because it can rarely # be implemented much more efficiently than by running through the vector. @@ -188,13 +282,41 @@ DeclareOperation( "Unpack", [IsVectorObj] ); # Note that \{\}\:\= is left out here since it tempts the programmer # to use constructions like A{[1..3]} := B{[4,5,6]} which produces # an intermediate object. Use CopySubVector instead! + + # The list operations Position and so on seem to be unnecessary for # vectors and matrices and thus are left out to simplify the interface. +# TODO: actually -- why not allow `Position` anyway? What's the harm? + # Note that since Concatenation is a function using Append, it will # not work for vectors and it cannot be overloaded! # Thus we need: + +## <#GAPDoc Label="MatObj_ConcatenationOfVectors"> +## +## +## +## a vector object +## +## Returns a new vector containing the entries of V1, +## V2, etc. As prototype V1 is used. +## +## +## <#/GAPDoc> DeclareGlobalFunction( "ConcatenationOfVectors" ); +## <#GAPDoc Label="MatObj_ExtractSubVector"> +## +## +## a vector object +## +## Returns a new vector containing the entries of V +## at the positions in l. +## +## +## <#/GAPDoc> DeclareOperation( "ExtractSubVector", [IsVectorObj,IsList] ); # Does the same as slicing v{l} but is here to be similar to # ExtractSubMatrix. @@ -236,28 +358,71 @@ DeclareOperation( "ExtractSubVector", [IsVectorObj,IsList] ); # Note2: If sorting is not done lexicographically then the objects # in that representation cannot be lists! +# TODO: rename AddRowVector to AddVector; but keep in mind that +# historically there already was AddRowVector, so be careful to not break that + # The following "in place" operations exist with the same restrictions: DeclareOperation( "AddRowVector", [ IsVectorObj and IsMutable, IsVectorObj ] ); + +# vec = vec2 * scal DeclareOperation( "AddRowVector", - [ IsVectorObj and IsMutable, IsVectorObj, IsObject ] ); + [ IsVectorObj and IsMutable, IsVectorObj, IsObject ] ); +# vec = scal * vec2 DeclareOperation( "AddRowVector", - [ IsVectorObj and IsMutable,IsVectorObj,IsObject,IsPosInt,IsPosInt ] ); + [ IsVectorObj and IsMutable, IsObject, IsVectorObj ] ); + +# vec := vec2{[to..from]} * scal +DeclareOperation( "AddRowVector", + [ IsVectorObj and IsMutable, IsVectorObj, IsObject, IsPosInt, IsPosInt ] ); + +# vec := scal * vec2{[to..from]} +DeclareOperation( "AddRowVector", + [ IsVectorObj and IsMutable, IsObject, IsVectorObj, IsPosInt, IsPosInt ] ); + + + +# TODO: rename MultRowVector to MultVector; but keep in mind that +# historically there already was MultRowVector, so be careful to not break that DeclareOperation( "MultRowVector", [ IsVectorObj and IsMutable, IsObject ] ); + +# +# Also, make it explicit from which side we multiply +# DeclareOperation( "MultRowVectorFromLeft", +# [ IsVectorObj and IsMutable, IsObject ] ); +# DeclareOperation( "MultRowVectorFromRight", +# [ IsVectorObj and IsMutable, IsObject ] ); +#DeclareSynonym( "MultRowVector", MultRowVectorFromRight ); + +# do we really need the following? for what? is any code using this right now? +# ( a, pa, b, pb, s ) -> a{pa} := b{pb} * s; DeclareOperation( "MultRowVector", [ IsVectorObj and IsMutable, IsList, IsVectorObj, IsList, IsObject ] ); -# The following operations for scalars and vectors are possible of course -# only for scalars in the BaseDomain: -# *, / (of course only vector/scalar) +# maybe have this: vec := vec{[from..to]} * scal ?? cvec has it -# The following unary arithmetical operations are possible for vectors: + +# The following operations for scalars and vectors are possible for scalars in the BaseDomain +# (and often also for more, e.g. usually the scalar is allowed to be an integer regardless of the +# base domain): +# *, / (for /, we do not define /) + +# The following unary arithmetical operations are possible for vectors, assuming +# they are possible in the base domain (so all of them in fields, but e.g. in +# a proper semiring, there is in general no additive inverse): # AdditiveInverseImmutable, AdditiveInverseMutable, # AdditiveInverseSameMutability, ZeroImmutable, ZeroMutable, # ZeroSameMutability, IsZero, Characteristic -DeclareOperation( "ScalarProduct", [ IsVectorObj, IsVectorObj ] ); + +# ScalarProduct is already overloaded a lot, so perhaps we don't need to define +# it here, and just expect people to write vec1*vec2. +# +# TODO: document very explicitly that * for vectors does not care whether +# either vector is a row or column vector; it just always is the scalar product. +# +##DeclareOperation( "ScalarProduct", [ IsVectorObj, IsVectorObj ] ); # This is defined for two vectors of equal length, it returns the standard # scalar product. @@ -265,6 +430,16 @@ DeclareOperation( "ScalarProduct", [ IsVectorObj, IsVectorObj ] ); # The "representation-preserving" contructor methods: ############################################################################ +## <#GAPDoc Label="MatObj_ZeroVector"> +## +## +## a vector object +## +## Returns a new vector of length l in the same representation +## as V containing only zeros. +## +## +## <#/GAPDoc> DeclareOperation( "ZeroVector", [IsInt,IsVectorObj] ); # Returns a new mutable zero vector in the same rep as the given one with # a possible different length. @@ -273,7 +448,27 @@ DeclareOperation( "ZeroVector", [IsInt,IsMatrixObj] ); # Returns a new mutable zero vector in a rep that is compatible with # the matrix but of possibly different length. -DeclareOperation( "Vector", [IsList,IsVectorObj]); +# Operation to create vector objects. +# The first just delegate to NewVector: +# TODO: replace IsOperation by something nicer, like IsFilter (but sadly that +# isn't a filter...) +# TODO: actually, should these be global functions instead of operations? +DeclareOperation( "Vector", [IsOperation, IsSemiring, IsList]); +DeclareOperation( "Vector", [IsOperation, IsSemiring, IsVectorObj]); + +# Here we implement default choices for the representation, depending +# in base domain: +DeclareOperation( "Vector", [IsSemiring, IsList]); +DeclareOperation( "Vector", [IsSemiring, IsVectorObj]); + +# And here are the variants with example object (as last argument): +DeclareOperation( "Vector", [IsList, IsVectorObj]); +DeclareOperation( "Vector", [IsVectorObj, IsVectorObj]); + +# And here guess everything: +DeclareOperation( "Vector", [IsList]); + + # Creates a new vector in the same representation but with entries from list. # The length is given by the length of the first argument. # It is *not* guaranteed that the list is copied! @@ -284,6 +479,17 @@ DeclareOperation( "Vector", [IsList,IsVectorObj]); ## the matrix but of possibly different length given by the first ## argument. It is *not* guaranteed that the list is copied! +## <#GAPDoc Label="MatObj_ConstructingFilter_Vector"> +## +## +## a filter +## +## Returns a filter f such that if is +## called with f a vector in the same representation as V +## is produced. +## +## +## <#/GAPDoc> DeclareOperation( "ConstructingFilter", [IsVectorObj] ); DeclareConstructor( "NewVector", [IsVectorObj,IsSemiring,IsList] ); @@ -305,36 +511,63 @@ DeclareOperation( "ChangedBaseDomain", [IsVectorObj,IsSemiring] ); # Changes the base domain. A copy of the row vector in the first argument is # created, which comes in a "similar" representation but over the new # base domain that is given in the second argument. +# example: given a vector over GF(2), create a new vector over GF(4) with "identical" content +# so it's kind of a type conversion / coercion +# TODO: better name, e.g. VectorWithChangedBasedDomain +# or maybe just turn this into a constructor resp. a new constructor special case : +# like +# DeclareConstructor( "NewVector", [IsVectorObj,IsSemiring,IsVectorObj] ); + DeclareGlobalFunction( "MakeVector" ); # A convenience function for users to choose some appropriate representation # and guess the base domain if not supplied as second argument. # This is not guaranteed to be efficient and should never be used # in library or package code. +# usage: MakeVector( [,] ) +# TODO: explain that this is not something you should use a lot; +# instead you should use NewVector, ZeroVector, Vector, ... +# explain a typical use case / scenario for this function.. +# also: do we *really* need it? Max expects we'll find out as we convert the library + ############################################################################ # Some things that fit nowhere else: ############################################################################ +## <#GAPDoc Label="MatObj_Randomize_Vectors"> +## +## +## +## +## Replaces every entry in V with a random one from the base +## domain. If given, the random source Rs is used to compute the +## random elements. Note that in this case, the random function for the +## base domain must support the random source argument. +## +## +## <#/GAPDoc> DeclareOperation( "Randomize", [IsVectorObj and IsMutable] ); +DeclareOperation( "Randomize", [IsVectorObj and IsMutable,IsRandomSource] ); # Changes the mutable argument in place, every entry is replaced # by a random element from BaseDomain. -# The argument is also returned by the function. - -DeclareOperation( "Randomize", [IsVectorObj and IsMutable,IsRandomSource] ); -# The same, use the second argument to provide "randomness". +# The second argument is used to provide "randomness". # The vector argument is also returned by the function. +# TODO: change this to use InstallMethodWithRandomSource; for this, we'll have +# to change the argument order (a method for the old order, to ensure backwards +# compatibility, could remain). + ############################################################################# ## #O CopySubVector( , , , ) ## ## <#GAPDoc Label="CopySubVector"> ## -## -## +## ## -## returns nothing. Does dst{dcols} := src{scols} +## nothing +## Does dst{dcols} := src{scols} ## without creating an intermediate object and thus - at least in ## special cases - much more efficiently. For certain objects like ## compressed vectors this might be significantly more efficient if @@ -343,17 +576,46 @@ DeclareOperation( "Randomize", [IsVectorObj and IsMutable,IsRandomSource] ); ## ## <#/GAPDoc> ## +## TODO: link to ExtractSubVector +## +## TODO: In AddRowVector and MultRowVector, the destination is always the first argument; +## it would be better if we were consistent... +## +## TODO: Maybe also have a version of this as follows: +## CopySubVector( dst, dst_from, dst_to, src, src_form, src_to ); +DeclareOperation( "CopySubVector", + [IsVectorObj and IsMutable, IsList, IsVectorObj, IsList] ); +## TODO: the following declaration is deprecated and only kept for compatibility DeclareOperation( "CopySubVector", [IsVectorObj,IsVectorObj and IsMutable, IsList,IsList] ); +## <#GAPDoc Label="MatObj_WeightOfVector"> +## +## +## an integer +## +## Computes the Hamming weight of the vector V, i.e., the number of +## nonzero entries. +## +## +## <#/GAPDoc> DeclareOperation( "WeightOfVector", [IsVectorObj] ); -# This computes the Hamming weight of a vector, i.e. the number of -# nonzero entries. + +## <#GAPDoc Label="MatObj_DistanceOfVectors"> +## +## +## an integer +## +## Computes the Hamming distance of the vectors V1 and V2, +## i.e., the number of entries in which the vectors differ. The vectors +## must be of equal length. +## +## +## <#/GAPDoc> DeclareOperation( "DistanceOfVectors", [IsVectorObj, IsVectorObj] ); -# This computes the Hamming distance of two vectors, i.e. the number -# of positions, in which the vectors differ. The vectors must have the -# same length. + ############################################################################ @@ -367,7 +629,25 @@ DeclareOperation( "DistanceOfVectors", [IsVectorObj, IsVectorObj] ); # Attributes of matrices: ############################################################################ -# The following are guaranteed to be always set or cheaply calculable: +############################################################################# +## +#A BaseDomain( ) +## +## <#GAPDoc Label="MatObj_BaseDomain_IsMatrixObj"> +## +## +## +## +## The object matrix is defined over. +## Note that not all entries of matrix necessarily lie in +## BaseDomain( matrix ) with respect to in. +## E.g. matrix can be defined over a polynomial ring but also contain +## integers +## (see section ). +## +## +## <#/GAPDoc> +## DeclareAttribute( "BaseDomain", IsMatrixObj ); # Typically, the base domain will be a ring, it need not be commutative # nor associative. For non-associative base domains, the behavior of @@ -378,16 +658,27 @@ DeclareAttribute( "NumberColumns", IsMatrixObj ); DeclareSynonymAttr( "NrRows", NumberRows ); DeclareSynonymAttr( "NrCols", NumberColumns ); +DeclareAttribute( "OneOfBaseDomain", IsMatrixObj ); +DeclareAttribute( "ZeroOfBaseDomain", IsMatrixObj ); + + # HACK: Length and RowLength were in the old version of MatrixObj. We want # to get rid of them, but for now, keep thm to allow the few packages # already using this to continue working. DeclareAttribute( "Length", IsMatrixObj ); # ???? DeclareSynonymAttr( "RowLength", NumberColumns ); -DeclareAttribute( "DimensionsMat", IsMatrixObj ); # returns [rows,cols] + +# WARNING: the following attributes should not be stored if the matrix is mutable... +DeclareAttribute( "DimensionsMat", IsMatrixObj ); +# returns [NrRows(mat),NrCols(mat)] +# for backwards compatibility with existing code DeclareAttribute( "RankMat", IsMatrixObj ); DeclareOperation( "RankMatDestructive", [ IsMatrixObj ] ); +# TODO: danger: RankMat should not be stored for mutable matrices... +# + ############################################################################ # In the following sense matrices behave like lists: @@ -402,24 +693,18 @@ DeclareOperation( "[]", [IsMatrixObj,IsPosInt] ); # , # for matrices which are not row-lists. Efficient code will have to use MatElm and # SetMatElm instead. -# These should probably only be defined for RowListMatrices??? +# TODO: ... resp. it will use use M[i,j] +# TODO: provide a default method which creates a proxy object for the given row +# and translates accesses to it to corresponding MatElm / SetMatElm calls; +# creating such a proxy object prints an InfoWarning; +# but for the method for plist matrices, no warning is shown, as it is efficient +# anyway -DeclareOperation( "PositionNonZero", [IsMatrixObj] ); -DeclareOperation( "PositionNonZero", [IsMatrixObj, IsInt] ); +# TODO: maybe also add GetRow(mat, i) and GetColumn(mat, i) ??? +# these return IsVectorObj objects. +# these again must be objects which are "linked" to the original matrix, as above... +# TODO: perhaps also have ExtractRow(mat, i) and ExtractColumn(mat, i) -DeclareOperation( "PositionLastNonZero", [IsMatrixObj] ); -DeclareOperation( "PositionLastNonZero", [IsMatrixObj, IsInt] ); - -DeclareOperation( "Position", [IsMatrixObj, IsVectorObj] ); -DeclareOperation( "Position", [IsMatrixObj, IsVectorObj, IsInt] ); - -# This allows for usage of PositionSorted: -DeclareOperation( "PositionSortedOp", [IsMatrixObj, IsVectorObj] ); -DeclareOperation( "PositionSortedOp",[IsMatrixObj,IsVectorObj,IsFunction]); - -# I intentionally left out "PositionNot" here. - -# Note that "PositionSet" is a function only for lists. End of story. # Note that arbitrary matrices need not behave like lists with respect to the # following operations: @@ -452,8 +737,24 @@ DeclareOperation( "PositionSortedOp",[IsMatrixObj,IsVectorObj,IsFunction]); ## DeclareOperation( "ExtractSubMatrix", [IsMatrixObj,IsList,IsList] ); +# TODO: perhaps also add ExtractSubMatrix( mat, row_from, row_to, col_from, col_to ) ??? +# probably not needed... one can use ranges + TryNextMethod ... +# but let's look at some places where this function is used + # Creates a fully mutable copy of the matrix. DeclareOperation( "MutableCopyMat", [IsMatrixObj] ); +# so this amounts roughly to `mat -> List( mat, ShallowCopy )` +# except that it produce an object in the same representation + + +# TODO: perhaps also add a recursive version of ShallowCopy with a parameter limiting the depth.. + +# TODO: can we also have a version of ShallowCopy which has a name starting with "Mutable" +# to make it easier to discover??? Like MutableCopyVec (which might just be an alias +# for ShallowCopy). +# Or at least mention ShallowCopy in the MutableCopyMat + + ############################################################################# ## @@ -464,7 +765,9 @@ DeclareOperation( "MutableCopyMat", [IsMatrixObj] ); ## ## ## -## returns nothing. Does dst{drows}{dcols} := src{srows}{scols} +## nothing +## Does dst{drows}{dcols} := +## src{srows}{scols} ## without creating an intermediate object and thus - at least in ## special cases - much more efficiently. For certain objects like ## compressed vectors this might be significantly more efficient if @@ -476,6 +779,9 @@ DeclareOperation( "MutableCopyMat", [IsMatrixObj] ); DeclareOperation( "CopySubMatrix", [IsMatrixObj,IsMatrixObj, IsList,IsList,IsList,IsList] ); +# TODO: order of arguments? dst before src??? perhaps: +# dst, drows, dcols, src, srows, scols + ############################################################################ # New element access for matrices ############################################################################ @@ -496,9 +802,15 @@ DeclareOperation( "[]:=", [IsMatrixObj,IsPosInt,IsPosInt,IsObject] ); # The following are implicitly there for all objects, we mention them here # to have a complete interface description in one place: -# ShallowCopy is missing here since its behaviour depends on the matrix -# being in IsRowListMatrix or in IsFlatMatrix! +# ShallowCopy is missing here since its behaviour depends on the internal +# representation of the matrix objects (e.g. list-of-lists resp. list-of-vectors +# versus a "flat" matrix, or a shallow matrix, or ...)! +# TODO: for mutable matrices, a difference between StructuralCopy and MutableCopyMat +# occurs if the matrix is a row list, and one row is immutable, another mutable; +# then StructuralCopy will return a new list of vectors, which again contains the +# same (identical) immutable vectors, and mutable copies of the other vectors. +# Whereas with MutableCopyMat, all new rows will be mutable # DeclareGlobalFunction( "StructuralCopy", [IsMatrixObj] ); # DeclareOperation( "ViewObj", [IsMatrixObj] ); @@ -547,11 +859,33 @@ DeclareOperation( "[]:=", [IsMatrixObj,IsPosInt,IsPosInt,IsObject] ); # Problem: How about inverses of integer matrices that exist as # elements of rationals matrix? +# TODO: how to deal with this? If you have e.g. the integer matrix [[2]], +# then our generic inversion code produces the inverse [[1/2]] -- but that +# is not defined over the base domain (the integers). +# +# So perhaps the default method for inversion has to be careful... +# there is no problem over fields. +# For non-commutative or non-assoc domains, we don't provid a default. +# For commutative rings, we could compute the determinant alongside +# the inverse ("almost" for free), and then only need to check that +# the determinant is a unit in the base domain +# But be careful regarding zero divisors -- so perhaps use a different +# default method for non-fields, which avoids division +# +# Also, we can have other methods for e.g. subdomains of the cyclotomics. +# DeclareOperation( "TraceMat", [IsMatrixObj] ); # The sum of the diagonal entries. Error for a non-square matrix. + +#DeclareOperation( "DeterminantMat", [IsMatrixObj] ); +# TODO: this only makes sense over commutative domains; +# be careful regarding default implementation (base domain = field vs. +# base domain = any commutative ring, possibly with zero divisors) + + ############################################################################ # Rule: # Operations not sensibly defined return fail and do not trigger an error: @@ -559,6 +893,13 @@ DeclareOperation( "TraceMat", [IsMatrixObj] ); # One for non-square matrices. # Inverse for non-square matrices # Inverse for square, non-invertible matrices. + +# FIXME: what is the rationale for the above? OK, inverting a square matrix: +# it is non-trivial to decide if the matrix is invertible, so it might be +# useful to be able to just try and invert it, and do something else if that "fail"s +# But it is cheap to verify that a matrix is non-square, so why not error for that? +# Same for One, ... + # # An exception are properties: # IsOne for non-square matrices returns false. @@ -573,15 +914,21 @@ DeclareOperation( "TraceMat", [IsMatrixObj] ); ############################################################################ DeclareOperation( "ZeroMatrix", [IsInt,IsInt,IsMatrixObj] ); +DeclareOperation( "ZeroMatrix", [IsSemiring, IsInt,IsInt] ); # warning: NullMat has ring arg last +DeclareOperation( "ZeroMatrix", [IsOperation, IsSemiring, IsInt,IsInt] ); # Returns a new fully mutable zero matrix in the same rep as the given one with # possibly different dimensions. First argument is number of rows, second # is number of columns. DeclareConstructor( "NewZeroMatrix", [IsMatrixObj,IsSemiring,IsInt,IsInt]); +# constructor -> first argument must be a filter, like e.g. IsPlistMatrixRep +# # Returns a new fully mutable zero matrix over the base domain in the # 2nd argument. The integers are the number of rows and columns. DeclareOperation( "IdentityMatrix", [IsInt,IsMatrixObj] ); +DeclareOperation( "IdentityMatrix", [IsSemiring, IsInt] ); # warning: IdentityMat has ring arg last +DeclareOperation( "IdentityMatrix", [IsOperation, IsSemiring, IsInt] ); # Returns a new mutable identity matrix in the same rep as the given one with # possibly different dimensions. @@ -589,6 +936,15 @@ DeclareConstructor( "NewIdentityMatrix", [IsMatrixObj,IsSemiring,IsInt]); # Returns a new fully mutable identity matrix over the base domain in the # 2nd argument. The integer is the number of rows and columns. +# TODO: perhaps imitate what we do for e.g. group constructors, and allow +# user to omit the filter; in that case, try to choose a "good" default +# representation ???? + +# TODO: perhaps add DiagonalMatrix? + +# TODO: convert (New)IdentityMatrix and (New)ZeroMatrix to be more similar to Matrix() + + DeclareOperation( "CompanionMatrix", [IsUnivariatePolynomial,IsMatrixObj] ); # Returns the companion matrix of the first argument in the representation # of the second argument. Uses row-convention. The polynomial must be @@ -597,54 +953,105 @@ DeclareOperation( "CompanionMatrix", [IsUnivariatePolynomial,IsMatrixObj] ); DeclareConstructor( "NewCompanionMatrix", [IsMatrixObj, IsUnivariatePolynomial, IsSemiring] ); # The constructor variant of . +# TODO: get rid of NewCompanionMatrix, at least for now -- if somebody *REALLY* +# needs it, we can still reconsider... Instead allow this: +# DeclareOperation( "CompanionMatrix", [IsFilter, IsUnivariatePolynomial, IsSemiring] ); +# which roughly does this: +# InstallMethod( CompanionMatrix, ... +# function(filt, f, R) +# n := Degree(f); +# mat := NewZeroMatrix(filt, R, n, n); +# ... set entries of mat ... +# end); + + # The following are already declared in the library: # Eventually here will be the right place to do this. -DeclareOperation( "Matrix", [IsList,IsInt,IsMatrixObj]); -# Creates a new matrix in the same representation as the fourth argument -# but with entries from list, the second argument is the number of -# columns. The first argument can be: -# - a plain list of vectors of the correct row length in a representation -# fitting to the matrix rep. -# - a plain list of plain lists where each sublist has the length of the rows -# - a plain list with length rows*cols with matrix entries given row-wise -# If the first argument is empty, then the number of rows is zero. -# Otherwise the first entry decides which case is given. -# The outer list is guaranteed to be copied, however, the entries of that -# list (the rows) need not be copied. -# The following convenience versions exist: -# With two arguments the first must not be empty and must not be a flat -# list. Then the number of rows is deduced from the length of the first -# argument and the number of columns is deduced from the length of the -# element of the first argument (done with a generic method): -DeclareOperation( "Matrix", [IsList,IsMatrixObj] ); +# variant with new filter + base domain (dispatches to NewMatrix) +DeclareOperation( "Matrix", [IsOperation, IsSemiring, IsList, IsInt]); +DeclareOperation( "Matrix", [IsOperation, IsSemiring, IsList]); +DeclareOperation( "Matrix", [IsOperation, IsSemiring, IsMatrixObj]); + +# variant with new base domain -> "guesses" good rep, then dispatches to NewMatrix +DeclareOperation( "Matrix", [IsSemiring, IsList, IsInt]); +DeclareOperation( "Matrix", [IsSemiring, IsList]); +DeclareOperation( "Matrix", [IsSemiring, IsMatrixObj]); + +# the following two operations use DefaultFieldOfMatrix to "guess" the base domain +DeclareOperation( "Matrix", [IsList, IsInt]); +DeclareAttribute( "Matrix", IsList, "mutable"); # HACK: because there already is an attribute Matrix + +# variant with example object at end (input is first) +DeclareOperation( "Matrix", [IsList, IsInt, IsMatrixObj]); +DeclareOperation( "Matrix", [IsList, IsMatrixObj]); # <- no need to overload this one +DeclareOperation( "Matrix", [IsMatrixObj, IsMatrixObj]); + +# perhaps also (or instead?) have this: +#DeclareOperation( "MatrixWithRows", [IsList (of vectors),IsMatrixObj]); ?? +#DeclareOperation( "MatrixWithColumns", [IsList (of vectors),IsMatrixObj]); ?? + -# Note that it is not possible to generate a matrix via "Matrix" without -# a template matrix object. Use the constructor methods instead: DeclareConstructor( "NewMatrix", [IsMatrixObj, IsSemiring, IsInt, IsList] ); # Constructs a new fully mutable matrix. The first argument has to be a filter -# indicating the representation. The second the base domain, the third -# the row length and the last a list containing either row vectors -# of the right length or lists with base domain elements. +# indicating the representation, the second the base domain, the third +# the row length. The last argument can be: +# - a plain list of vector objects of correct length +# - a plain list of plain lists of correct length +# - a flat plain list with rows*cols entries in row major order +# (FoldList turns a flat list into a list of lists) +# where the corresponding entries must be in or compatible with the base domain. +# If the last argument already contains vector objects, they are copied. # The last argument is guaranteed not to be changed! -# If the last argument already contains row vectors, they are copied. +# TODO: what does "flat" above mean??? +# TODO: from an old comment on Matrix / NewMatrix, wrt to the last argument +# The outer list is guaranteed to be copied, however, the entries of that +# list (the rows) need not be copied. +# TODO: Isn't it inconsistent to copy the rows if they are vector objects, but otherwise not? Also: matobjplist.gi uses NewVector, which copies its list-argument + +# FIXME: why is IsInt,IsList reversed compared to Matrix(), where it is IsList,IsInt + + + + +# given a matrix , produce a filter such that NewMatrix called with this filter +# will produce a matrix in the same representation as DeclareOperation( "ConstructingFilter", [IsMatrixObj] ); +# TODO: what does this do? +# Implementation: given an n x m matrix , create a new zero vector of +# length n (= NrRows) in a representation "compatible" with that of , i.e. +# there "should be" a fast action * +# FIXME: is this really useful? Compare to, say, `ExtractRow(mat,1)` etc. ??? DeclareOperation( "CompatibleVector", [IsMatrixObj] ); DeclareOperation( "ChangedBaseDomain", [IsMatrixObj,IsSemiring] ); # Changes the base domain. A copy of the matrix in the first argument is # created, which comes in a "similar" representation but over the new # base domain that is given in the second argument. +# TODO: better name, e.g. MatrixWithChangedBasedDomain +# or maybe just turn this into a constructor resp. a new constructor special case : +# like +# DeclareConstructor( "NewMatrix", [IsMatrixObj,IsSemiring,IsMatrixObj] ); + + + -DeclareGlobalFunction( "MakeMatrix" ); +# usage: MakeVector( [,] ) # A convenience function for users to choose some appropriate representation # and guess the base domain if not supplied as second argument. # This is not guaranteed to be efficient and should never be used # in library or package code. +# It is mainly useful to help migrate existing code incrementally to use +# the new MatrixObj interface. + + +# TODO: how useful are all these constructors in practice? +# Ideally we should try to limit their number, focusing on a few useful ones.. +# best to see what actual code needs ############################################################################ @@ -657,17 +1064,54 @@ DeclareOperation( "Randomize", [IsMatrixObj and IsMutable,IsRandomSource] ); # by a random element from BaseDomain. # The second version will come when we have random sources. +# TODO: only keep the first operation and suggest using InstallMethodWithRandomSource + + DeclareAttribute( "TransposedMatImmutable", IsMatrixObj ); DeclareOperation( "TransposedMatMutable", [IsMatrixObj] ); +# TODO: one problem with DeclareAttribute: it only really makes sense if the +# matrix is immutable; but we have no way of declaring this; so it is very +# easy to end up having a mutable, attribute storing matrix, which stores +# its transpose "by accident", and then gets modified later on +# -DeclareOperation( "IsDiagonalMat", [IsMatrixObj] ); +DeclareOperation( "IsDiagonalMat", [IsMatrixObj] ); DeclareOperation( "IsUpperTriangularMat", [IsMatrixObj] ); DeclareOperation( "IsLowerTriangularMat", [IsMatrixObj] ); +# TODO: if we allow attributes, we might just as well do the above to be +# declared as properties, so that this information is stored; but once +# again, we would only want to allow this for immutable matrix objects. +# ... + +# TODO: what about the following (and also note the names...): +# - IsScalarMat, IsSquareMat, ... ? + +# TODO: Why do we call them RankMat, IsDiagonalMat, etc. +# and not RankMatrix, IsDiagonalMatrix, etc. ? +# in contrast to Matrix(), NewMatrix(), ExtractSubMatrix(), ... +# +# One reasons is backwards compatibility of course, but this could +# also be achieved with synonyms. +# Another argument for "Mat" is brevity, but is that really so important/useful +# versus "clarity" +# Still, we could unify these names, always using "Matrix". +# +# Of course we could also introduce a rule when to use "Matrix" and when to use "Mat", +# but this rule will invariably be ignored or forgotten and inconsistency will +# occur. +# +# So new plan: Always use "Matrix", but provide "Mat" aliases for +# backwards compatibility (if the operation already exists with that name), +# and perhaps (????) for "convenience" +# DeclareOperation( "KroneckerProduct", [IsMatrixObj,IsMatrixObj] ); # The result is fully mutable. +# FIXME: what is the purpose of Unfold and Fold? +# Maybe remove them, and wait if somebody asks for / needs this + DeclareOperation( "Unfold", [IsMatrixObj, IsVectorObj] ); # Concatenates all rows of a matrix to one single vector in the same # representation as the given template vector. Usually this must @@ -682,6 +1126,7 @@ DeclareOperation( "Fold", [IsVectorObj, IsPosInt, IsMatrixObj] ); DeclareOperation( "Unpack", [IsRowListMatrix] ); # It guarantees to copy, that is changing the returned object does # not change the original object. +# TODO: this is already used by the fining package ############################################################################ @@ -695,6 +1140,17 @@ DeclareOperation( "Unpack", [IsRowListMatrix] ); # List operations with some restrictions: ############################################################################ +# TODO: what is this good for? Theory: it would probably help when migrating +# existing code to MatrixObj, as you could change your lists-of-list matrices +# to MatrixObjs, and then incrementally adapt your code to *not* use the +# following operations. +# +# In that case, we should make it very clear that using these functions is +# discouraged.... + +# TODO: let's see if this is really useful when e.g. adapting the library. +# If not, then we might not even need `IsRowListMatrix` + DeclareOperation( "[]:=", [IsRowListMatrix,IsPosInt,IsObject] ); # Only guaranteed to work for the position in [1..Length(VECTOR)] and only # for elements in a suitable vector type. @@ -745,7 +1201,11 @@ DeclareOperation( "ListOp", [IsRowListMatrix, IsFunction] ); ############################################################################ # DeclareOperation( "*", [IsVectorObj, IsMatrixObj] ); +# DeclareOperation( "*", [IsMatrixObj, IsVectorObj] ); +# TODO: the following does the same as "*", but is useful +# as convenience for Orbit-ish operations. +# We otherwise discourage its use in "naked" code -- use * instead # DeclareOperation( "^", [IsVectorObj, IsMatrixObj] ); @@ -756,3 +1216,22 @@ DeclareOperation( "ListOp", [IsRowListMatrix, IsFunction] ); # AsList # AddCoeffs + + + +# while we are at it also make the naming of following more uniform ??? + +# TriangulizeIntegerMat +# TriangulizedIntegerMat(Transform) +# HermiteNormalFormIntegerMat(Transform) +# SmithNormalFormIntegerMat(Transform) +# and contrast them to these: +# BaseIntMat +# BaseIntersectionIntMats +# ComplementIntMat +# DeterminantIntMat +# DiagonalizeIntMat +# NormalFormIntMat +# NullspaceIntMat +# SolutionIntMat +# SolutionNullspaceIntMat diff --git a/lib/matobjplist.gd b/lib/matobjplist.gd index 079eb36c73..646f0af260 100644 --- a/lib/matobjplist.gd +++ b/lib/matobjplist.gd @@ -41,10 +41,6 @@ BindGlobal( "ELSPOS", 2 ); DeclareFilter( "IsIntVector" ); DeclareFilter( "IsFFEVector" ); -# Another pair of filters that slow down things: -DeclareFilter( "IsCheckingVector" ); -DeclareFilter( "IsCheckingMatrix" ); - ############################################################################ # Constructors: ############################################################################ diff --git a/lib/matobjplist.gi b/lib/matobjplist.gi index 1baecceec8..435cca51ac 100644 --- a/lib/matobjplist.gi +++ b/lib/matobjplist.gi @@ -17,9 +17,10 @@ InstallGlobalFunction( MakePlistVectorType, function( basedomain, filter ) local T, filter2; + filter2 := filter and IsMutable; if HasCanEasilyCompareElements(Representative(basedomain)) and CanEasilyCompareElements(Representative(basedomain)) then - filter2 := filter and IsMutable and CanEasilyCompareElements; + filter2 := filter2 and CanEasilyCompareElements; fi; if IsIdenticalObj(basedomain,Integers) then T := NewType(FamilyObj(basedomain), @@ -35,20 +36,20 @@ InstallGlobalFunction( MakePlistVectorType, end); InstallMethod( NewVector, "for IsPlistVectorRep, a ring, and a list", - [ IsPlistVectorRep and IsCheckingVector, IsRing, IsList ], + [ IsPlistVectorRep, IsRing, IsList ], function( filter, basedomain, l ) local typ, v; - typ := MakePlistVectorType(basedomain,filter); + typ := MakePlistVectorType(basedomain,IsPlistVectorRep); v := [basedomain,ShallowCopy(l)]; Objectify(typ,v); return v; end ); InstallMethod( NewZeroVector, "for IsPlistVectorRep, a ring, and an int", - [ IsPlistVectorRep and IsCheckingVector, IsRing, IsInt ], + [ IsPlistVectorRep, IsRing, IsInt ], function( filter, basedomain, l ) local typ, v; - typ := MakePlistVectorType(basedomain,filter); + typ := MakePlistVectorType(basedomain,IsPlistVectorRep); v := [basedomain,Zero(basedomain)*[1..l]]; Objectify(typ,v); return v; @@ -56,25 +57,30 @@ InstallMethod( NewZeroVector, "for IsPlistVectorRep, a ring, and an int", InstallMethod( NewMatrix, "for IsPlistMatrixRep, a ring, an int, and a list", - [ IsPlistMatrixRep and IsCheckingMatrix, IsRing, IsInt, IsList ], + [ IsPlistMatrixRep, IsRing, IsInt, IsList ], function( filter, basedomain, rl, l ) - local m,i,e,filter2; - if IsIdenticalObj(IsPlistMatrixRep,filter) then - filter2 := IsPlistVectorRep; - else - filter2 := IsPlistVectorRep and IsCheckingVector; - fi; + local nd, filterVectors, m, e, filter2, i; + # check if l is flat list + if Length(l) > 0 and not IsVectorObj(l[1]) then + nd := NestingDepthA(l); + if nd < 2 or nd mod 2 = 1 then + l := List([0,rl..Length(l)-rl], i -> l{[i+1..i+rl]}); + elif Length(l) mod rl <> 0 then + Error( "NewMatrix: Length of l is not a multiple of rl" ); + fi; + fi; + filterVectors := IsPlistVectorRep; m := 0*[1..Length(l)]; - e := NewVector(filter2, basedomain, []); for i in [1..Length(l)] do if IsVectorObj(l[i]) and IsPlistVectorRep(l[i]) then m[i] := ShallowCopy(l[i]); else - m[i] := Vector( l[i], e ); + m[i] := NewVector( filterVectors, basedomain, l[i] ); fi; od; + e := NewVector(filterVectors, basedomain, []); m := [basedomain,e,rl,m]; - filter2 := filter and IsMutable; + filter2 := IsPlistMatrixRep and IsMutable; if HasCanEasilyCompareElements(Representative(basedomain)) and CanEasilyCompareElements(Representative(basedomain)) then filter2 := filter2 and CanEasilyCompareElements; @@ -86,14 +92,10 @@ InstallMethod( NewMatrix, InstallMethod( NewZeroMatrix, "for IsPlistMatrixRep, a ring, and two ints", - [ IsPlistMatrixRep and IsCheckingMatrix, IsRing, IsInt, IsInt ], + [ IsPlistMatrixRep, IsRing, IsInt, IsInt ], function( filter, basedomain, rows, cols ) local m,i,e,filter2; - if IsIdenticalObj(IsPlistMatrixRep,filter) then - filter2 := IsPlistVectorRep; - else - filter2 := IsPlistVectorRep and IsCheckingVector; - fi; + filter2 := IsPlistVectorRep; m := 0*[1..rows]; e := NewVector(filter2, basedomain, []); for i in [1..rows] do @@ -107,16 +109,12 @@ InstallMethod( NewZeroMatrix, InstallMethod( NewIdentityMatrix, "for IsPlistMatrixRep, a ring, and an int", - [ IsPlistMatrixRep and IsCheckingMatrix, IsRing, IsInt ], + [ IsPlistMatrixRep, IsRing, IsInt ], function( filter, basedomain, rows ) - local m,i,e,filter2; - if IsIdenticalObj(IsPlistMatrixRep,filter) then - filter2 := IsPlistVectorRep; - else - filter2 := IsPlistVectorRep and IsCheckingVector; - fi; + local filterVectors, m, e, i; + filterVectors := IsPlistVectorRep; m := 0*[1..rows]; - e := NewVector(IsPlistVectorRep, basedomain, []); + e := NewVector(filterVectors, basedomain, []); for i in [1..rows] do m[i] := ZeroVector( rows, e ); m[i][i] := One(basedomain); @@ -139,16 +137,12 @@ InstallMethod( ViewObj, "for a plist vector", [ IsPlistVectorRep ], else Print("<"); fi; - if IsCheckingVector(v) then Print("checking "); fi; Print("plist vector over ",v![BDPOS]," of length ",Length(v![ELSPOS]),">"); end ); InstallMethod( PrintObj, "for a plist vector", [ IsPlistVectorRep ], function( v ) Print("NewVector(IsPlistVectorRep"); - if IsCheckingVector(v) then - Print(" and IsCheckingVector"); - fi; if IsFinite(v![BDPOS]) and IsField(v![BDPOS]) then Print(",GF(",Size(v![BDPOS]),"),",v![ELSPOS],")"); else @@ -160,9 +154,6 @@ InstallMethod( String, "for a plist vector", [ IsPlistVectorRep ], function( v ) local st; st := "NewVector(IsPlistVectorRep"; - if IsCheckingVector(v) then - Append(st," and IsCheckingVector"); - fi; if IsFinite(v![BDPOS]) and IsField(v![BDPOS]) then Append(st,Concatenation( ",GF(",String(Size(v![BDPOS])),"),", String(v![ELSPOS]),")" )); @@ -176,7 +167,6 @@ InstallMethod( String, "for a plist vector", [ IsPlistVectorRep ], InstallMethod( Display, "for a plist vector", [ IsPlistVectorRep ], function( v ) Print( "\n"); end ); @@ -262,19 +252,6 @@ InstallMethod( \[\], "for a plist vector and a positive integer", return v![ELSPOS][p]; end ); -InstallMethod( \[\]\:\=, - "for a checking plist vector, a positive integer, and an obj", - [ IsPlistVectorRep and IsCheckingVector, IsPosInt, IsObject ], - function( v, p, ob ) - if p > Length(v![ELSPOS]) then - ErrorNoReturn("\\[\\]\\:\\=: Assignment out of bounds!"); - fi; - if not ob in v![BDPOS] then - ErrorNoReturn("\\[\\]\\:\\=: Object to be assigned is not in base domain"); - fi; - v![ELSPOS][p] := ob; - end ); - InstallMethod( \[\]\:\=, "for a plist vector, a positive integer, and an obj", [ IsPlistVectorRep, IsPosInt, IsObject ], function( v, p, ob ) @@ -360,28 +337,6 @@ InstallMethod( \+, "for two plist vectors", [a![BDPOS],SUM_LIST_LIST_DEFAULT(a![ELSPOS],b![ELSPOS])]); end ); -InstallMethod( \+, "for two checking plist vectors", - [ IsPlistVectorRep and IsCheckingVector, - IsPlistVectorRep and IsCheckingVector], - function( a, b ) - local ty; - if Length(a![ELSPOS]) <> Length(b![ELSPOS]) then - Error("\\+: Cannot add vectors of different length"); - return fail; - fi; - if not IsIdenticalObj(a![BDPOS],b![BDPOS]) then - Error("\\+: Cannot add vectors over different base domains"); - return fail; - fi; - if not IsMutable(a) and IsMutable(b) then - ty := TypeObj(b); - else - ty := TypeObj(a); - fi; - return Objectify(ty, - [a![BDPOS],SUM_LIST_LIST_DEFAULT(a![ELSPOS],b![ELSPOS])]); - end ); - InstallMethod( \-, "for two plist vectors", [ IsPlistVectorRep, IsPlistVectorRep ], function( a, b ) @@ -395,28 +350,6 @@ InstallMethod( \-, "for two plist vectors", [a![BDPOS],DIFF_LIST_LIST_DEFAULT(a![ELSPOS],b![ELSPOS])]); end ); -InstallMethod( \-, "for two checking plist vectors", - [ IsPlistVectorRep and IsCheckingVector, - IsPlistVectorRep and IsCheckingVector], - function( a, b ) - local ty; - if Length(a![ELSPOS]) <> Length(b![ELSPOS]) then - Error("\\-: Cannot subtract vectors of different length"); - return fail; - fi; - if not IsIdenticalObj(a![BDPOS],b![BDPOS]) then - Error("\\-: Cannot subtract vectors over different base domains"); - return fail; - fi; - if not IsMutable(a) and IsMutable(b) then - ty := TypeObj(b); - else - ty := TypeObj(a); - fi; - return Objectify(ty, - [a![BDPOS],DIFF_LIST_LIST_DEFAULT(a![ELSPOS],b![ELSPOS])]); - end ); - InstallMethod( \=, "for two plist vectors", [ IsPlistVectorRep, IsPlistVectorRep ], function( a, b ) @@ -435,19 +368,6 @@ InstallMethod( AddRowVector, "for two plist vectors", ADD_ROW_VECTOR_2( a![ELSPOS], b![ELSPOS] ); end ); -InstallMethod( AddRowVector, "for two checking plist vectors", - [ IsPlistVectorRep and IsCheckingVector and IsMutable, - IsPlistVectorRep and IsCheckingVector ], - function( a, b ) - if Length(a![ELSPOS]) <> Length(b![ELSPOS]) then - ErrorNoReturn("AddRowVector: Cannot add vectors of different length"); - fi; - if not IsIdenticalObj(a![BDPOS],b![BDPOS]) then - ErrorNoReturn("AddRowVector: Cannot add vectors over different base domains"); - fi; - ADD_ROW_VECTOR_2( a![ELSPOS], b![ELSPOS] ); - end ); - # Better method for integer vectors: InstallMethod( AddRowVector, "for two plist vectors", [ IsPlistVectorRep and IsMutable and IsIntVector, @@ -462,19 +382,6 @@ InstallMethod( AddRowVector, "for two plist vectors, and a scalar", ADD_ROW_VECTOR_3( a![ELSPOS], b![ELSPOS], s ); end ); -InstallMethod( AddRowVector, "for two checking plist vectors, and a scalar", - [ IsPlistVectorRep and IsCheckingVector and IsMutable, - IsPlistVectorRep and IsCheckingVector, IsObject ], - function( a, b, s ) - if Length(a![ELSPOS]) <> Length(b![ELSPOS]) then - ErrorNoReturn("AddRowVector: Cannot subtract vectors of different length"); - fi; - if not IsIdenticalObj(a![BDPOS],b![BDPOS]) then - ErrorNoReturn("AddRowVector: Cannot add vectors over different base domains"); - fi; - ADD_ROW_VECTOR_3( a![ELSPOS], b![ELSPOS], s ); - end ); - # Better method for integer vectors: InstallMethod( AddRowVector, "for two plist vectors, and a scalar", [ IsPlistVectorRep and IsIntVector and IsMutable, @@ -495,20 +402,6 @@ InstallMethod( AddRowVector, ADD_ROW_VECTOR_5( a![ELSPOS], b![ELSPOS], s, from, to ); end ); -InstallMethod( AddRowVector, - "for two checking plist vectors, a scalar, and two positions", - [ IsPlistVectorRep and IsCheckingVector and IsMutable, - IsPlistVectorRep and IsCheckingVector, IsObject, IsPosInt, IsPosInt ], - function( a, b, s, from, to ) - if Length(a![ELSPOS]) <> Length(b![ELSPOS]) then - ErrorNoReturn("AddRowVector: Cannot add vectors of different length"); - fi; - if not IsIdenticalObj(a![BDPOS],b![BDPOS]) then - ErrorNoReturn("AddRowVector: Cannot add vectors over different base domains"); - fi; - ADD_ROW_VECTOR_5( a![ELSPOS], b![ELSPOS], s, from, to ); - end ); - # Better method for integer vectors: InstallMethod( AddRowVector, "for two integer plist vectors, a scalar, and two positions", @@ -528,15 +421,6 @@ InstallMethod( MultRowVector, "for a plist vector, and a scalar", MULT_ROW_VECTOR_2(v![ELSPOS],s); end ); -InstallMethod( MultRowVector, "for a checking plist vector, and a scalar", - [ IsPlistVectorRep and IsCheckingVector and IsMutable, IsObject ], - function( v, s ) - if not s in v![BDPOS] then - ErrorNoReturn("MultRowVector: Scalar is not in base domain"); - fi; - MULT_ROW_VECTOR_2(v![ELSPOS],s); - end ); - InstallMethod( MultRowVector, "for a plist vector, and a scalar", [ IsPlistVectorRep and IsIntVector and IsMutable, IsObject ], function( v, s ) @@ -555,25 +439,6 @@ InstallMethod( MultRowVector, a![ELSPOS]{pa} := b![ELSPOS]{pb} * s; end ); -InstallMethod( MultRowVector, - "for a checking plist vector, a list, a ch. plist vector, a list, a scalar", - [ IsPlistVectorRep and IsCheckingVector and IsMutable, IsList, - IsPlistVectorRep and IsCheckingVector, IsList, IsObject ], - function( a, pa, b, pb, s ) - local l; - l := Length(a![ELSPOS]); - if not ForAll(pa,x->x >= 1 and x <= l) then - ErrorNoReturn("MultRowVector: Positions pa must lie in first vector"); - fi; - if not IsIdenticalObj(a![BDPOS],b![BDPOS]) then - ErrorNoReturn("MultRowVector: Cannot add vectors over different base domains"); - fi; - if not s in a![BDPOS] then - ErrorNoReturn("MultRowVector: Scalar not in base domain"); - fi; - a![ELSPOS]{pa} := b![ELSPOS]{pb} * s; - end ); - InstallMethod( \*, "for a plist vector and a scalar", [ IsPlistVectorRep, IsScalar ], function( v, s ) @@ -581,16 +446,6 @@ InstallMethod( \*, "for a plist vector and a scalar", [v![BDPOS],PROD_LIST_SCL_DEFAULT(v![ELSPOS],s)] ); end ); -InstallMethod( \*, "for a checking plist vector and a scalar", - [ IsPlistVectorRep, IsScalar ], - function( v, s ) - if not s in v![BDPOS] then - ErrorNoReturn("\\*: Scalar not in base domain"); - fi; - return Objectify( TypeObj(v), - [v![BDPOS],PROD_LIST_SCL_DEFAULT(v![ELSPOS],s)] ); - end ); - InstallMethod( \*, "for a scalar and a plist vector", [ IsScalar, IsPlistVectorRep ], function( s, v ) @@ -598,16 +453,6 @@ InstallMethod( \*, "for a scalar and a plist vector", [v![BDPOS],PROD_SCL_LIST_DEFAULT(s,v![ELSPOS])] ); end ); -InstallMethod( \*, "for a scalar and a checking plist vector", - [ IsScalar, IsPlistVectorRep ], - function( s, v ) - if not s in v![BDPOS] then - ErrorNoReturn("\\*: Scalar not in base domain"); - fi; - return Objectify( TypeObj(v), - [v![BDPOS],PROD_SCL_LIST_DEFAULT(s,v![ELSPOS])] ); - end ); - InstallMethod( \/, "for a plist vector and a scalar", [ IsPlistVectorRep, IsScalar ], function( v, s ) @@ -615,21 +460,6 @@ InstallMethod( \/, "for a plist vector and a scalar", [v![BDPOS],PROD_LIST_SCL_DEFAULT(v![ELSPOS],s^-1)] ); end ); -InstallMethod( \/, "for a checking plist vector and a scalar", - [ IsPlistVectorRep and IsCheckingVector, IsScalar ], - function( v, s ) - local res; - if not s in v![BDPOS] then - ErrorNoReturn("\\/: Scalar not in base domain"); - fi; - res := Objectify( TypeObj(v), - [v![BDPOS],PROD_LIST_SCL_DEFAULT(v![ELSPOS],s^-1)] ); - if IsIntVector(v) and not ForAll(res![ELSPOS],IsInt) then - ErrorNoReturn("\\/: Result is not an integer vector"); - fi; - return res; - end ); - InstallMethod( AdditiveInverseSameMutability, "for a plist vector", [ IsPlistVectorRep ], function( v ) @@ -701,16 +531,6 @@ InstallMethod( CopySubVector, "for two plist vectors and two lists", b![ELSPOS]{pb} := a![ELSPOS]{pa}; end ); -InstallMethod( CopySubVector, "for two plist vectors and two lists", - [ IsPlistVectorRep, IsPlistVectorRep and IsMutable and IsCheckingVector, - IsList, IsList ], - function( a,b,pa,pb ) - if not(ForAll(pb,x->x <= Length(b![ELSPOS]))) then - ErrorNoReturn("CopySubVector: Positions in pb must lie in vector b"); - fi; - b![ELSPOS]{pb} := a![ELSPOS]{pa}; - end ); - ############################################################################ ############################################################################ @@ -833,23 +653,6 @@ InstallMethod( \[\]\:\=, m![ROWSPOS][p] := v; end ); -InstallMethod( \[\]\:\=, - "for a checking plist matrix, a positive integer, and a plist vector", - [ IsPlistMatrixRep and IsCheckingMatrix and IsMutable, - IsPosInt, IsPlistVectorRep ], - function( m, p, v ) - if p > Length(m![ROWSPOS])+1 then - ErrorNoReturn("\\[\\]\\:\\=: Matrices must be dense, you cannot assign here"); - fi; - if not IsIdenticalObj(m![BDPOS],v![BDPOS]) then - ErrorNoReturn("\\[\\]\\:\\=: Vector and matrix must be over the same base domain"); - fi; - if not m![RLPOS] = Length(v![ELSPOS]) then - ErrorNoReturn("\\[\\]\\:\\=: Vector does not have the right length"); - fi; - m![ROWSPOS][p] := v; - end ); - InstallMethod( \{\}, "for a plist matrix and a list", [ IsPlistMatrixRep, IsList ], function( m, p ) @@ -864,41 +667,12 @@ InstallMethod( Add, "for a plist matrix and a plist vector", Add(m![ROWSPOS],v); end ); -InstallMethod( Add, "for a checking plist matrix and a plist vector", - [ IsPlistMatrixRep and IsCheckingMatrix and IsMutable, - IsPlistVectorRep ], - function( m, v ) - if not IsIdenticalObj(m![BDPOS],v![BDPOS]) then - ErrorNoReturn("Add: Vector and matrix must be over the same base domain"); - fi; - if not m![RLPOS] = Length(v![ELSPOS]) then - ErrorNoReturn("Add: Vector does not have the right length"); - fi; - Add(m![ROWSPOS],v); - end ); - InstallMethod( Add, "for a plist matrix, a plist vector, and a pos. int", [ IsPlistMatrixRep and IsMutable, IsPlistVectorRep, IsPosInt ], function( m, v, p ) Add(m![ROWSPOS],v,p); end ); -InstallMethod( Add, "for a checking plist matrix, a plist vector, and a pos", - [ IsPlistMatrixRep and IsMutable and IsCheckingMatrix, - IsPlistVectorRep, IsPosInt ], - function( m, v, p ) - if p > Length(m![ROWSPOS])+1 then - ErrorNoReturn("Add: Matrices must be dense, you cannot assign here"); - fi; - if not IsIdenticalObj(m![BDPOS],v![BDPOS]) then - ErrorNoReturn("Add: Vector and matrix must be over the same base domain"); - fi; - if not m![RLPOS] = Length(v![ELSPOS]) then - ErrorNoReturn("Add: Vector does not have the right length"); - fi; - Add(m![ROWSPOS],v,p); - end ); - InstallMethod( Remove, "for a plist matrix", [ IsPlistMatrixRep and IsMutable ], function( m ) @@ -933,42 +707,12 @@ InstallMethod( \{\}\:\=, "for a plist matrix, a list, and a plist matrix", m![ROWSPOS]{pp} := n![ROWSPOS]; end ); -InstallMethod( \{\}\:\=, - "for a checking plist matrix, a list, and a plist matrix", - [ IsPlistMatrixRep and IsMutable and IsCheckingMatrix, IsList, - IsPlistMatrixRep ], - function( m, pp, n ) - if not(ForAll(pp,x->x >= 1 and x <= Length(m![ROWSPOS])+1)) then - ErrorNoReturn("\\{\\}\\:\\=: Positions must be in the matrix"); - fi; - if not IsIdenticalObj(m![BDPOS],n![BDPOS]) then - ErrorNoReturn("\\{\\}\\:\\=: Both matrices must be over the same base domain"); - fi; - if not m![RLPOS] = n![RLPOS] then - ErrorNoReturn("\\{\\}\\:\\=: Row lengths are not equal"); - fi; - m![ROWSPOS]{pp} := n![ROWSPOS]; - end ); - InstallMethod( Append, "for two plist matrices", [ IsPlistMatrixRep and IsMutable, IsPlistMatrixRep ], function( m, n ) Append(m![ROWSPOS],n![ROWSPOS]); end ); -InstallMethod( Append, "for a checking plist matrix, and a plist matrix", - [ IsPlistMatrixRep and IsMutable and IsCheckingMatrix, - IsPlistMatrixRep ], - function( m, n ) - if not IsIdenticalObj(m![BDPOS],n![BDPOS]) then - ErrorNoReturn("Append: Both matrices must be over the same base domain"); - fi; - if not m![RLPOS] = n![RLPOS] then - ErrorNoReturn("Append: Row lengths are not equal"); - fi; - Append(m![ROWSPOS],n![ROWSPOS]); - end ); - InstallMethod( ShallowCopy, "for a plist matrix", [ IsPlistMatrixRep ], function( m ) @@ -1005,100 +749,6 @@ InstallMethod( Unpack, "for a plist matrix", return List(m![ROWSPOS],v->ShallowCopy(v![ELSPOS])); end ); -InstallMethod( PositionNonZero, "for a plist matrix", - [ IsPlistMatrixRep ], - function( m ) - local i,l,le; - l := m![ROWSPOS]; - le := Length(l); - i := 1; - while i <= le and IsZero(l[i]) do i := i + 1; od; - return i; - end ); - -InstallMethod( PositionNonZero, "for a plist matrix, and a position", - [ IsPlistMatrixRep, IsInt ], - function( m, p ) - local i,l,le; - l := m![ROWSPOS]; - le := Length(l); - i := p+1; - while i <= le and IsZero(l[i]) do i := i + 1; od; - if i > le+1 then - return le+1; - else - return i; - fi; - end ); - -InstallMethod( PositionLastNonZero, "for a plist matrix", - [ IsPlistMatrixRep ], - function( m ) - local i,l,le; - l := m![ROWSPOS]; - le := Length(l); - i := le; - while i >= 1 and IsZero(l[i]) do i := i - 1; od; - return i; - end ); - -InstallMethod( PositionLastNonZero, "for a plist matrix, and a position", - [ IsPlistMatrixRep, IsInt ], - function( m, p ) - local i,l,le; - l := m![ROWSPOS]; - le := Length(l); - i := p-1; - while i >= 1 and IsZero(l[i]) do i := i - 1; od; - if i < 0 then - return 0; - else - return i; - fi; - end ); - -InstallMethod( Position, "for a plist matrix, and a plist vector", - [ IsPlistMatrixRep, IsPlistVectorRep ], - function( m, v ) - local i,l,le; - l := m![ROWSPOS]; - le := Length(l); - i := 1; - while i <= le and l[i] <> v do i := i + 1; od; - if i > le then - return fail; - else - return i; - fi; - end ); - -InstallMethod( Position, "for a plist matrix, and a plist vector", - [ IsPlistMatrixRep, IsPlistVectorRep, IsInt ], - function( m, v, p ) - local i,l,le; - l := m![ROWSPOS]; - le := Length(l); - i := p+1; - while i <= le and l[i] <> v do i := i + 1; od; - if i > le then - return fail; - else - return i; - fi; - end ); - -InstallMethod( PositionSortedOp, "for a plist matrix, and a plist vector", - [ IsPlistMatrixRep, IsPlistVectorRep ], - function( m, v ) - return POSITION_SORTED_LIST(m![ROWSPOS],v); - end ); - -InstallMethod( PositionSortedOp, "for a plist matrix, and a plist vector", - [ IsPlistMatrixRep, IsPlistVectorRep, IsFunction ], - function( m, v, f ) - return POSITION_SORTED_LIST_COMP(m![ROWSPOS],v); - end ); - InstallMethod( MutableCopyMat, "for a plist matrix", [ IsPlistMatrixRep ], function( m ) @@ -1149,25 +799,6 @@ InstallOtherMethod( CopySubMatrix, od; end ); -InstallMethod( CopySubMatrix, - "for a plist matrix and a checking plist matrix and four lists", - [ IsPlistMatrixRep, - IsPlistMatrixRep and IsCheckingMatrix and IsMutable, - IsList, IsList, IsList, IsList ], - function( m, n, srcrows, dstrows, srccols, dstcols ) - local i; - if not ForAll(dstcols,x->x <= n![RLPOS]) then - ErrorNoReturn("CopySubMatrix: Destination column positions out of range"); - fi; - if not(ForAll(dstcols,x->x <= Length(n![ROWSPOS])+1)) then - ErrorNoReturn("CopySubMatrix: Destination row positions out of range"); - fi; - for i in [1..Length(srcrows)] do - n![ROWSPOS][dstrows[i]]![ELSPOS]{dstcols} := - m![ROWSPOS][srcrows[i]]![ELSPOS]{srccols}; - od; - end ); - InstallMethod( MatElm, "for a plist matrix and two positions", [ IsPlistMatrixRep, IsPosInt, IsPosInt ], function( m, row, col ) @@ -1180,21 +811,20 @@ InstallMethod( SetMatElm, "for a plist matrix, two positions, and an object", m![ROWSPOS][row]![ELSPOS][col] := ob; end ); -InstallMethod( SetMatElm, - "for a checking plist matrix, two positions, and an object", - [ IsPlistMatrixRep and IsCheckingMatrix and IsMutable, - IsPosInt, IsPosInt, IsObject ], +InstallMethod( \[\], "for a plist matrix and two positions", + [ IsPlistMatrixRep, IsPosInt, IsPosInt ], + function( m, row, col ) + return m![ROWSPOS][row]![ELSPOS][col]; + end ); + +InstallMethod( \[\]\:\=, "for a plist matrix, two positions, and an object", + [ IsPlistMatrixRep and IsMutable, IsPosInt, IsPosInt, IsObject ], function( m, row, col, ob ) - if row > Length( m![ROWSPOS] ) or col > m![RLPOS] then - ErrorNoReturn("SetMatElm: Row or column number out of bounds"); - fi; - if not ob in m![BDPOS] then - ErrorNoReturn("SetMatElm: Value to be assigned is not in base domain"); - fi; m![ROWSPOS][row]![ELSPOS][col] := ob; end ); + ############################################################################ # Printing and viewing methods: ############################################################################ @@ -1202,7 +832,6 @@ InstallMethod( SetMatElm, InstallMethod( ViewObj, "for a plist matrix", [ IsPlistMatrixRep ], function( m ) Print("<"); - if IsCheckingMatrix(m) then Print("checking "); fi; if not IsMutable(m) then Print("immutable "); fi; Print(Length(m![ROWSPOS]),"x",m![RLPOS],"-matrix over ",m![BDPOS],">"); end ); @@ -1210,9 +839,6 @@ InstallMethod( ViewObj, "for a plist matrix", [ IsPlistMatrixRep ], InstallMethod( PrintObj, "for a plist matrix", [ IsPlistMatrixRep ], function( m ) Print("NewMatrix(IsPlistMatrixRep"); - if IsCheckingMatrix(m) then - Print(" and IsCheckingMatrix"); - fi; if IsFinite(m![BDPOS]) and IsField(m![BDPOS]) then Print(",GF(",Size(m![BDPOS]),"),"); else @@ -1225,7 +851,6 @@ InstallMethod( Display, "for a plist matrix", [ IsPlistMatrixRep ], function( m ) local i; Print("<"); - if IsCheckingMatrix(m) then Print("checking "); fi; if not IsMutable(m) then Print("immutable "); fi; Print(Length(m![ROWSPOS]),"x",m![RLPOS],"-matrix over ",m![BDPOS],":\n"); for i in [1..Length(m![ROWSPOS])] do @@ -1243,9 +868,6 @@ InstallMethod( String, "for plist matrix", [ IsPlistMatrixRep ], function( m ) local st; st := "NewMatrix(IsPlistMatrixRep"; - if IsCheckingMatrix(m) then - Append(st," and IsCheckingMatrix"); - fi; Add(st,','); if IsFinite(m![BDPOS]) and IsField(m![BDPOS]) then Append(st,"GF("); @@ -1279,27 +901,6 @@ InstallMethod( \+, "for two plist matrices", SUM_LIST_LIST_DEFAULT(a![ROWSPOS],b![ROWSPOS])]); end ); -InstallMethod( \+, "for two checking plist matrices", - [ IsPlistMatrixRep and IsCheckingMatrix, - IsPlistMatrixRep and IsCheckingMatrix ], - function( a, b ) - local ty; - if Length(a![ROWSPOS]) <> Length(b![ROWSPOS]) or - a![RLPOS] <> b![RLPOS] then - ErrorNoReturn("\\+: Dimensions do not fit together"); - fi; - if not IsIdenticalObj(a![BDPOS],b![BDPOS]) then - ErrorNoReturn("\\+: BaseDomains are not the same"); - fi; - if not IsMutable(a) and IsMutable(b) then - ty := TypeObj(b); - else - ty := TypeObj(a); - fi; - return Objectify(ty,[a![BDPOS],a![EMPOS],a![RLPOS], - SUM_LIST_LIST_DEFAULT(a![ROWSPOS],b![ROWSPOS])]); - end ); - InstallMethod( \-, "for two plist matrices", [ IsPlistMatrixRep, IsPlistMatrixRep ], function( a, b ) @@ -1313,27 +914,6 @@ InstallMethod( \-, "for two plist matrices", DIFF_LIST_LIST_DEFAULT(a![ROWSPOS],b![ROWSPOS])]); end ); -InstallMethod( \-, "for two checking plist matrices", - [ IsPlistMatrixRep and IsCheckingMatrix, - IsPlistMatrixRep and IsCheckingMatrix ], - function( a, b ) - local ty; - if Length(a![ROWSPOS]) <> Length(b![ROWSPOS]) or - a![RLPOS] <> b![RLPOS] then - ErrorNoReturn("\\-: Dimensions do not fit together"); - fi; - if not IsIdenticalObj(a![BDPOS],b![BDPOS]) then - ErrorNoReturn("\\-: BaseDomains are not the same"); - fi; - if not IsMutable(a) and IsMutable(b) then - ty := TypeObj(b); - else - ty := TypeObj(a); - fi; - return Objectify(ty,[a![BDPOS],a![EMPOS],a![RLPOS], - DIFF_LIST_LIST_DEFAULT(a![ROWSPOS],b![ROWSPOS])]); - end ); - InstallMethod( \*, "for two plist matrices", [ IsPlistMatrixRep, IsPlistMatrixRep ], function( a, b ) @@ -1744,18 +1324,6 @@ InstallMethod( Fold, "for a plist vector, a positive int, and a plist matrix", # fi; # end ); -InstallMethod( ConstructingFilter, "for a checking plist vector", - [ IsPlistVectorRep and IsCheckingVector ], - function( v ) - return IsPlistVectorRep and IsCheckingVector; - end ); - -InstallMethod( ConstructingFilter, "for a checking plist matrix", - [ IsPlistMatrixRep and IsCheckingMatrix ], - function( m ) - return IsPlistMatrixRep and IsCheckingMatrix; - end ); - InstallMethod( ConstructingFilter, "for a plist vector", [ IsPlistVectorRep ], function( v ) @@ -1774,12 +1342,6 @@ InstallMethod( ChangedBaseDomain, "for a plist vector, and a domain", return NewVector( IsPlistVectorRep, r, v![ELSPOS] ); end ); -InstallMethod( ChangedBaseDomain, "for a checking plist vector, and a domain", - [ IsPlistVectorRep and IsCheckingVector, IsRing ], - function( v, r ) - return NewVector(IsPlistVectorRep and IsCheckingVector, r, v![ELSPOS]); - end ); - InstallMethod( ChangedBaseDomain, "for a plist matrix, and a domain", [ IsPlistMatrixRep, IsRing ], function( m, r ) @@ -1787,13 +1349,6 @@ InstallMethod( ChangedBaseDomain, "for a plist matrix, and a domain", List(m![ROWSPOS], x-> x![ELSPOS])); end ); -InstallMethod( ChangedBaseDomain, "for a checking plist matrix, and a domain", - [ IsPlistMatrixRep and IsCheckingMatrix, IsRing ], - function( m, r ) - return NewMatrix(IsPlistMatrixRep and IsCheckingMatrix, r, NumberColumns(m), - List(m![ROWSPOS], x-> x![ELSPOS])); - end ); - InstallMethod( CompatibleVector, "for a plist matrix", [ IsPlistMatrixRep ], function( v ) @@ -1803,7 +1358,7 @@ InstallMethod( CompatibleVector, "for a plist matrix", InstallMethod( NewCompanionMatrix, "for IsPlistMatrixRep, a polynomial and a ring", [ IsPlistMatrixRep, IsUnivariatePolynomial, IsRing ], - function( ty, po, bd ) + function( filter, po, bd ) local i,l,ll,n,one; one := One(bd); l := CoefficientsOfUnivariatePolynomial(po); @@ -1812,7 +1367,7 @@ InstallMethod( NewCompanionMatrix, Error("CompanionMatrix: polynomial is not monic"); return fail; fi; - ll := NewMatrix(ty,bd,n,[]); + ll := NewMatrix(IsPlistMatrixRep,bd,n,[]); l := Vector(-l{[1..n]},CompatibleVector(ll)); for i in [1..n-1] do Add(ll,ZeroMutable(l)); @@ -1822,24 +1377,3 @@ InstallMethod( NewCompanionMatrix, return ll; end ); -InstallMethod( NewCompanionMatrix, - "for IsPlistMatrixRep and IsCheckingMatrix, a polynomial and a ring", - [ IsPlistMatrixRep and IsCheckingMatrix, IsUnivariatePolynomial, IsRing ], - function( ty, po, bd ) - local i,l,ll,n,one; - one := One(bd); - l := CoefficientsOfUnivariatePolynomial(po); - n := Length(l)-1; - if not IsOne(l[n+1]) then - Error("CompanionMatrix: polynomial is not monic"); - return fail; - fi; - ll := NewMatrix(ty,bd,n,[]); - l := Vector(-l{[1..n]},CompatibleVector(ll)); - for i in [1..n-1] do - Add(ll,ZeroMutable(l)); - ll[i][i+1] := one; - od; - Add(ll,l); - return ll; - end ); diff --git a/lib/matrix.gi b/lib/matrix.gi index 7e3c7e7523..336d33051c 100644 --- a/lib/matrix.gi +++ b/lib/matrix.gi @@ -23,15 +23,19 @@ InstallMethod(Zero, [IsRectangularTable and IsAdditiveElementWithZeroCollColl and IsInternalRep], ZERO_ATTR_MAT); -# Fallback methods for matrix entry access. We lower the rank -# so that these are only used as a last resort. +# Fallback methods for matrix entry access. InstallOtherMethod( \[\], [ IsMatrix, IsPosInt, IsPosInt ], - -RankFilter(IsMatrix), + 1-RankFilter(IsMatrix), function (m,i,j) return m[i][j]; end ); InstallOtherMethod( \[\]\:\=, [ IsMatrix and IsMutable, IsPosInt, IsPosInt, IsObject ], - -RankFilter(IsMatrix), + 1-RankFilter(IsMatrix), function (m,i,j,o) m[i][j]:=o; end ); +InstallOtherMethod( Unpack, + [ IsMatrix ], + ShallowCopy ); + + ############################################################################# ## #F PrintArray( ) . . . . . . . . . . . . . . . . pretty print matrix @@ -159,7 +163,7 @@ InstallMethod( IsGeneralizedCartanMatrix, ## #M IsDiagonalMat() ## -InstallMethod( IsDiagonalMat, +InstallOtherMethod( IsDiagonalMat, "for a matrix", [ IsMatrix ], function( mat ) @@ -183,7 +187,7 @@ InstallOtherMethod( IsDiagonalMat, [ IsEmpty ], ReturnTrue ); ## #M IsUpperTriangularMat() ## -InstallMethod( IsUpperTriangularMat, +InstallOtherMethod( IsUpperTriangularMat, "for a matrix", [ IsMatrix ], function( mat ) @@ -204,7 +208,7 @@ InstallMethod( IsUpperTriangularMat, ## #M IsLowerTriangularMat() ## -InstallMethod( IsLowerTriangularMat, +InstallOtherMethod( IsLowerTriangularMat, "for a matrix", [ IsMatrix ], function( mat ) @@ -251,7 +255,7 @@ DeclareRepresentation( "IsNullMapMatrix", IsMatrix, [ ] ); BindGlobal( "NullMapMatrix", Objectify( NewType( ListsFamily, IsNullMapMatrix ), MakeImmutable([ ]) ) ); -InstallMethod( Length, +InstallOtherMethod( Length, "for null map matrix", [ IsNullMapMatrix ], null -> 0 ); @@ -534,7 +538,7 @@ BindGlobal( "Matrix_MinimalPolynomialSameField", function( fld, mat, ind ) piv := j; repeat AddRowVector(w,base[piv],-w[piv], piv, n); - piv := PositionNot(w, zero, piv); + piv := PositionNonZero(w, piv); until piv > n or not IsBound(base[piv]); if piv <= n then MultRowVector(w,Inverse(w[piv])); @@ -1127,13 +1131,11 @@ InstallMethod( IsZero, [ IsMatrix ], function( mat ) local ncols, # number of columns - zero, # zero coefficient row; # loop over rows in 'obj' ncols:= DimensionsMat( mat )[2]; - zero:= Zero( mat[1,1] ); for row in mat do - if PositionNot( row, zero ) <= ncols then + if PositionNonZero( row ) <= ncols then return false; fi; od; @@ -1203,6 +1205,51 @@ function( mat ) end ); +############################################################################# +## +#M BaseDomain( ) +#M OneOfBaseDomain( ) +#M ZeroOfBaseDomain( ) +#M BaseDomain( ) +#M OneOfBaseDomain( ) +#M ZeroOfBaseDomain( ) +## +## Note that a matrix that is a plain list is necessarily nonempty +## and dense, hence it is safe to access the first row. +## Since the rows are homogeneous and nonempty, +## one can also access the first entry in the first row. +## +InstallOtherMethod( BaseDomain, + "generic method for a row vector", + [ IsRowVector ], + DefaultRing ); + +InstallOtherMethod( OneOfBaseDomain, + "generic method for a row vector", + [ IsRowVector ], + v -> One( v[1] ) ); + +InstallOtherMethod( ZeroOfBaseDomain, + "generic method for a row vector", + [ IsRowVector ], + v -> Zero( v[1] ) ); + +InstallOtherMethod( BaseDomain, + "generic method for a matrix that is a plain list", + [ IsMatrix and IsPlistRep ], + mat -> BaseDomain( mat[1] ) ); + +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] ) ); + + ############################################################################# ## #M DepthOfUpperTriangularMatrix( ) @@ -1470,18 +1517,19 @@ InstallMethod( DeterminantMatDivFree, V; ## graph # check that the argument is a square matrix and set the size - n := Length(M); - if not n = Length(M[1]) or not IsRectangularTable(M) then + n := NumberRows( M ); + if not n = NumberColumns( M ) or not IsRectangularTable(M) then Error("DeterminantMatDivFree: must be a square matrix"); fi; ## initialze the final sum, the vertex set, initial parity ## and level indexes ## - zero := Zero(M[1,1]); + zero := ZeroOfBaseDomain( M ); Vs := zero; V := []; - pmone := (-One(M[1,1]))^((n mod 2)+1); + 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 @@ -1517,14 +1565,17 @@ InstallMethod( DeterminantMatDivFree, 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[a][c][clevel] * M[ v, w ]; V[b][b][nlevel]:= V[b][b][nlevel]- - V[a][c][clevel]*M[v][u]; + V[a][c][clevel] * M[ v, u ]; od; od; od; @@ -1546,7 +1597,7 @@ InstallMethod( DeterminantMatDivFree, ## 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]; + Vs := Vs + V[n-u+1][n-v+1][clevel] * M[ v, u ]; od; od; @@ -1560,7 +1611,7 @@ InstallMethod( DeterminantMatDivFree, ## #M DimensionsMat( ) ## -InstallMethod( DimensionsMat, +InstallOtherMethod( DimensionsMat, [ IsMatrix ], function( A ) if IsRectangularTable(A) then @@ -1882,7 +1933,7 @@ InstallMethod( MutableCopyMat, "generic method", [IsList], ## #M MutableTransposedMat( ) . . . . . . . . . . transposed of a matrix ## -InstallMethod( MutableTransposedMat, +InstallOtherMethod( MutableTransposedMat, "generic method", [ IsRectangularTable and IsMatrix ], function( mat ) @@ -2192,7 +2243,7 @@ end ); ## #M RankMat( ) . . . . . . . . . . . . . . . . . . . rank of a matrix ## -InstallMethod( RankMatDestructive, +InstallOtherMethod( RankMatDestructive, "generic method for mutable matrices", [ IsMatrix and IsMutable ], function( mat ) @@ -2203,7 +2254,7 @@ InstallMethod( RankMatDestructive, return mat; end ); -InstallMethod( RankMat, +InstallOtherMethod( RankMat, "generic method for matrices", [ IsMatrix ], mat -> RankMatDestructive( MutableCopyMat( mat ) ) ); @@ -2249,7 +2300,7 @@ InstallMethod( SemiEchelonMatDestructive, fi; od; - j := PositionNot( row, zero ); + j := PositionNonZero( row ); if j <= ncols then # We found a new basis vector. @@ -2354,7 +2405,7 @@ InstallMethod( SemiEchelonMatTransformationDestructive, fi; od; - j:= PositionNot( row, zero ); + j:= PositionNonZero( row ); if j <= ncols then # We found a new basis vector. @@ -2428,10 +2479,10 @@ InstallGlobalFunction( SemiEchelonMatsNoCo, function( mats ) # Get the first nonzero column. i:= 1; - j:= PositionNot( mat[1], zero ); + j:= PositionNonZero( mat[1] ); while n < j and i < m do i:= i + 1; - j:= PositionNot( mat[i], zero ); + j:= PositionNonZero( mat[i] ); od; if j <= n then @@ -2515,24 +2566,22 @@ InstallMethod( IsMonomialMatrix, "generic method for matrices", [ IsMatrix ], function( M ) - local zero, # zero of the base ring - len, # length of rows + local len, # length of rows found, # store positions of nonzero elements row, # loop over rows j; # position of first non-zero element - zero:= Zero(M[1,1]); len:= Length( M[1] ); if Length( M ) <> len then return false; fi; found:= BlistList( M, [] ); for row in M do - j := PositionNot( row, zero ); + j := PositionNonZero( row ); if len < j or found[j] then return false; fi; - if PositionNot( row, zero, j ) <= len then + if PositionNonZero( row, j ) <= len then return false; fi; found[j] := true; @@ -2597,7 +2646,7 @@ end ); ## #M KroneckerProduct( , ) ## -InstallMethod( KroneckerProduct, +InstallOtherMethod( KroneckerProduct, "generic method for matrices", IsIdenticalObj, [ IsMatrix, IsMatrix ], diff --git a/lib/memory.gi b/lib/memory.gi index fab85b36fb..cee094cd65 100644 --- a/lib/memory.gi +++ b/lib/memory.gi @@ -364,17 +364,17 @@ InstallOtherMethod(CycleStructurePerm, # Matrix methods: -InstallMethod( DimensionsMat, "for a matrix with memory",true, - [ IsMatrix and IsObjWithMemory ], 0, +InstallOtherMethod( DimensionsMat, "for a matrix with memory", + [ IsMatrix and IsObjWithMemory ], function(M) return DimensionsMat(M!.el); end); -InstallMethod( Length, "for a matrix with memory", +InstallOtherMethod( Length, "for a matrix with memory", [ IsMatrix and IsObjWithMemory ], M -> Length(M!.el) ) ; -InstallMethod( ELM_LIST, "for a matrix with memory",true, - [ IsMatrix and IsObjWithMemory, IsPosInt ], 0, +InstallOtherMethod( ELM_LIST, "for a matrix with memory", + [ IsMatrix and IsObjWithMemory, IsPosInt ], function(M,i) return M!.el[i]; end); diff --git a/lib/padics.gi b/lib/padics.gi index 6ac4757eed..f2e20008eb 100644 --- a/lib/padics.gi +++ b/lib/padics.gi @@ -950,7 +950,7 @@ InstallMethod( IsZero, 0, function ( x ) - if PositionNot( x![2], 0 ) > Length( x![2] ) then + if PositionNonZero( x![2] ) > Length( x![2] ) then return true; fi; return false; diff --git a/lib/pcgs.gi b/lib/pcgs.gi index 6aeba23bf7..83b91e189e 100644 --- a/lib/pcgs.gi +++ b/lib/pcgs.gi @@ -354,7 +354,7 @@ InstallMethod( DepthOfPcElement, 0, function( pcgs, elm ) - return PositionNot( ExponentsOfPcElement( pcgs, elm ), 0 ); + return PositionNonZero( ExponentsOfPcElement( pcgs, elm ) ); end ); ############################################################################# @@ -367,7 +367,7 @@ InstallMethod( DepthAndLeadingExponentOfPcElement, function( pcgs, elm ) local e,p; e:=ExponentsOfPcElement( pcgs, elm ); - p:=PositionNot( e, 0 ); + p:=PositionNonZero( e ); if p>Length(e) then return [p,0]; else @@ -528,7 +528,7 @@ function( pcgs, elm ) local exp, dep; exp := ExponentsOfPcElement( pcgs, elm ); - dep := PositionNot( exp, 0 ); + dep := PositionNonZero( exp ); if Length(exp) < dep then return fail; else diff --git a/lib/pcgsmodu.gi b/lib/pcgsmodu.gi index 10f5e5f553..3c9cc3bcef 100644 --- a/lib/pcgsmodu.gi +++ b/lib/pcgsmodu.gi @@ -586,7 +586,7 @@ function( pcgs, elm ) local exp, dep; exp := ExponentsOfPcElement( pcgs, elm ); - dep := PositionNot( exp, 0 ); + dep := PositionNonZero( exp ); if Length(exp) < dep then return fail; else @@ -730,7 +730,7 @@ function( pcgs, elm ) if d > Length(num) then return Length(pcgs)+1; elif d in pcgs!.moduloDepths then - return PositionNot( ExponentsOfPcElement( pcgs, elm ), 0 ); + return PositionNonZero( ExponentsOfPcElement( pcgs, elm ) ); else return pcgs!.depthMap[d]; fi; diff --git a/lib/ringsc.gi b/lib/ringsc.gi index 85d7445c8c..639783437b 100644 --- a/lib/ringsc.gi +++ b/lib/ringsc.gi @@ -80,7 +80,7 @@ InstallMethod( PrintObj, return; fi; - depth := PositionNot( elm, 0 ); + depth := PositionNonZero( elm ); if len < depth then @@ -137,7 +137,7 @@ function( elm ) return ""; fi; - depth := PositionNot( elm, 0 ); + depth := PositionNonZero( elm ); if len < depth then @@ -556,7 +556,7 @@ local p, j, f, fj, g, q, gj, m, k, i; e:=e[1]; fi; - #p:=PositionNot(e,0); + #p:=PositionNonZero(e); # find the position of largest order p:=-1; f:=1; @@ -1132,7 +1132,7 @@ local ones,s, moduli, orders, offsets, o, p, newmod, t, nams, gens, e, f, D, i, for i in [1..Length(list)] do o:=[]; for j in [1..Length(s[i][1])] do - p:=PositionNot(s[i][1][j],0); + p:=PositionNonZero(s[i][1][j]); if moduli[i][p]=0 then Add(o,0); else diff --git a/lib/vecmat.gi b/lib/vecmat.gi index 8ee0d41408..68a685744d 100644 --- a/lib/vecmat.gi +++ b/lib/vecmat.gi @@ -1655,18 +1655,6 @@ InstallMethod( Append, "for GF2 vectors", IsGF2VectorRep and IsList], 0, APPEND_GF2VEC); -############################################################################# -## -#M PositionCanonical( , ) . . for GF2 matrices -## -InstallMethod( PositionCanonical, - "for internally represented lists, fall back on `Position'", - true, # the list may be non-homogeneous. - [ IsList and IsGF2MatrixRep, IsObject ], 0, - function( list, obj ) - return Position( list, obj, 0 ); -end ); - ############################################################################# ## #M ShallowCopy( ) . . . for GF2 vectors @@ -2122,9 +2110,8 @@ InstallMethod(DeterminantMatDestructive, ## -InstallMethod(RankMatDestructive, +InstallOtherMethod(RankMatDestructive, "kernel method for plain list of GF2 vectors", - true, [IsMatrix and IsPlistRep and IsFFECollColl and IsMutable], GF2_AHEAD_OF_8BIT_RANK, RANK_LIST_GF2VECS); @@ -2233,27 +2220,15 @@ InstallMethod( Matrix, "for a list of vecs, an integer, and a gf2 mat", ConvertToMatrixRep(li,2); return li; end ); -BindGlobal( "PositionLastNonZeroFunc", + +InstallMethod( PositionLastNonZero, "for a row vector obj", + [IsVectorObj], function(l) local i; i := Length(l); while i >= 1 and IsZero(l[i]) do i := i - 1; od; return i; end ); -BindGlobal( "PositionLastNonZeroFunc2", - function(l,pos) - local i; - i := pos-1; - while i >= 1 and IsZero(l[i]) do i := i - 1; od; - return i; - end ); - -InstallMethod( PositionLastNonZero, "for a row vector obj", - [IsVectorObj], PositionLastNonZeroFunc ); -InstallMethod( PositionLastNonZero, "for a matrix obj", - [IsMatrixObj], PositionLastNonZeroFunc ); -InstallMethod( PositionLastNonZero, "for a matrix obj, and an index", - [IsMatrixObj, IsPosInt], PositionLastNonZeroFunc2 ); InstallMethod( ExtractSubMatrix, "for a gf2 matrix, and two lists", [IsGF2MatrixRep, IsList, IsList], diff --git a/lib/vspcmat.gi b/lib/vspcmat.gi index c1043e5882..f4a5f99021 100644 --- a/lib/vspcmat.gi +++ b/lib/vspcmat.gi @@ -338,10 +338,10 @@ BindGlobal( "HeadsInfoOfSemiEchelonizedMats", function( mats, dims ) # Get the pivot. row:= 1; - j:= PositionNot( mats[i][row], zero ); + j:= PositionNonZero( mats[i][row] ); while dimcol < j and row < dimrow do row:= row + 1; - j:= PositionNot( mats[i][row], zero ); + j:= PositionNonZero( mats[i][row] ); od; if dimrow < row or mats[i][ row ][j] <> one then @@ -1068,7 +1068,7 @@ InstallMethod( CloseMutableBasis, # If necessary add the sifted vector, and update the basis info. for i in [ 1 .. m ] do - j := PositionNot( v[i], zero ); + j := PositionNonZero( v[i] ); if j <= n then scalar:= Inverse( v[i][j] ); for k in [ 1 .. m ] do diff --git a/lib/vspcrow.gi b/lib/vspcrow.gi index bd9a9fa91a..077bca45ff 100644 --- a/lib/vspcrow.gi +++ b/lib/vspcrow.gi @@ -318,7 +318,7 @@ InstallMethod( Coefficients, # Compute the coefficients of the base vectors. v:= ShallowCopy( v ); - i:= PositionNot( v, zero ); + i:= PositionNonZero( v ); while i <= len do pos:= heads[i]; if pos <> 0 then @@ -327,7 +327,7 @@ InstallMethod( Coefficients, else return fail; fi; - i:= PositionNot( v, zero ); + i:= PositionNonZero( v ); od; # Return the coefficients. @@ -414,7 +414,7 @@ HeadsInfoOfSemiEchelonizedMat := function( mat, dim ) # Loop over the columns. for i in [ 1 .. nrows ] do - j:= PositionNot( mat[i], zero ); + j:= PositionNonZero( mat[i] ); if dim < j or mat[i][j] <> one then return fail; fi; @@ -776,7 +776,7 @@ end); InstallMethod( IsZero, "for a row vector", [ IsRowVector ], - v -> IsEmpty( v ) or Length( v ) < PositionNot( v, Zero( v[1] ) ) ); + v -> IsEmpty( v ) or Length( v ) < PositionNonZero( v ) ); ############################################################################# @@ -1575,13 +1575,13 @@ InstallMethod( CloseMutableBasis, for j in [ 1 .. ncols ] do if zero <> v[j] and heads[j] <> 0 then -#T better loop with `PositionNot'? +#T better loop with `PositionNonZero'? AddRowVector( v, basisvectors[ heads[j] ], - v[j] ); fi; od; # If necessary add the sifted vector, and update the basis info. - j := PositionNot( v, zero ); + j := PositionNonZero( v ); if j <= ncols then MultRowVector( v, Inverse( v[j] ) ); Add( basisvectors, v ); @@ -1683,8 +1683,7 @@ InstallOtherMethod( SiftedVector, InstallGlobalFunction( OnLines, function( vec, g ) local c; vec:= OnPoints( vec, g ); - c:= PositionNonZero( vec ); # better than PositionNot, as no `Zero' - # element needs to be created + c:= PositionNonZero( vec ); if c <= Length( vec ) then # Normalize from the *left* if the matrices act from the right! @@ -1706,8 +1705,7 @@ function( v ) local depth; if 0 < Length(v) then - depth:=PositionNonZero( v ); # better than PositionNot, as no `Zero' - # element needs to be created + depth:=PositionNonZero( v ); if depth <= Length(v) then return Inverse(v[depth]) * v; else diff --git a/lib/zlattice.gi b/lib/zlattice.gi index a7b5c4c63d..2cc3ccd2bb 100644 --- a/lib/zlattice.gi +++ b/lib/zlattice.gi @@ -132,7 +132,7 @@ InstallGlobalFunction( PadicCoefficients, coeff:= []; step:= 0; p2:= ( prime - 1 ) / 2; - while PositionNot( b, 0 ) <= n and step < depth do + while PositionNonZero( b ) <= n and step < depth do step:= step + 1; coeff[ step ]:= ShallowCopy( b * Amodpinv ); for i in [ 1 .. n ] do diff --git a/tst/testinstall/MatrixObj/ConcatenationOfVectors.tst b/tst/testinstall/MatrixObj/ConcatenationOfVectors.tst new file mode 100644 index 0000000000..e3d42c6561 --- /dev/null +++ b/tst/testinstall/MatrixObj/ConcatenationOfVectors.tst @@ -0,0 +1,23 @@ +gap> START_TEST("ConcatenationOfVectors.tst"); +gap> l1 := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> l2 := [6,2,7,4,5,6]; +[ 6, 2, 7, 4, 5, 6 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, l1); + +gap> v2 := Vector(IsPlistVectorRep, Rationals, l2); + +gap> vv := ConcatenationOfVectors( v1, v2, v2 ); + +gap> Unpack( vv ); +[ 1, 2, 3, 4, 5, 6, 6, 2, 7, 4, 5, 6, 6, 2, 7, 4, 5, 6 ] +gap> v3 := Vector(GF(5), l1*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> v4 := Vector(GF(5), l2*One(GF(5))); +[ Z(5)^0, Z(5), Z(5), Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> vv := ConcatenationOfVectors( [ v3, v4, v3 ] ); +< mutable compressed vector length 18 over GF(5) > +gap> Unpack( vv ); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0, Z(5)^0, Z(5), Z(5), Z(5)^2, + 0*Z(5), Z(5)^0, Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> STOP_TEST("ConcatenationOfVectors.tst",1); diff --git a/tst/testinstall/MatrixObj/CopySubVector.tst b/tst/testinstall/MatrixObj/CopySubVector.tst new file mode 100644 index 0000000000..29b6c08ab8 --- /dev/null +++ b/tst/testinstall/MatrixObj/CopySubVector.tst @@ -0,0 +1,20 @@ +gap> START_TEST("CopySubVector.tst"); +gap> l1 := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, l1); + +gap> l2 := [1,1,1,2,2,2,3,3,3]; +[ 1, 1, 1, 2, 2, 2, 3, 3, 3 ] +gap> v2 := Vector(IsPlistVectorRep, Rationals, l2); + +gap> CopySubVector( v1, [2,4,6], v2, [1,2,4] ); +gap> Unpack(v1); +[ 1, 1, 3, 1, 5, 2 ] +gap> v3 := Vector(GF(5), l1*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> v4 := Vector(GF(5), l2*One(GF(5))); +[ Z(5)^0, Z(5)^0, Z(5)^0, Z(5), Z(5), Z(5), Z(5)^3, Z(5)^3, Z(5)^3 ] +gap> CopySubVector( v4, [2,4,6], v3, [1,2,3] ); +gap> v4; +[ Z(5)^0, Z(5)^0, Z(5)^0, Z(5), Z(5), Z(5)^3, Z(5)^3, Z(5)^3, Z(5)^3 ] +gap> STOP_TEST("CopySubVector.tst",1); diff --git a/tst/testinstall/MatrixObj/DistanceOfVectors.tst b/tst/testinstall/MatrixObj/DistanceOfVectors.tst new file mode 100644 index 0000000000..0e746944da --- /dev/null +++ b/tst/testinstall/MatrixObj/DistanceOfVectors.tst @@ -0,0 +1,32 @@ +gap> START_TEST( "DistanceOfVectors.tst" ); +gap> l1 := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, l1); + +gap> l2 := [0,2,3,0,5,6]; +[ 0, 2, 3, 0, 5, 6 ] +gap> v2 := Vector(IsPlistVectorRep, Rationals, l2); + +gap> l3 := [0,0,0,0,0,0]; +[ 0, 0, 0, 0, 0, 0 ] +gap> v3 := Vector(IsPlistVectorRep, Rationals, l3); + +gap> DistanceOfVectors( v1, v2 ); +2 +gap> DistanceOfVectors( v2, v3 ); +4 +gap> DistanceOfVectors( v1, v3 ); +6 +gap> v4 := Vector(GF(5), l1*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> v5 := Vector(GF(5), l2*One(GF(5))); +[ 0*Z(5), Z(5), Z(5)^3, 0*Z(5), 0*Z(5), Z(5)^0 ] +gap> v6 := Vector(GF(5), l3*One(GF(5))); +[ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ] +gap> DistanceOfVectors( v4, v5 ); +2 +gap> DistanceOfVectors( v5, v6 ); +3 +gap> DistanceOfVectors( v4, v6 ); +5 +gap> STOP_TEST( "DistanceOfVectors.tst", 1); diff --git a/tst/testinstall/MatrixObj/ExtractSubvector.tst b/tst/testinstall/MatrixObj/ExtractSubvector.tst new file mode 100644 index 0000000000..8b3edff6b0 --- /dev/null +++ b/tst/testinstall/MatrixObj/ExtractSubvector.tst @@ -0,0 +1,14 @@ +gap> START_TEST("ExtractSubVector.tst"); +gap> l1 := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> v3 := Vector(GF(5), l1*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> ExtractSubVector( v3, [1,2,4] ); +[ Z(5)^0, Z(5), Z(5)^2 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, l1); + +gap> ExtractSubVector( v1, [1,2,4] ); + +gap> Unpack( last ); +[ 1, 2, 4 ] +gap> STOP_TEST("ExtractSubVector.tst",1); diff --git a/tst/testinstall/MatrixObj/IdentityMatrix.tst b/tst/testinstall/MatrixObj/IdentityMatrix.tst new file mode 100644 index 0000000000..99666855da --- /dev/null +++ b/tst/testinstall/MatrixObj/IdentityMatrix.tst @@ -0,0 +1,74 @@ +gap> START_TEST("IdentityMatrix.tst"); + +# +gap> m:=IdentityMatrix( GF(2), 5 ); Display(m); + + 1 . . . . + . 1 . . . + . . 1 . . + . . . 1 . + . . . . 1 +gap> m:=IdentityMatrix( GF(3), 5 ); Display(m); +[ [ Z(3)^0, 0*Z(3), 0*Z(3), 0*Z(3), 0*Z(3) ], + [ 0*Z(3), Z(3)^0, 0*Z(3), 0*Z(3), 0*Z(3) ], + [ 0*Z(3), 0*Z(3), Z(3)^0, 0*Z(3), 0*Z(3) ], + [ 0*Z(3), 0*Z(3), 0*Z(3), Z(3)^0, 0*Z(3) ], + [ 0*Z(3), 0*Z(3), 0*Z(3), 0*Z(3), Z(3)^0 ] ] + 1 . . . . + . 1 . . . + . . 1 . . + . . . 1 . + . . . . 1 +gap> m:=IdentityMatrix( GF(4), 5 ); Display(m); +[ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], + [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ], + [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ], + [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ], + [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ] + 1 . . . . + . 1 . . . + . . 1 . . + . . . 1 . + . . . . 1 +gap> m:=IdentityMatrix( Integers, 5 ); Display(m); +<5x5-matrix over Integers> +<5x5-matrix over Integers: +[[ 1, 0, 0, 0, 0 ] + [ 0, 1, 0, 0, 0 ] + [ 0, 0, 1, 0, 0 ] + [ 0, 0, 0, 1, 0 ] + [ 0, 0, 0, 0, 1 ] +]> +gap> m:=IdentityMatrix( Integers mod 4, 5 ); Display(m); +<5x5-matrix over (Integers mod 4)> +<5x5-matrix over (Integers mod 4): +[[ ZmodnZObj( 1, 4 ), ZmodnZObj( 0, 4 ), ZmodnZObj( 0, 4 ), ZmodnZObj( 0, 4 ), + ZmodnZObj( 0, 4 ) ] + [ ZmodnZObj( 0, 4 ), ZmodnZObj( 1, 4 ), ZmodnZObj( 0, 4 ), ZmodnZObj( 0, 4 ), + ZmodnZObj( 0, 4 ) ] + [ ZmodnZObj( 0, 4 ), ZmodnZObj( 0, 4 ), ZmodnZObj( 1, 4 ), ZmodnZObj( 0, 4 ), + ZmodnZObj( 0, 4 ) ] + [ ZmodnZObj( 0, 4 ), ZmodnZObj( 0, 4 ), ZmodnZObj( 0, 4 ), ZmodnZObj( 1, 4 ), + ZmodnZObj( 0, 4 ) ] + [ ZmodnZObj( 0, 4 ), ZmodnZObj( 0, 4 ), ZmodnZObj( 0, 4 ), ZmodnZObj( 0, 4 ), + ZmodnZObj( 1, 4 ) ] +]> + +# +gap> m:=IdentityMatrix( Integers, 0 ); Display(m); +<0x0-matrix over Integers> +<0x0-matrix over Integers: +]> + +# some error checking +gap> m:=IdentityMatrix( GF(2), -1 ); +Error, must be a non-negative integer (not a integer) +gap> m:=IdentityMatrix( GF(3), -1 ); +Error, must be a non-negative integer (not a integer) +gap> m:=IdentityMatrix( GF(4), -1 ); +Error, must be a non-negative integer (not a integer) +gap> m:=IdentityMatrix( Integers mod 4, -1 ); +Error, must be a non-negative integer (not a integer) + +# +gap> STOP_TEST("IdentityMatrix.tst"); diff --git a/tst/testinstall/MatrixObj/ListOp.tst b/tst/testinstall/MatrixObj/ListOp.tst new file mode 100644 index 0000000000..f47377f3c7 --- /dev/null +++ b/tst/testinstall/MatrixObj/ListOp.tst @@ -0,0 +1,14 @@ +gap> START_TEST("ListOp.tst"); +gap> ll := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, ll); + +gap> List( v1 ); +[ 1, 2, 3, 4, 5, 6 ] +gap> List( v1, x->x^2 ); +[ 1, 4, 9, 16, 25, 36 ] +gap> v5 := Vector(GF(5), ll*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> List( v5, x->x^2 ); +[ Z(5)^0, Z(5)^2, Z(5)^2, Z(5)^0, 0*Z(5), Z(5)^0 ] +gap> STOP_TEST( "ListOp.tst", 1); diff --git a/tst/testinstall/MatrixObj/Matrix.tst b/tst/testinstall/MatrixObj/Matrix.tst new file mode 100644 index 0000000000..dcf299f0e8 --- /dev/null +++ b/tst/testinstall/MatrixObj/Matrix.tst @@ -0,0 +1,43 @@ +gap> START_TEST("Matrix.tst"); + +# +gap> m := Matrix( [[1,2],[3,4]] ); +<2x2-matrix over Rationals> +gap> Display(m); +<2x2-matrix over Rationals: +[[ 1, 2 ] + [ 3, 4 ] +]> + +# +gap> m := Matrix( [[1,2],[3,4]] * Z(2) ); + +gap> Display(m); + 1 . + 1 . + +# +gap> m := Matrix( IsGF2MatrixRep, GF(2), [[1,2],[3,4]] * Z(2) ); + +gap> Display(m); + 1 . + 1 . + +# +gap> m := Matrix( Is8BitMatrixRep, GF(4), [[1,2],[3,4]] * Z(2) ); +[ [ Z(2)^0, 0*Z(2) ], [ Z(2)^0, 0*Z(2) ] ] +gap> Display(m); + 1 . + 1 . + +# +gap> m := Matrix( IsPlistMatrixRep, GF(2), [[1,2],[3,4]] * Z(2) ); +<2x2-matrix over GF(2)> +gap> Display(m); +<2x2-matrix over GF(2): +[[ Z(2)^0, 0*Z(2) ] + [ Z(2)^0, 0*Z(2) ] +]> + +# +gap> STOP_TEST("Matrix.tst"); diff --git a/tst/testinstall/MatrixObj/PositionNonZero.tst b/tst/testinstall/MatrixObj/PositionNonZero.tst new file mode 100644 index 0000000000..c32ec3ccec --- /dev/null +++ b/tst/testinstall/MatrixObj/PositionNonZero.tst @@ -0,0 +1,32 @@ +gap> START_TEST( "PositionNonZero.tst" ); +gap> l1 := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, l1); + +gap> PositionNonZero( v1 ); +1 +gap> l2 := [0,2,3,0,5,6]; +[ 0, 2, 3, 0, 5, 6 ] +gap> v2 := Vector(IsPlistVectorRep, Rationals, l2); + +gap> PositionNonZero( v2 ); +2 +gap> l3 := [0,0,0,0,0,0]; +[ 0, 0, 0, 0, 0, 0 ] +gap> v3 := Vector(IsPlistVectorRep, Rationals, l3); + +gap> PositionNonZero( v3 ); +7 +gap> v4 := Vector(GF(5), l1*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> PositionNonZero( v4 ); +1 +gap> v5 := Vector(GF(5), l2*One(GF(5))); +[ 0*Z(5), Z(5), Z(5)^3, 0*Z(5), 0*Z(5), Z(5)^0 ] +gap> PositionNonZero( v5 ); +2 +gap> v6 := Vector(GF(5), l3*One(GF(5))); +[ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ] +gap> PositionNonZero( v6 ); +7 +gap> STOP_TEST( "PositionNonZero.tst", 1); diff --git a/tst/testinstall/MatrixObj/Randomize.tst b/tst/testinstall/MatrixObj/Randomize.tst new file mode 100644 index 0000000000..8680965c7f --- /dev/null +++ b/tst/testinstall/MatrixObj/Randomize.tst @@ -0,0 +1,22 @@ +gap> START_TEST("Randomize.tst"); +gap> ll := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, ll); + +gap> Randomize( v1 ); +gap> Unpack( v1 ); +[ -2/3, 2, 1, -4, 0, 1 ] +gap> Randomize( v1 ); +gap> Unpack( v1 ); +[ -1, -1, 1/2, 0, -2, 1/2 ] +gap> v2 := Vector(GF(5), ll*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> Randomize( v2 ); +[ 0*Z(5), Z(5)^2, Z(5)^0, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> v2; +[ 0*Z(5), Z(5)^2, Z(5)^0, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> Randomize( v2 ); +[ Z(5), Z(5)^0, 0*Z(5), Z(5)^3, Z(5)^3, Z(5)^3 ] +gap> v2; +[ Z(5), Z(5)^0, 0*Z(5), Z(5)^3, Z(5)^3, Z(5)^3 ] +gap> STOP_TEST( "Randomize.tst", 1); diff --git a/tst/testinstall/MatrixObj/TraceMat.tst b/tst/testinstall/MatrixObj/TraceMat.tst new file mode 100644 index 0000000000..fc916a8abf --- /dev/null +++ b/tst/testinstall/MatrixObj/TraceMat.tst @@ -0,0 +1,16 @@ +gap> START_TEST( "TraceMat.tst" ); +gap> l := [[1,2],[3,4]]; +[ [ 1, 2 ], [ 3, 4 ] ] +gap> m1 := Matrix(l); +<2x2-matrix over Rationals> +gap> m2 := Matrix(Integers,l); +<2x2-matrix over Integers> +gap> m3 := Matrix(GF(7),l*One(GF(7))); +[ [ Z(7)^0, Z(7)^2 ], [ Z(7), Z(7)^4 ] ] +gap> TraceMat( m1 ); +5 +gap> TraceMat( m2 ); +5 +gap> TraceMat( m3 ); +Z(7)^5 +gap> STOP_TEST( "TraceMat.tst", 1); diff --git a/tst/testinstall/MatrixObj/Unpack.tst b/tst/testinstall/MatrixObj/Unpack.tst new file mode 100644 index 0000000000..7b52ca20bd --- /dev/null +++ b/tst/testinstall/MatrixObj/Unpack.tst @@ -0,0 +1,12 @@ +gap> START_TEST("Unpack.tst"); +gap> ll := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, ll); + +gap> Unpack( v1 ); +[ 1, 2, 3, 4, 5, 6 ] +gap> v5 := Vector(GF(5), ll*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> Unpack( v5 ); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> STOP_TEST( "Unpack.tst", 1); diff --git a/tst/testinstall/MatrixObj/WeightOfVector.tst b/tst/testinstall/MatrixObj/WeightOfVector.tst new file mode 100644 index 0000000000..f97ff59a71 --- /dev/null +++ b/tst/testinstall/MatrixObj/WeightOfVector.tst @@ -0,0 +1,32 @@ +gap> START_TEST( "WeightOfVector.tst" ); +gap> l1 := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, l1); + +gap> WeightOfVector( v1 ); +6 +gap> l2 := [0,2,3,0,5,6]; +[ 0, 2, 3, 0, 5, 6 ] +gap> v2 := Vector(IsPlistVectorRep, Rationals, l2); + +gap> WeightOfVector( v2 ); +4 +gap> l3 := [0,0,0,0,0,0]; +[ 0, 0, 0, 0, 0, 0 ] +gap> v3 := Vector(IsPlistVectorRep, Rationals, l3); + +gap> WeightOfVector( v3 ); +0 +gap> v4 := Vector(GF(5), l1*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> WeightOfVector( v4 ); +5 +gap> v5 := Vector(GF(5), l2*One(GF(5))); +[ 0*Z(5), Z(5), Z(5)^3, 0*Z(5), 0*Z(5), Z(5)^0 ] +gap> WeightOfVector( v5 ); +3 +gap> v6 := Vector(GF(5), l3*One(GF(5))); +[ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ] +gap> WeightOfVector( v6 ); +0 +gap> STOP_TEST( "WeightOfVector.tst", 1); diff --git a/tst/testinstall/MatrixObj/ZeroVector.tst b/tst/testinstall/MatrixObj/ZeroVector.tst new file mode 100644 index 0000000000..00ffbad7b9 --- /dev/null +++ b/tst/testinstall/MatrixObj/ZeroVector.tst @@ -0,0 +1,17 @@ +gap> START_TEST("ZeroVector.tst"); +gap> l1 := [1,2,3,4,5,6]; +[ 1, 2, 3, 4, 5, 6 ] +gap> v1 := Vector(IsPlistVectorRep, Rationals, l1); + +gap> v0 := ZeroVector( 15, v1 ); + +gap> Unpack( v0 ); +[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ] +gap> v3 := Vector(GF(5), l1*One(GF(5))); +[ Z(5)^0, Z(5), Z(5)^3, Z(5)^2, 0*Z(5), Z(5)^0 ] +gap> v0 := ZeroVector( 12, v3 ); +< mutable compressed vector length 12 over GF(5) > +gap> Unpack( v0 ); +[ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), + 0*Z(5), 0*Z(5), 0*Z(5) ] +gap> STOP_TEST("ZeroVector.tst",1);