gh-39733: Make rational matrix rref default to flint_multimodular, add suboptions for flint algorithm
    
Because one of the algorithms used by flint is multimodular, it ought to
be faster than the implementation in Python.

At least after we upgrade to a version after
flintlib/flint#2129 .

(p/s: if someone uses the old version, the current choice of flint might
be slower in some cases, see the linked issue. An alternative which is
likely always faster is to explicitly use the multimodular algorithm in
flint. Do you think the current implementation is fine, or should we
provide an explicit `flint_multimodular` option instead?)

Fixes #39197

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [x] The title is concise and informative.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [ ] I have created tests covering the changes.
- [ ] I have updated the documentation and checked the documentation
preview.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - #12345: short description why this is a dependency -->
<!-- - #34567: ... -->


#39204
    
URL: #39733
Reported by: user202729
Reviewer(s): Travis Scrimshaw

Author: Release Manager

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

@@ -205,16 +205,16 @@ def matrix_rational_echelon_form_multimodular(Matrix self, height_guess=None, pr
         sage: A.pivots()
         (0, 1, 2)
 
-    A small benchmark, showing that the multimodular algorithm is sometimes faster
-    and sometimes slower than the flint algorithm::
+    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", "padic", "multimodular"):
+        ....:     for algorithm in ("flint:fflu", "flint:multimodular", "padic", "multimodular"):
         ....:         # classical is too slow
         ....:         B = copy.copy(A)
         ....:         t = walltime()
@@ -227,10 +227,15 @@ def matrix_rational_echelon_form_multimodular(Matrix self, height_guess=None, pr
         ....:             cancel_alarm()
         ....:         data.append((round(walltime(t), 4), algorithm))
         ....:     return sorted(data)
-        sage: benchmark(20, 20, 10000)  # long time (multimodular wins)
-        [...'multimodular'...'flint'...]
-        sage: benchmark(39, 40, 200)  # long time (flint wins)
-        [...'flint'...'multimodular'...]
+        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