Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
| Download
GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Views: 418346#############################################################################
##
## GaussDense.gi Gauss package Simon Goertzen
##
## Copyright 2007-2008 Lehrstuhl B für Mathematik, RWTH Aachen
##
## Implementation stuff for Gauss algorithms on dense (IsMatrix) matrices.
##
#############################################################################
##
InstallMethod( EchelonMatTransformationDestructive,
"generic method for matrices",
[ IsMatrix and IsMutable ],
function( mat )
local zero, # zero of the ring of <mat>
nrows, # number of rows in <mat>
ncols, # number of columns in <mat>
vectors, # list of basis vectors
field,
heads, # list of pivot positions in 'vectors'
i, # loop over rows
j, # loop over columns
T, # transformation matrix
coeffs, # list of coefficient vectors for 'vectors'
relations, # basis vectors of the null space of 'mat'
row,
head,
x,
row2,
rank,
list,
a;
nrows := Length( mat );
ncols := Length( mat[1] );
zero := Zero( mat[1][1] );
heads := ListWithIdenticalEntries( ncols, 0 );
vectors := [];
if Characteristic( zero ) mod 2 = 0 then
if IsGF2VectorRep( mat[1] ) then
field := GF( 2 );
else
field := GF( 2^Q_VEC8BIT( mat[1] ) );
fi;
else
field := Field( zero );
fi;
T := IdentityMat( nrows, field );
coeffs := [];
relations := [];
for i in [ 1 .. nrows ] do
row := mat[i];
row2 := T[i];
# Reduce the row with the known basis vectors.
for j in [ 1 .. ncols ] do
head := heads[j];
if head <> 0 then
x := - row[j];
if x <> zero then
AddRowVector( row2, coeffs[ head ], x );
AddRowVector( row, vectors[ head ], x );
fi;
fi;
od;
j:= PositionNot( row, zero );
if j <= ncols then
# We found a new basis vector.
x := Inverse( row[j] );
if x = fail then
TryNextMethod();
fi;
Add( coeffs, row2 * x );
Add( vectors, row * x );
heads[j]:= Length( vectors );
else
Add( relations, row2 );
fi;
od;
# gauss upwards:
list := Filtered( heads, x->x<>0 );
rank := Length( list );
for j in [ncols,ncols-1..1] do
head := heads[j];
if head <> 0 then
a := Difference( [1..head-1], heads{[j+1..ncols]} );
for i in a do
row := vectors[i];
row2 := coeffs[i];
x := - row[j];
if x <> zero then
AddRowVector( row2, coeffs[head], x );
AddRowVector( row, vectors[head], x );
fi;
od;
fi;
od;
#order rows:
vectors := vectors{list};
coeffs := coeffs{list};
list := Filtered( [1..ncols], j -> heads[j] <> 0 );
heads{list} := [1..rank]; #just for compatibilty, vectors are ordered already
return rec( heads := heads,
vectors := vectors,
coeffs := coeffs,
relations := relations );
end );
##
InstallMethod( EchelonMatTransformation,
"generic method for matrices",
[ IsMatrix ],
function( mat )
local copymat, v, vc, f;
copymat := [];
f := DefaultFieldOfMatrix(mat);
for v in mat do
vc := ShallowCopy(v);
ConvertToVectorRepNC(vc,f);
Add(copymat, vc);
od;
return EchelonMatTransformationDestructive( copymat );
end);
##
InstallMethod( EchelonMatDestructive,
"generic method for matrices",
[ IsMatrix and IsMutable ],
function( mat )
local zero, # zero of the ring of <mat>
nrows, # number of rows in <mat>
ncols, # number of columns in <mat>
vectors, # list of basis vectors
heads, # list of pivot positions in 'vectors'
i, # loop over rows
j, # loop over columns
row,
head,
x,
row2,
rank,
list,
a;
nrows := Length( mat );
ncols := Length( mat[1] );
zero := Zero( mat[1][1] );
heads := ListWithIdenticalEntries( ncols, 0 );
vectors := [];
for i in [ 1 .. nrows ] do
row := mat[i];
# Reduce the row with the known basis vectors.
for j in [ 1 .. ncols ] do
head := heads[j];
if head <> 0 then
x := - row[j];
if x <> zero then
AddRowVector( row, vectors[ head ], x );
fi;
fi;
od;
j:= PositionNot( row, zero );
if j <= ncols then
# We found a new basis vector.
x := Inverse( row[j] );
if x = fail then
TryNextMethod();
fi;
Add( vectors, row * x );
heads[j]:= Length( vectors );
fi;
od;
# gauss upwards:
list := Filtered( heads, x->x<>0 );
rank := Length( list );
for j in [ncols,ncols-1..1] do
head := heads[j];
if head <> 0 then
a := Difference( [1..head-1], heads{[j+1..ncols]} );
for i in a do
row := vectors[i];
x := - row[j];
if x <> zero then
AddRowVector( row, vectors[head], x );
fi;
od;
fi;
od;
#order rows:
vectors := vectors{list};
#ConvertToMatrixRepNC( vectors ); FIXME: Is this important or neccessary?
list := Filtered( [1..ncols], j -> heads[j] <> 0 );
heads{list} := [1..rank]; #just for compatibility, vectors are ordered already
return rec( heads := heads,
vectors := vectors );
end );
##
InstallMethod( EchelonMat,
"generic method for matrices",
[ IsMatrix ],
function( mat )
local copymat, v, vc, f;
copymat := [];
f := DefaultFieldOfMatrix(mat);
for v in mat do
vc := ShallowCopy(v);
ConvertToVectorRepNC(vc,f);
Add(copymat, vc);
od;
return EchelonMatDestructive( copymat );
end);
##
InstallMethod( ReduceMatWithEchelonMat,
"for general matrices over a ring, second argument must be in REF",
[ IsMatrix, IsMatrix ],
function( mat, N )
local nrows1,
ncols,
nrows2,
M,
f,
v,
vc,
zero,
i,
row2,
j,
k,
row1,
x;
nrows1 := Length( mat );
nrows2 := Length( N );
ncols := Length( mat[1] );
if ncols <> Length( N[1] ) then
Error( "column dimensions do not match! aborting ReduceMat!" );
fi;
M := [];
f := DefaultFieldOfMatrix( mat );
for v in mat do
vc := ShallowCopy( v );
ConvertToVectorRepNC( vc, f );
Add( M, vc );
od;
zero := Zero( M[1][1] );
for i in [1 .. nrows2] do
row2 := N[i];
j := PositionNot( row2, zero );
if j <= ncols then
for k in [1 .. nrows1] do
row1 := M[k];
x := - row1[j];
if x <> zero then
AddRowVector( row1, row2, x );
fi;
od;
fi;
od;
return rec( reduced_matrix := M );
end );
##
InstallMethod( ReduceMat,
[ IsMatrix, IsMatrix ],
function( mat, N )
return ReduceMatWithEchelonMat( mat, N );
end
);
##
InstallMethod( ReduceMatWithEchelonMatTransformation,
"for general matrices over a ring, second argument must be in REF",
[ IsMatrix, IsMatrix ],
function( mat, N )
local nrows1,
ncols,
nrows2,
M,
f,
v,
vc,
T,
zero,
i,
row2,
j,
k,
row1,
x;
nrows1 := Length( mat );
nrows2 := Length( N );
ncols := Length( mat[1] );
if ncols <> Length( N[1] ) then
Error( "column dimensions do not match! aborting ReduceMatTransformation!" );
fi;
M := [];
f := DefaultFieldOfMatrix( mat );
for v in mat do
vc := ShallowCopy( v );
ConvertToVectorRepNC( vc, f );
Add( M, vc );
od;
T := NullMat( nrows1, nrows2, f );
zero := Zero( M[1][1] );
for i in [1 .. nrows2] do
row2 := N[i];
j := PositionNot( row2, zero );
if j <= ncols then
for k in [1 .. nrows1] do
row1 := M[k];
x := - row1[j];
if x <> zero then
AddRowVector( row1, row2, x );
T[k][i] := x;
fi;
od;
fi;
od;
return rec( reduced_matrix := M, transformation := T);
end
);
##
InstallMethod( ReduceMatTransformation,
[ IsMatrix, IsMatrix ],
function( mat, N )
return ReduceMatWithEchelonMatTransformation( mat, N );
end
);
##
InstallMethod( KernelEchelonMatDestructive,
"generic method for matrices",
[ IsMatrix and IsMutable, IsList ],
function( mat, L )
local zero, # zero of the ring of <mat>
nrows, # number of rows in <mat>
ncols, # number of columns in <mat>
vectors, # list of basis vectors
heads, # list of pivot positions in 'vectors'
i, # loop over rows
j, # loop over columns
T, # transformation matrix
coeffs, # list of coefficient vectors for 'vectors'
relations, # basis vectors of the null space of 'mat'
row,
head,
x,
row2;
nrows := Length( mat );
ncols := Length( mat[1] );
zero := Zero( mat[1][1] );
heads := ListWithIdenticalEntries( ncols, 0 );
vectors := [];
T := IdentityMat( nrows, zero ){[1..nrows]}{L};
coeffs := [];
relations := [];
for i in [ 1 .. nrows ] do
row := mat[i];
row2 := T[i];
# Reduce the row with the known basis vectors.
for j in [ 1 .. ncols ] do
head := heads[j];
if head <> 0 then
x := - row[j];
if x <> zero then
AddRowVector( row2, coeffs[ head ], x );
AddRowVector( row, vectors[ head ], x );
fi;
fi;
od;
j:= PositionNot( row, zero );
if j <= ncols then
# We found a new basis vector.
x := Inverse( row[j] );
if x = fail then
TryNextMethod();
fi;
Add( coeffs, row2 * x );
Add( vectors, row * x );
heads[j]:= Length( vectors );
else
Add( relations, row2 );
fi;
od;
return rec( relations := relations );
end );
##
InstallGlobalFunction( KernelMat,
function( arg )
local copymat,
f,
v,
vc;
if IsSparseMatrix( arg[1] ) then
return CallFuncList( KernelMatSparse, arg );
fi;
copymat := [];
f := DefaultFieldOfMatrix( arg[1] );
if f = fail then
Error( "No KernelMat for dense matrices over non-fields!" );
fi;
for v in arg[1] do
vc := ShallowCopy(v);
ConvertToVectorRepNC(vc,f);
Add(copymat, vc);
od;
if Length( arg ) = 1 then
return KernelEchelonMatDestructive( copymat, [1..Length( arg[1] )] );
elif Length( arg ) > 1 then
return KernelEchelonMatDestructive( copymat, arg[2] );
fi;
end );