From b1bd112a4c4b01a03ac1f63e3a33d57121b15e5e Mon Sep 17 00:00:00 2001 From: Alexander Hulpke Date: Wed, 16 Feb 2022 15:38:28 -0700 Subject: [PATCH 1/6] PositionNonZero for vector objects with start value --- lib/matobjplist.gi | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/lib/matobjplist.gi b/lib/matobjplist.gi index b484cb273d..861799927c 100644 --- a/lib/matobjplist.gi +++ b/lib/matobjplist.gi @@ -280,6 +280,12 @@ InstallMethod( PositionNonZero, "for a plist vector", [ IsPlistVectorRep ], return PositionNonZero( v![ELSPOS] ); end ); +InstallOtherMethod( PositionNonZero, "for a plist vector and start", + [ IsPlistVectorRep,IsInt ], + function( v,s ) + return PositionNonZero( v![ELSPOS],s ); + end ); + InstallMethod( PositionLastNonZero, "for a plist vector", [ IsPlistVectorRep ], function( v ) local els,i; From f65e003a34206a3a49471c7bdc0d512baa21fa23 Mon Sep 17 00:00:00 2001 From: Alexander Hulpke Date: Wed, 16 Feb 2022 15:51:13 -0700 Subject: [PATCH 2/6] ENHANCE Methods for matrix objects Relaxed some method installations, respectively almost duplicated some methods, so that basic linear algebra becomes available for matrix objects (at least over fields). This might be far from optimal and is just a first approximation. --- lib/matrix.gi | 216 ++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 199 insertions(+), 17 deletions(-) diff --git a/lib/matrix.gi b/lib/matrix.gi index 011774208d..4aef74f6e0 100644 --- a/lib/matrix.gi +++ b/lib/matrix.gi @@ -2153,6 +2153,22 @@ InstallMethod( NullspaceMat, [ IsOrdinaryMatrix ], mat -> SemiEchelonMatTransformation(mat).relations ); +InstallOtherMethod(NullspaceMat,"matrix objects",[IsMatrixObj], +function(mat) +local r,nr,nc,n,i,j; + r:=SemiEchelonMatTransformation(mat).relations; + if Length(r)=0 then return r;fi; + nr:=Length(r); + nc:=Length(r[1]); + n:=ZeroMatrix(BaseDomain(mat),nr,nc); + for i in [1..nr] do + for j in [1..nc] do + n[i,j]:=r[i][j]; + od; + od; + return n; +end); + InstallMethod( NullspaceMatDestructive, "generic method for ordinary matrices", [ IsOrdinaryMatrix and IsMutable], @@ -2348,9 +2364,9 @@ InstallOtherMethod( RankMat, ## #M SemiEchelonMat( ) ## -InstallMethod( SemiEchelonMatDestructive, +InstallOtherMethod( SemiEchelonMatDestructive, "generic method for matrices", - [ IsMatrix and IsMutable ], + [ IsMatrixOrMatrixObj and IsMutable ], function( mat ) local zero, # zero of the field of nrows, # number of rows in @@ -2440,10 +2456,14 @@ InstallMethod( SemiEchelonMatTransformation, return SemiEchelonMatTransformationDestructive( copymat ); end); -InstallMethod( SemiEchelonMatTransformationDestructive, - "generic method for matrices", - [ IsMatrix and IsMutable], - function( mat ) +InstallOtherMethod(SemiEchelonMatTransformation,"matrix objects",[IsMatrixObj], +function( mat ) + local copymat; + copymat:=MutableCopyMatrix(mat); + return SemiEchelonMatTransformationDestructive( copymat ); +end); + +BindGlobal("DoSemiEchelonMatTransformationDestructive", function(f, mat ) local zero, # zero of the field of nrows, # number of rows in ncols, # number of columns in @@ -2454,12 +2474,11 @@ InstallMethod( SemiEchelonMatTransformationDestructive, T, # transformation matrix coeffs, # list of coefficient vectors for 'vectors' relations, # basis vectors of the null space of 'mat' - row, head, x, row2,f; + row, head, x, row2; nrows := NrRows( mat ); ncols := NrCols( mat ); - f := DefaultFieldOfMatrix(mat); if f = fail then f := mat[1,1]; fi; @@ -2514,6 +2533,20 @@ InstallMethod( SemiEchelonMatTransformationDestructive, end ); +InstallOtherMethod( SemiEchelonMatTransformationDestructive, + "generic method for matrices", + [ IsMatrix and IsMutable], + mat->DoSemiEchelonMatTransformationDestructive( + DefaultFieldOfMatrix(mat),mat)); + +InstallOtherMethod( SemiEchelonMatTransformationDestructive, + "generic method for matrix objects over fields", + [ IsMatrixObj and IsRowListMatrix and IsMutable],function(mat) +local f; + f:=BaseDomain(mat); + if not IsField(f) then TryNextMethod();fi; + return DoSemiEchelonMatTransformationDestructive(f,mat); +end); ############################################################################# @@ -2783,11 +2816,11 @@ end); ## ## One solution of * = or `fail'. ## -InstallMethod( SolutionMatDestructive, +InstallOtherMethod( SolutionMatDestructive, "generic method", IsCollsElms, - [ IsOrdinaryMatrix and IsMutable, - IsRowVector and IsMutable], + [ IsMatrixOrMatrixObj and IsMutable, + IsListOrCollection and IsMutable], function( mat, vec ) local i,ncols,sem, vno, z,x, row, sol; ncols := Length(vec); @@ -2890,11 +2923,10 @@ end); #end ); # -InstallMethod( SolutionMat, +InstallOtherMethod( SolutionMat, "generic method for ordinary matrix and vector", IsCollsElms, - [ IsOrdinaryMatrix, - IsRowVector ], + [ IsMatrixOrMatrixObj,IsListOrCollection ], function ( mat, vec ) return SolutionMatDestructive( MutableCopyMatrix( mat ), ShallowCopy( vec ) ); @@ -2921,7 +2953,7 @@ end ); ## second position a base of the intersection. Both bases are in ## semi-echelon form. ## -InstallMethod( SumIntersectionMat, +InstallMethod( SumIntersectionMat,"ordinary matrices", IsIdenticalObj, [ IsMatrix, IsMatrix ], function( M1, M2 ) @@ -2988,9 +3020,79 @@ function( M1, M2 ) return [ sum, int ]; end ); +InstallOtherMethod( SumIntersectionMat,"MatricObject", IsIdenticalObj, + [ IsMatrixObj, IsMatrixObj ], +function( M1, M2 ) + local n, # number of columns + mat, # matrix for Zassenhaus algorithm + zero, # zero vector + v, # loop over 'M1' and 'M2' + heads, # list of leading positions + sum, # base of the sum + i, # loop over rows of 'mat' + int; # base of the intersection + + if NrRows( M1 ) = 0 then + return [ M2, M1 ]; + elif NrRows( M2 ) = 0 then + return [ M1, M2 ]; + elif NrCols( M1 ) <> NrCols( M2 ) then + Error( "dimensions of matrices are not compatible" ); + elif ZeroOfBaseDomain( M1 ) <> ZeroOfBaseDomain( M2 ) then + Error( "fields of matrices are not compatible" ); + fi; + + n:= NrCols( M1 ); + mat:= []; + zero:= ZeroVector(n,M1); + + # Set up the matrix for Zassenhaus' algorithm. + mat:= ZeroMatrix(BaseDomain(M1),NrRows(M1)+NrRows(M2),2*n); + + i:=1; + for v in List(M1) do + CopySubVector(v,mat[i],[1..n],[1..n]); + CopySubVector(v,mat[i],[1..n],[n+1..2*n]); + #v:= ShallowCopy( v ); + #Append( v, v ); + i:=i+1; + od; + for v in List(M2) do + CopySubVector(v,mat[i],[1..n],[1..n]); + #mat[i]{[n+1..2*n]}:=zero; # not needed, as initially 0 + #v:= ShallowCopy( v ); + #Append( v, zero ); + i:=i+1; + od; + + # Transform `mat' into semi-echelon form. + mat := SemiEchelonMatDestructive( mat ); + heads := mat.heads; + mat := mat.vectors; + + # Extract the bases for the sum \ldots + sum:= []; + for i in [ 1 .. n ] do + if heads[i] <> 0 then + Add( sum, mat[ heads[i] ]{ [ 1 .. n ] } ); + fi; + od; + + # \ldots and the intersection. + int:= []; + for i in [ n+1 .. Length( heads ) ] do + if heads[i] <> 0 then + Add( int, mat[ heads[i] ]{ [ n+1 .. 2*n ] } ); + fi; + od; + + # return the result + return [ sum, int ]; +end ); + ############################################################################# -## +#end #M TriangulizeMat( ) . . . . . bring a matrix in upper triangular form ## InstallMethod( TriangulizeMat, @@ -3065,6 +3167,78 @@ InstallMethod( TriangulizeMat, Info( InfoMatrix, 1, "TriangulizeMat returns" ); end ); +InstallOtherMethod( TriangulizeMat, + "generic method for mutable matrix obj", + [ IsMatrixObj and IsMutable ], + function ( mat ) + local m, n, i, j, k, row, zero, x, row2; + + Info( InfoMatrix, 1, "TriangulizeMat called" ); + + m := NrRows(mat); + if m>0 then + + # get the size of the matrix + n := NrCols(mat); + zero := ZeroOfBaseDomain( mat ); + + # make sure that the rows are mutable + for i in [ 1 .. m ] do + if not IsMutable( mat[i] ) then + mat[i]:= ShallowCopy( mat[i] ); + fi; + od; + + # run through all columns of the matrix + i := 0; + for k in [1..n] do + + # find a nonzero entry in this column + j := i + 1; + while j <= m and mat[j,k] = zero do j := j + 1; od; + + # if there is a nonzero entry + if j <= m then + + # increment the rank + Info( InfoMatrix, 2, " nonzero columns: ", k ); + i := i + 1; + + # make its row the current row and normalize it + row := mat[j]; + mat[j] := mat[i]; + x:= Inverse( row[k] ); + if x = fail then + TryNextMethod(); + fi; + MultVector( row, x ); + mat[i] := row; + + # clear all entries in this column + for j in [1..i-1] do + row2 := mat[j]; + x := row2[k]; + if x <> zero then + AddRowVector( row2, row, - x ); + fi; + od; + for j in [i+1..m] do + row2 := mat[j]; + x := row2[k]; + if x <> zero then + AddRowVector( row2, row, - x ); + fi; + od; + + fi; + + od; + + fi; + + Info( InfoMatrix, 1, "TriangulizeMat returns" ); +end ); + InstallOtherMethod( TriangulizeMat, "for an empty list", @@ -3079,6 +3253,14 @@ local m; return m; end); +InstallOtherMethod( TriangulizedMat," matrix objects", [IsMatrixObj], +function ( mat ) +local m; + m:=MutableCopyMatrix(mat); + TriangulizeMat(m); + return m; +end); + InstallOtherMethod( TriangulizedMat, "for an empty list", [ IsList and IsEmpty ], mat -> []); @@ -4202,7 +4384,7 @@ BindGlobal("POW_MAT_INT", function(mat, n) local d, addb, trafo, value, t, ti, mm, pol, ind; d := NrRows(mat); # finding a better break even point probably also depends on q - if n < 2^QuoInt(3*d,4) then + if n < 2^QuoInt(3*d,4) or not IsField(DefaultFieldOfMatrix(mat)) then return POW_OBJ_INT(mat, n); fi; # helper function to build up a semi-echelon basis From ce880cfdb9a2d445475d670b8d0b6ea8e3c82906 Mon Sep 17 00:00:00 2001 From: Alexander Hulpke Date: Wed, 16 Feb 2022 15:53:34 -0700 Subject: [PATCH 3/6] ENHANCE: Matrix objects for ZmodnZ Create matrix/vector objects for (Integers mod m) (as ZmodnZ) that store the entries as list of integers and modulo-reduce in the arithmetic. This is significantly faster than to keep lists of lists of ZmodnZObj, since then there is a permanent packing/unpacking of ZmodnZObj. --- lib/matobj.gi | 6 + lib/matobjnz.gd | 19 + lib/matobjnz.gi | 1543 +++++++++++++++++++++++++++++++++++++++++++++++ lib/read3.g | 1 + lib/read5.g | 1 + 5 files changed, 1570 insertions(+) create mode 100644 lib/matobjnz.gd create mode 100644 lib/matobjnz.gi diff --git a/lib/matobj.gi b/lib/matobj.gi index 0dc81b0b20..3dd194e6eb 100644 --- a/lib/matobj.gi +++ b/lib/matobj.gi @@ -147,6 +147,8 @@ function( basedomain ) return IsGF2VectorRep; elif IsFinite(basedomain) and IsField(basedomain) and Size(basedomain) <= 256 then return Is8BitVectorRep; + elif IsFinite(basedomain) and IsZmodnZObj(Zero(basedomain)) then + return IsZmodnZVectorRep; fi; return IsPlistVectorRep; end); @@ -156,6 +158,8 @@ function( basedomain ) return IsGF2MatrixRep; elif IsFinite(basedomain) and IsField(basedomain) and Size(basedomain) <= 256 then return Is8BitMatrixRep; + elif IsFinite(basedomain) and IsZmodnZObj(Zero(basedomain)) then + return IsZmodnZMatrixRep; fi; return IsPlistMatrixRep; end); @@ -1345,6 +1349,8 @@ InstallMethod( Characteristic, [ IsMatrixOrMatrixObj ], M -> Characteristic( BaseDomain( M ) ) ); +InstallOtherMethod(DefaultFieldOfMatrix,[IsMatrixOrMatrixObj],BaseDomain); + ############################################################################# ## diff --git a/lib/matobjnz.gd b/lib/matobjnz.gd new file mode 100644 index 0000000000..520d8a19a4 --- /dev/null +++ b/lib/matobjnz.gd @@ -0,0 +1,19 @@ +############################################################################# +## +## This file is part of GAP, a system for computational discrete algebra. +## +## SPDX-License-Identifier: GPL-2.0-or-later +## +## Copyright of GAP belongs to its developers, whose names are too numerous +## to list here. Please refer to the COPYRIGHT file for details. +## + +# represent vectors/matrices over Z/nZ by nonnegative integer lists +# in the range [0..n-1], but reduce after +# arithmetic. This way avoid always wrapping all entries separately + +DeclareRepresentation( "IsZmodnZVectorRep", + IsVectorObj and IsPositionalObjectRep and HasBaseDomain, [] ); + +DeclareRepresentation( "IsZmodnZMatrixRep", + IsRowListMatrix and IsPositionalObjectRep and HasBaseDomain, [] ); diff --git a/lib/matobjnz.gi b/lib/matobjnz.gi new file mode 100644 index 0000000000..227154fe8a --- /dev/null +++ b/lib/matobjnz.gi @@ -0,0 +1,1543 @@ +############################################################################# +## +## This file is part of GAP, a system for computational discrete algebra. +## +## SPDX-License-Identifier: GPL-2.0-or-later +## +## Copyright of GAP belongs to its developers, whose names are too numerous +## to list here. Please refer to the COPYRIGHT file for details. +## + +# represent vectors/matrices over Z/nZ by nonnegative integer lists +# in the range [0..n-1], but reduce after +# arithmetic. This way avoid always wrapping all entries separately + +BindGlobal("ZNZVECREDUCE",function(v,l,m) +local i; + for i in [1..l] do + if v[i]<0 or v[i]>=m then v[i]:=v[i] mod m;fi; + od; +end); + +InstallMethod( ConstructingFilter, "for a zmodnz vector", + [ IsZmodnZVectorRep ], + function( v ) + return IsZmodnZVectorRep; + end ); + +InstallOtherMethod( ConstructingFilter, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + return IsZmodnZMatrixRep; + end ); + +InstallMethod( CompatibleVectorFilter, "zmodnz", + [ IsZmodnZMatrixRep ], + M -> IsZmodnZVectorRep ); + +############################################################################ +# Vectors +############################################################################ + +InstallMethod( NewVector, "for IsZmodnZVectorRep, a ring, and a list", + [ IsZmodnZVectorRep, IsRing, IsList ], + function( filter, basedomain, l ) + local typ, v; + typ:=NewType(FamilyObj(basedomain),IsZmodnZVectorRep and IsMutable and + CanEasilyCompareElements); + # force list of integers + if FamilyObj(basedomain)=FamilyObj(l) then + l:=List(l,Int); + else + l:=ShallowCopy(l); + fi; + v := [basedomain,l]; + Objectify(typ,v); + return v; + end ); + +InstallMethod( NewZeroVector, "for IsZmodnZVectorRep, a ring, and an int", + [ IsZmodnZVectorRep, IsRing, IsInt ], + function( filter, basedomain, l ) + local typ, v; + typ:=NewType(FamilyObj(basedomain),IsZmodnZVectorRep and IsMutable and + CanEasilyCompareElements); + # represent list as integers + v := [basedomain,0*[1..l]]; + Objectify(typ,v); + return v; + end ); + +InstallMethod( ViewObj, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + if not IsMutable(v) then + Print(""); + end ); + +InstallMethod( PrintObj, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + Print("NewVector(IsZmodnZVectorRep"); + if IsField(v![BDPOS]) then + Print(",GF(",Size(v![BDPOS]),"),",v![ELSPOS],")"); + else + Print(",",String(v![BDPOS]),",",v![ELSPOS],")"); + fi; + end ); + +InstallMethod( String, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + local st; + st := "NewVector(IsZmodnZVectorRep"; + if IsField(v![BDPOS]) then + Append(st,Concatenation( ",GF(",String(Size(v![BDPOS])),"),", + String(v![ELSPOS]),")" )); + else + Append(st,Concatenation( ",",String(v![BDPOS]),",", + String(v![ELSPOS]),")" )); + fi; + return st; + end ); + +InstallMethod( Display, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + Print( "\n"); + end ); + +InstallMethod( BaseDomain, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + return v![BDPOS]; + end ); + +InstallMethod( Length, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + return Length(v![ELSPOS]); + end ); + +InstallMethod( ShallowCopy, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + local res; + res := Objectify(TypeObj(v),[v![BDPOS],ShallowCopy(v![ELSPOS])]); + if not IsMutable(v) then SetFilterObj(res,IsMutable); fi; + return res; + end ); + +# StructuralCopy works automatically + +InstallMethod( PostMakeImmutable, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + MakeImmutable( v![ELSPOS] ); + end ); + +############################################################################ +# Representation preserving constructors: +############################################################################ + +# not needed according to MH +# InstallMethod( ZeroVector, "for an integer and a zmodnz vector", +# [ IsInt, IsZmodnZVectorRep ], +# function( l, t ) +# local v; +# v := Objectify(TypeObj(t), +# [t![BDPOS],ListWithIdenticalEntries(l,0)]); +# if not IsMutable(v) then SetFilterObj(v,IsMutable); fi; +# return v; +# end ); +# +# InstallMethod( ZeroVector, "for an integer and a zmodnz matrix", +# [ IsInt, IsZmodnZMatrixRep ], +# function( l, m ) +# local v; +# v := Objectify(TypeObj(m![EMPOS]), +# [m![BDPOS],ListWithIdenticalEntries(l,0)]); +# if not IsMutable(v) then SetFilterObj(v,IsMutable); fi; +# return v; +# end ); + +InstallMethod( Vector, "for a plain list and a zmodnz vector",IsIdenticalObj, + [ IsList and IsPlistRep, IsZmodnZVectorRep ], + function( l, t ) + local v; + # force list of integers + if FamilyObj(t![BDPOS])=FamilyObj(l) then l:=List(l,Int); fi; + v := Objectify(TypeObj(t),[t![BDPOS],l]); + if not IsMutable(v) then SetFilterObj(v,IsMutable); fi; + return v; + end ); + +InstallMethod( Vector, "for a list and a zmodnz vector", + [ IsList, IsZmodnZVectorRep ], + function( l, t ) + local v; + v := ShallowCopy(l); + if IsGF2VectorRep(l) then + PLAIN_GF2VEC(v); + elif Is8BitVectorRep(l) then + PLAIN_VEC8BIT(v); + fi; + v := Objectify(TypeObj(t),[t![BDPOS],v]); + if not IsMutable(v) then SetFilterObj(v,IsMutable); fi; + return v; + end ); + + +############################################################################ +# A selection of list operations: +############################################################################ + +InstallMethod( \[\], "for a zmodnz vector and a positive integer", + [ IsZmodnZVectorRep, IsPosInt ], + function( v, p ) + return ZmodnZObj(ElementsFamily(FamilyObj(v)),v![ELSPOS][p]); + end ); + +InstallMethod( \[\]\:\=, "for a zmodnz vector, a positive integer, and an obj", + [ IsZmodnZVectorRep, IsPosInt, IsObject ], + function( v, p, ob ) + v![ELSPOS][p] := Int(ob); + end ); + +InstallMethod( \{\}, "for a zmodnz vector and a list", + [ IsZmodnZVectorRep, IsList ], + function( v, l ) + return Objectify(TypeObj(v),[v![BDPOS],v![ELSPOS]{l}]); + end ); + +InstallMethod( PositionNonZero, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + return PositionNonZero( v![ELSPOS] ); + end ); + +InstallMethod( PositionNonZero, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + return PositionNonZero( v![ELSPOS] ); + end ); + +InstallOtherMethod( PositionNonZero, "for a zmodnz vector and start", + [ IsZmodnZVectorRep,IsInt ], + function( v,s ) + return PositionNonZero( v![ELSPOS],s ); + end ); + +InstallMethod( PositionLastNonZero, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + local els,i; + els := v![ELSPOS]; + i := Length(els); + while i > 0 and IsZero(els[i]) do i := i - 1; od; + return i; + end ); + +InstallMethod( ListOp, "for a zmodnz vector", [ IsZmodnZVectorRep ], +function( v ) +local fam; + fam:=ElementsFamily(FamilyObj(v)); + return List([1..Length(v![ELSPOS])],x->ZmodnZObj(fam,v![ELSPOS][x])); +end ); + +InstallMethod( ListOp, "for a zmodnz vector and a function", + [ IsZmodnZVectorRep, IsFunction ], +function( v, f ) +local fam; + fam:=ElementsFamily(FamilyObj(v)); + return List(List([1..Length(v![ELSPOS])],x->ZmodnZObj(fam,v![ELSPOS][x])),f); +end ); + +InstallMethod( Unpack, "for a zmodnz vector", + [ IsZmodnZVectorRep ], +function( v ) +local fam; + fam:=ElementsFamily(FamilyObj(v)); + return List([1..Length(v![ELSPOS])],x->ZmodnZObj(fam,v![ELSPOS][x])); +end ); + +############################################################################ +# Arithmetical operations: +############################################################################ + +InstallMethod( \+, "for two zmodnz vectors",IsIdenticalObj, + [ IsZmodnZVectorRep, IsZmodnZVectorRep ], +function( a, b ) +local ty,i,m,mu; + if not IsMutable(a) and IsMutable(b) then + ty := TypeObj(b); + else + ty := TypeObj(a); + fi; + m:=Size(a![BDPOS]); + b:=SUM_LIST_LIST_DEFAULT(a![ELSPOS],b![ELSPOS]); + if not IsMutable(b) then mu:=true;b:=ShallowCopy(b); + else mu:=false;fi; + for i in [1..Length(b)] do if b[i]>=m then b[i]:=b[i] mod m;fi;od; + if mu then MakeImmutable(b);fi; + return Objectify(ty,[a![BDPOS],b]); +end ); + +InstallOtherMethod( \+, "for zmodnz vector and plist",IsIdenticalObj, + [ IsZmodnZVectorRep, IsList ], +function( a, b ) + return a+Vector(BaseDomain(a),b); +end ); + +InstallOtherMethod( \+, "for plist and zmodnz vector",IsIdenticalObj, + [ IsList,IsZmodnZVectorRep ], +function( a, b ) + return Vector(BaseDomain(b),a)+b; +end ); + +InstallMethod( \-, "for two zmodnz vectors",IsIdenticalObj, + [ IsZmodnZVectorRep, IsZmodnZVectorRep ], +function( a, b ) +local ty,i,m,mu; + if not IsMutable(a) and IsMutable(b) then + ty := TypeObj(b); + else + ty := TypeObj(a); + fi; + m:=Size(a![BDPOS]); + b:=a![ELSPOS] - b![ELSPOS]; + if not IsMutable(b) then mu:=true;b:=ShallowCopy(b); + else mu:=false;fi; + for i in [1..Length(b)] do if b[i]<0 then b[i]:=b[i] mod m;fi;od; + if mu then MakeImmutable(b);fi; + return Objectify(ty,[a![BDPOS],b]); +end ); + +InstallOtherMethod( \-, "for zmodnz vector and plist",IsIdenticalObj, + [ IsZmodnZVectorRep, IsList ], +function( a, b ) + return a-Vector(BaseDomain(a),b); +end ); + +InstallOtherMethod( \-, "for plist and zmodnz vector",IsIdenticalObj, + [ IsList,IsZmodnZVectorRep ], +function( a, b ) + return Vector(BaseDomain(b),a)-b; +end ); + +InstallMethod( \=, "for two zmodnz vectors",IsIdenticalObj, + [ IsZmodnZVectorRep, IsZmodnZVectorRep ], + function( a, b ) + return EQ_LIST_LIST_DEFAULT(a![ELSPOS],b![ELSPOS]); + end ); + +InstallMethod( \=, "for zmodnz vector and plist",IsIdenticalObj, + [ IsZmodnZVectorRep, IsPlistRep ], +function( a, b ) + return a![ELSPOS]=List(b,x->x![1]); +end ); + +InstallMethod( \=, "for plist an zmodnz vector",IsIdenticalObj, + [ IsPlistRep,IsZmodnZVectorRep], +function(b,a) + return a![ELSPOS]=List(b,x->x![1]); +end ); + +InstallMethod( \<, "for two zmodnz vectors",IsIdenticalObj, + [ IsZmodnZVectorRep, IsZmodnZVectorRep ], + function( a, b ) + return LT_LIST_LIST_DEFAULT(a![ELSPOS],b![ELSPOS]); + end ); + +InstallMethod( AddRowVector, "for two zmodnz vectors", + [ IsZmodnZVectorRep and IsMutable, IsZmodnZVectorRep ], +function( a, b ) +local i,m; + a:=a![ELSPOS]; + ADD_ROW_VECTOR_2_FAST( a, b![ELSPOS] ); + m:=Size(b![BDPOS]); + for i in [1..Length(a)] do if a[i]>=m then a[i]:=a[i] mod m;fi;od; +end ); + +InstallMethod( AddRowVector, "for two zmodnz vectors, and a scalar", + [ IsZmodnZVectorRep and IsMutable, IsZmodnZVectorRep, IsObject ], +function( a, b, s ) +local i,m; + if IsZmodnZObj(s) then s:=Int(s);fi; + a:=a![ELSPOS]; + if IsSmallIntRep(s) then + ADD_ROW_VECTOR_3_FAST( a, b![ELSPOS], s ); + else + ADD_ROW_VECTOR_3( a, b![ELSPOS], s ); + fi; + m:=Size(b![BDPOS]); + if s>=0 then + for i in [1..Length(a)] do if a[i]>=m then a[i]:=a[i] mod m;fi;od; + else + for i in [1..Length(a)] do if a[i]<0 then a[i]:=a[i] mod m;fi;od; + fi; +end ); + +InstallOtherMethod( AddRowVector, "for zmodnz vector, plist, and a scalar", + [ IsZmodnZVectorRep and IsMutable, IsPlistRep, IsObject ], +function( a, b, s ) +local i,m; + if not ForAll(b,IsModulusRep) then TryNextMethod();fi; + if IsZmodnZObj(s) then s:=Int(s);fi; + m:=Size(a![BDPOS]); + a:=a![ELSPOS]; + b:=List(b,x->x![1]); + + if IsSmallIntRep(s) then + ADD_ROW_VECTOR_3_FAST( a, b, s ); + else + ADD_ROW_VECTOR_3( a, b, s ); + fi; + if s>=0 then + for i in [1..Length(a)] do if a[i]>=m then a[i]:=a[i] mod m;fi;od; + else + for i in [1..Length(a)] do if a[i]<0 then a[i]:=a[i] mod m;fi;od; + fi; +end ); + +InstallMethod( AddRowVector, + "for two zmodnz vectors, a scalar, and two positions", + [ IsZmodnZVectorRep and IsMutable, IsZmodnZVectorRep, + IsObject, IsPosInt, IsPosInt ], +function( a, b, s, from, to ) +local i,m; + if IsZmodnZObj(s) then s:=Int(s);fi; + a:=a![ELSPOS]; + if IsSmallIntRep(s) then + ADD_ROW_VECTOR_5_FAST( a, b![ELSPOS], s, from, to ); + else + ADD_ROW_VECTOR_5( a, b![ELSPOS], s, from, to ); + fi; + m:=Size(b![BDPOS]); + if s>=0 then + for i in [1..Length(a)] do if a[i]>=m then a[i]:=a[i] mod m;fi;od; + else + for i in [1..Length(a)] do if a[i]<0 then a[i]:=a[i] mod m;fi;od; + fi; +end ); + +InstallMethod( MultVectorLeft, + "for a zmodnz vector, and an object", + [ IsZmodnZVectorRep and IsMutable, IsObject ], +function( v, s ) +local i,m; + m:=Size(v![BDPOS]); + if IsZmodnZObj(s) then s:=Int(s);fi; + v:=v![ELSPOS]; + MULT_VECTOR_2_FAST(v,s); + if s>=0 then + for i in [1..Length(v)] do if v[i]>=m then v[i]:=v[i] mod m;fi;od; + else + for i in [1..Length(v)] do if v[i]<0 then v[i]:=v[i] mod m;fi;od; + fi; +end ); + +# The four argument version of MultVectorLeft / ..Right uses the generic +# implementation in matobj.gi + +BindGlobal("ZMODNZVECSCAMULT", +function( w, s ) +local i,m,t,b,v; + t:=TypeObj(w); + b:=w![BDPOS]; + m:=Size(b); + if IsZmodnZObj(s) then s:=Int(s);fi; + v:=PROD_LIST_SCL_DEFAULT(w![ELSPOS],s); + if not IsMutable(v) then + v:=ShallowCopy(v); + fi; + if s>=0 then + for i in [1..Length(v)] do if v[i]>=m then v[i]:=v[i] mod m;fi;od; + else + for i in [1..Length(v)] do if v[i]<0 then v[i]:=v[i] mod m;fi;od; + fi; + if not IsMutable(w![ELSPOS]) then MakeImmutable(v);fi; + return Objectify(t,[b,v]); +end ); + +InstallMethod( \*, "for a zmodnz vector and a scalar", + [ IsZmodnZVectorRep, IsScalar ],ZMODNZVECSCAMULT); + +InstallMethod( \*, "for a scalar and a zmodnz vector", + [ IsScalar, IsZmodnZVectorRep ], +function( s, v ) + return ZMODNZVECSCAMULT(v,s); +end ); + +InstallMethod( \/, "for a zmodnz vector and a scalar", + [ IsZmodnZVectorRep, IsScalar ], +function( v, s ) + return ZMODNZVECSCAMULT(v,s^-1); +end ); + +BindGlobal("ZMODNZVECADDINVCLEANUP",function(m,l) +local i; + if IsMutable(l) then + for i in [1..Length(l)] do if l[i]<0 then l[i]:=l[i] mod m;fi;od; + else + l:=ShallowCopy(l); + for i in [1..Length(l)] do if l[i]<0 then l[i]:=l[i] mod m;fi;od; + MakeImmutable(l); + fi; + return l; +end); + +InstallMethod( AdditiveInverseSameMutability, "for a zmodnz vector", + [ IsZmodnZVectorRep ], + function( v ) + return Objectify( TypeObj(v), + [v![BDPOS],ZMODNZVECADDINVCLEANUP(Size(v![BDPOS]), + AdditiveInverseSameMutability(v![ELSPOS]))] ); + end ); + +InstallMethod( AdditiveInverseImmutable, "for a zmodnz vector", + [ IsZmodnZVectorRep ], + function( v ) + local res; + res := Objectify( TypeObj(v), + [v![BDPOS],ZMODNZVECADDINVCLEANUP(Size(v![BDPOS]), + AdditiveInverseSameMutability(v![ELSPOS]))] ); + MakeImmutable(res); + return res; + end ); + +InstallMethod( AdditiveInverseMutable, "for a zmodnz vector", + [ IsZmodnZVectorRep ], + function( v ) + local res; + res := Objectify(TypeObj(v), + [v![BDPOS],ZMODNZVECADDINVCLEANUP(Size(v![BDPOS]), + AdditiveInverseMutable(v![ELSPOS]))]); + if not IsMutable(v) then SetFilterObj(res,IsMutable); fi; + return res; + end ); + +# redundant according to MH +# InstallMethod( ZeroSameMutability, "for a zmodnz vector", [ IsZmodnZVectorRep ], +# function( v ) +# return Objectify(TypeObj(v),[v![BDPOS],ZeroSameMutability(v![ELSPOS])]); +# end ); +# +# InstallMethod( ZeroImmutable, "for a zmodnz vector", [ IsZmodnZVectorRep ], +# function( v ) +# local res; +# res := Objectify(TypeObj(v),[v![BDPOS],ZeroImmutable(v![ELSPOS])]); +# MakeImmutable(res); +# return res; +# end ); + +InstallMethod( ZeroMutable, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + local res; + res := Objectify(TypeObj(v), + [v![BDPOS],ZeroMutable(v![ELSPOS])]); + if not IsMutable(v) then SetFilterObj(res,IsMutable); fi; + return res; + end ); + +InstallMethod( IsZero, "for a zmodnz vector", [ IsZmodnZVectorRep ], + function( v ) + return IsZero( v![ELSPOS] ); + end ); + +#InstallMethodWithRandomSource( Randomize, +# "for a random source and a mutable zmodnz vector", +# [ IsRandomSource, IsZmodnZVectorRep and IsMutable ], +# function( rs, v ) +# local bd,i; +# bd := v![BDPOS]; +# for i in [1..Length(v![ELSPOS])] do +# v![ELSPOS][i] := Random( rs, bd ); +# od; +# return v; +# end ); + +InstallMethod( CopySubVector, "for two zmodnz vectors and two lists", + [ IsZmodnZVectorRep, IsZmodnZVectorRep and IsMutable, IsList, IsList ], + function( a,b,pa,pb ) + # The following should eventually go into the kernel: + b![ELSPOS]{pb} := a![ELSPOS]{pa}; + end ); + +############################################################################ +# Matrices +############################################################################ + +InstallMethod( NewMatrix, + "for IsZmodnZMatrixRep, a ring, an int, and a list", + [ IsZmodnZMatrixRep, IsRing, IsInt, IsList ], + function( filter, basedomain, rl, l ) + local nd, filterVectors, m, e, filter2, i; + + # If applicable then replace a flat list 'l' by a nested list + # of lists of length 'rl'. + if Length(l) > 0 and not IsVectorObj(l[1]) then + nd := NestingDepthA(l); + if nd < 2 or nd mod 2 = 1 then + if Length(l) mod rl <> 0 then + Error( "NewMatrix: Length of l is not a multiple of rl" ); + fi; + l := List([0,rl..Length(l)-rl], i -> l{[i+1..i+rl]}); + fi; + fi; + + filterVectors := IsZmodnZVectorRep; + m := 0*[1..Length(l)]; + for i in [1..Length(l)] do + if IsVectorObj(l[i]) and IsZmodnZVectorRep(l[i]) then + m[i] := ShallowCopy(l[i]); + else + m[i] := NewVector( filterVectors, basedomain, l[i] ); + fi; + od; + e := NewVector(filterVectors, basedomain, []); + m := [basedomain,e,rl,m]; + filter2 := IsZmodnZMatrixRep and IsMutable; + if HasCanEasilyCompareElements(Representative(basedomain)) and + CanEasilyCompareElements(Representative(basedomain)) then + filter2 := filter2 and CanEasilyCompareElements; + fi; + Objectify( NewType(CollectionsFamily(FamilyObj(basedomain)), + filter2), m ); + return m; + end ); + +InstallMethod( NewZeroMatrix, + "for IsZmodnZMatrixRep, a ring, and two ints", + [ IsZmodnZMatrixRep, IsRing, IsInt, IsInt ], + function( filter, basedomain, rows, cols ) + local m,i,e,filter2; + filter2 := IsZmodnZVectorRep; + m := 0*[1..rows]; + e := NewVector(filter2, basedomain, []); + for i in [1..rows] do + m[i] := ZeroVector( cols, e ); + od; + m := [basedomain,e,cols,m]; + Objectify( NewType(CollectionsFamily(FamilyObj(basedomain)), + filter and IsMutable), m ); + return m; + end ); + +InstallMethod( NewIdentityMatrix, + "for IsZmodnZMatrixRep, a ring, and an int", + [ IsZmodnZMatrixRep, IsRing, IsInt ], + function( filter, basedomain, dim ) + local mat, one, i; + mat := NewZeroMatrix(filter, basedomain, dim, dim); +# one := One(basedomain); + for i in [1..dim] do + mat[i,i] := 1; + od; + return mat; + end ); + +InstallOtherMethod( BaseDomain, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + return m![BDPOS]; + end ); + +InstallMethod( NumberRows, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + return Length(m![ROWSPOS]); + end ); + +InstallMethod( NumberColumns, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + return m![RLPOS]; + end ); + +# InstallMethod( DimensionsMat, "for a zmodnz matrix", +# [ IsZmodnZMatrixRep ], +# function( m ) +# return [Length(m![ROWSPOS]),m![RLPOS]]; +# end ); + + +############################################################################ +# Representation preserving constructors: +############################################################################ + +# redundant according to MH +# InstallMethod( ZeroMatrix, "for two integers and a zmodnz matrix", +# [ IsInt, IsInt, IsZmodnZMatrixRep ], +# function( rows,cols,m ) +# local l,t,res; +# t := m![EMPOS]; +# l := List([1..rows],i->ZeroVector(cols,t)); +# res := Objectify( TypeObj(m), [m![BDPOS],t,cols,l] ); +# if not IsMutable(m) then +# SetFilterObj(res,IsMutable); +# fi; +# return res; +# end ); + +InstallMethod( IdentityMatrix, "for an integer and a zmodnz matrix", + [ IsInt, IsZmodnZMatrixRep ], + function( rows,m ) + local i,l,o,t,res; + t := m![EMPOS]; + l := List([1..rows],i->ZeroVector(rows,t)); + o := One(m![BDPOS]); + for i in [1..rows] do + l[i][i] := o; + od; + res := Objectify( TypeObj(m), [m![BDPOS],t,rows,l] ); + if not IsMutable(m) then + SetFilterObj(res,IsMutable); + fi; + return res; + end ); + +InstallMethod( Matrix, "for a list and a zmodnz matrix", + [ IsList, IsInt, IsZmodnZMatrixRep ], + function( rows,rowlen,m ) + local i,l,nrrows,res,t; + t := m![EMPOS]; + if Length(rows) > 0 then + if IsVectorObj(rows[1]) and IsZmodnZVectorRep(rows[1]) then + nrrows := Length(rows); + l := rows; + elif IsList(rows[1]) then + nrrows := Length(rows); + l := ListWithIdenticalEntries(Length(rows),0); + for i in [1..Length(rows)] do + l[i] := Vector(rows[i],t); + od; + else # a flat initializer: + nrrows := Length(rows)/rowlen; + l := ListWithIdenticalEntries(nrrows,0); + for i in [1..nrrows] do + l[i] := Vector(rows{[(i-1)*rowlen+1..i*rowlen]},t); + od; + fi; + else + l := []; + nrrows := 0; + fi; + res := Objectify( TypeObj(m), [m![BDPOS],t,rowlen,l] ); + if not IsMutable(m) then + SetFilterObj(res,IsMutable); + fi; + return res; + end ); + +############################################################################ +# Printing and viewing methods: +############################################################################ + +InstallMethod( ViewObj, "for a zmodnz matrix", [ IsZmodnZMatrixRep ], + function( m ) + Print("<"); + if not IsMutable(m) then Print("immutable "); fi; + Print(Length(m![ROWSPOS]),"x",m![RLPOS],"-matrix over ",m![BDPOS],">"); + end ); + +InstallMethod( PrintObj, "for a zmodnz matrix", [ IsZmodnZMatrixRep ], + function( m ) + Print("NewMatrix(IsZmodnZMatrixRep"); + if IsFinite(m![BDPOS]) and IsField(m![BDPOS]) then + Print(",GF(",Size(m![BDPOS]),"),"); + else + Print(",",String(m![BDPOS]),","); + fi; + Print(NumberColumns(m),",",Unpack(m),")"); + end ); + +InstallMethod( Display, "for a zmodnz matrix", [ IsZmodnZMatrixRep ], + function( m ) + local i; + Print("<"); + if not IsMutable(m) then Print("immutable "); fi; + Print(Length(m![ROWSPOS]),"x",m![RLPOS],"-matrix over ",m![BDPOS],":\n"); + Display(List(m![ROWSPOS],x->x![ELSPOS])); +# for i in [1..Length(m![ROWSPOS])] do +# if i = 1 then +# Print("["); +# else +# Print(" "); +# fi; +# Print(m![ROWSPOS][i]![ELSPOS],"\n"); +# od; + Print("]>\n"); + end ); + +InstallMethod( String, "for zmodnz matrix", [ IsZmodnZMatrixRep ], + function( m ) + local st; + st := "NewMatrix(IsZmodnZMatrixRep"; + Add(st,','); + if IsFinite(m![BDPOS]) and IsField(m![BDPOS]) then + Append(st,"GF("); + Append(st,String(Size(m![BDPOS]))); + Append(st,"),"); + else + Append(st,String(m![BDPOS])); + Append(st,","); + fi; + Append(st,String(NumberColumns(m))); + Add(st,','); + Append(st,String(Unpack(m))); + Add(st,')'); + return st; + end ); + + +############################################################################ +# A selection of list operations: +############################################################################ + +InstallOtherMethod( \[\], "for a zmodnz matrix and a positive integer", +#T Once the declaration of '\[\]' for 'IsMatrixObj' disappears, +#T we can use 'InstallMethod'. + [ IsZmodnZMatrixRep, IsPosInt ], + function( m, p ) + return m![ROWSPOS][p]; + end ); + +InstallOtherMethod( \[\]\:\=, + "for a zmodnz matrix, a positive integer, and a zmodnz vector", + [ IsZmodnZMatrixRep and IsMutable, IsPosInt, IsZmodnZVectorRep ], + function( m, p, v ) + m![ROWSPOS][p] := v; + end ); + +InstallOtherMethod( \{\}, "for a zmodnz matrix and a list", + [ IsZmodnZMatrixRep, IsList ], + function( m, p ) + local l; + l := m![ROWSPOS]{p}; + return Objectify(TypeObj(m),[m![BDPOS],m![EMPOS],m![RLPOS],l]); + end ); + +InstallMethod( Add, "for a zmodnz matrix and a zmodnz vector", + [ IsZmodnZMatrixRep and IsMutable, IsZmodnZVectorRep ], + function( m, v ) + Add(m![ROWSPOS],v); + end ); + +InstallMethod( Add, "for a zmodnz matrix, a zmodnz vector, and a pos. int", + [ IsZmodnZMatrixRep and IsMutable, IsZmodnZVectorRep, IsPosInt ], + function( m, v, p ) + Add(m![ROWSPOS],v,p); + end ); + +InstallMethod( Remove, "for a zmodnz matrix", + [ IsZmodnZMatrixRep and IsMutable ], + m -> Remove( m![ROWSPOS] ) ); + +InstallMethod( Remove, "for a zmodnz matrix, and a position", + [ IsZmodnZMatrixRep and IsMutable, IsPosInt ], + function( m, p ) + Remove( m![ROWSPOS],p ); + end ); +#T must return the removed row if it was bound + +InstallMethod( IsBound\[\], "for a zmodnz matrix, and a position", + [ IsZmodnZMatrixRep, IsPosInt ], + function( m, p ) + return p <= Length(m![ROWSPOS]); + end ); + +InstallMethod( Unbind\[\], "for a zmodnz matrix, and a position", + [ IsZmodnZMatrixRep and IsMutable, IsPosInt ], + function( m, p ) + if p <> Length(m![ROWSPOS]) then + ErrorNoReturn("Unbind\\[\\]: Matrices must stay dense, you cannot Unbind here"); + fi; + Unbind( m![ROWSPOS][p] ); + end ); + +InstallMethod( \{\}\:\=, "for a zmodnz matrix, a list, and a zmodnz matrix", + [ IsZmodnZMatrixRep and IsMutable, IsList, + IsZmodnZMatrixRep ], + function( m, pp, n ) + m![ROWSPOS]{pp} := n![ROWSPOS]; + end ); + +InstallMethod( Append, "for two zmodnz matrices", + [ IsZmodnZMatrixRep and IsMutable, IsZmodnZMatrixRep ], + function( m, n ) + Append(m![ROWSPOS],n![ROWSPOS]); + end ); + +InstallMethod( ShallowCopy, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local res; + res := Objectify(TypeObj(m),[m![BDPOS],m![EMPOS],m![RLPOS], + ShallowCopy(m![ROWSPOS])]); + if not IsMutable(m) then + SetFilterObj(res,IsMutable); + fi; +#T 'ShallowCopy' MUST return a mutable object +#T if such an object exists at all! + return res; + end ); + +InstallMethod( PostMakeImmutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + MakeImmutable( m![ROWSPOS] ); + end ); + +InstallOtherMethod( ListOp, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + return List(m![ROWSPOS]); + end ); + +InstallOtherMethod( ListOp, "for a zmodnz matrix and a function", + [ IsZmodnZMatrixRep, IsFunction ], + function( m, f ) + return List(m![ROWSPOS],f); + end ); + +InstallOtherMethod( Unpack, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], +function( m ) +local fam; + fam:=ElementsFamily(FamilyObj(BaseDomain(m))); + return List(m![ROWSPOS],v-> + List([1..Length(v![ELSPOS])],x->ZmodnZObj(fam,v![ELSPOS][x]))); + end ); + + +InstallMethod( MutableCopyMat, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local l,res; + l := List(m![ROWSPOS],ShallowCopy); + res := Objectify(TypeObj(m),[m![BDPOS],m![EMPOS],m![RLPOS],l]); + if not IsMutable(m) then + SetFilterObj(res,IsMutable); + fi; + return res; + end); + +InstallMethod( ExtractSubMatrix, "for a zmodnz matrix, and two lists", + [ IsZmodnZMatrixRep, IsList, IsList ], + function( m, p, q ) + local i,l; + l := m![ROWSPOS]{p}; + for i in [1..Length(l)] do + l[i] := Objectify(TypeObj(l[i]),[l[i]![BDPOS],l[i]![ELSPOS]{q}]); + od; + return Objectify(TypeObj(m),[m![BDPOS],m![EMPOS],Length(q),l]); + end ); + +InstallMethod( CopySubMatrix, "for two zmodnz matrices and four lists", + [ IsZmodnZMatrixRep, IsZmodnZMatrixRep and IsMutable, + IsList, IsList, IsList, IsList ], + function( m, n, srcrows, dstrows, srccols, dstcols ) + local i; + # This eventually should go into the kernel without creating + # a intermediate objects: + for i in [1..Length(srcrows)] do + n![ROWSPOS][dstrows[i]]![ELSPOS]{dstcols} := + m![ROWSPOS][srcrows[i]]![ELSPOS]{srccols}; + od; + end ); + +# InstallOtherMethod( CopySubMatrix, +# "for two zmodnzs -- fallback in case of bad rep.", +# [ IsZmodnZRep, IsZmodnZRep and IsMutable, +# IsList, IsList, IsList, IsList ], +# function( m, n, srcrows, dstrows, srccols, dstcols ) +# local i; +# # in this representation all access probably has to go through the +# # generic method selection, so it is not clear whether there is an +# # improvement in moving this into the kernel. +# for i in [1..Length(srcrows)] do +# n[dstrows[i]]{dstcols}:=m[srcrows[i]]{srccols}; +# od; +# end ); + +InstallMethod( MatElm, "for a zmodnz matrix and two positions", + [ IsZmodnZMatrixRep, IsPosInt, IsPosInt ], + function( m, row, col ) + return ZmodnZObj(ElementsFamily(FamilyObj(m![BDPOS])), + m![ROWSPOS][row]![ELSPOS][col]); + end ); + +InstallMethod( SetMatElm, "for a zmodnz matrix, two positions, and an object", + [ IsZmodnZMatrixRep and IsMutable, IsPosInt, IsPosInt, IsObject ], + function( m, row, col, ob ) + m![ROWSPOS][row]![ELSPOS][col] := Int(ob); + end ); + + +############################################################################ +# Arithmetical operations: +############################################################################ + +InstallMethod( \+, "for two zmodnz matrices", + [ IsZmodnZMatrixRep, IsZmodnZMatrixRep ], + function( a, b ) + local ty; + 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 zmodnz matrices", + [ IsZmodnZMatrixRep, IsZmodnZMatrixRep ], + function( a, b ) + local ty; + 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 zmodnz matrices",IsIdenticalObj, + [ IsZmodnZMatrixRep, IsZmodnZMatrixRep ], + function( a, b ) + # Here we do full checking since it is rather cheap! + local i,j,l,ty,v,w,m,r; + if not IsMutable(a) and IsMutable(b) then + ty := TypeObj(b); + else + ty := TypeObj(a); + fi; + if not a![RLPOS] = Length(b![ROWSPOS]) then + ErrorNoReturn("\\*: Matrices do not fit together"); + fi; + if not IsIdenticalObj(a![BDPOS],b![BDPOS]) then + ErrorNoReturn("\\*: Matrices not over same base domain"); + fi; + r:=BaseDomain(a); + m:=Size(r); + l := ListWithIdenticalEntries(Length(a![ROWSPOS]),0); + for i in [1..Length(l)] do + if b![RLPOS] = 0 then + l[i] := b![EMPOS]; + else + v := a![ROWSPOS][i]; + + # do arithmetic over Z first and reduce afterwards + w:=ListWithIdenticalEntries(b![RLPOS],0); + v:=v![ELSPOS]; + for j in [1..a![RLPOS]] do + AddRowVector(w,b![ROWSPOS][j]![ELSPOS],v[j]); + #if (j mod 1000=0) and not ForAll(w,IsSmallIntRep) then + # ZNZVECREDUCE(w,b![RLPOS],m); + #fi; + od; + ZNZVECREDUCE(w,b![RLPOS],m); + w:=Vector(r,w); + + l[i] := w; + fi; + od; + if not IsMutable(a) and not IsMutable(b) then + MakeImmutable(l); + fi; + return Objectify( ty, [a![BDPOS],a![EMPOS],b![RLPOS],l] ); + end ); + +InstallMethod( \=, "for two zmodnz matrices",IsIdenticalObj, + [ IsZmodnZMatrixRep, IsZmodnZMatrixRep ], + function( a, b ) + return EQ_LIST_LIST_DEFAULT(a![ROWSPOS],b![ROWSPOS]); + end ); + +InstallMethod( \<, "for two zmodnz matrices",IsIdenticalObj, + [ IsZmodnZMatrixRep, IsZmodnZMatrixRep ], + function( a, b ) + return LT_LIST_LIST_DEFAULT(a![ROWSPOS],b![ROWSPOS]); + end ); + +InstallMethod( AdditiveInverseSameMutability, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local l; + l := List(m![ROWSPOS],AdditiveInverseSameMutability); + if not IsMutable(m) then + MakeImmutable(l); + fi; + return Objectify( TypeObj(m), [m![BDPOS],m![EMPOS],m![RLPOS],l] ); + end ); + +InstallMethod( AdditiveInverseImmutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local l,res; + l := List(m![ROWSPOS],AdditiveInverseImmutable); + res := Objectify( TypeObj(m), [m![BDPOS],m![EMPOS],m![RLPOS],l] ); + MakeImmutable(res); + return res; + end ); + +InstallMethod( AdditiveInverseMutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local l,res; + l := List(m![ROWSPOS],AdditiveInverseMutable); + res := Objectify( TypeObj(m), [m![BDPOS],m![EMPOS],m![RLPOS],l] ); + if not IsMutable(m) then + SetFilterObj(res,IsMutable); + fi; + return res; + end ); + +InstallMethod( ZeroSameMutability, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local l; + l := List(m![ROWSPOS],ZeroSameMutability); + if not IsMutable(m) then + MakeImmutable(l); + fi; + return Objectify( TypeObj(m), [m![BDPOS],m![EMPOS],m![RLPOS],l] ); + end ); + +InstallMethod( ZeroImmutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local l,res; + l := List(m![ROWSPOS],ZeroImmutable); + res := Objectify( TypeObj(m), [m![BDPOS],m![EMPOS],m![RLPOS],l] ); + MakeImmutable(res); + return res; + end ); + +InstallMethod( ZeroMutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local l,res; + l := List(m![ROWSPOS],ZeroMutable); + res := Objectify( TypeObj(m), [m![BDPOS],m![EMPOS],m![RLPOS],l] ); + if not IsMutable(m) then + SetFilterObj(res,IsMutable); + fi; + return res; + end ); + +InstallMethod( IsZero, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local i; + for i in [1..Length(m![ROWSPOS])] do + if not IsZero(m![ROWSPOS][i]) then + return false; + fi; + od; + return true; + end ); + +InstallMethod( IsOne, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local i,j,n; + if Length(m![ROWSPOS]) <> m![RLPOS] then + #Error("IsOne: Matrix must be square"); + return false; + fi; + n := m![RLPOS]; + for i in [1..n] do + if not IsOne(m![ROWSPOS][i]![ELSPOS][i]) then return false; fi; + for j in [1..i-1] do + if not IsZero(m![ROWSPOS][i]![ELSPOS][j]) then return false; fi; + od; + for j in [i+1..n] do + if not IsZero(m![ROWSPOS][i]![ELSPOS][j]) then return false; fi; + od; + od; + return true; + end ); + +InstallMethod( OneSameMutability, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local o; + if m![RLPOS] <> Length(m![ROWSPOS]) then + #Error("OneSameMutability: Matrix is not square"); + #return; + return fail; + fi; + o := IdentityMatrix(m![RLPOS],m); + if not IsMutable(m) then + MakeImmutable(o); + fi; + return o; + end ); + +InstallMethod( OneMutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + if m![RLPOS] <> Length(m![ROWSPOS]) then + #Error("OneMutable: Matrix is not square"); + #return; + return fail; + fi; + return IdentityMatrix(m![RLPOS],m); + end ); + +InstallMethod( OneImmutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local o; + if m![RLPOS] <> Length(m![ROWSPOS]) then + #Error("OneImmutable: Matrix is not square"); + #return; + return fail; + fi; + o := IdentityMatrix(m![RLPOS],m); + MakeImmutable(o); + return o; + end ); + +# For the moment we delegate to the fast kernel arithmetic for plain +# lists of plain lists: + +InstallMethod( InverseMutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local n,modulus; + modulus:=Size(BaseDomain(m)); + if m![RLPOS] <> Length(m![ROWSPOS]) then + #Error("InverseMutable: Matrix is not square"); + #return; + return fail; + fi; + # Make a plain list of lists: + n := List(m![ROWSPOS],x->x![ELSPOS]); + n := InverseMutable(n); # Invert! + if n = fail then return fail; fi; + n:=List(n,x->List(x,y->y mod modulus)); + return Matrix(n,Length(n),m); + end ); + +InstallMethod( InverseImmutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local n,modulus; + modulus:=Size(BaseDomain(m)); + if m![RLPOS] <> Length(m![ROWSPOS]) then + #Error("InverseMutable: Matrix is not square"); + #return; + return fail; + fi; + # Make a plain list of lists: + n := List(m![ROWSPOS],x->x![ELSPOS]); + n := InverseMutable(n); # Invert! + if n = fail then return fail; fi; + n:=List(n,x->List(x,y->y mod modulus)); + n := Matrix(n,Length(n),m); + MakeImmutable(n); + return n; + end ); + +InstallMethod( InverseSameMutability, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local n,modulus; + modulus:=Size(BaseDomain(m)); + if m![RLPOS] <> Length(m![ROWSPOS]) then + #Error("InverseMutable: Matrix is not square"); + #return; + return fail; + fi; + # Make a plain list of lists: + n := List(m![ROWSPOS],x->x![ELSPOS]); + n := InverseMutable(n); # Invert! + if n = fail then return fail; fi; + n:=List(n,x->List(x,y->y mod modulus)); + n := Matrix(n,Length(n),m); + if not IsMutable(m) then + MakeImmutable(n); + fi; + return n; + end ); + +InstallMethod( RankMat, "for a zmodnz matrix", [ IsZmodnZMatrixRep ], +function( m ) + m:=MutableCopyMat(m); + m:=SemiEchelonMatDestructive(m); + if m<>fail then m:=Length(m.vectors);fi; + return m; +end); + + +#InstallMethodWithRandomSource( Randomize, +# "for a random source and a mutable zmodnz matrix", +# [ IsRandomSource, IsZmodnZMatrixRep and IsMutable ], +# function( rs, m ) +# local v; +# for v in m![ROWSPOS] do +# Randomize( rs, v ); +# od; +# return m; +# end ); + +InstallMethod( TransposedMatMutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local i,n,v; + n := ListWithIdenticalEntries(m![RLPOS],0); + for i in [1..m![RLPOS]] do + v := Vector(List(m![ROWSPOS],v->v![ELSPOS][i]),m![EMPOS]); + n[i] := v; + od; + return Objectify(TypeObj(m),[m![BDPOS],m![EMPOS],Length(m![ROWSPOS]),n]); + end ); + +InstallMethod( TransposedMatImmutable, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( m ) + local n; + n := TransposedMatMutable(m); + MakeImmutable(n); + return n; + end ); + +InstallMethod( \*, "for a zmodnz vector and a zmodnz matrix", + IsElmsColls, [ IsZmodnZVectorRep, IsZmodnZMatrixRep ], + function( v, m ) + local i,res,s,r; + r:=BaseDomain(v); + # do arithmetic over Z first so that we reduce only once + res:=ListWithIdenticalEntries(m![RLPOS],0); + for i in [1..Length(v![ELSPOS])] do + s := v![ELSPOS][i]; + if not IsZero(s) then + AddRowVector(res,m![ROWSPOS][i]![ELSPOS],s); + #if (i mod 100=0) and not ForAll(res,IsSmallIntRep) then + # ZNZVECREDUCE(res,Length(res),Size(r)); + #fi; + fi; + od; + ZNZVECREDUCE(res,Length(res),Size(r)); + res:=Vector(r,res); + + if not IsMutable(v) and not IsMutable(m) then + MakeImmutable(res); + fi; + return res; + end ); + +InstallOtherMethod( \*, "for a plist vector and a zmodnz matrix", + IsElmsColls, [ IsList, IsZmodnZMatrixRep ], + function( v, m ) + local i,res,s,r; + r:=BaseDomain(m); + # do arithmetic over Z first so that we reduce only once + res:=ListWithIdenticalEntries(m![RLPOS],0); + for i in [1..Length(v)] do + s := v[i]; + if not IsZero(s) then + AddRowVector(res,m![ROWSPOS][i]![ELSPOS],Int(s)); + #if (i mod 100=0) and not ForAll(res,IsSmallIntRep) then + # ZNZVECREDUCE(res,Length(res),Size(r)); + #fi; + fi; + od; + ZNZVECREDUCE(res,Length(res),Size(r)); + res:=Vector(r,res); + + if not IsMutable(v) and not IsMutable(m) then + MakeImmutable(res); + fi; + return res; + end ); + +InstallMethod( CompatibleVector, "for a zmodnz matrix", + [ IsZmodnZMatrixRep ], + function( v ) + return NewZeroVector(IsZmodnZVectorRep,BaseDomain(v),NumberRows(v)); + end ); + +InstallMethod( DeterminantMat, "for a zmodnz matrix", [ IsZmodnZMatrixRep ], +function( a ) +local m; + m:=Size(BaseDomain(a)); + a:=List(a![ROWSPOS],x->x![ELSPOS]); + return ZmodnZObj(DeterminantMat(a),m); +end ); + + +# Minimal/Characteristic Polynomial stuff +############################################################################# +## +## Variant of +#F Matrix_OrderPolynomialInner( , , , ) +## +BindGlobal( "ZModnZMOPI",function( fld, mat, vec, vecs) + local d, w, p, one, zero, zeroes, piv, pols, x, t,i; + Info(InfoMatrix,3,"Order Polynomial Inner on ",NrRows(mat), + " x ",NrCols(mat)," matrix over ",fld," with ", + Number(vecs)," basis vectors already given"); + d := Length(vec); + pols := []; + one := One(fld); + zero := Zero(fld); + zeroes := []; + + # this loop runs images of under powers of + # trying to reduce them with smaller powers (and tracking the polynomial) + # or with vectors from as passed in + # when we succeed, we know the order polynomial + + repeat + w := ShallowCopy(vec); + p := ShallowCopy(zeroes); + Add(p,one); + #p:=ZmodnZVec(fam,p); + p:=Vector(fld,p); + piv := PositionNonZero(w,0); + + # + # Strip as far as we can + # + + while piv <= d and IsBound(vecs[piv]) do + x := -w[piv]; + if IsBound(pols[piv]) then + #AddCoeffs(p, pols[piv], x); + #p:=p+pols[piv]*x; + t:=pols[piv]*x; + for i in [1..Length(t)] do + p[i]:=p[i]+t[i]; + od; + fi; + AddRowVector(w, vecs[piv], x, piv, d); + #w:=w+vecs[piv]*x; + piv := PositionNonZero(w,piv); + od; + + # + # if something is left then we don't have the order poly yet + # update tables etc. + # + + if piv <=d then + x := Inverse(w[piv]); + MultVector(p, x); + #p:=p*x; + #MakeImmutable(p); + pols[piv] := p; + MultVector(w, x ); + #w:=w*x; + #MakeImmutable(w); + vecs[piv] := w; + vec := vec*mat; + Add(zeroes,zero); + fi; + until piv > d; + MakeImmutable(p); + Info(InfoMatrix,3,"Order Polynomial returns ",p); + return p; +end ); + +InstallOtherMethod( MinimalPolynomial, "ZModnZ, spinning over field", + IsElmsCollsX, + [ IsField and IsFinite, IsMatrixObj, IsPosInt ], +function( fld, mat, ind ) + local i, n, ords, base, vec, one,cp,zero, fam, + processVec, mp, dim, span,op,w, piv,j; + + Info(InfoMatrix,1,"Minimal Polynomial called on ", + NrRows(mat)," x ",NrCols(mat)," matrix over ",fld); + n := NrRows(mat); + base := []; + dim := 0; # should be number of bound positions in base + one := One(fld); + zero := Zero(fld); + fam := ElementsFamily(FamilyObj(fld)); + mp:=[one]; + #keep coeffs + #mp := UnivariatePolynomialByCoefficients( fam, mp,ind); + while dim < n do + vec:=ZeroVector(n,mat[1]); + for i in [1..n] do + if (not IsBound(base[i])) and Random([0,1])=1 then vec[i]:=one;fi; + od; + if IsZero(vec) then + vec[Random(1,n)] := one; # make sure it's not zero + fi; + span := []; + op := ZModnZMOPI( fld, mat, vec, span); + op:=List(op); + mp:=QUOTREM_LAURPOLS_LISTS(ProductCoeffs(mp,op),GcdCoeffs(mp,AsList(op)))[1]; + mp:=mp/mp[Length(mp)]; + Info(InfoMatrix,2,"So Far ",dim,", Span=",Length(span)); + + for j in [1..Length(span)] do + if IsBound(span[j]) then + if dim < n then + if not IsBound(base[j]) then + base[j] := span[j]; + dim := dim+1; + else + w := ShallowCopy(span[j]); + piv := j; + repeat + AddRowVector(w,base[piv],-w[piv], piv, n); + piv := PositionNonZero(w, piv); + until piv > n or not IsBound(base[piv]); + if piv <= n then + #MultVector(w,Inverse(w[piv])); + w:=w*Inverse(w[piv]); + #MakeImmutable(w); + base[piv] := w; + dim := dim+1; + fi; + fi; + fi; + fi; + od; + od; + mp := UnivariatePolynomialByCoefficients( fam, mp,ind); + Assert(3, IsZero(Value(mp,mat))); + Info(InfoMatrix,1,"Minimal Polynomial returns ", mp); + return mp; +end); + + +InstallOtherMethod( CharacteristicPolynomialMatrixNC, "zmodnz spinning over field", + IsElmsCollsX, + [ IsField, IsMatrixObj, IsPosInt ], function( fld, mat, ind) +local i, n, ords, base, imat, vec, one,cp,op,zero,fam; + Info(InfoMatrix,1,"Characteristic Polynomial called on ", + NrRows(mat)," x ",NrCols(mat)," matrix over ",fld); + imat := ImmutableMatrix(fld,mat); + n := NrRows(mat); + base := []; + vec := ZeroOp(mat[1]); + one := One(fld); + zero := Zero(fld); + fam := ElementsFamily(FamilyObj(fld)); + cp:=[one]; + cp := UnivariatePolynomialByCoefficients(fam,cp,ind); + for i in [1..n] do + if not IsBound(base[i]) then + vec[i] := one; + op := Unpack(ZModnZMOPI( fld, imat, vec, base)); + cp := cp * UnivariatePolynomialByCoefficients( fam,op,ind); + vec[i] := zero; + fi; + od; + Assert(2, Length(CoefficientsOfUnivariatePolynomial(cp)) = n+1); + if AssertionLevel()>=3 then + # cannot use Value(cp,imat), as this uses characteristic polynomial + n:=Zero(imat); + one:=One(imat); + for i in Reversed(CoefficientsOfUnivariatePolynomial(cp)) do + n:=n*imat+(i*one); + od; + Assert(3,IsZero(n)); + fi; + Info(InfoMatrix,1,"Characteristic Polynomial returns ", cp); + return cp; +end ); diff --git a/lib/read3.g b/lib/read3.g index 5e296e7201..76b4e05c0b 100644 --- a/lib/read3.g +++ b/lib/read3.g @@ -106,6 +106,7 @@ ReadLib( "wordass.gd" ); ReadLib( "matobj2.gd" ); ReadLib( "matobjplist.gd" ); +ReadLib( "matobjnz.gd" ); # files dealing with rewriting systems ReadLib( "rws.gd" ); diff --git a/lib/read5.g b/lib/read5.g index 2460261efa..07963b6d65 100644 --- a/lib/read5.g +++ b/lib/read5.g @@ -113,6 +113,7 @@ ReadLib( "vecmat.gi" ); ReadLib( "vec8bit.gi" ); ReadLib( "mat8bit.gi" ); ReadLib( "matobjplist.gi" ); +ReadLib( "matobjnz.gi" ); ReadLib( "meataxe.gi" ); ReadLib( "meatauto.gi" ); From 5eaa5a424b71ded617415758a3dc9f7a21ec9657 Mon Sep 17 00:00:00 2001 From: Alexander Hulpke Date: Wed, 16 Feb 2022 15:56:39 -0700 Subject: [PATCH 4/6] ENHANCE: Use matrix objects in Dixon-Schneider By using matrix objects (in particular for GF(p), p>65536), matrix arithmetic becomes significantly faster and can speed up calculation cost for larger groups substantially. This resolves #4716 --- lib/ctblgrp.gi | 129 +++++++++++++++++++++---------------------------- 1 file changed, 56 insertions(+), 73 deletions(-) diff --git a/lib/ctblgrp.gi b/lib/ctblgrp.gi index cc8a058ddd..e274136fe7 100644 --- a/lib/ctblgrp.gi +++ b/lib/ctblgrp.gi @@ -91,6 +91,7 @@ local D,C,cl,pl; SortParallel(cl,pl,ClassComparison); D.perm:=PermList(pl); D.permlist:=pl; + D.invpermlist:=ListPerm(D.perm^-1,Length(D.permlist)); D.currentInverseClassNo:=0; D.characterTable:=C; @@ -104,32 +105,6 @@ local D,C,cl,pl; return D; end); -# ############################################################################# -# ## -# #F DxPowerClass(,,) . . . . . . . . . . . . . evaluate powermap -# ## -# DxPowerClass := function(D,nu,power) -# local p,primes,cl; -# cl:=nu; -# power:=power mod D.characterTable.orders[cl]; -# if power=0 then -# return 1; -# elif power=1 then -# return cl; -# else -# primes:=Factors(power); -# for p in primes do -# if not IsBound(D.characterTable.powermap[p]) then -# return D.ClassElement(D, -# Representative(D.classes[nu])^power); -# else -# cl:=D.characterTable.powermap[p][cl]; -# fi; -# od; -# return cl; -# fi; -# end; - ############################################################################# ## #F DxCalcAllPowerMaps() . . . . . . . calculate power maps for char.table @@ -270,13 +245,11 @@ end; ############################################################################# ## -#F DxNiceBasis(v,space) +#F DxIsInSpace(v,space) ## DxIsInSpace := function(v,r) - if not IsBound(r.space) then - r.space:=VectorSpace(Field(r.base[1]),r.base); - fi; - return v in r.space; +local nr,nc,m,i,j; + return SolutionMat(Matrix(BaseDomain(r.base[1]),r.base),v)<>fail; end; ############################################################################# @@ -286,9 +259,9 @@ end; DxNiceBasis := function(d,r) local b; if not IsBound(r.niceBasis) then - b:=List(r.base,i->i{d.permlist}); # copied and permuted according to order + b:=Matrix(BaseDomain(r.base[1]),List(r.base,i->i{d.permlist})); # copied and permuted according to order TriangulizeMat(b); - b:=List(b,i->Permuted(i,d.perm)); # permuted back + b:=List(b,i->i{d.invpermlist}); # permuted back r.niceBasis:=Immutable(b); fi; Assert(1,Length(r.niceBasis)=Length(r.base)); @@ -336,6 +309,7 @@ DxRegisterModularChar := function(D,c) local d,p; # it may happen,that an irreducible character will be registered twice: # 2-dim space,1 Orbit,combinatoric split. Avoid this! + if IsPlistRep(c) then c:=Vector(c); fi; if not(c in D.modulars) then Add(D.modulars,c); d:=Int(c[1]); @@ -412,25 +386,28 @@ local i,pcomp,m,r,D,neue,tm,news,opt; od; if opt<>true then - pcomp:=NullspaceMat(D.projectionMat*TransposedMat(D.modulars)); - for i in [1..Length(D.raeume)] do - r:=D.raeume[i]; - if r.dim = Length(r.base[1]) then - # trivial case: Intersection with full space in the beginning - r:=rec(base:=pcomp); - else - r:=rec(base:=SumIntersectionMat(pcomp,r.base)[2]); - fi; - r.dim:=Length(r.base); - # note stabilizer - if IsBound(D.raeume[i].stabilizer) then - r.stabilizer:=D.raeume[i].stabilizer; - fi; - if r.dim>0 then - DxActiveCols(D,r); - fi; - D.raeume[i]:=r; - od; + pcomp:=NullspaceMat(D.projectionMat*TransposedMatImmutable(Matrix(D.field,D.modulars))); + if Length(pcomp)>0 then + for i in [1..Length(D.raeume)] do + r:=D.raeume[i]; + if r.dim = Length(r.base[1]) then + # trivial case: Intersection with full space in the beginning + r:=rec(base:=pcomp); + else + r:=rec(base:=SumIntersectionMat(pcomp, + Matrix(BaseDomain(r.base[1]),r.base))[2]); + fi; + r.dim:=Length(r.base); + # note stabilizer + if IsBound(D.raeume[i].stabilizer) then + r.stabilizer:=D.raeume[i].stabilizer; + fi; + if r.dim>0 then + DxActiveCols(D,r); + fi; + D.raeume[i]:=r; + od; + fi; fi; D.raeume:=Filtered(D.raeume,i->i.dim>0); end ); @@ -627,7 +604,7 @@ local ml,nu,ret,r,p,v,alo,ofs,orb,i,j,inv,b; alo:=Length(b); od; od; - inv:=b^(-1); + inv:=Matrix(BaseDomain(b[1]),b)^(-1); for i in ml do v:=i*inv; for r in [1..Length(D.raeume)] do @@ -636,7 +613,7 @@ local ml,nu,ret,r,p,v,alo,ofs,orb,i,j,inv,b; p[j]:=v[j]; od; p:=p*b; - if p<>nu then + if not IsZero(p) then AddSet(ret,DxLiftCharacter(D,p)); fi; od; @@ -653,7 +630,7 @@ DxEigenbase := function(M,f) local dim,i,k,eigenvalues,base,minpol,bases; k:=Length(M); - minpol:=MinimalPolynomial(f,M); + minpol:=MinimalPolynomial(BaseDomain(M),M); Assert(2,IsDuplicateFree(RootsOfUPol(minpol))); eigenvalues:=Set(RootsOfUPol(minpol)); @@ -664,7 +641,6 @@ DxEigenbase := function(M,f) if base=[] then Error("This can`t happen: Wrong Eigenvalue ???"); else - #TriangulizeMat(base); dim:=dim+Length(base); Add(bases,base); fi; @@ -705,7 +681,7 @@ InstallGlobalFunction(SplitStep,function(D,bestMat) if ForAny(raeume,i->i.dim>1) then bestMatCol:=D.requiredCols[bestMat]; bestMatSplit:=D.splitBases[bestMat]; - M:= NullMat( k, k, 0 ); + M:= ZeroMatrix(Integers,k,k); Info(InfoCharacterTable,1,"Matrix ",bestMat,",Representative of Order ", Order(D.classreps[bestMat]), ",Centralizer: ",D.centralizers[bestMat]); @@ -729,15 +705,21 @@ InstallGlobalFunction(SplitStep,function(D,bestMat) dim:=raum.dim; # cut out the 'active part' for computation of an eigenbase activeCols:=DxActiveCols(D,raum); - N:= NullMat( dim, dim, o ); + N:= ZeroMatrix(f,dim,dim); for row in [1..dim] do - Row:=base[row]*M; + #Row:=base[row]*M; + Row:=CompatibleVector(M); + for col in [1..Length(Row)] do + Row[col]:=base[row,col]; + od; + Row:=Row*M; for col in [1..dim] do - N[row][col]:=Row[activeCols[col]]; + N[row,col]:=Row[activeCols[col]]; od; od; eigen:=DxEigenbase(N,f); # Base umrechnen + base:=Matrix(BaseDomain(base[1]),base); eigenbase:=List(eigen.base,i->List(i,j->j*base)); #eigenvalues:=List(eigen.values,i->i/D.classiz[bestMat]); @@ -1092,7 +1074,7 @@ local i,newRaeume,raum,neuer,j,ch,irrs,mods,incirrs,incmods,nb,rt,neuc; Info(InfoCharacterTable,1,"TwoDimSpace image"); nb:=D.asCharacterMorphism(raum.base[1],j); neuc:=List(irrs,i->D.asCharacterMorphism(i,j)); - if not ForAny([i+1..Length(D.raeume)],i->nb in D.raeume[i].space) then + if not ForAny([i+1..Length(D.raeume)],i->DxIsInSpace(nb,D.raeume[i])) then incirrs:=Union(incirrs,neuc); incmods:=Union(incmods,List(mods,i->D.asCharacterMorphism(i,j))); else @@ -1747,7 +1729,7 @@ end; StandardClassMatrixColumn := function(D,M,r,t) local c,gt,s,z,i,T,w,e,j,p,orb; if t=1 then - M[D.inversemap[r]][t]:=D.classiz[r]; + M[D.inversemap[r],t]:=D.classiz[r]; else orb:=DxGaloisOrbits(D,r); z:=D.classreps[t]; @@ -1756,9 +1738,9 @@ StandardClassMatrixColumn := function(D,M,r,t) p:=RepresentativeAction(Stabilizer(orb.group,r),c,t); if p<>fail then # was the first column of the galois class active? - if ForAny(M,i->i[c]>0) then + if ForAny([1..NrRows(M)],i->M[i,c]>0) then for i in D.classrange do - M[i^p][t]:=M[i][c]; + M[i^p,t]:=M[i,c]; od; Info(InfoCharacterTable,2,"by GaloisImage"); return; @@ -1787,19 +1769,19 @@ StandardClassMatrixColumn := function(D,M,r,t) else # only strong test possible s:=D.ClassElement(D,e); fi; - M[s][t]:=M[s][t]+T[2][i]; + M[s,t]:=M[s,t]+T[2][i]; od; if w then # weak discrimination possible ? gt:=Set(Filtered(orb.orbits,i->Length(i)>1)); for i in gt do if i[1] in orb.identifees then # were these classes detected weakly ? - e:=M[i[1]][t]; + e:=M[i[1],t]; if e>0 then Info(InfoCharacterTable,2,"GaloisIdentification ",i,": ",e); fi; for j in i do - M[j][t]:=e/Length(i); + M[j,t]:=e/Length(i); od; fi; od; @@ -1808,7 +1790,7 @@ StandardClassMatrixColumn := function(D,M,r,t) for i in [1..Length(T[1])] do s:=D.ClassElement(D,T[1][i] * z); Unbind(T[1][i]); - M[s][t]:=M[s][t]+T[2][i]; + M[s,t]:=M[s,t]+T[2][i]; od; fi; fi; @@ -1997,7 +1979,7 @@ local G, # group D.field:=f; D.one:=One(f); D.z:=z; - r:=rec(base:=Immutable( IdentityMat(k,D.one) ),dim:=k); + r:=rec(base:=List(IdentityMat(k,D.one),Vector),dim:=k); D.raeume:=[r]; # Galois group operating on the columns @@ -2011,11 +1993,12 @@ local G, # group D.galOp:=[]; D.irreducibles:=[]; - M:= IdentityMat(k); + fk:=D.one/(Size(G) mod prime); + M:= ZeroMatrix(D.field,k,k); for i in [1..k] do - M[i][i]:=D.classiz[i] mod prime; + M[i,i]:=(D.classiz[i] mod prime)*D.one*fk; od; - D.projectionMat:=M*(D.one/(Size(G) mod prime)); + D.projectionMat:=M; #if (USECTPGROUP or Size(G)<2000 or k*10>=Size(G)) # and IsBound(G.isAgGroup) and G.isAgGroup @@ -2043,7 +2026,7 @@ local G, # group # CharacterMorphisms. D.raeume[1].stabilizer:=CharacterMorphismGroup(D); m:=First(D.classes,i->Size(i)>1); - if Size(m)>8 then + if m<>fail and Size(m)>8 then D.maycent:=true; fi; fi; From eec16caa7f200c4c43a0f73c10e3946bd9358fbe Mon Sep 17 00:00:00 2001 From: Alexander Hulpke Date: Wed, 16 Feb 2022 16:27:56 -0700 Subject: [PATCH 5/6] Test and changes in test output and examples --- doc/ref/ctbl.xml | 33 ------------------------ tst/testinstall/MatrixObj/ZeroVector.tst | 4 +-- tst/teststandard/ctblfuns.tst | 3 +++ 3 files changed, 5 insertions(+), 35 deletions(-) diff --git a/doc/ref/ctbl.xml b/doc/ref/ctbl.xml index ff51a29929..3dcf1c4a1c 100644 --- a/doc/ref/ctbl.xml +++ b/doc/ref/ctbl.xml @@ -515,39 +515,6 @@ gap> Length( ro.irreducibles ); gap> DxIncludeIrreducibles( d, ro.irreducibles ); ]]>

-The tensor products of the nonlinear characters among each other are reduced -with the irreducible characters. -The result is split according to the spaces found, which yields characters -of smaller norms, but no new irreducibles. -

- nlc:= Filtered( d.irreducibles, i -> i[1] > 1 );; -gap> t:= Tensored( nlc, nlc );; -gap> ro:= ReducedCharacters( c, d.irreducibles, t );; ro.irreducibles; -[ ] -gap> List( ro.remainders, i -> ScalarProduct( c, i, i) ); -[ 2, 2, 4, 4, 4, 4, 13, 13, 18, 18, 19, 21, 21, 36, 36, 29, 34, 34, - 42, 34, 48, 54, 62, 68, 68, 78, 84, 84, 88, 90, 159, 169, 169, 172, - 172, 266, 271, 271, 268, 274, 274, 280, 328, 373, 373, 456, 532, - 576, 679, 683, 683, 754, 768, 768, 890, 912, 962, 1453, 1453, 1601, - 1601, 1728, 1739, 1739, 1802, 2058, 2379, 2414, 2543, 2744, 2744, - 2920, 3078, 3078, 4275, 4275, 4494, 4760, 5112, 5115, 5115, 5414, - 6080, 6318, 7100, 7369, 7369, 7798, 8644, 10392, 12373, 12922, - 14122, 14122, 18948, 21886, 24641, 24641, 25056, 38942, 44950, - 78778 ] -gap> t:= SplitCharacters( d, ro.remainders );; -gap> List( t, i -> ScalarProduct( c, i, i ) ); -[ 2, 2, 4, 2, 2, 4, 4, 3, 6, 5, 5, 9, 9, 4, 12, 13, 18, 18, 20, 18, - 20, 24, 26, 32, 32, 16, 42, 59, 69, 69, 72, 72, 36, 72, 78, 78, 84, - 122, 117, 127, 117, 127, 64, 132, 100, 144, 196, 256, 456, 532, - 576, 679, 683, 683, 754, 768, 768, 890, 912, 962, 1453, 1453, 1601, - 1601, 1728, 1739, 1739, 1802, 2058, 2379, 2414, 2543, 2744, 2744, - 2920, 3078, 3078, 4275, 4275, 4494, 4760, 5112, 5115, 5115, 5414, - 6080, 6318, 7100, 7369, 7369, 7798, 8644, 10392, 12373, 12922, - 14122, 14122, 18948, 21886, 24641, 24641, 25056, 38942, 44950, - 78778 ] -]]> -

Finally we calculate the characters induced from all cyclic subgroups and obtain the missing irreducibles by applying the LLL-algorithm to them.

diff --git a/tst/testinstall/MatrixObj/ZeroVector.tst b/tst/testinstall/MatrixObj/ZeroVector.tst index f50f424ead..de4aa2b3d1 100644 --- a/tst/testinstall/MatrixObj/ZeroVector.tst +++ b/tst/testinstall/MatrixObj/ZeroVector.tst @@ -82,9 +82,9 @@ Error, ZeroVector: length must be non-negative # gap> ZeroVector(Integers mod 4, 2); - + gap> ZeroVector(Integers mod 4, 0); - + gap> ZeroVector(Integers mod 4, -1); Error, ZeroVector: length must be non-negative diff --git a/tst/teststandard/ctblfuns.tst b/tst/teststandard/ctblfuns.tst index a1121d32d8..8b94ed3d96 100644 --- a/tst/teststandard/ctblfuns.tst +++ b/tst/teststandard/ctblfuns.tst @@ -52,5 +52,8 @@ gap> FrobeniusCharacterValue( 82*E(16)+E(16)^5, 269 ); gap> FrobeniusCharacterValue( E(16), 269 ); 162+256z+143z2+219z3 +# Dixon-Schneider test that also exercises MatrixObjects over Z/nZ +gap> Irr(MathieuGroup(24));; + # gap> STOP_TEST( "ctblfuns.tst" ); From e6898446bc6a89e48c9ed37c3be18cb637ac005e Mon Sep 17 00:00:00 2001 From: Alexander Hulpke Date: Wed, 16 Feb 2022 19:07:59 -0700 Subject: [PATCH 6/6] HACK: Rank method higher than in package The semigroups package has a method for `Matrix` that causes problems. --- lib/matobj.gi | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lib/matobj.gi b/lib/matobj.gi index 3dd194e6eb..2bc5fd5987 100644 --- a/lib/matobj.gi +++ b/lib/matobj.gi @@ -297,6 +297,10 @@ InstallMethod( Matrix, InstallMethod( Matrix, [ IsSemiring, IsList ], +# rank higher than a method in the semigroups package, which otherwise jumps +# in and causes an error when testing +# line 318 of semigroups-3.4.0/gap/elements/semiringmat.gi +20, function( R, list ) if Length(list) = 0 then Error( "list must be not empty" );