[GAP Forum] Bug in modular determinant
Alexander Konovalov
alexander.konovalov at gmail.com
Tue Oct 28 21:57:51 GMT 2008
Dear Ángel,
Thank you for reporting this. This is already fixed and will appear in
the next release GAP 4.5.
In the meantime you may use the following workaround, extracted from
the GAP development version:
InstallMethod( DeterminantMatDestructive,"nonprime residue rings",
[ IsOrdinaryMatrix and
CategoryCollections(CategoryCollections(IsZmodnZObjNonprime)) and
IsMutable],
DeterminantMatDivFree);
#############################################################################
##
#M DeterminantMatDivFree( <M> )
##
## Division free method. This is an alternative to the fraction free
method
## when division of matrix entries is expensive or not possible.
##
## This method implements a division free algorithm found in
## Mahajan and Vinay \cite{MV97}.
##
## The run time is $O(n^4)$
## Auxillary storage size $n^2+n + C$
##
## Our implementation has two runtime optimizations (both noted
## by Mahajan and Vinay)
## 1. Partial monomial sums, subtractions, and products are done at
## each level.
## 2. Prefix property is maintained allowing for a pruning of many
## vertices at each level
##
## and two auxillary storage size optimizations
## 1. only the upper triangular and diagonal portion of the
## auxillary storage is used.
## 2. Level information storage is reused (2 levels).
##
## This code was implemented by:
## Timothy DeBaillie
## Robert Morse
## Marcus Wassmer
##
InstallMethod( DeterminantMatDivFree,
"Division-free method",
[ IsMatrix ],
function ( M )
local u,v,w,i, ## indices
a,b,c,x,y, ## temp indices
temp, ## temp variable
nlevel, ## next level
clevel, ## current level
pmone, ## plus or minus one
zero, ## zero of the ring
n, ## size of the matrix
Vs, ## final sum
V; ## graph
# check that the argument is a square matrix and set the size
n := Length(M);
if not n = Length(M[1]) or not IsRectangularTable(M) then
Error("DeterminantMatDivFree: <mat> must be a square
matrix");
fi;
## initialze the final sum, the vertex set, initial parity
## and level indexes
##
zero := Zero(M[1][1]);
Vs := zero;
V := [];
pmone := (-One(M[1][1]))^((n mod 2)+1);
clevel := 1; nlevel := 2;
## Vertices are indexed [u,v,i] holding the (partial) monomials
## whose sums will form the determinant
## where i = depth in the tree (current and next reusing
## level storage)
## u,v indices in the matrix
##
## Only the upper triangular portion of the storage space is
## needed. It is easier to create lower triangular data type
## which we do here and index via index arithmetic.
##
for u in [1..n] do
Add(V,[]);
for v in [1..u] do
Add(V[u],[zero,zero]);
od;
## Initialize the level 0 nodes with +/- one, depending on
## the initial parity determined by the size of the matrix
##
V[u][u][clevel] := pmone;
od;
## Here are the $O(n^4)$ edges labeled by the elements of
## the matrix $M$. We build up products of the labels which form
## the monomials which make up the determinant.
##
## 1. Parity of monomials are maintained implicitly.
## 2. Partial sums for some vertices are not part of the final
## answer and can be pruned.
##
for i in [0..n-2] do
for u in [1..i+2] do ## <---- pruning of vertices
for v in [u..n] do ## (maintains the prefix
property)
for w in [u+1..n] do
## translate indices to lower triangluar
coordinates
##
a := n-u+1; b := n-w+1; c := n-v+1;
V[a][b][nlevel]:= V[a][b][nlevel]+
V[a][c][clevel]*M[v][w];
V[b][b][nlevel]:= V[b][b][nlevel]-
V[a][c][clevel]*M[v][u];
od;
od;
od;
## set the new current and next level. The new next level
## is intialized to zero
##
temp := nlevel; nlevel := clevel; clevel := temp;
for x in [1..n] do
for y in [1..x] do
V[x][y][nlevel] := zero;
od;
od;
od;
## with the final level, we form the last monomial product and
then
## sum these monomials (parity has been accounted for)
## to find the determinant.
##
for u in [1..n] do
for v in [u..n] do
Vs := Vs + V[n-u+1][n-v+1][clevel]*M[v][u];
od;
od;
## Return the final sum
##
return Vs;
end);
Then after reading this code into GAP your example works:
gap> a := ZmodnZObj(1,27)*[[21,4,11,1],[0,25,11,1],[0,2,15,1],
[13,19,4,1]];
[ [ ZmodnZObj( 21, 27 ), ZmodnZObj( 4, 27 ), ZmodnZObj( 11, 27 ),
ZmodnZObj( 1, 27 ) ],
[ ZmodnZObj( 0, 27 ), ZmodnZObj( 25, 27 ), ZmodnZObj( 11, 27 ),
ZmodnZObj( 1, 27 ) ],
[ ZmodnZObj( 0, 27 ), ZmodnZObj( 2, 27 ), ZmodnZObj( 15, 27 ),
ZmodnZObj( 1, 27 ) ],
[ ZmodnZObj( 13, 27 ), ZmodnZObj( 19, 27 ), ZmodnZObj( 4, 27 ),
ZmodnZObj( 1, 27 ) ] ]
gap> DeterminantMat(a);
ZmodnZObj( 12, 27 )
Best wishes,
Alexander
On 23 Oct 2008, at 16:27, Decanato Fac. Matemáticas wrote:
> Dear GAP Forum,
>
> I encountered the following problem when I tried to calculate the
> determinant of a matrix of integers modulo 27:
>
> gap> a := ZmodnZObj(1,27)*[[21,4,11,1],[0,25,11,1],[0,2,15,1],
> [13,19,4,1]];
> [ [ ZmodnZObj( 21, 27 ), ZmodnZObj( 4, 27 ), ZmodnZObj( 11, 27 ),
> ZmodnZObj( 1, 27 ) ],
> [ ZmodnZObj( 0, 27 ), ZmodnZObj( 25, 27 ), ZmodnZObj( 11, 27 ),
> ZmodnZObj( 1, 27 ) ],
> [ ZmodnZObj( 0, 27 ), ZmodnZObj( 2, 27 ), ZmodnZObj( 15, 27 ),
> ZmodnZObj( 1, 27 ) ],
> [ ZmodnZObj( 13, 27 ), ZmodnZObj( 19, 27 ), ZmodnZObj( 4, 27 ),
> ZmodnZObj( 1, 27 ) ] ]
> gap> Determinant(a);
> Error, no method found! For debugging hints type ?Recovery from
> NoMethodFound
> Error, no 2nd choice method found for `MultRowVector' on 2 arguments
> called fr\
> om
> MultRowVector( row2, Inverse( det ) ); called from
> DeterminantMatDestructive( MutableCopyMat( mat ) ) called from
> <function>( <arguments> ) called from read-eval-loop
> Entering break read-eval-print loop ...
> you can 'quit;' to quit to outer loop, or
> you can 'return;' to continue
> brk>
>
> Ángel
More information about the Forum
mailing list