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/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; diff --git a/lib/matobj.gi b/lib/matobj.gi index 0dc81b0b20..2bc5fd5987 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); @@ -293,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" ); @@ -1345,6 +1353,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/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; 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 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" ); 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" );