<p><em>(Elementary Divisors and Integer Matrices, by Frank Lübeck)</em></p>
<p>This chapter describes the functions defined in the <strong class="pkg">GAP</strong>4 package <strong class="pkg">EDIM</strong>. The main functions implement variants of an algorithm for computing for a given prime <span class="SimpleMath">p</span> the <span class="SimpleMath">p</span>-parts of the elementary divisors of an integer matrix. These algorithms use a <span class="SimpleMath">p</span>-adic method and are described by the author in <a href="chapBib.html#biBL98">[Lüb02]</a> (see <code class="func">ElementaryDivisorsPPartRk</code> (<a href="chap1.html#X813B0D73868CD751"><span class="RefLink">1.2-1</span></a>)).</p>
<p>These functions were already applied to integer matrices of dimension greater than <span class="SimpleMath">11000</span> (which had many non-trivial elementary divisors which were products of small primes).</p>
<p>Furthermore there are functions for finding the biggest elementary divisor of an invertible integer matrix and the inverse of a rational invertible matrix (see <code class="func">ExponentSquareIntMatFullRank</code> (<a href="chap1.html#X86DA61D978D2D889"><span class="RefLink">1.3-3</span></a>) and <code class="func">InverseRatMat</code> (<a href="chap1.html#X7A9656D47C4D2D16"><span class="RefLink">1.3-1</span></a>)). These algorithms use <span class="SimpleMath">p</span>-adic approximations, explained in <a href="chap1.html#X83FD1AAB7EB2E934"><span class="RefLink">1.7</span></a>.</p>
<p>Finally we distribute implementations of some other algorithms for finding elementary divisors or normal forms of integer matrices: A <span class="SimpleMath">p</span>-modular algorithm by Havas and Sterling from <a href="chapBib.html#biBHS79">[HS79]</a> (see <code class="func">ElementaryDivisorsPPartHavasSterling</code> (<a href="chap1.html#X82EC4724865F4DF9"><span class="RefLink">1.2-2</span></a>)) and LLL-based algorithms for extended greatest common divisors of integers (see <code class="func">GcdexIntLLL</code> (<a href="chap1.html#X799B1A5285D00859"><span class="RefLink">1.5-1</span></a>)) and for Hermite normal forms of integer matrices with (very nice) transforming matrices (see <code class="func">HermiteIntMatLLL</code> (<a href="chap1.html#X7C6AE6777B72F9D2"><span class="RefLink">1.5-2</span></a>)).</p>
<p>Please, send me an e-mail (<span class="URL"><a href="mailto:Frank.Luebeck@Math.RWTH-Aachen.De">Frank.Luebeck@Math.RWTH-Aachen.De</a></span>) if you have any questions, remarks, suggestions, etc. concerning this mini-package. Also, I would like to hear about applications of this package.</p>
<h4>1.1 <span class="Heading">Installation of the <strong class="pkg">EDIM</strong> package</span></h4>
<p>To install this package first unpack it inside some GAP root directory in the subdirectory <code class="file">pkg</code> (see <a href="../../../doc/ref/chap76.html#X82473E4B8756C6CD"><span class="RefLink">Reference: Installing a GAP Package</span></a>). Then the <strong class="pkg">EDIM</strong> package can already be loaded and used. But we strongly recommend to compile a kernel function as well during installation, otherwise the function <code class="func">ElementaryDivisorsPPartRkExpSmall</code> (<a href="chap1.html#X813B0D73868CD751"><span class="RefLink">1.2-1</span></a>) will not be available.</p>
<p>To install the kernel function go to the directory <code class="file">pkg/EDIM-...</code> to which the package was extracted and call</p>
<p>where <code class="code">path</code> is a path to the main <strong class="pkg">GAP</strong> root directory (if not given, the default <code class="file">../..</code> is assumed). Afterwards call <code class="code">make</code> to compile a binary file.</p>
<p>If you have installed several GAP kernels repeat these two steps for each of them. You can run a test of the installation by typing <code class="code">make test</code>.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ InfoEDIM</code></td><td class="tdright">( info class )</td></tr></table></div>
<p>This is an <code class="keyw">Info</code> class for the <strong class="pkg">EDIM</strong>-package. By <code class="code">SetInfoLevel(InfoEDIM, 1);</code> you can switch on the printing of some information during the computations of certain <strong class="pkg">EDIM</strong>-functions.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ ElementaryDivisorsPPartRk</code>( <var class="Arg">A</var>, <var class="Arg">p</var>[, <var class="Arg">rk</var>] )</td><td class="tdright">( function )</td></tr></table></div>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ ElementaryDivisorsPPartRkI</code>( <var class="Arg">A</var>, <var class="Arg">p</var>, <var class="Arg">rk</var> )</td><td class="tdright">( function )</td></tr></table></div>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ ElementaryDivisorsPPartRkII</code>( <var class="Arg">A</var>, <var class="Arg">p</var>, <var class="Arg">rk</var> )</td><td class="tdright">( function )</td></tr></table></div>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ ElementaryDivisorsPPartRkExp</code>( <var class="Arg">A</var>, <var class="Arg">p</var>, <var class="Arg">rk</var>, <var class="Arg">exp</var> )</td><td class="tdright">( function )</td></tr></table></div>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ ElementaryDivisorsPPartRkExpSmall</code>( <var class="Arg">A</var>, <var class="Arg">p</var>, <var class="Arg">rk</var>, <var class="Arg">exp</var>, <var class="Arg">il</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>These functions return a list <span class="SimpleMath">[m_1, m_2, ..., m_r]</span> where <span class="SimpleMath">m_i</span> is the number of nonzero elementary divisors of <var class="Arg">A</var> divisible by <span class="SimpleMath"><var class="Arg">p</var>^i</span> (see <code class="func">ElementaryDivisorsMat</code> (<a href="../../../doc/ref/chap24.html#X7AC4D74F81908109"><span class="RefLink">Reference: ElementaryDivisorsMat</span></a>) for a definition of the elementary divisors).</p>
<p>The algorithms for these functions are described in <a href="chapBib.html#biBL98">[Lüb02]</a>.</p>
<p><var class="Arg">A</var> must be a matrix with integer entries, <var class="Arg">p</var> a prime, and <var class="Arg">rk</var> the rank of <var class="Arg">A</var> (as rational matrix). In the first version of the command <var class="Arg">rk</var> is computed, if it is not given.</p>
<p>The first version of the command delegates its job to the fourth version by trying growing values for <var class="Arg">exp</var>, see below.</p>
<p>The second and third versions implement the main algorithm described in <a href="chapBib.html#biBL98">[Lüb02]</a> and a variation. Here <code class="func">ElementaryDivisorsPPartRkII</code> has a bit more overhead, but can be advantageous because the intermediate entries during the computation can be much smaller.</p>
<p>In the fourth form <var class="Arg">exp</var> must be an upper bound for the highest power of <var class="Arg">p</var> appearing in an elementary divisor of <var class="Arg">A</var>. This information allows reduction of matrix entries modulo <span class="SimpleMath"><var class="Arg">p</var>^{<var class="Arg">exp+1</var>}</span> during the computation.</p>
<p>If <var class="Arg">exp</var> is too small or the given <var class="Arg">rk</var> is too large the function returns <code class="keyw">fail</code>.</p>
<p>If <span class="SimpleMath"><var class="Arg">p</var>^<var class="Arg">exp</var></span> is small enough we use internally a kernel function which can also be used directly in the fifth form of the command. There <var class="Arg">il</var> can be <span class="SimpleMath">0</span> or <span class="SimpleMath">1</span> where in the second case some information is printed during the computation.</p>
<p>More precisely, <code class="func">ElementaryDivisorsPPartRkExpSmall</code> is applicable when <span class="SimpleMath"><var class="Arg">p</var>^{<var class="Arg">exp</var>+1} < 2^{b-4}</span> and <span class="SimpleMath">(<var class="Arg">p</var>^{<var class="Arg">exp</var>+1} - 1)(p-1) < 2^b</span> where <span class="SimpleMath">b=32</span> in a 32-bit version of <strong class="pkg">GAP</strong> and <span class="SimpleMath">b=64</span> in a 64-bit version of <strong class="pkg">GAP</strong>.</p>
<p>This last form of the function was already succesfully applied to dense matrices of rank up to <span class="SimpleMath">11000</span>.</p>
<p>Note that you have to compile a file (see <a href="chap1.html#X84D3F5E77E3BF046"><span class="RefLink">1.1</span></a>) while installing this package, if you want to have this kernel function available.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ ElementaryDivisorsPPartHavasSterling</code>( <var class="Arg">A</var>, <var class="Arg">p</var>, <var class="Arg">d</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>For an integer matrix <var class="Arg">A</var> and a prime <var class="Arg">p</var> this function returns a list <span class="SimpleMath">[m_1, m_2, ..., m_r]</span> where <span class="SimpleMath">m_i</span> is the number of nonzero elementary divisors of <var class="Arg">A</var> divisible by <span class="SimpleMath"><var class="Arg">p</var>^i</span>.</p>
<p>An upper bound <var class="Arg">d</var> for the highest power of <var class="Arg">p</var> appearing in an elementary divisor of <var class="Arg">A</var> must be given. Smaller <var class="Arg">d</var> improve the performance of the algorithm considerably.</p>
<p>This is an implementation of the modular algorithm described in <a href="chapBib.html#biBHS79">[HS79]</a>.</p>
<p>We added a slight improvement: we divide the considered submatrices by the <var class="Arg">p</var>-part of the greatest common divisor of all entries (and lower the <var class="Arg">d</var> appropriately). This reduces the size of the entries and often shortens the pivot search.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ InverseRatMat</code>( <var class="Arg">A</var>[, <var class="Arg">p</var>] )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns the inverse of an invertible matrix over the rational numbers.</p>
<p>It first computes the inverse modulo some prime <var class="Arg">p</var>, computes from this a <var class="Arg">p</var>-adic approximation to the inverse and finally constructs the rational entries from their <var class="Arg">p</var>-adic approximations. See section <a href="chap1.html#X83FD1AAB7EB2E934"><span class="RefLink">1.7</span></a> for more details.</p>
<p>This seems to be better than <strong class="pkg">GAP</strong>'s standard Gauß algorithm (A^-1) already for small matrices. (Try, e.g., RandomMat(20,20,[-10000..10000]) or RandomMat(100,100).)
<p>The optional argument <var class="Arg">p</var> should be a prime such that <var class="Arg">A</var> modulo <var class="Arg">p</var> is invertible (default is <span class="SimpleMath"><var class="Arg">p</var>=251</span>). If <var class="Arg">A</var> is not invertible modulo <var class="Arg">p</var> then <var class="Arg">p</var> is automatically replaced by the next prime.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ RationalSolutionIntMat</code>( <var class="Arg">A</var>, <var class="Arg">v</var>[, <var class="Arg">p</var>[, <var class="Arg">invA</var>]] )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns the solution <span class="SimpleMath">x</span> of the system of linear equations <span class="SimpleMath">x <var class="Arg">A</var> = <var class="Arg">v</var></span>.</p>
<p>Here, <var class="Arg">A</var> must be a matrix with integer entries which is invertible over the rationals and <var class="Arg">v</var> must be a vector with integer entries of the appropriate length.</p>
<p>The optional arguments are a prime <var class="Arg">p</var> such that <span class="SimpleMath"><var class="Arg">A</var> mod p</span> is invertible (if not given, <span class="SimpleMath">p = 251</span> is assumed) and the inverse <var class="Arg">invA</var> of <span class="SimpleMath"><var class="Arg">A</var> mod p</span>.</p>
<p>The solution is computed via <span class="SimpleMath">p</span>-adic approximation as explained in <a href="chap1.html#X83FD1AAB7EB2E934"><span class="RefLink">1.7</span></a>.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ ExponentSquareIntMatFullRank</code>( <var class="Arg">A</var>[, <var class="Arg">p</var>[, <var class="Arg">nr</var>]] )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns the biggest elementary divisor of a square integer matrix <var class="Arg">A</var> of full rank.</p>
<p>For such a matrix <var class="Arg">A</var> the least common multiple of the denominators of all entries of the inverse matrix <span class="SimpleMath"><var class="Arg">A</var>^-1</span> is exactly the biggest elementary divisor of <var class="Arg">A</var>.</p>
<p>This function is implemented by a slight modification of <code class="func">InverseRatMat</code> (<a href="chap1.html#X7A9656D47C4D2D16"><span class="RefLink">1.3-1</span></a>). The third argument <var class="Arg">nr</var> tells the function to return the least common multiple of the first <var class="Arg">nr</var> rows of the rational inverse matrix only. Very often the function will already return the biggest elementary divisor with <span class="SimpleMath"><var class="Arg">nr</var>=2</span> or <span class="SimpleMath">3</span> (and the command without this argument would spend most time in checking, that this is correct).</p>
<p>The optional argument <var class="Arg">p</var> should be a prime such that <var class="Arg">A</var> modulo <var class="Arg">p</var> is invertible (default is <span class="SimpleMath"><var class="Arg">p</var>=251</span>). If <var class="Arg">A</var> is not invertible modulo <var class="Arg">p</var> then <var class="Arg">p</var> is automatically replaced by the next prime.</p>
<div class="example"><pre>
<span class="GAPprompt">gap></span> <span class="GAPinput">ReadPackage("edim", "tst/mat");</span>
Reading 242x242 integer matrix 'mat' with elementary divisors 'eldiv'.
true
<span class="GAPprompt">gap></span> <span class="GAPinput">inv := InverseRatMat(mat);; time; </span>
840
<span class="GAPprompt">gap></span> <span class="GAPinput">ExponentSquareIntMatFullRank(mat, 101, 3); # same as without the `3'
115200
</pre></div>
<h4>1.4 <span class="Heading">All Elementary Divisors Using p-adic Method</span></h4>
<p>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.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ ElementaryDivisorsSquareIntMatFullRank</code>( <var class="Arg">A</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns a list of nonzero elementary divisors of an integer matrix <var class="Arg">A</var>.</p>
<p>Here we start with computing the biggest elementary divisor via <code class="func">ExponentSquareIntMatFullRank</code> (<a href="chap1.html#X86DA61D978D2D889"><span class="RefLink">1.3-3</span></a>). If it runs into a problem because <var class="Arg">A</var> is singular modulo a choosen prime (it starts by default with 251) then the prime is automatically replaced by the next one.</p>
<p>The rest is done using <code class="func">ElementaryDivisorsPPartRkExp</code> (<a href="chap1.html#X813B0D73868CD751"><span class="RefLink">1.2-1</span></a>) and <code class="func">RankMod</code> (<a href="chap1.html#X8312EDA78209B4EA"><span class="RefLink">1.6-5</span></a>).</p>
<p>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.</p>
<p>Note that this function may for many matrices not be the best choice for computing all elementary divisors. You may first try the standard <strong class="pkg">GAP</strong> library routines for Smith normal form instead of this function. Nevertheless remember <code class="func">ElementaryDivisorsSquareIntMatFullRank</code> for hard and big examples. It is particularly good when the largest elementary divisor is a very small factor of the determinant.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ ElementaryDivisorsIntMatDeterminant</code>( <var class="Arg">A</var>, <var class="Arg">det</var>[, <var class="Arg">rk</var>] )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns a list of nonzero elementary divisors of an integer matrix <var class="Arg">A</var>.</p>
<p>Here <var class="Arg">det</var> must be an integer which is a multiple of the biggest determinant divisor of <var class="Arg">A</var>. If the matrix does not have full rank then its rank <var class="Arg">rk</var> must be given, too.</p>
<p>The argument <var class="Arg">det</var> can be given in the form of <code class="code">Collected(FactorsInt(<var class="Arg">det</var>))</code>.</p>
<p>This function handles prime divisors of <var class="Arg">det</var> with multiplicity smaller than 4 specially, for the other prime divisors <span class="SimpleMath">p</span> it delegates to <code class="func">ElementaryDivisorsPPartRkExp</code> (<a href="chap1.html#X813B0D73868CD751"><span class="RefLink">1.2-1</span></a>) where the <var class="Arg">exp</var> argument is the multiplicity of the <span class="SimpleMath">p</span> in <var class="Arg">det</var>. (Note that this is not very good when <span class="SimpleMath">p</span> has actually a much smaller multiplicity in the largest elementary divisor.)</p>
<h4>1.5 <span class="Heading">Gcd and Normal Forms Using LLL</span></h4>
<p>The <strong class="pkg">EDIM</strong>-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 <a href="chapBib.html#biBHMM98">[HMM98]</a> by Havas, Majewski and Matthews.</p>
<p>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).</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ GcdexIntLLL</code>( <var class="Arg">n1</var>, <var class="Arg">n2</var>, <var class="Arg">...</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns for integers <span class="SimpleMath"><var class="Arg">n1</var>, <var class="Arg">n2</var>, ...</span> a list <span class="SimpleMath">[g, [c_1, c_2, ...]]</span>, where <span class="SimpleMath">g = c_1<var class="Arg">n1</var> + c_2<var class="Arg">n2</var> + ...</span> is the greatest common divisor of the <var class="Arg">ni</var>. Here all the <span class="SimpleMath">c_i</span> are usually very small.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ HermiteIntMatLLL</code>( <var class="Arg">A</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>This returns the Hermite normal form of an integer matrix <var class="Arg">A</var> and uses the LLL-algorithm to avoid entry explosion.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ HermiteIntMatLLLTrans</code>( <var class="Arg">A</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns a pair of matrices <span class="SimpleMath">[H, L]</span> where <span class="SimpleMath">H = L <var class="Arg">A</var></span> is the Hermite normal form of an integer matrix <var class="Arg">A</var>. The transforming matrix <span class="SimpleMath">L</span> can have surprisingly small entries.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ SmithIntMatLLL</code>( <var class="Arg">A</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns the Smith normal form of an integer matrix <var class="Arg">A</var> using the LLL-algorithm to avoid entry explosion.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ SmithIntMatLLLTrans</code>( <var class="Arg">A</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns <span class="SimpleMath">[S, L, R]</span> where <span class="SimpleMath">S = L <var class="Arg">A</var> R</span> is the Smith normal form of an integer matrix <var class="Arg">A</var>.</p>
<p>We apply the algorithm for Hermite normal form several times to get the Smith normal form, that is not in the paper <a href="chapBib.html#biBHMM98">[HMM98]</a>. The transforming matrices need not be as nice as for the Hermite normal form.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ RatNumberFromModular</code>( <var class="Arg">n</var>, <var class="Arg">k</var>, <var class="Arg">l</var>, <var class="Arg">x</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns <span class="SimpleMath">r/s = <var class="Arg">x</var> mod <var class="Arg">n</var></span>, if it exists. More precisely:</p>
<p><var class="Arg">n</var>, <var class="Arg">k</var>, <var class="Arg">l</var> must be positive integers with <span class="SimpleMath">2<var class="Arg">k</var><var class="Arg">l</var> ≤ <var class="Arg">n</var></span> and <var class="Arg">x</var> an integer with <span class="SimpleMath">-<var class="Arg">n</var>/2 < <var class="Arg">x</var> ≤ <var class="Arg">n</var>/2</span>. If it exists this function returns a rational number <span class="SimpleMath">r/s</span> with <span class="SimpleMath">0 < s < <var class="Arg">l</var></span>, <span class="SimpleMath">gcd(s, <var class="Arg">n</var>) = 1</span>, <span class="SimpleMath">-<var class="Arg">k</var> < r < <var class="Arg">k</var></span> and <span class="SimpleMath">r/s</span> congruent to <span class="SimpleMath"><var class="Arg">x</var> mod <var class="Arg">n</var></span> (i.e., <span class="SimpleMath"><var class="Arg">n</var> ∣ r - s <var class="Arg">x</var></span>). Such an <span class="SimpleMath">r/s</span> is unique. The function returns <code class="keyw">fail</code> if such a number does not exist.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ InverseIntMatMod</code>( <var class="Arg">A</var>, <var class="Arg">p</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns an inverse matrix modulo a prime <var class="Arg">p</var> or <code class="keyw">fail</code>. More precisely:</p>
<p><var class="Arg">A</var> must be an integer matrix and <var class="Arg">p</var> a prime such that <varclass="Arg">A</var> is invertible modulo <var class="Arg">p</var>. This function returns an integer matrix <var class="Arg">inv</var> with entries in the range <span class="SimpleMath">]-<var class="Arg">p</var>/2 ... <var class="Arg">p</var>/2]</span> such that <var class="Arg">inv</var><var class="Arg">A</var> reduced modulo p is the identity matrix.</p>
<p>It returns <code class="keyw">fail</code> if the inverse modulo <var class="Arg">p</var> does not exist. This function is particularly fast for primes smaller 256.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ HadamardBoundIntMat</code>( <var class="Arg">A</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>The Hadamard bound for a square integer matrix <var class="Arg">A</var> is the product of Euclidean norms of the nonzero rows (or columns) of <var class="Arg">A</var>. It is an upper bound for the absolute value of the determinant of <var class="Arg">A</var>.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ CheapFactorsInt</code>( <var class="Arg">n</var>[, <var class="Arg">nr</var>] )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns a list of factors of an integer <var class="Arg">n</var>, including "small" prime factors - here the optional argument <var class="Arg">nr</var> is the number of iterations for `FactorsRho' (default is 2000).
<p>This is only a slight modification of the library function <code class="func">FactorsInt</code> (<a href="../../../doc/ref/chap14.html#X82C989DB84744B36"><span class="RefLink">Reference: FactorsInt</span></a>) which avoids an error message when the number is not completely factored.</p>
<div class="func"><table class="func" width="100%"><tr><td class="tdleft"><code class="func">‣ RankMod</code>( <var class="Arg">A</var>, <var class="Arg">p</var> )</td><td class="tdright">( function )</td></tr></table></div>
<p>This function returns the rank of an integer matrix <var class="Arg">A</var> modulo <var class="Arg">p</var>. Here <var class="Arg">p</var> must not necessarily be a prime. If it is not and this function returns an integer, then this is the rank of <var class="Arg">A</var> for all prime divisors of <varclass="Arg">p</var>.</p>
<p>If during the computation a factorisation of <var class="Arg">p</var> is found (because some pivot entry has nontrivial greatest common divisor with <var class="Arg">p</var>) then the function is recursively applied to the found factors <code class="code">f_i</code> of <var class="Arg">p</var>. The result is then given in the form <code class="code">[[f_1, rk_1], [f_2, rk_2], ...]</code>.</p>
<p>The idea to make this function useful for non primes was to use it with large factors of the biggest elementary divisor of <var class="Arg">A</var> whose prime factorization cannot be found easily.</p>
<h4>1.7 <span class="Heading">InverseRatMat - the Algorithm</span></h4>
<p>The idea is to recover a rational matrix from an <span class="SimpleMath">l</span>-adic approximation for some prime <span class="SimpleMath">l</span>. 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 <a href="chapBib.html#biBD82">[Dix82]</a> of Dixon.</p>
<p>Let <span class="SimpleMath">A</span> be an invertible matrix over the rational numbers. By multiplying with a constant we may assume that its entries are in fact integers.</p>
<p>(1) We first describe how to find an <span class="SimpleMath">l</span>-adic approximation of <span class="SimpleMath">A^-1</span>. Find a prime <span class="SimpleMath">l</span> such that <span class="SimpleMath">A</span> is invertible modulo <span class="SimpleMath">l</span> and let <span class="SimpleMath">B</span> be the integer matrix with entries in the range <span class="SimpleMath">]-l/2,l/2]</span> such that <span class="SimpleMath">BA</span> is congruent to the identity matrix modulo <span class="SimpleMath">l</span>. (This can be computed fast by usual Gauß elimination.)</p>
<p>Now let <span class="SimpleMath">v ∈ ℤ^r</span> be a row vector. Define two sequences <span class="SimpleMath">v_i</span> and <span class="SimpleMath">x_i</span> of row vectors in <span class="SimpleMath">ℤ^r</span> by: <span class="SimpleMath">x_0 := 0 ∈ ℤ^r</span>, <span class="SimpleMath">v_0 := -v</span> and for <span class="SimpleMath">i > 0</span> set <span class="SimpleMath">x_i</span> to the vector congruent to <span class="SimpleMath">-v_i-1 B</span> modulo <span class="SimpleMath">l</span> having entries in the range <span class="SimpleMath">]-l/2, l/2]</span>. Then all entries of <span class="SimpleMath">x_i A + v_i-1</span> are divisible by <span class="SimpleMath">l</span> and we set <span class="SimpleMath">v_i := (1/l) ⋅ (x_i A + v_i-1)</span>.</p>
<p>Induction shows that for <span class="SimpleMath">y_i := ∑_k=1^i l^k-1 x_k</span> we have <span class="SimpleMath">y_i A = v + l^i v_i</span> for all <span class="SimpleMath">i ≥ 0</span>. Hence the sequence <span class="SimpleMath">y_i</span>, <span class="SimpleMath">i ≥ 0</span>, gives an <span class="SimpleMath">l</span>-adic approximation to the vector <span class="SimpleMath">y ∈ ℚ^r</span> with <span class="SimpleMath">y A = v</span>.</p>
<p>(2) The second point is to show how we can get the vector <span class="SimpleMath">y</span> from a sufficiently good approximation <span class="SimpleMath">y_i</span>. Note that the sequence of <span class="SimpleMath">y_i</span> becomes constant for <span class="SimpleMath">i ≥ i_0</span> if all entries of <span class="SimpleMath">y</span> are integers of absolute value smaller than <span class="SimpleMath">l^i_0 / 2</span> because of our choice of representatives of residue classes modulo <span class="SimpleMath">l</span> in the interval <span class="SimpleMath">]-l/2, l/2]</span>.</p>
<p>More generally consider <span class="SimpleMath">a / b ∈ ℚ</span> with <span class="SimpleMath">b > 0</span> and <span class="SimpleMath">a, b</span> coprime. Then there is for each <span class="SimpleMath">n ∈ ℕ</span> which is coprime to <span class="SimpleMath">b</span> a unique <span class="SimpleMath">c ∈ ℤ</span> with <span class="SimpleMath">-n / 2 < c ≤ n / 2</span> and <span class="SimpleMath">a ≡ c b mod n</span>. This <span class="SimpleMath">c</span> can be computed via the extended Euclidean algorithm applied to <span class="SimpleMath">b</span> and <span class="SimpleMath">n</span>.</p>
<p>Now let <span class="SimpleMath">n, α, β ∈ ℕ</span> with <span class="SimpleMath">2 α β ≤ n</span>. Then the map <span class="SimpleMath">{a/b ∈ ℚ ∣ -α ≤ a ≤ α, 1 ≤ b < β } → ]-n/2, n/2]</span>, <span class="SimpleMath">a/b ↦ c</span> (defined as above) is injective (since for <span class="SimpleMath">a/b</span>, <span class="SimpleMath">a'/b'</span> in the above set we have <span class="SimpleMath">a b' - a' b ≡ 0 mod n</span> if and only if <span class="SimpleMath">a b' - a' b = 0</span>).</p>
<p>In practice we can use for any <span class="SimpleMath">c ∈ ]-n/2, n/2]</span> a certain extended Euclidean algorithm applied to <span class="SimpleMath">n</span> and <span class="SimpleMath">c</span> to decide if <span class="SimpleMath">c</span> is in the image of the above map and to find the corresponding <span class="SimpleMath">a/b</span> if it exists.</p>
<p>(3) To put things together we apply (2) to the entries of the vectors <span class="SimpleMath">y_i</span> constructed in (1), choosing <span class="SimpleMath">n = l^i</span>, <span class="SimpleMath">α = sqrtn/2</span> and <span class="SimpleMath">β = sqrtn</span>. If we have found this way a candidate for <span class="SimpleMath">y</span> we can easily check if it is correct by computing <span class="SimpleMath">y A</span>. If <span class="SimpleMath">μ</span> is the maximal absolute value of all numerators and denominators of the entries of <span class="SimpleMath">y</span> it is clear from (2) that we will find <span class="SimpleMath">y</span> from <span class="SimpleMath">y_i</span> if <span class="SimpleMath">l^i > 2 μ^2</span>.</p>
<p>(4) If we take as <span class="SimpleMath">v</span> in (1) to(3) all standard unit vectors we clearly get the rows of <span class="SimpleMath">A^-1</span>. But we can do it better. Namely we can take as <span class="SimpleMath">v</span> the standard unit vectors multiplied by the least common multiple <span class="SimpleMath">ϵ</span> of the denominators of the already computed entries of <span class="SimpleMath">A^-1</span>. In many examples this <span class="SimpleMath">ϵ</span> actually equals <span class="SimpleMath">ϵ_r</span> after the computation of the first or first few rows. Therefore we will often find quickly the next row of <span class="SimpleMath">A^-1</span> already in (1), because we find a <span class="SimpleMath">v_i = 0</span> such that the sequence of <span class="SimpleMath">y_i</span> becomes constant (<span class="SimpleMath">=y</span>).</p>
<h5>1.7-1 <span class="Heading">Rank of Integer Matrix</span></h5>
<p>The following strategy has shown to be useful in proving that some very big integer matrix is not invertible.</p>
<ul>
<li><p>Check the rank modulo some small primes, say with <code class="func">RankMod</code> (<a href="chap1.html#X8312EDA78209B4EA"><span class="RefLink">1.6-5</span></a>).</p>
</li>
<li><p>If the rank seems less than the number of rows choose a prime <span class="SimpleMath">p</span>, a collection of lines which is linearly independent modulo <span class="SimpleMath">p</span>, and another line linearly dependend on these. Guess that this last line is also linearly dependend on the chosen collection over the rational numbers (maybe check modulo several small primes).</p>
</li>
<li><p>Find columns of the collection of lines which give an invertible matrix modulo some prime.</p>
</li>
<li><p>Then use <code class="func">RationalSolutionIntMat</code> (<a href="chap1.html#X8302B31E86B3AFDB"><span class="RefLink">1.3-2</span></a>) with the invertible submatrix and corresponding entries of the linearly dependend row to prove this.</p>
</li>
</ul>
<p>Guessing the rank of a matrix from the rank modulo several primes, chosing a maximal set of lines which are linearly independent modulo some primes, and using <code class="func">RationalSolutionIntMat</code> (<a href="chap1.html#X8302B31E86B3AFDB"><span class="RefLink">1.3-2</span></a>) with the remaining lines, one may also find the exact rank of a huge integer matrix.</p>
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
