This file contains a bugfix for a performance problem in GAP 3.4.3.
You should apply this bugfix if you are computing with integer matrices.
The problem is that the function 'DiagonalizeMatrix', which is used by
'ElementaryDivisorsMat' and by all the function 'AbelianInvariants...',
is not careful enough to keep the entries small.
HOW TO APPLY
The problem is a performance problem, because it may cause a computation
to take much longer that it should. Thus the bugfix has low priority.Go to the GAP directory (the directory with the 'lib/' subdirectory),
name this file 'bugfix08.dif', and issue the command:patch -p0 < bugfix08.difThis workaround changes only the library. Thus you need not recompile
the GAP kernel.
VERSION
GAP 3.4.3.0
DESCRIPTION
'DiagonalizeMatrix', which is used by 'ElementaryDivisorsMat' and by all
the functions 'AbelianInvariants...', is not careful enough to keep the
entries small. Thus the size of intermediate entries may grow
exponentially.
CORRECT BEHAVIOUR
gap> M := [ > [ 63851, -4582,-14295, 36409, 31705,-15403,-15221,-26632, > -1419,-35034,-33355, 14971,-27398, -64, 66589 ], > [ 12285,-51674, -3928, 21474,-57991, 45051, 14076,-74043, > 56389, 38873,-61766,-32221, 50255, 4823,-28850 ], > [ 64,-48464,-29455,-13125,-44592,-10902,-10389,-13507, > 19844, 33041, 36473, 42632,-52949, 7623,-42507 ], > [ -4823, 42496,-15407,-45251,-52043, 10367, 83895,-28792, > -20969, 48444, 16651,-27093,-22652,-25856, -3285 ], > [ -7623,-46432, -8892, -3614, 5233,-25527, 5832, -9893, > -21233, 3393,-20748, 1381, 25463,-24033, 21948 ], > [ 25856,-12008, 14231,-16133,-25740, 14048,-15403,-12659, > -30926, 14295,-19395, 22910, -9225, 2553,-51311 ], > [ 24033,-32034, 8751,-55273, 41195, 6515, 45051, 45380, > -29069, 3928, 20419, 21667, 1108,-30892, 14177 ], > [ -2553, 19910, 37078,-39838,-16195, 23123,-10902, 27179, > 22077, 29455,-18686, 5617,-48979, -4509, 30080 ], > [ 30892, 40630,-10449, 31051, 15972, -5480, 10367, 14329, > -27746, 15407,-24381,-54678,-18553, 5699, 21319 ], > [ 4509, -1972,-15141, -1309, -5331, 28327,-25527, 28488, > 12373, 8892,-35329, -4543,-25774,-14010,-32265 ], > [ -5699, -4262,-11678, 1252, 6003,-47527, 14048, 13077, > 3295,-14231, 25352,-18861, 16635, -6913, 1566 ], > [ 14010, 3112,-39643, -655,-17194, -4692, 6515, 29143, > -5442, -8751, 12701, 2368, 183, 14267, 12667 ], > [ 6913, -7678,-41587, 23471,-23397, 3463, 23123,-10394, > 27929,-37078, 8165,-20649, 2236, 20832,-16395 ], > [-14267, 13974, 16148, 32402,-20719,-27965, -5480, -3259, > -10137, 10449,-14922, 33553, 13955,-31317, 1456 ], > [-20832, -5168, 1131,-27917, -3088, -1944, 28327, 17523, > 10660, 15141, 2639, 9860, -4969, 4095,-11097 ] ];; gap> InfoMatrix1 := Print;; InfoMatrix2 := Print;; InfoMatrix3 := Print;; gap> Collected( ElementaryDivisorsMat( M ) ); [ [ 1, 1 ], [ 78901, 1 ], [ 157802, 13 ] ] gap> # this should not produce intermediate entries with more than 10 digits gap> # the old code produced intermediate entries with more than 100 digits
COMMENT
This patch consists of a new implementation of 'DiagonalizeMatrix'.
It tries to keep the entries small through careful selection of pivots.
The selected pivot is the entry for which the product of row and column
norm is minimal (this need not be the entry with minimal absolute value).
It is based upon ideas by George Havas and code by Bohdan Majewski.
PATCH
Prereq: 3.33.1.1
--- lib/matrix.g Mon May 6 14:10:17 1996
+++ lib/matrix.g Mon May 6 14:38:57 1996
@@ -2,13 +2,16 @@
##
#A matrix.g GAP library Martin Schoenert
##
-#A @(#)$Id: 1.html,v 1.2 2004/04/21 15:06:18 felsch Exp $
+#A @(#)$Id: 1.html,v 1.2 2004/04/21 15:06:18 felsch Exp $
##
#Y Copyright 1990-1992, Lehrstuhl D fuer Mathematik, RWTH Aachen, Germany
##
## This file contains those functions that mainly deal with matrices.
##
#H $Log: 1.html,v $
#H Revision 1.2 2004/04/21 15:06:18 felsch
#H Corrected links in the Forum Archive pages. VF
#H
#H Revision 1.1.1.1 2004/04/20 13:39:30 felsch
#H The final GAP-Forum archive until 2003.
#H
#H Revision 1.5 2003/06/12 19:20:33 gap
#H Further update. AH
#H
#H Revision 1.4 2003/06/12 17:28:25 gap
#H Address updates by JN. AH
#H
#H Revision 1.3 1997/08/15 11:19:33 gap
#H New forum setup. AH
#H
#H Revision 1.2 1997/04/24 15:32:47 gap
#H These files were replaced by the versions in WWW. The content is basically the
#H same but the formatting has been much more friendly towards the HTML-Converter.
#H AH
#H
#H Revision 1.1 1996/10/30 13:07:04 gap
#H added forum archive and translation files.
#H
+#H Revision 3.33.1.2 1996/05/06 12:20:27 mschoene
+#H changed 'DiagonalizeMat' to use norms to select the pivot
+#H
#H Revision 3.33.1.1 1994/08/24 15:22:00 sam
#H improved 'TransposedMat'
#H
@@ -126,6 +129,7 @@
##
if not IsBound(InfoMatrix1) then InfoMatrix1 := Ignore; fi;
if not IsBound(InfoMatrix2) then InfoMatrix2 := Ignore; fi;
+if not IsBound(InfoMatrix3) then InfoMatrix3 := Ignore; fi;
#############################################################################
@@ -1351,149 +1355,363 @@
fi;
end;
+
#############################################################################
##
-#F DiagonalizeMat( <mat> ) . . . . . . . . . . diagonalize an integer matrix
+#F BestQuoInt(<n>,<m>)
##
-DiagonalizeMat := function ( mat )
- local i, # current position
- k, l, # row and column loop variables
- r, c, # row and column index of minimal element
- e, f, # entries of the matrix
- g, # (extended) gcd of <e> and <f>
- v, w, # rows of the matrix or row and column norms
- m, # maximal entry, only for information
- isClearCol; #
+## 'BestQuoInt' returns the best quotient <q> of the integers <n> and <m>.
+## This is the quotient such that '<n>-<q>\*<m>' has minimal absolute value.
+## If there are two quotients whose remainders have the same absolute value,
+## then the quotient with the smaller absolute value is choosen.
+##
+BestQuoInt := function ( n, m )
+ if 0 <= m and 0 <= n then
+ return QuoInt( n + QuoInt( m - 1, 2 ), m );
+ elif 0 <= m then
+ return QuoInt( n - QuoInt( m - 1, 2 ), m );
+ elif 0 <= n then
+ return QuoInt( n - QuoInt( m + 1, 2 ), m );
+ else
+ return QuoInt( n + QuoInt( m + 1, 2 ), m );
+ fi;
+end;
- InfoMatrix1("#I DiagonalizeMat called\n");
- if mat <> [] then
+#############################################################################
+##
+#F DiagonalizeIntMatNormDriven(<mat>) . . . . diagonalize an integer matrix
+##
+## 'DiagonalizeIntMatNormDriven' diagonalizes the integer matrix <mat>.
+##
+## It tries to keep the entries small through careful selection of pivots.
+##
+## First it selects a nonzero entry for which the product of row and column
+## norm is minimal (this need not be the entry with minimal absolute value).
+## Then it brings this pivot to the upper left corner and makes it positive.
+##
+## Next it subtracts multiples of the first row from the other rows, so that
+## the new entries in the first column have absolute value at most pivot/2.
+## Likewise it subtracts multiples of the 1st column from the other columns.
+##
+## If afterwards not all new entries in the first column and row are zero,
+## then it selects a new pivot from those entries (again driven by product
+## of norms) and reduces the first column and row again.
+##
+## If finally all offdiagonal entries in the first column and row are zero,
+## then it starts all over again with the submatrix '<mat>{[2..]}{[2..]}'.
+##
+## It is based upon ideas by George Havas and code by Bohdan Majewski.
+## G. Havas and B. Majewski, Integer Matrix Diagonalization, JSC, to appear
+##
+DiagonalizeIntMatNormDriven := function ( mat )
+ local nrrows, # number of rows (length of <mat>)
+ nrcols, # number of columns (length of <mat>[1])
+ rownorms, # norms of rows
+ colnorms, # norms of columns
+ d, # diagonal position
+ pivk, pivl, # position of a pivot
+ norm, # product of row and column norms of the pivot
+ clear, # are the row and column cleared
+ row, # one row
+ col, # one column
+ ent, # one entry of matrix
+ quo, # quotient
+ h, # gap width in shell sort
+ k, l, # loop variables
+ max, omax; # maximal entry and overall maximal entry
- InfoMatrix2("#I divisors: \c");
- m := 0;
-
- # loop over all rows respectively columns of the matrix
- for i in [1..Minimum( Length(mat), Length(mat[1]) )] do
-
- # compute the row and column norms
- v := 0 * [1..Length(mat)];
- w := 0 * [1..Length(mat[1])];
- for k in [i..Length(mat)] do
- for l in [i..Length(mat[1])] do
- e := mat[k][l];
- if 0 < e then
- v[k] := v[k] + e;
- w[l] := w[l] + e;
- else
- v[k] := v[k] - e;
- w[l] := w[l] - e;
- fi;
- od;
- od;
-
- # find the element with the smallest absolute value in the matrix
- f := 0;
- for k in [i..Length(mat)] do
- for l in [i..Length(mat[1])] do
- e := mat[k][l];
- if 0 < e and (f=0 or e<f or e=f and v[k]*w[l]<g) then
- f := e; r := k; c := l; g := v[k]*w[l];
- elif e < 0 and (f=0 or -e<f or -e=f and v[k]*w[l]<g) then
- f := -e; r := k; c := l; g := v[k]*w[l];
- fi;
- if 0 < e and m < e then
- m := e;
- elif e < 0 and m < -e then
- m := -e;
- fi;
- od;
- od;
-
- # if there is no nonzero entry we are done
- if f = 0 then
- InfoMatrix2("\n");
- InfoMatrix1("#I DiagonalizeMat returns\n");
- return;
- fi;
-
- # move the minimal element to position 'mat[i][i]' make it positive
- if i <> r then
- v := mat[i]; mat[i] := mat[r]; mat[r] := v;
- fi;
- if i <> c then
- for k in [i..Length(mat)] do
- e := mat[k][i]; mat[k][i] := mat[k][c]; mat[k][c] := e;
- od;
- fi;
- if mat[i][i] < 0 then
- mat[i] := - mat[i];
- fi;
-
- # now clear the column i and the row i
- isClearCol := false;
- while not isClearCol do
-
- # clear the column i using unimodular row operations such that
- # mat[i][i] becomes gcd(mat[i][i],mat[r][i]) and mat[r][i] = 0
- k := i + 1;
- while k <= Length(mat) do
- e := mat[i][i]; f := mat[k][i];
- if f mod e = 0 then
- mat[k] := mat[k] - f/e * mat[i];
- elif f <> 0 then
- g := Gcdex( e, f );
- v := mat[i]; w := mat[k];
- mat[i] := g.coeff1 * v + g.coeff2 * w;
- mat[k] := g.coeff3 * v + g.coeff4 * w;
- fi;
- k := k + 1;
- od;
-
- isClearCol := true;
-
- # clear the row i using unimodular column operations such that
- # mat[i][i] becomes gcd(mat[i][i],mat[i][c]) and mat[i][c] = 0
- # after such an operation we may have to clear column i again
- l := i + 1;
- while l <= Length(mat[1]) and isClearCol do
- e := mat[i][i]; f := mat[i][l];
- if f mod e = 0 then
- mat[i][l] := 0;
- elif f <> 0 then
- g := Gcdex( e, f );
- for k in [i..Length(mat)] do
- v := mat[k][i]; w := mat[k][l];
- mat[k][i] := g.coeff1 * v + g.coeff2 * w;
- mat[k][l] := g.coeff3 * v + g.coeff4 * w;
- od;
- isClearCol := false;
- fi;
- l := l + 1;
- od;
-
- od;
-
- InfoMatrix2(mat[i][i]," \c");
- od;
-
- InfoMatrix2("\n");
- InfoMatrix2("#I maximal entry ",m,"\n");
+ # give some information
+ InfoMatrix1("#I DiagonalizeMat called\n");
+ InfoMatrix2("#I divisors: \c");
+ omax := 0;
+
+ # get the number of rows and columns
+ nrrows := Length( mat );
+ if nrrows <> 0 then
+ nrcols := Length( mat[1] );
+ else
+ nrcols := 0;
fi;
+ rownorms := [];
+ colnorms := [];
+
+ # loop over the diagonal positions
+ d := 1;
+ while d <= nrrows and d <= nrcols do
+ InfoMatrix3("\n#I d=",d," \c");
+
+ # find the maximal entry
+ if InfoMatrix3 <> Ignore then
+ max := 0;
+ for k in [ d .. nrrows ] do
+ for l in [ d .. nrcols ] do
+ ent := mat[k][l];
+ if 0 < ent and max < ent then
+ max := ent;
+ elif ent < 0 and max < -ent then
+ max := -ent;
+ fi;
+ od;
+ od;
+ InfoMatrix3("max=",max," \c");
+ if omax < max then omax := max; fi;
+ fi;
+
+ # compute the Euclidean norms of the rows and columns
+ for k in [ d .. nrrows ] do
+ row := mat[k];
+ rownorms[k] := row * row;
+ od;
+ for l in [ d .. nrcols ] do
+ col := mat{[d..nrrows]}[l];
+ colnorms[l] := col * col;
+ od;
+ InfoMatrix3("n\c");
+
+ # push rows containing only zeroes down and forget about them
+ for k in [ nrrows, nrrows-1 .. d ] do
+ if k < nrrows and rownorms[k] = 0 then
+ row := mat[k];
+ mat[k] := mat[nrrows];
+ mat[nrrows] := row;
+ norm := rownorms[k];
+ rownorms[k] := rownorms[nrrows];
+ rownorms[nrrows] := norm;
+ fi;
+ if rownorms[nrrows] = 0 then
+ nrrows := nrrows - 1;
+ fi;
+ od;
+
+ # quit if there are no more nonzero entries
+ if nrrows < d then
+ #N 1996/04/30 mschoene should 'break'
+ InfoMatrix3("\n#I overall maximal entry ",omax);
+ InfoMatrix1("\n#I DiagonalizeMat returns\n");
+ return;
+ fi;
- InfoMatrix1("#I DiagonalizeMat returns\n");
+ # push columns containing only zeroes right and forget about them
+ for l in [ nrcols, nrcols-1 .. d ] do
+ if l < nrcols and colnorms[l] = 0 then
+ col := mat{[d..nrrows]}[l];
+ mat{[d..nrrows]}[l] := mat{[d..nrrows]}[nrcols];
+ mat{[d..nrrows]}[nrcols] := col;
+ norm := colnorms[l];
+ colnorms[l] := colnorms[nrcols];
+ colnorms[nrcols] := norm;
+ fi;
+ if colnorms[nrcols] = 0 then
+ nrcols := nrcols - 1;
+ fi;
+ od;
+
+ # sort the rows with respect to their norms
+ h := 1; while 9 * h + 4 < nrrows-(d-1) do h := 3 * h + 1; od;
+ while 0 < h do
+ for l in [ h+1 .. nrrows-(d-1) ] do
+ norm := rownorms[l+(d-1)];
+ row := mat[l+(d-1)];
+ k := l;
+ while h+1 <= k and norm < rownorms[k-h+(d-1)] do
+ rownorms[k+(d-1)] := rownorms[k-h+(d-1)];
+ mat[k+(d-1)] := mat[k-h+(d-1)];
+ k := k - h;
+ od;
+ rownorms[k+(d-1)] := norm;
+ mat[k+(d-1)] := row;
+ od;
+ h := QuoInt( h, 3 );
+ od;
+
+ # choose a pivot in the '<mat>{[<d>..]}{[<d>..]}' submatrix
+ # the pivot must be the topmost nonzero entry in its column,
+ # now that the rows are sorted with respect to their norm
+ pivk := 0; pivl := 0;
+ norm := Maximum(rownorms) * Maximum(colnorms) + 1;
+ for l in [ d .. nrcols ] do
+ k := d;
+ while mat[k][l] = 0 do
+ k := k + 1;
+ od;
+ if rownorms[k] * colnorms[l] < norm then
+ pivk := k; pivl := l;
+ norm := rownorms[k] * colnorms[l];
+ fi;
+ od;
+ InfoMatrix3("p\c");
+
+ # move the pivot to the diagonal and make it positive
+ if d <> pivk then
+ row := mat[d];
+ mat[d] := mat[pivk];
+ mat[pivk] := row;
+ fi;
+ if d <> pivl then
+ col := mat{[d..nrrows]}[d];
+ mat{[d..nrrows]}[d] := mat{[d..nrrows]}[pivl];
+ mat{[d..nrrows]}[pivl] := col;
+ fi;
+ if mat[d][d] < 0 then
+ mat[d] := - mat[d];
+ fi;
+
+ # now perform row operations so that the entries in the
+ # <d>-th column have absolute value at most pivot/2
+ clear := true;
+ row := mat[d];
+ for k in [ d+1 .. nrrows ] do
+ quo := BestQuoInt( mat[k][d], mat[d][d] );
+ if quo = 1 then
+ mat[k] := mat[k] - row;
+ elif quo = -1 then
+ mat[k] := mat[k] + row;
+ elif quo <> 0 then
+ mat[k] := mat[k] - quo * row;
+ fi;
+ clear := clear and mat[k][d] = 0;
+ od;
+ InfoMatrix3("c\c");
+
+ # now perform column operations so that the entries in
+ # the <d>-th row have absolute value at most pivot/2
+ col := mat{[d..nrrows]}[d];
+ for l in [ d+1 .. nrcols ] do
+ quo := BestQuoInt( mat[d][l], mat[d][d] );
+ if quo = 1 then
+ mat{[d..nrrows]}[l] := mat{[d..nrrows]}[l] - col;
+ elif quo = -1 then
+ mat{[d..nrrows]}[l] := mat{[d..nrrows]}[l] + col;
+ elif quo <> 0 then
+ mat{[d..nrrows]}[l] := mat{[d..nrrows]}[l] - quo * col;
+ fi;
+ clear := clear and mat[d][l] = 0;
+ od;
+ InfoMatrix3("r \c");
+
+ # repeat until the <d>-th row and column are totally cleared
+ while not clear do
+
+ # compute the Euclidean norms of the rows and columns
+ # that have a nonzero entry in the <d>-th column resp. row
+ for k in [ d .. nrrows ] do
+ if mat[k][d] <> 0 then
+ row := mat[k];
+ rownorms[k] := row * row;
+ fi;
+ od;
+ for l in [ d .. nrcols ] do
+ if mat[d][l] <> 0 then
+ col := mat{[d..nrrows]}[l];
+ colnorms[l] := col * col;
+ fi;
+ od;
+ InfoMatrix3("n\c");
+
+ # choose a pivot in the <d>-th row or <d>-th column
+ pivk := 0; pivl := 0;
+ norm := Maximum(rownorms) * Maximum(colnorms) + 1;
+ for l in [ d+1 .. nrcols ] do
+ if 0 <> mat[d][l] and rownorms[d] * colnorms[l] < norm then
+ pivk := d; pivl := l;
+ norm := rownorms[d] * colnorms[l];
+ fi;
+ od;
+ for k in [ d+1 .. nrrows ] do
+ if 0 <> mat[k][d] and rownorms[k] * colnorms[d] < norm then
+ pivk := k; pivl := d;
+ norm := rownorms[k] * colnorms[d];
+ fi;
+ od;
+ InfoMatrix3("p\c");
+
+ # move the pivot to the diagonal and make it positive
+ if d <> pivk then
+ row := mat[d];
+ mat[d] := mat[pivk];
+ mat[pivk] := row;
+ fi;
+ if d <> pivl then
+ col := mat{[d..nrrows]}[d];
+ mat{[d..nrrows]}[d] := mat{[d..nrrows]}[pivl];
+ mat{[d..nrrows]}[pivl] := col;
+ fi;
+ if mat[d][d] < 0 then
+ mat[d] := - mat[d];
+ fi;
+
+ # now perform row operations so that the entries in the
+ # <d>-th column have absolute value at most pivot/2
+ clear := true;
+ row := mat[d];
+ for k in [ d+1 .. nrrows ] do
+ quo := BestQuoInt( mat[k][d], mat[d][d] );
+ if quo = 1 then
+ mat[k] := mat[k] - row;
+ elif quo = -1 then
+ mat[k] := mat[k] + row;
+ elif quo <> 0 then
+ mat[k] := mat[k] - quo * row;
+ fi;
+ clear := clear and mat[k][d] = 0;
+ od;
+ InfoMatrix3("c\c");
+
+ # now perform column operations so that the entries in
+ # the <d>-th row have absolute value at most pivot/2
+ col := mat{[d..nrrows]}[d];
+ for l in [ d+1.. nrcols ] do
+ quo := BestQuoInt( mat[d][l], mat[d][d] );
+ if quo = 1 then
+ mat{[d..nrrows]}[l] := mat{[d..nrrows]}[l] - col;
+ elif quo = -1 then
+ mat{[d..nrrows]}[l] := mat{[d..nrrows]}[l] + col;
+ elif quo <> 0 then
+ mat{[d..nrrows]}[l] := mat{[d..nrrows]}[l] - quo * col;
+ fi;
+ clear := clear and mat[d][l] = 0;
+ od;
+ InfoMatrix3("r \c");
+
+ od;
+
+ # print the diagonal entry (for information only)
+ InfoMatrix3("div=");
+ InfoMatrix2(mat[d][d]," \c");
+
+ # go on to the next diagonal position
+ d := d + 1;
+
+ od;
+
+ # close with some more information
+ InfoMatrix3("\n#I overall maximal entry ",omax);
+ InfoMatrix1("\n#I DiagonalizeMat returns\n");
end;
+DiagonalizeIntMat := DiagonalizeIntMatNormDriven;
+
+
+#############################################################################
+##
+#F DiagonalizeMat(<mat>) . . . . . . . . . . . . . . . diagonalize a matrix
+##
+#N 1996/05/06 mschoene should be extended for other rings
+##
+DiagonalizeMat := DiagonalizeIntMat;
+
#############################################################################
##
-#F ElementaryDivisorsMat( <mat> ) elementary divisors of an integer matrix
+#F ElementaryDivisorsMat(<mat>) . . . . . . elementary divisors of a matrix
##
## 'ElementaryDivisors' returns a list of the elementary divisors, i.e., the
## unique <d> with '<d>[<i>]' divides '<d>[<i>+1]' and <mat> is equivalent
## to a diagonal matrix with the elements '<d>[<i>]' on the diagonal.
##
ElementaryDivisorsMat := function ( mat )
- local divs, gcd, m, n, i, k;
+ local divs, gcd, zero, m, n, i, k;
# make a copy to avoid changing the original argument
mat := Copy( mat );
@@ -1507,16 +1725,22 @@
for i in [1..Minimum(m,n)] do
divs[i] := mat[i][i];
od;
+ if divs <> [] then zero := divs[1] - divs[1]; fi;
# transform the divisors so that every divisor divides the next
for i in [1..Length(divs)-1] do
for k in [i+1..Length(divs)] do
- if divs[i] <> 0 and divs[k] mod divs[i] <> 0 then
- gcd := GcdInt( divs[i], divs[k] );
+ if divs[i] = zero and divs[k] <> zero then
+ divs[i] := divs[k];
+ divs[k] := zero;
+ elif divs[i] <> zero
+ and EuclideanRemainder( divs[k], divs[i] ) <> zero then
+ gcd := Gcd( divs[i], divs[k] );
divs[k] := divs[k] / gcd * divs[i];
divs[i] := gcd;
fi;
od;
+ divs[i] := StandardAssociate( divs[i] );
od;
return divs;
END OF bugfix08.dif ________________________________________________________