Allow singular matrices in saturate (#1102)

* Allow singular matrices in `saturate`

* Evade computing the hnf transform

* `solve` solves it

src/Misc/Matrix.jl

@@ -115,30 +115,48 @@ end
     saturate(A::ZZMatrix) -> ZZMatrix
 
 Computes the saturation of $A$, that is, a basis for $\mathbf{Q}\otimes M \cap
-\mathbf{Z}^n$, where $M$ is the row span of $A$ and $n$ the number of rows of
-$A$.
+\mathbf{Z}^n$, where $M$ is the row span of $A$ and $n$ the number of columns
+of $A$.
 
-Equivalently, return $TA$ for an invertible rational matrix $T$ such that $TA$
-is integral and the elementary divisors of $TA$ are all trivial.
+Equivalently, return $TA$ (minus zero rows) for an invertible rational matrix
+$T$ such that $TA$ is integral and the elementary divisors of $TA$ are all
+trivial.
+
+# Examples
+
+```jldoctest
+julia> saturate(ZZ[1 2 3;4 5 6;7 0 7])
+[1    2    3]
+[1    1    1]
+[0   -1   -1]
+
+julia> saturate(ZZ[1 2 3;4 5 6;7 8 9])
+[1   2   3]
+[1   1   1]
+
+julia> saturate(ZZ[1 2 3;4 5 6])
+[1   2   3]
+[1   1   1]
+
+julia> saturate(ZZ[1 2;3 4;5 6])
+[1   2]
+[1   1]
+
+```
 """
-function saturate(A::ZZMatrix)
-  #row saturation: want
-  #  TA in Z, T in Q and elem_div TA = [1]
-  #
-  #  AT = H (in HNF), then A = HT^-1 and H^-1A = T^-1
-  # since T is uni-mod, H^-1 A is in Z with triv. elm. div
-  H = transpose(A)
-  H = hnf!(H)
-  H = sub(H, 1:ncols(H), 1:ncols(H))
-    ccall((:fmpz_mat_transpose, libflint), Nothing,
-           (Ref{ZZMatrix}, Ref{ZZMatrix}), H, H)
-  Hi, d = pseudo_inv(H)
-  S = Hi*A
-  divexact!(S, S, d)
-  #@hassert :HNF 1  d*Sd == S
+function saturate(A::ZZMatrix) :: ZZMatrix
+  # Let AU = [H 0matrix] in HNF and HS = A = [H 0matrix]U⁻¹
+  # We have S == U⁻¹[1:rank(H), :] in ZZ with trivial elementary divisors.
+  # For any invertible H' with H'H = 1, also S = H'HS = H'A.
+  H = hnf!(transpose(A))
+  H = transpose(H[1:rank(H), :])
+  S = solve(H, A)
+  @assert rank(S) == rank(H)
   return S
 end
 
+transpose!(A::Union{ZZMatrix, QQMatrix}) = transpose!(A, A)
+
 function transpose!(A::ZZMatrix, B::ZZMatrix)
   ccall((:fmpz_mat_transpose, libflint), Nothing,
     (Ref{ZZMatrix}, Ref{ZZMatrix}), A, B)