gh-40856: Workaround for linbox charpoly/minpoly issues
    
As in the title. Fix #37068. Can
be removed once the upstream issue https://github.com/linbox-
team/linbox/issues/53 is fixed.

It's possible to do a more expensive check (e.g. check if the matrix
satisfy the minpoly with something like Freivalds' algorithm), but I
guess this suffice (and it's very cheap).

### 📝 Checklist

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

- [ ] The title is concise and informative.
- [ ] The description explains in detail what this PR is about.
- [ ] 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: ... -->
    
URL: #40856
Reported by: user202729
Reviewer(s):

Author: Release Manager

src/sage/matrix/matrix_integer_dense.pyx

@@ -1411,10 +1411,13 @@ cdef class Matrix_integer_dense(Matrix_dense):
             fmpz_mat_charpoly(g._poly, self._matrix)
             sig_off()
         elif algorithm == 'linbox':
-            g = (<Polynomial_integer_dense_flint> PolynomialRing(ZZ, names=var).gen())._new()
-            sig_on()
-            linbox_fmpz_mat_charpoly(g._poly, self._matrix)
-            sig_off()
+            while True:  # linbox is unreliable, see :issue:`37068`
+                g = (<Polynomial_integer_dense_flint> PolynomialRing(ZZ, names=var).gen())._new()
+                sig_on()
+                linbox_fmpz_mat_charpoly(g._poly, self._matrix)
+                sig_off()
+                if g.lc() == 1:
+                    break
         elif algorithm == 'generic':
             g = Matrix_dense.charpoly(self, var)
         else:
@@ -1494,10 +1497,13 @@ cdef class Matrix_integer_dense(Matrix_dense):
             return g.change_variable_name(var)
 
         if algorithm == 'linbox':
-            g = (<Polynomial_integer_dense_flint> PolynomialRing(ZZ, names=var).gen())._new()
-            sig_on()
-            linbox_fmpz_mat_minpoly(g._poly, self._matrix)
-            sig_off()
+            while True:  # linbox is unreliable, see :issue:`37068`
+                g = (<Polynomial_integer_dense_flint> PolynomialRing(ZZ, names=var).gen())._new()
+                sig_on()
+                linbox_fmpz_mat_minpoly(g._poly, self._matrix)
+                sig_off()
+                if g.lc() == 1:
+                    break
         elif algorithm == 'generic':
             g = Matrix_dense.minpoly(self, var)
         else:

src/sage/matrix/matrix_integer_sparse.pyx

@@ -901,16 +901,19 @@ cdef class Matrix_integer_sparse(Matrix_sparse):
         cdef linbox.DensePolynomial_integer * p = new linbox.DensePolynomial_integer(givZZ, <size_t> self._nrows)
         cdef Polynomial_integer_dense_flint g = (<Polynomial_integer_dense_flint> R.gen())._new()
 
-        sig_on()
-        linbox.charpoly(p[0], M[0])
-        sig_off()
+        while True:  # linbox is unreliable, see :issue:`37068`
+            sig_on()
+            linbox.charpoly(p[0], M[0])
+            sig_off()
 
-        cdef size_t i
-        fmpz_poly_fit_length(g._poly, p.size())
-        for i in range(p.size()):
-            tmp = p[0][i].get_mpz_const()
-            fmpz_poly_set_coeff_mpz(g._poly, i, tmp)
-        _fmpz_poly_set_length(g._poly, p.size())
+            fmpz_poly_fit_length(g._poly, p.size())
+            for i in range(p.size()):
+                tmp = p[0][i].get_mpz_const()
+                fmpz_poly_set_coeff_mpz(g._poly, i, tmp)
+            _fmpz_poly_set_length(g._poly, p.size())
+
+            if g.lc() == 1:
+                break
 
         del M
         del p
@@ -1000,18 +1003,21 @@ cdef class Matrix_integer_sparse(Matrix_sparse):
         cdef linbox.SparseMatrix_integer * M = new_linbox_matrix_integer_sparse(givZZ, self)
         cdef linbox.DensePolynomial_integer * p = new linbox.DensePolynomial_integer(givZZ, <size_t> self._nrows)
         cdef Polynomial_integer_dense_flint g = (<Polynomial_integer_dense_flint> R.gen())._new()
+        cdef mpz_t tmp
 
-        sig_on()
-        linbox.minpoly(p[0], M[0])
-        sig_off()
+        while True:  # linbox is unreliable, see :issue:`37068`
+            sig_on()
+            linbox.minpoly(p[0], M[0])
+            sig_off()
 
-        cdef size_t i
-        cdef mpz_t tmp
-        fmpz_poly_fit_length(g._poly, p.size())
-        for i in range(p.size()):
-            tmp = p[0][i].get_mpz_const()
-            fmpz_poly_set_coeff_mpz(g._poly, i, tmp)
-        _fmpz_poly_set_length(g._poly, p.size())
+            fmpz_poly_fit_length(g._poly, p.size())
+            for i in range(p.size()):
+                tmp = p[0][i].get_mpz_const()
+                fmpz_poly_set_coeff_mpz(g._poly, i, tmp)
+            _fmpz_poly_set_length(g._poly, p.size())
+
+            if g.lc() == 1:
+                break
 
         del M
         del p