  
  [1m[4m[31m1. The [1mEDIM[1m[4m[31m-Package[0m
  
  [22m[36m(Elementary Divisors and Integer Matrices, by Frank Lübeck)[0m
  
  This  chapter describes the functions defined in the mini [1mGAP[0m4 package [1mEDIM[0m.
  The  main  functions  implement variants of an algorithm for computing for a
  given  prime  p the p-parts of the elementary divisors of an integer matrix.
  These  algorithms  use  a  p-adic  method and are described by the author in
  [L01] (see [1m[34mElementaryDivisorsPPartRk[0m ([1m1.2-1[0m)).
  
  These  functions  were  already  applied  to  integer  matrices of dimension
  greater  than  11000  (which  had many non-trivial elementary divisors which
  were products of small primes).
  
  Furthermore  there  are functions for finding the biggest elementary divisor
  of  an  invertible  integer  matrix and the inverse of a rational invertible
  matrix (see [1m[34mExponentSquareIntMatFullRank[0m ([1m1.3-2[0m) and [1m[34mInverseRatMat[0m ([1m1.3-1[0m)).
  These algorithms use p-adic approximations, explained in [1m1.7[0m.
  
  Finally  we  distribute implementations of some other algorithms for finding
  elementary  divisors  or  normal  forms  of  integer  matrices:  A p-modular
  algorithm     by     Havas     and     Sterling     from     [HS79]     (see
  [1m[34mElementaryDivisorsPPartHavasSterling[0m  ([1m1.2-2[0m))  and LLL-based algorithms for
  extended  greatest common divisors of integers (see [1m[34mGcdexIntLLL[0m ([1m1.5-1[0m)) and
  for  Hermite  normal forms of integer matrices with (very nice) transforming
  matrices (see [1m[34mHermiteIntMatLLL[0m ([1m1.5-2[0m)).
  
  By  default the [1mEDIM[0m is automatically loaded by [1mGAP[0m when it is installed. If
  the  automatic  loading  is  disabled in your installation you must load the
  package with [22m[32mRequirePackage("edim");[0m before its functions become available.
  
  Please,  send me an e-mail ([34mmailto:Frank.Luebeck@Math.RWTH-Aachen.De[0m) if you
  have any questions, remarks, suggestions, etc. concerning this mini-package.
  Also, I would like to hear about applications of this package.
  
  Frank Lübeck
  
  
  [1m[4m[31m1.1 Installation of the [1mEDIM[1m[4m[31m-package[0m
  
  To  install  this package (after extracting the packages archive file to the
  GAP  home  directory) go to the directory [1mpkg/edim[0m (the directory containing
  this README file) and call
  
  [22m[32m/bin/sh ./configure [path][0m
  
  where  [22m[32mpath[0m  is  a  path  to  the main [1mGAP[0m root directory (if not given, the
  default [1m../..[0m is assumed).
  
  Afterwards call [22m[32mmake[0m to compile a binary file.
  
  If  you  installed  GAP  on  several  architectures,  you  must execute this
  [22m[32mconfigure/make[0m   step   on  each  of  the  architectures  immediately  after
  configuring  GAP  itself  on  this  architecture. The package will also work
  without      this      step     but     then     the     kernel     function
  [1m[34mElementaryDivisorsPPartRkExpSmall[0m                ([1m1.2-1[0m)                (see
  [1m[34mElementaryDivisorsPPartRkExpSmall[0m ([1m1.2-1[0m)) will not be available.
  
  The  [1mINSTALL[0m  in the main package directory also explains how to compile the
  kernel  function  into  a  static  module.  You  can  also run a test of the
  installation by typing [22m[32mmake test[0m.
  
  [1m[4m[31m1.1-1 InfoEDIM[0m
  
  [1m[34m> InfoEDIM________________________________________________________[0minfo class
  
  This  is  an  [22m[32mInfo[0m class for the [1mEDIM[0m-package. By [22m[32mSetInfoLevel(InfoEDIM, 1);[0m
  you  can  switch on the printing of some information during the computations
  of certain [1mEDIM[0m-functions.
  
  
  [1m[4m[31m1.2 p-Parts of Elementary Divisors[0m
  
  Here we explain the main functions of the package.
  
  [1m[4m[31m1.2-1 ElementaryDivisorsPPartRk[0m
  
  [1m[34m> ElementaryDivisorsPPartRk( [0m[22m[34mA, p[, rk][0m[1m[34m ) __________________________[0mfunction
  [1m[34m> ElementaryDivisorsPPartRkI( [0m[22m[34mA, p, rk[0m[1m[34m ) ___________________________[0mfunction
  [1m[34m> ElementaryDivisorsPPartRkII( [0m[22m[34mA, p, rk[0m[1m[34m ) __________________________[0mfunction
  [1m[34m> ElementaryDivisorsPPartRkExp( [0m[22m[34mA, p, rk, exp[0m[1m[34m ) ____________________[0mfunction
  [1m[34m> ElementaryDivisorsPPartRkExpSmall( [0m[22m[34mA, p, rk, exp, il[0m[1m[34m ) ___________[0mfunction
  
  These  functions  return a list [m_1, m_2, ..., m_r] where m_i is the number
  of    nonzero   elementary   divisors   of   [22m[34mA[0m   divisible   by   [22m[34mp[0m^i   (see
  [1m[34mElementaryDivisorsMat[0m ([1mReference: ElementaryDivisorsMat[0m) for a definition of
  the elementary divisors).
  
  The algorithms for these functions are described in [L01].
  
  [22m[34mA[0m must be a matrix with integer entries, [22m[34mp[0m a prime, and [22m[34mrk[0m the rank of [22m[34mA[0m (as
  rational  matrix). In the first version of the command [22m[34mrk[0m is computed, if it
  is not given.
  
  The  first version of the command delegates its job to the fourth version by
  trying growing values for [22m[34mexp[0m, see below.
  
  The  second  and  third  versions  implement the main algorithm described in
  [L01]  and  a  variation.  Here  [1m[34mElementaryDivisorsPPartRkII[0m  has a bit more
  overhead,  but  can  be advantageous because the intermediate entries during
  the computation can be much smaller.
  
  In  the  fourth  form  [22m[34mexp[0m must be an upper bound for the highest power of [22m[34mp[0m
  appearing  in  an elementary divisor of [22m[34mA[0m. This information allows reduction
  of matrix entries modulo [22m[34mp[0m^[22m[34mexp[0m during the computation.
  
  If  [22m[34mexp[0m  is  too  small  or the given [22m[34mrk[0m is not correct the function returns
  `fail'.
  
  As  long as [22m[34mp[0m^[22m[34mexp[0m is smaller than 2^28 and [22m[34mp[0m^[22m[34mexp[0m + 2 is smaller than 2^31 we
  use  internally  a  kernel  function  which can also be used directly in the
  fifth  form  of the command. There [22m[34mil[0m can be 0 or 1 where in the second case
  some information is printed during the computation.
  
  This  last  form  of  the  function was already succesfully applied to dense
  matrices of rank up to 11000.
  
  Note  that  you  have  to  compile  a  file  (see [1m1.1[0m) while installing this
  package, if you want to have this kernel function available.
  
  [22m[35m---------------------------  Example  ----------------------------[0m
    [22m[35mgap> ReadPkg("edim",  "tst/mat");[0m
    [22m[35mReading 242x242 integer matrix 'mat' and its elementary divisors 'eldiv'.[0m
    [22m[35mgap> ElementaryDivisorsPPartRkI(mat, 2, 242); time; # mat has full rank[0m
    [22m[35m[ 94, 78, 69, 57, 23, 23, 9, 2, 2, 0 ][0m
    [22m[35m9170[0m
    [22m[35mgap> ElementaryDivisorsPPartRkExpSmall(mat, 2, 242, 10, 0); time;[0m
    [22m[35m[ 94, 78, 69, 57, 23, 23, 9, 2, 2, 0 ][0m
    [22m[35m560[0m
  [22m[35m------------------------------------------------------------------[0m
  
  [1m[4m[31m1.2-2 ElementaryDivisorsPPartHavasSterling[0m
  
  [1m[34m> ElementaryDivisorsPPartHavasSterling( [0m[22m[34mA, p, d[0m[1m[34m ) __________________[0mfunction
  
  For  an  integer  matrix  [22m[34mA[0m and a prime [22m[34mp[0m this function returns a list [m_1,
  m_2,  ...,  m_r] where m_i is the number of nonzero elementary divisors of [22m[34mA[0m
  divisible by [22m[34mp[0m^i.
  
  An  upper  bound  [22m[34md[0m  for  the  highest power of [22m[34mp[0m appearing in an elementary
  divisor  of  [22m[34mA[0m  must  be  given.  Smaller  [22m[34md[0m  improve the performance of the
  algorithm considerably.
  
  This is an implementation of the modular algorithm described in [HS79].
  
  We  added  a slight improvement: we divide the considered submatrices by the
  [22m[34mp[0m-part  of  the  greatest  common  divisor  of  all entries (and lower the [22m[34md[0m
  appropriately).  This reduces the size of the entries and often shortens the
  pivot search.
  
  [22m[35m---------------------------  Example  ----------------------------[0m
    [22m[35mgap> ReadPkg("edim",  "tst/mat");[0m
    [22m[35mReading 242x242 integer matrix 'mat' and its elementary divisors 'eldiv'.[0m
    [22m[35mgap> ElementaryDivisorsPPartHavasSterling(mat, 2, 10); time;[0m
    [22m[35m[ 94, 78, 69, 57, 23, 23, 9, 2, 2 ][0m
    [22m[35m25390[0m
  [22m[35m------------------------------------------------------------------[0m
  
  
  [1m[4m[31m1.3 Inverse of Rational Matrices[0m
  
  [1m[4m[31m1.3-1 InverseRatMat[0m
  
  [1m[34m> InverseRatMat( [0m[22m[34mA[, p][0m[1m[34m ) __________________________________________[0mfunction
  
  This  function returns the inverse of an invertible matrix over the rational
  numbers.
  
  It  first  computes  the  inverse  modulo some prime [22m[34mp[0m, computes from this a
  [22m[34mp[0m-adic  approximation  to  the  inverse  and finally constructs the rational
  entries from their [22m[34mp[0m-adic approximations. See section [1m1.7[0m for more details.
  
  This seems to be better than [1mGAP[0m's standard Gauß algorithm ([22m[32mA\^{}-1[0m) already
  for   small   matrices.   (Try,  e.g.,  [22m[32mRandomMat(20,20,[-10000..10000])[0m  or
  [22m[32mRandomMat(100,100)[0m.)
  
  The optional argument [22m[34mp[0m should be a prime such that [22m[34mA[0m modulo [22m[34mp[0m is invertible
  (default  is [22m[34mp[0m=251). If [22m[34mA[0m is not invertible modulo [22m[34mp[0m then [22m[34mp[0m is automatically
  replaced by the next prime.
  
  [1m[4m[31m1.3-2 ExponentSquareIntMatFullRank[0m
  
  [1m[34m> ExponentSquareIntMatFullRank( [0m[22m[34mA[, p[, nr]][0m[1m[34m ) _____________________[0mfunction
  
  This  function  returns  the  biggest elementary divisor of a square integer
  matrix [22m[34mA[0m of full rank.
  
  For  such  a  matrix  [22m[34mA[0m the least common multiple of the denominators of all
  entries of the inverse matrix [22m[34mA[0m^-1 is exactly the biggest elementary divisor
  of [22m[34mA[0m.
  
  This  function  is  implemented  by  a  slight modification of [1m[34mInverseRatMat[0m
  ([1m1.3-1[0m). The third argument [22m[34mnr[0m tells the function to return the least common
  multiple  of  the  first  [22m[34mnr[0m  rows of the rational inverse matrix only. Very
  often  the  function will already return the biggest elementary divisor with
  [22m[34mnr[0m=2  or  3  (and the command without this argument would spend most time in
  checking, that this is correct).
  
  The optional argument [22m[34mp[0m should be a prime such that [22m[34mA[0m modulo [22m[34mp[0m is invertible
  (default  is [22m[34mp[0m=251). If [22m[34mA[0m is not invertible modulo [22m[34mp[0m then [22m[34mp[0m is automatically
  replaced by the next prime.
  
  [22m[35m---------------------------  Example  ----------------------------[0m
    [22m[35mgap> ReadPkg("edim",  "tst/mat");[0m
    [22m[35mReading 242x242 integer matrix 'mat' and its elementary divisors 'eldiv'.[0m
    [22m[35mgap> inv := InverseRatMat(mat);; time;                      [0m
    [22m[35m14960[0m
    [22m[35mgap> ExponentSquareIntMatFullRank(mat, 101, 3); # same as without the `3'[0m
    [22m[35m115200[0m
  [22m[35m------------------------------------------------------------------[0m
  
  
  [1m[4m[31m1.4 All Elementary Divisors Using p-adic Method[0m
  
  In  the  following  two  functions  we put things together. In particular we
  handle  the  prime  parts  of the elementary divisors efficiently for primes
  appearing  with  low  powers  in the highest elementary divisor respectively
  determinant divisor.
  
  [1m[4m[31m1.4-1 ElementaryDivisorsSquareIntMatFullRank[0m
  
  [1m[34m> ElementaryDivisorsSquareIntMatFullRank( [0m[22m[34mA[0m[1m[34m ) ______________________[0mfunction
  
  This  function  returns  a list of nonzero elementary divisors of an integer
  matrix [22m[34mA[0m.
  
  Here   we   start   with   computing  the  biggest  elementary  divisor  via
  [1m[34mExponentSquareIntMatFullRank[0m ([1m1.3-2[0m). If it runs into a problem because [22m[34mA[0m is
  singular  modulo  a  choosen  prime (it starts by default with 251) then the
  prime is automatically replaced by the next one.
  
  The  rest  is  done  using  [1m[34mElementaryDivisorsPPartRkExp[0m ([1m1.2-1[0m) and [1m[34mRankMod[0m
  ([1m1.6-5[0m).
  
  The  function  fails  if the biggest elementary divisor cannot be completely
  factored  and  the  non-factored  part  is  not  a  divisor  of  the biggest
  elementary divisor only.
  
  Note  that  this  function  may for many matrices not be the best choice for
  computing  all  elementary  divisors.  You  may  first  try the standard [1mGAP[0m
  library   routines   for   Smith  normal  form  instead  of  this  function.
  Nevertheless  remember  [1m[34mElementaryDivisorsSquareIntMatFullRank[0m  for hard and
  big examples. It is particularly good when the largest elementary divisor is
  a very small factor of the determinant.
  
  [22m[35m---------------------------  Example  ----------------------------[0m
    [22m[35mgap> Collected(ElementaryDivisorsSquareIntMatFullRank(mat));      [0m
    [22m[35m[ [ 1, 49 ], [ 3, 99 ], [ 6, 7 ], [ 30, 9 ], [ 60, 9 ], [ 120, 2 ], [0m
    [22m[35m  [ 360, 10 ], [ 720, 22 ], [ 3600, 12 ], [ 14400, 14 ], [ 28800, 7 ], [0m
    [22m[35m  [ 115200, 2 ] ][0m
    [22m[35mgap> time;[0m
    [22m[35m10180[0m
    [22m[35mgap> last2 = Collected(DiagonalOfMat(NormalFormIntMat(mat, 1).normal));[0m
    [22m[35mtrue[0m
    [22m[35mgap> time;[0m
    [22m[35m188490[0m
  [22m[35m------------------------------------------------------------------[0m
  
  [1m[4m[31m1.4-2 ElementaryDivisorsIntMatDeterminant[0m
  
  [1m[34m> ElementaryDivisorsIntMatDeterminant( [0m[22m[34mA, det[, rk][0m[1m[34m ) ______________[0mfunction
  
  This  function  returns  a list of nonzero elementary divisors of an integer
  matrix [22m[34mA[0m.
  
  Here  [22m[34mdet[0m  must be an integer which is a multiple of the biggest determinant
  divisor of [22m[34mA[0m. If the matrix does not have full rank then its rank [22m[34mrk[0m must be
  given, too.
  
  The argument [22m[34mdet[0m can be given in the form of [22m[32mCollected(FactorsInt([22m[34mdet[0m))[0m.
  
  This function handles prime divisors of [22m[34mdet[0m with multiplicity smaller than 4
  specially,    for   the   other   prime   divisors   p   it   delegates   to
  [1m[34mElementaryDivisorsPPartRkExp[0m   ([1m1.2-1[0m)   where   the  [22m[34mexp[0m  argument  is  the
  multiplicity  of  the  p in [22m[34mdet[0m. (Note that this is not very good when p has
  actually a much smaller multiplicity in the largest elementary divisor.)
  
  [22m[35m---------------------------  Example  ----------------------------[0m
    [22m[35mgap> ReadPkg("edim",  "tst/mat");[0m
    [22m[35mReading 242x242 integer matrix 'mat' and its elementary divisors 'eldiv'.[0m
    [22m[35mgap> # not so good:[0m
    [22m[35mgap> ElementaryDivisorsIntMatDeterminant(mat,Product(eldiv)) = [0m
    [22m[35mConcatenation([1..49]*0+1, eldiv); time;[0m
    [22m[35mtrue[0m
    [22m[35m185770[0m
  [22m[35m------------------------------------------------------------------[0m
  
  
  [1m[4m[31m1.5 Gcd and Normal Forms Using LLL[0m
  
  The   [1mEDIM[0m-mini   package  also  contains  implementations  of  an  extended
  Gcd-algorithm for integers and a Hermite and Smith normal form algorithm for
  integer  matrices  using LLL-techiques. They are well described in the paper
  [HMM98] by Havas, Majewski and Matthews.
  
  They  are particularly useful if one wants to have the normal forms together
  with  transforming  matrices. These transforming matrices have spectacularly
  nice (i.e., "small") entries in cases of input matrices which are non-square
  or not of full rank (otherwise the transformation to the Hermite normal form
  is unique).
  
  In detail:
  
  [1m[4m[31m1.5-1 GcdexIntLLL[0m
  
  [1m[34m> GcdexIntLLL( [0m[22m[34mn1, n2, ...[0m[1m[34m ) _______________________________________[0mfunction
  
  This  function returns for integers [22m[34mn1[0m, [22m[34mn2[0m, ... a list [g, [c_1, c_2, ...]],
  where g = c_1[22m[34mn1[0m + c_2[22m[34mn2[0m + ... is the greatest common divisor of the [22m[34mni[0m. Here
  all the c_i are usually very small.
  
  [22m[35m---------------------------  Example  ----------------------------[0m
    [22m[35mgap> GcdexIntLLL( 517, 244, -304, -872, -286, 854, 866, 224, -765, -38);[0m
    [22m[35m[ 1, [ 0, 0, 0, 0, 1, 0, 1, 1, 1, 1 ] ][0m
  [22m[35m------------------------------------------------------------------[0m
  
  [1m[4m[31m1.5-2 HermiteIntMatLLL[0m
  
  [1m[34m> HermiteIntMatLLL( [0m[22m[34mA[0m[1m[34m ) ____________________________________________[0mfunction
  
  This  returns  the  Hermite  normal form of an integer matrix [22m[34mA[0m and uses the
  LLL-algorithm to avoid entry explosion.
  
  [1m[4m[31m1.5-3 HermiteIntMatLLLTrans[0m
  
  [1m[34m> HermiteIntMatLLLTrans( [0m[22m[34mA[0m[1m[34m ) _______________________________________[0mfunction
  
  This function returns a pair of matrices [H, L] where H = L [22m[34mA[0m is the Hermite
  normal  form  of  an  integer  matrix  [22m[34mA[0m. The transforming matrix L can have
  surprisingly small entries.
  
  [22m[35m---------------------------  Example  ----------------------------[0m
    [22m[35mgap> ReadPkg("edim",  "tst/mat2");[0m
    [22m[35mReading 34x34 integer matrix 'mat2' and its elementary divisors 'eldiv2'.[0m
    [22m[35mgap> tr := HermiteIntMatLLLTrans(mat2);; Maximum(List(Flat(tr[2]), AbsInt));[0m
    [22m[35m606[0m
    [22m[35mgap> tr[2]*mat2 = tr[1];                                                [0m
    [22m[35mtrue[0m
  [22m[35m------------------------------------------------------------------[0m
  
  [1m[4m[31m1.5-4 SmithIntMatLLL[0m
  
  [1m[34m> SmithIntMatLLL( [0m[22m[34mA[0m[1m[34m ) ______________________________________________[0mfunction
  
  This function returns the Smith normal form of an integer matrix [22m[34mA[0m using the
  LLL-algorithm to avoid entry explosion.
  
  [1m[4m[31m1.5-5 SmithIntMatLLLTrans[0m
  
  [1m[34m> SmithIntMatLLLTrans( [0m[22m[34mA[0m[1m[34m ) _________________________________________[0mfunction
  
  This  function returns [S, L, R] where S = L [22m[34mA[0m R is the Smith normal form of
  an integer matrix [22m[34mA[0m.
  
  We  apply  the  algorithm  for  Hermite normal form several times to get the
  Smith  normal  form,  that  is  not  in  the paper [HMM98]. The transforming
  matrices need not be as nice as for the Hermite normal form.
  
  [22m[35m---------------------------  Example  ----------------------------[0m
    [22m[35mgap> ReadPkg("edim",  "tst/mat2");[0m
    [22m[35mReading 34x34 integer matrix 'mat2' and its elementary divisors 'eldiv2'.[0m
    [22m[35mgap> tr := SmithIntMatLLLTrans(mat2);;[0m
    [22m[35mgap> tr[2] * mat2 * tr[3] = tr[1];    [0m
    [22m[35mtrue[0m
  [22m[35m------------------------------------------------------------------[0m
  
  
  [1m[4m[31m1.6 Utility Functions from the [1mEDIM[1m[4m[31m-package[0m
  
  [1m[4m[31m1.6-1 RatNumberFromModular[0m
  
  [1m[34m> RatNumberFromModular( [0m[22m[34mn, k, l, x[0m[1m[34m ) _______________________________[0mfunction
  
  This function returns r/s = [22m[34mx[0m mod [22m[34mn[0m, if it exists. More precisely:
  
  [22m[34mn[0m, [22m[34mk[0m, [22m[34ml[0m must be positive integers with 2[22m[34mk[0m[22m[34ml[0m <= [22m[34mn[0m and [22m[34mx[0m an integer with -[22m[34mn[0m/2 <
  [22m[34mx[0m  <= [22m[34mn[0m/2. If it exists this function returns a rational number r/s with 0 <
  s < [22m[34ml[0m, gcd(s, [22m[34mn[0m) = 1, -[22m[34mk[0m < r < [22m[34mk[0m and r/s congruent to [22m[34mx[0m mod [22m[34mn[0m (i.e., [22m[34mn[0m | r -
  s [22m[34mx[0m). Such an r/s is unique. The function returns [22m[32mfail[0m if such a number does
  not exist.
  
  [1m[4m[31m1.6-2 InverseIntMatMod[0m
  
  [1m[34m> InverseIntMatMod( [0m[22m[34mA, p[0m[1m[34m ) _________________________________________[0mfunction
  
  This  function  returns  an  inverse  matrix  modulo a prime [22m[34mp[0m or [22m[32mfail[0m. More
  precisely:
  
  [22m[34mA[0m  must  be an integer matrix and [22m[34mp[0m a prime such that [22m[34mA[0m is invertible modulo
  [22m[34mp[0m.  This  function  returns  an integer matrix [22m[34minv[0m with entries in the range
  ]-[22m[34mp[0m/2 ... [22m[34mp[0m/2] such that [22m[34minv[0m[22m[34mA[0m reduced modulo p is the identity matrix.
  
  It  returns  [22m[32mfail[0m  if  the inverse modulo [22m[34mp[0m does not exist. This function is
  particularly fast for primes smaller 256.
  
  [1m[4m[31m1.6-3 HadamardBoundIntMat[0m
  
  [1m[34m> HadamardBoundIntMat( [0m[22m[34mA[0m[1m[34m ) _________________________________________[0mfunction
  
  The Hadamard bound for a square integer matrix [22m[34mA[0m is the product of Euclidean
  norms  of  the  nonzero rows (or columns) of [22m[34mA[0m. It is an upper bound for the
  absolute value of the determinant of [22m[34mA[0m.
  
  [1m[4m[31m1.6-4 CheapFactorsInt[0m
  
  [1m[34m> CheapFactorsInt( [0m[22m[34mn[, nr][0m[1m[34m ) _______________________________________[0mfunction
  
  This  function  returns a list of factors of an integer [22m[34mn[0m, including "small"
  prime  factors  -  here the optional argument [22m[34mnr[0m is the number of iterations
  for `FactorsRho' (default is 2000).
  
  This  is  only  a  slight  modification  of  the library function [1m[34mFactorsInt[0m
  ([1mReference: FactorsInt[0m) which avoids an error message when the number is not
  completely factored.
  
  [1m[4m[31m1.6-5 RankMod[0m
  
  [1m[34m> RankMod( [0m[22m[34mA, p[0m[1m[34m ) __________________________________________________[0mfunction
  
  This  function returns the rank of an integer matrix [22m[34mA[0m modulo [22m[34mp[0m. Here [22m[34mp[0m must
  not  necessarily  be  a  prime.  If  it  is not and this function returns an
  integer, then this is the rank of [22m[34mA[0m for all prime divisors of [22m[34mp[0m.
  
  If  during the computation a factorisation of [22m[34mp[0m is found (because some pivot
  entry  has  nontrivial  greatest common divisor with [22m[34mp[0m) then the function is
  recursively applied to the found factors [22m[32mf\_i[0m of [22m[34mp[0m. The result is then given
  in the form [22m[32m[[f\_1, rk\_1], [f\_2, rk\_2], ...][0m.
  
  The  idea  to  make  this  function useful for non primes was to use it with
  large   factors   of  the  biggest  elementary  divisor  of  [22m[34mA[0m  whose  prime
  factorization cannot be found easily.
  
  [22m[35m---------------------------  Example  ----------------------------[0m
    [22m[35mgap> ReadPkg("edim",  "tst/mat");[0m
    [22m[35mReading 242x242 integer matrix 'mat' and its elementary divisors 'eldiv'.[0m
    [22m[35mgap> RankMod(mat, 5);[0m
    [22m[35m155[0m
    [22m[35mgap> RankMod(mat, (2*79*4001));[0m
    [22m[35m[ [ 2, 148 ], [ 79, 242 ], [ 4001, 242 ] ][0m
  [22m[35m------------------------------------------------------------------[0m
  
  
  [1m[4m[31m1.7 InverseRatMat - the Algorithm[0m
  
  The  idea  is  to recover a rational matrix from an l-adic approximation for
  some prime l. This description came out of discussions with Jürgen Müller. I
  thank  John  Cannon for pointing out that the basic idea already appeared in
  the paper [D82] of Dixon.
  
  Let A be an invertible matrix over the rational numbers. By multiplying with
  a constant we may assume that its entries are in fact integers.
  
  (1)  We  first  describe how to find an l-adic approximation of A^-1. Find a
  prime  l  such that A is invertible modulo l and let B be the integer matrix
  with  entries  in  the  range  ]-l/2,l/2]  such  that BA is congruent to the
  identity  matrix  modulo  l.  (This  can  be  computed  fast  by  usual Gauß
  elimination.)
  
  Now  let  v  in Z^r be a row vector. Define two sequences v_i and x_i of row
  vectors  in  Z^r by: x_0 := 0 in Z^r, v_0 := -v and for i > 0 set x_i to the
  vector  congruent  to  -v_i-1  B modulo l having entries in the range ]-l/2,
  l/2]. Then all entries of x_i A + v_i-1 are divisible by l and we set v_i :=
  (1/l) * (x_i A + v_i-1).
  
  Induction  shows that for y_i := sum_k=1^i l^k-1 x_k we have y_i A = v + l^i
  v_i  for  all  i  >=  0.  Hence  the  sequence  y_i, i >= 0, gives an l-adic
  approximation to the vector y in Q^r with y A = v.
  
  (2)  The  second  point  is  to  show  how  we  can  get the vector y from a
  sufficiently  good  approximation y_i. Note that the sequence of y_i becomes
  constant  for  i  >=  i_0 if all entries of y are integers of absolute value
  smaller  than  l^i_0 / 2 because of our choice of representatives of residue
  classes modulo l in the interval ]-l/2, l/2].
  
  More  generally  consider a / b in Q with b > 0 and a, b coprime. Then there
  is  for each n in N which is coprime to b a unique c in Z with -n / 2 < c <=
  n  /  2 and a = c b mod n. This c can be computed via the extended Euclidean
  algorithm applied to b and n.
  
  Now  let n, alpha, beta in N with 2 alpha beta <= n. Then the map a/b in Q |
  -alpha  <=  a  <=  alpha, 1 <= b < beta -> ]-n/2, n/2], a/b -> c (defined as
  above) is injective (since for a/b, a'/b' in the above set we have a b' - a'
  b = 0 mod n if and only if a b' - a' b = 0).
  
  In practice we can use for any c in ]-n/2, n/2] a certain extended Euclidean
  algorithm applied to n and c to decide if c is in the image of the above map
  and to find the corresponding a/b if it exists.
  
  (3)  To  put  things together we apply (2) to the entries of the vectors y_i
  constructed  in  (1), choosing n = l^i, alpha = sqrtn/2 and beta = sqrtn. If
  we  have  found  this  way  a  candidate  for y we can easily check if it is
  correct  by  computing  y  A.  If  mu  is  the maximal absolute value of all
  numerators and denominators of the entries of y it is clear from (2) that we
  will find y from y_i if l^i > 2 mu^2.
  
  (4)  If  we  take as v in (1) to(3) all standard unit vectors we clearly get
  the  rows  of  A^-1.  But  we  can do it better. Namely we can take as v the
  standard unit vectors multiplied by the least common multiple epsilon of the
  denominators  of the already computed entries of A^-1. In many examples this
  epsilon  actually  equals  epsilon_r  after  the computation of the first or
  first  few  rows.  Therefore we will often find quickly the next row of A^-1
  already  in  (1),  because  we  find a v_i = 0 such that the sequence of y_i
  becomes constant (=y).
  
