> < ^ Date: Fri, 10 May 1996 16:42:00 +0100 (MET)
> < ^ From: Martin Schoenert <martin.schoenert@math.rwth-aachen.de >
^ Subject: BUGFIX 08 (Performance problem in 'DiagonalizeMatrix')

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.dif

This 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 ________________________________________________________

> < [top]