Exact small-entry matrix algebra in quadratic logarithmic workspace

Every $d\times d$ minor, $0\le d\le n$, has absolute determinant at most $d!\le (n+2{)}^n$ by [F1], interpreting the empty determinant as $1$. Put $B=(n+2{)}^n$, $M=16(n+2{)}^4$, $L_M=\mathrm{lcm}(1,\dots ,M)$ and $P_M={\prod }_{p\le M,p\text{ prime}}p$. All these are integers; only $B$ and the small bound $M$, not the complete integers $L_M,P_M$, will be stored.

The bound $L_M\ge 2^{M-1}\mathrm{/}M$ is elementary. For $N=M-1$ and any $0\le j\le N$, counting multiples of $p,p^2,\dots$ in factorials gives $v_p\left(\genfrac{}{}{0px}{}{N}{j}\right)={\sum }_{a\ge 1}(\lfloor \frac{N}{p^a}\rfloor -\lfloor \frac{j}{p^a}\rfloor -\lfloor \frac{N-j}{p^a}\rfloor )\mathrm{.}$ Each summand is zero or one and vanishes for $p^a>N$. Thus this valuation is at most $\lfloor {\log }_{p}N\rfloor$, at most the exponent of $p$ in $L_M$. Unique factorization makes every $\left(\genfrac{}{}{0px}{}{N}{j}\right)$ divide $L_M$. The sum of the $M$ binomial coefficients is $2^{M-1}$, so at least one is at least $2^{M-1}\mathrm{/}M$, proving the bound.

The exponent of every prime in $L_M$ is at most ${\log }_{2}M$. Consequently $L_M\le P_M^{{\log }_{2}M}$ and ${\log }_2{P}_M\ge \frac{M-1-{\log }_{2}M}{{\log }_{2}M}\mathrm{.}$ For the chosen $M$, ${\log }_{2}M=4+4{\log }_2(n+2)\le 4(n+2)$ and $M-1-{\log }_{2}M\ge M\mathrm{/}2$. The displayed lower bound is therefore at least $2(n+2{)}^3$, which is greater than $n{\log }_2(n+2)+3$. Hence $P_M>8B$. Enumerate integers $2,\dots ,M$, testing each for primality by trial division up to its square root. There are polynomially many trials on $O(\log (n+2))$-bit integers. Stop the prime stream as soon as its running product $Q$ exceeds $2B$; the bound just proved guarantees that this happens by $M$. Immediately before the last multiplication $Q\le 2B$, so the final product is at most $2BM$. The product accumulator thus has $O(n\log (n+2))$ bits.

For a specified square submatrix $D$, compute its determinant modulo each successive prime by [F2]. Every nonzero residue $h$ modulo a prime $p$ is a unit, with an explicit small-space inverse: apply Euclidean descent to $(p,h)$, maintaining $r_i=s_{i}p+t_{i}h$ with initial coefficient pairs $(1,0),(0,1)$ and updates $(s_{i+1},t_{i+1})=(s_{i-1},t_{i-1})-q_i(s_i,t_i)$. The last nonzero remainder is $1$, since a common divisor of $p$ and $0<h<p$ is $1$ by primality. Its coefficient $t_i$ is an inverse modulo $p$. This also constructs Bézout identities for arbitrary integer pairs, and induction gives one for any finite list and its gcd. Thus elimination is over a field; if a pivot cannot be found in a square trailing block, the determinant residue is zero. Row interchanges change its sign, and otherwise the product of the unscaled pivots gives the determinant. This uses one matrix of $O(n^2)$ residues, each with $O(\log (n+2))$ bits, and $O(n^3)$ modular operations per prime. Combine each residue $r$ with the previous representative $a$ modulo $Q$ by $a_{\mathrm{n}\mathrm{e}\mathrm{w}}=a+Q((r-a)Q^{-1}<mrow><mi\; mathvariant="normal">m</mi><mi\; mathvariant="normal">o</mi><mi\; mathvariant="normal">d</mi></mrow>p),Q_{\mathrm{n}\mathrm{e}\mathrm{w}}=Qp\mathrm{.}$ The inverse exists by [F3] and is computed from $Q<mrow><mi\; mathvariant="normal">m</mi><mi\; mathvariant="normal">o</mi><mi\; mathvariant="normal">d</mi></mrow>p$, so its Euclidean inputs have only $O(\log (n+2))$ bits. After $Q>2B$, the unique representative in $[-B,B]$ is the integer determinant by step 1.1. Only the current residue matrix, the two reconstruction integers and small counters are retained; no list of primes, residues or previous matrices is stored.

For rank, perform prime-field elimination on the original rectangular matrix over the same prime stream, retaining the largest rank found and the original row indices of a basis attaining it. These indices can be selected during elimination without a transformation matrix: after the first $k$ pivots, the unselected rows have only had combinations of the $k$ selected original rows subtracted, so the next pivot row is independent of them exactly when its remainder is nonzero. The rank equals the largest order of a nonzero minor over either field. To see the needed direction directly, elimination selects $r$ original independent rows and $r$ pivot columns; their $r\times r$ submatrix becomes triangular with nonzero pivots after invertible row operations, so it has nonzero determinant. Conversely, $r+1$ rows are dependent and multilinearity makes every minor of that order zero. A nonzero modular minor is a nonzero integer minor, so every modular rank is at most the rational rank $r$. If $r>0$, choose any nonzero $r\times r$ integer minor $d$. If all sampled ranks were less than $r$, every sampled prime would divide $d$. Their product $Q>2B$ would divide $d$, contradicting $0<\mathrm{\mid }d\mathrm{\mid }\le B$. Therefore the maximum modular rank is exactly $r$, and its retained original rows form a rational row basis. For $r=0$ every entry and modular rank is zero.

If the number of columns is $b$ and the rank is $b-1$, let $R$ consist of the $b-1$ independent original rows from step 4.2. For $1\le j\le b$, compute $w_j=(-1{)}^{j+1}\det R_{\widehat{j}}$ by step 4.1, retaining the vector $w$. Some cofactor is nonzero because $R$ has full row rank. Expanding the determinant of the matrix obtained by adding any row of $R$ as a duplicate row gives $Rw=0$ by [F1]. Every row of $A$ lies in the rational row span of $R$, hence $Aw=0$. Divide $w$ by the positive gcd of its entries; the result is primitive. If a rational multiple $cu$ of a primitive integer vector $u$ is integral, a Bézout identity ${\sum }_j{a}_j{u}_j=1$ gives $c={\sum }_j{a}_j(c{u}_j)\in \mathbb{Z}$, so it generates all integer points on that line. Either overall sign may be selected. For $b=1$, $R$ is empty and its one empty cofactor is $1$, as required.

All counters bounded by $M$ have $O(\log (n+2))$ bits. The residue matrix uses $O(n^2\log (n+2))$ bits, the row-index list $O(n\log (n+2))$ bits, and each reconstructed determinant has $O(n\log (n+2))$ bits. The at most $n$ retained cofactors therefore also use $O(n^2\log (n+2))$ bits; the same residue matrix storage is reused for each one. Trial division, at most $M$ prime eliminations per determinant and at most $n$ cofactor determinants take polynomially many arithmetic operations on polynomial-length integers. Binary multiplication and long division take polynomial time; the Euclidean descent takes $O(\mathrm{\ell })$ divisions on $\mathrm{\ell }$-bit inputs because after two nonzero remainder steps the remainder is at most half its previous value. The coefficient-tracking algorithm is needed computationally only for inverses with inputs less than $p\le M$. It has $O(\log p)$ updates and quotients at most $p$, so even the crude coefficient bound obtained by multiplying $2p$ per update has $O((\log p{)}^2)$ bits. Thus these small inverses need only $O((\log (n+2){)}^2)$ extra workspace, which fits the asserted bound. Normalizing the cofactor vector needs only ordinary remainder gcds on $O(n\log (n+2))$-bit integers; the existential Bézout identity used to prove primitivity need not be computed or stored. All integer multiplications and divisions can be done by the digit-by-digit algorithms using a constant number of buffers of the operand/output lengths. Thus reconstruction, inversions, gcd normalization, rank and kernel calculations all have polynomial bit time and the asserted workspace. Empty matrices are handled directly by rank zero and empty determinant one. This proves all claims.