From 70809e993ef1cf0cb84ed231eb2781cf05750ab0 Mon Sep 17 00:00:00 2001
From: Travis Scrimshaw <tcscrims@gmail.com>
Date: Wed, 3 Apr 2024 14:56:31 +0900
Subject: [PATCH] Don't do inexact rings, better error messages, and cache the
 results

---
 src/sage/matrix/matrix2.pyx | 42 ++++++++++++++++++++++++++++++++++---
 1 file changed, 39 insertions(+), 3 deletions(-)

diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx
index 637b69c7d37..9a7fdd393c5 100644
--- a/src/sage/matrix/matrix2.pyx
+++ b/src/sage/matrix/matrix2.pyx
@@ -11729,16 +11729,52 @@ cdef class Matrix(Matrix1):
             Traceback (most recent call last):
             ...
             RuntimeError: Some eigenvalue does not exist in Rational Field.
+
+        TESTS::
+
+            sage: X = random_matrix(QQ, 4)
+            sage: S, N = X.jordan_decomposition()
+            sage: X == S + N
+            True
+            sage: S is X.jordan_decomposition()[0]  # result is cached
+            True
+            sage: N is X.jordan_decomposition()[1]  # result is cached
+            True
+            sage: A = matrix(ZZ, 5, 5, {(0,1): -1, (1,0): 1, (2,3): -1})
+            sage: A.jordan_decomposition()
+            Traceback (most recent call last):
+            ...
+            ValueError: unable to compute Jordan decomposition
+            sage: B = A.change_ring(RR)
+            sage: B.jordan_decomposition()
+            Traceback (most recent call last):
+            ...
+            NotImplementedError: Jordan decomposition not implemented over inexact rings
         """
+        if not self.base_ring().is_exact():
+            raise NotImplementedError("Jordan decomposition not implemented over inexact rings")
+        JD = self.fetch('jordan_decomposition')
+        if JD is not None:
+            return JD
         f = self.minpoly()
         h = f // f.gcd(f.diff())
         o, p, q = h.xgcd(h.diff())
-        assert o.is_one()
+        if not o.is_one():
+            raise ValueError("unable to compute Jordan decomposition")
         A = self
         hq = h * q
-        while h(A):
+        # very bad bound, but requires no extra computation
+        # a better bound is the maximum multiplicity in the minpoly,
+        #   but this requires factoring the minpoly
+        for _ in range(self.nrows()):
+            if not h(A):
+                ret = (A, self - A)
+                ret[0].set_immutable()
+                ret[1].set_immutable()
+                self.cache('jordan_decomposition', ret)
+                return ret
             A -= hq(A)
-        return (A, self - A)
+        raise ValueError("Jordan decomposition does not exist")
 
     def diagonalization(self, base_field=None):
         """