Make rational matrix rref default to flint_multimodular, add suboptions for flint algorithm

src/sage/matrix/matrix_rational_dense.pyx

@@ -1499,24 +1499,33 @@ cdef class Matrix_rational_dense(Matrix_dense):
 
         - ``algorithm`` -- an optional specification of an algorithm. One of
 
-          - ``None``: (default) uses flint for small dimension and multimodular otherwise
+          - ``None``: (default) try to pick the best choice,
 
-          - ``'flint'``: use the flint library,
+          - ``'flint'``: use flint library
+            `function <https://flintlib.org/doc/fmpq_mat.html#c.fmpq_mat_rref>`_,
+            which automatically chooses between
+            `classical algorithm <https://flintlib.org/doc/fmpq_mat.html#c.fmpq_mat_rref_classical>`_
+            (Gaussian elimination),
+            `fraction-free multimodular <https://flintlib.org/doc/fmpz_mat.html#c.fmpz_mat_rref_mul>`_,
+            and `fraction-free LU decomposition <https://flintlib.org/doc/fmpz_mat.html#c.fmpz_mat_rref_fflu>`_,
+
+          - ``'flint:classical'``, ``'flint:multimodular'``, ``'flint:fflu'``: use the
+            flint library as above, but select an algorithm explicitly,
 
           - ``'padic'``: an algorithm based on the IML `p`-adic solver,
 
-          - ``'multimodular'``: uses a multimodular algorithm the uses
-            linbox modulo many primes (likely to be faster when coefficients
-            are huge),
+          - ``'multimodular'``: uses a multimodular algorithm implemented in Cython
+            that uses linbox modulo many primes,
+            see :func:`~sage.matrix.misc.matrix_rational_echelon_form_multimodular`,
 
           - ``'classical'``: just clear each column using Gauss elimination.
 
         - ``height_guess``, ``**kwds`` -- all passed to the
-          multimodular algorithm; ignored by other algorithms
+          ``'multimodular'`` algorithm; ignored by other algorithms
 
         - ``proof`` -- boolean or ``None`` (default: None, see
           proof.linear_algebra or sage.structure.proof). Passed to the
-          multimodular algorithm. Note that the Sage global default is
+          ``'multimodular'`` algorithm. Note that the Sage global default is
           ``proof=True``.
 
         EXAMPLES::
@@ -1548,7 +1557,8 @@ cdef class Matrix_rational_dense(Matrix_dense):
         Echelonizing a matrix in place throws away the cache of
         the old matrix (:issue:`14506`)::
 
-            sage: for algo in ["flint", "padic", "multimodular", "classical"]:
+            sage: for algo in ["flint", "padic", "multimodular", "classical", "flint:classical",
+            ....:              "flint:multimodular", "flint:fflu"]:
             ....:      a = Matrix(QQ, [[1,2],[3,4]])
             ....:      _ = a.det()          # fills the cache
             ....:      _ = a._clear_denom() # fills the cache
@@ -1560,13 +1570,10 @@ cdef class Matrix_rational_dense(Matrix_dense):
         self.check_mutability()
 
         if algorithm is None:
-            if self._nrows <= 25 or self._ncols <= 25:
-                algorithm = 'flint'
-            else:
-                algorithm = 'multimodular'
+            algorithm = 'flint:multimodular'
 
-        if algorithm == 'flint':
-            pivots = self._echelonize_flint()
+        if algorithm in ('flint', 'flint:classical', 'flint:multimodular', 'flint:fflu'):
+            pivots = self._echelonize_flint(algorithm)
         elif algorithm == 'multimodular':
             pivots = self._echelonize_multimodular(height_guess, proof, **kwds)
         elif algorithm == 'classical':
@@ -1656,20 +1663,24 @@ cdef class Matrix_rational_dense(Matrix_dense):
         self.cache('rank', len(E.pivots()))
         return E
 
-    def _echelonize_flint(self):
+    def _echelonize_flint(self, algorithm: str):
         r"""
+        INPUT: See :meth:`echelonize` for the options.
+        Only options that use flint are allowed, passing other algorithms may
+        trigger undefined behavior.
+
         EXAMPLES::
 
             sage: m = matrix(QQ, 4, range(16))
-            sage: m._echelonize_flint()
+            sage: m._echelonize_flint("flint")
             (0, 1)
             sage: m
             [ 1  0 -1 -2]
             [ 0  1  2  3]
             [ 0  0  0  0]
             [ 0  0  0  0]
             sage: m = matrix(QQ, 4, 6, [-1,0,0,-2,-1,-2,-1,0,0,-2,-1,0,3,3,-2,0,0,3,-2,-3,1,1,-2,3])
-            sage: m._echelonize_flint()
+            sage: m._echelonize_flint("flint")
             (0, 1, 2, 5)
             sage: m
             [   1    0    0    2    1    0]
@@ -1678,10 +1689,40 @@ cdef class Matrix_rational_dense(Matrix_dense):
             [   0    0    0    0    0    1]
         """
         self.clear_cache()
+
+        if fmpq_mat_is_empty(self._matrix):
+            return ()
+
         cdef long r
+        cdef fmpz_mat_t Aclear
+        cdef fmpz_t den
 
         sig_on()
-        r = fmpq_mat_rref(self._matrix, self._matrix)
+
+        if algorithm == 'flint':
+            r = fmpq_mat_rref(self._matrix, self._matrix)
+        elif algorithm == 'flint:classical':
+            r = fmpq_mat_rref_classical(self._matrix, self._matrix)
+        else:
+            # copied from fmpq_mat_rref_fraction_free
+            fmpz_mat_init(Aclear, self._nrows, self._ncols)
+            fmpq_mat_get_fmpz_mat_rowwise(Aclear, NULL, self._matrix)
+            fmpz_init(den)
+
+            if algorithm == 'flint:fflu':
+                r = fmpz_mat_rref_fflu(Aclear, den, Aclear)
+            else:
+                assert algorithm == 'flint:multimodular'
+                r = fmpz_mat_rref_mul(Aclear, den, Aclear)
+
+            if r == 0:
+                fmpq_mat_zero(self._matrix)
+            else:
+                fmpq_mat_set_fmpz_mat_div_fmpz(self._matrix, Aclear, den)
+
+            fmpz_mat_clear(Aclear)
+            fmpz_clear(den)
+
         sig_off()
 
         # compute pivots

src/sage/matrix/misc.pyx

@@ -204,6 +204,49 @@ def matrix_rational_echelon_form_multimodular(Matrix self, height_guess=None, pr
         [                0                 0                 1                 2]
         sage: A.pivots()
         (0, 1, 2)
+
+    A small benchmark, showing that flint fraction-free multimodular algorithm
+    is always faster than the fraction-free multimodular algorithm implemented in Python::
+
+        sage: import copy
+        sage: def benchmark(num_row, num_col, entry_size, timeout=2, integer_coefficient=True):
+        ....:     A = matrix(QQ, [[
+        ....:         randint(1, 2^entry_size) if integer_coefficient else ZZ(randint(1, 2^entry_size))/randint(1, 2^entry_size)
+        ....:         for col in range(num_col)] for row in range(num_row)])
+        ....:     data=[]
+        ....:     for algorithm in ("flint:fflu", "flint:multimodular", "padic", "multimodular"):
+        ....:         # classical is too slow
+        ....:         B = copy.copy(A)
+        ....:         t = walltime()
+        ....:         alarm(timeout)
+        ....:         try:
+        ....:             B.echelonize(algorithm=algorithm)
+        ....:         except AlarmInterrupt:
+        ....:             pass
+        ....:         finally:
+        ....:             cancel_alarm()
+        ....:         data.append((round(walltime(t), 4), algorithm))
+        ....:     return sorted(data)
+        sage: benchmark(20, 20, 10000)  # long time
+        [...'flint:multimodular'...'multimodular'...'flint:fflu'...]
+        sage: benchmark(39, 40, 200)  # long time
+        [...'flint:multimodular'...'flint:fflu'...'multimodular'...]
+
+    In older versions of flint
+    before this `issue <https://github.com/flintlib/flint/issues/2129>`_
+    is fixed, ``algorithm='flint'`` (automatic choice) may be slower than
+    ``algorithm='flint:multimodular'``.
+
+    In this case, there are more columns than rows, which means the resulting
+    matrix has height much higher than the input matrix. We check that the function
+    does not take too long::
+
+        sage: A = matrix(QQ, [[randint(1, 2^500) for col in range(40)] for row in range(20)])
+        sage: t = walltime()
+        sage: A.echelonize(algorithm="multimodular")  # long time
+        sage: t = walltime(t)  # long time
+        sage: (t < 10, t)  # long time
+        (True, ...)
     """
     if proof is None:
         from sage.structure.proof.proof import get_flag
@@ -221,10 +264,12 @@ def matrix_rational_echelon_form_multimodular(Matrix self, height_guess=None, pr
         height_guess = 10000000*(height+100)
     tm = verbose("height_guess = %s" % height_guess, level=2, caller_name="multimod echelon")
 
+    cdef Integer M
+    from sage.arith.misc import integer_floor as floor
     if proof:
-        M = self._ncols * height_guess * height + 1
+        M = floor(max(1, self._ncols * height_guess * height + 1))
     else:
-        M = height_guess + 1
+        M = floor(max(1, height_guess + 1))
 
     if self.is_sparse():
         from sage.matrix.matrix_modn_sparse import MAX_MODULUS
@@ -309,7 +354,7 @@ def matrix_rational_echelon_form_multimodular(Matrix self, height_guess=None, pr
         except ValueError as msg:
             verbose(msg, level=2)
             verbose("Not enough primes to do CRT lift; redoing with several more primes.", level=2, caller_name="multimod echelon")
-            M = prod * p*p*p
+            M <<= M.bit_length() // 5 + 1
             continue
 
         if not proof:
@@ -321,7 +366,7 @@ def matrix_rational_echelon_form_multimodular(Matrix self, height_guess=None, pr
             verbose("Validity of result checked.", level=2, caller_name="multimod echelon")
             break
         verbose("Validity failed; trying again with more primes.", level=2, caller_name="multimod echelon")
-        M = prod * p*p*p
+        M <<= M.bit_length() // 5 + 1
     #end while
     verbose("total time",tm, level=2, caller_name="multimod echelon")
     return E, tuple(best_pivots)