]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_element.py
eja: more orthonormalization fixes.
[sage.d.git] / mjo / eja / eja_element.py
index d99a7d873ab69e511186b7f0bc3ba3127a45c85a..d3e9a33ceba6e1fc4e58b417f3b953e2e2d1c3d7 100644 (file)
@@ -2,10 +2,10 @@ from sage.matrix.constructor import matrix
 from sage.modules.free_module import VectorSpace
 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
 
-from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
+from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 from mjo.eja.eja_utils import _mat2vec
 
-class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
+class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
     """
     An element of a Euclidean Jordan algebra.
     """
@@ -444,8 +444,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ALGORITHM:
 
-        We appeal to the quadratic representation as in Koecher's
-        Theorem 12 in Chapter III, Section 5.
+        In general we appeal to the quadratic representation as in
+        Koecher's Theorem 12 in Chapter III, Section 5. But if the
+        parent algebra's "characteristic polynomial of" coefficients
+        happen to be cached, then we use Proposition II.2.4 in Faraut
+        and Korányi which gives a formula for the inverse based on the
+        characteristic polynomial and the Cayley-Hamilton theorem for
+        Euclidean Jordan algebras::
 
         SETUP::
 
@@ -515,22 +520,19 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             ....:    x.operator().inverse()(J.one()) == x.inverse() )
             True
 
-        Proposition II.2.4 in Faraut and Korányi gives a formula for
-        the inverse based on the characteristic polynomial and the
-        Cayley-Hamilton theorem for Euclidean Jordan algebras::
+        Check that the fast (cached) and slow algorithms give the same
+        answer::
 
-            sage: set_random_seed()
-            sage: J = ComplexHermitianEJA(3)
-            sage: x = J.random_element()
-            sage: while not x.is_invertible():
-            ....:     x = J.random_element()
-            sage: r = J.rank()
-            sage: a = x.characteristic_polynomial().coefficients(sparse=False)
-            sage: expected  = (-1)^(r+1)/x.det()
-            sage: expected *= sum( a[i+1]*x^i for i in range(r) )
-            sage: x.inverse() == expected
+            sage: set_random_seed()                              # long time
+            sage: J = random_eja(field=QQ, orthonormalize=False) # long time
+            sage: x = J.random_element()                         # long time
+            sage: while not x.is_invertible():                   # long time
+            ....:     x = J.random_element()                     # long time
+            sage: slow = x.inverse()                             # long time
+            sage: _ = J._charpoly_coefficients()                 # long time
+            sage: fast = x.inverse()                             # long time
+            sage: slow == fast                                   # long time
             True
-
         """
         if not self.is_invertible():
             raise ValueError("element is not invertible")
@@ -587,6 +589,18 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: (not J.is_trivial()) and J.zero().is_invertible()
             False
 
+        Test that the fast (cached) and slow algorithms give the same
+        answer::
+
+            sage: set_random_seed()                              # long time
+            sage: J = random_eja(field=QQ, orthonormalize=False) # long time
+            sage: x = J.random_element()                         # long time
+            sage: slow = x.is_invertible()                       # long time
+            sage: _ = J._charpoly_coefficients()                 # long time
+            sage: fast = x.is_invertible()                       # long time
+            sage: slow == fast                                   # long time
+            True
+
         """
         if self.is_zero():
             if self.parent().is_trivial():
@@ -816,10 +830,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ALGORITHM:
 
-        For now, we skip the messy minimal polynomial computation
-        and instead return the dimension of the vector space spanned
-        by the powers of this element. The latter is a bit more
-        straightforward to compute.
+        .........
 
         SETUP::
 
@@ -867,12 +878,59 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         """
-        if self.is_zero() and not self.parent().is_trivial():
+        n = self.parent().dimension()
+
+        if n == 0:
+            # The minimal polynomial is an empty product, i.e. the
+            # constant polynomial "1" having degree zero.
+            return 0
+        elif self.is_zero():
             # The minimal polynomial of zero in a nontrivial algebra
-            # is "t"; in a trivial algebra it's "1" by convention
-            # (it's an empty product).
+            # is "t", and is of degree one.
             return 1
-        return self.subalgebra_generated_by().dimension()
+        elif n == 1:
+            # If this is a nonzero element of a nontrivial algebra, it
+            # has degree at least one. It follows that, in an algebra
+            # of dimension one, the degree must be actually one.
+            return 1
+
+        # BEWARE: The subalgebra_generated_by() method uses the result
+        # of this method to construct a basis for the subalgebra. That
+        # means, in particular, that we cannot implement this method
+        # as ``self.subalgebra_generated_by().dimension()``.
+
+        # Algorithm: keep appending (vector representations of) powers
+        # self as rows to a matrix and echelonizing it. When its rank
+        # stops increasing, we've reached a redundancy.
+
+        # Given the special cases above, we can assume that "self" is
+        # nonzero, the algebra is nontrivial, and that its dimension
+        # is at least two.
+        M = matrix([(self.parent().one()).to_vector()])
+        old_rank = 1
+
+        # Specifying the row-reduction algorithm can e.g.  help over
+        # AA because it avoids the RecursionError that gets thrown
+        # when we have to look too hard for a root.
+        #
+        # Beware: QQ supports an entirely different set of "algorithm"
+        # keywords than do AA and RR.
+        algo = None
+        from sage.rings.all import QQ
+        if self.parent().base_ring() is not QQ:
+            algo = "scaled_partial_pivoting"
+
+        for d in range(1,n):
+            M = matrix(M.rows() + [(self**d).to_vector()])
+            M.echelonize(algo)
+            new_rank = M.rank()
+            if new_rank == old_rank:
+                return new_rank
+            else:
+                old_rank = new_rank
+
+        return n
+
 
 
     def left_matrix(self):
@@ -999,7 +1057,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
                 # in the "normal" case without us having to think about it.
                 return self.operator().minimal_polynomial()
 
-        A = self.subalgebra_generated_by(orthonormalize_basis=False)
+        A = self.subalgebra_generated_by(orthonormalize=False)
         return A(self).operator().minimal_polynomial()
 
 
@@ -1108,10 +1166,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         P = self.parent()
         left_mult_by_self = lambda y: self*y
         L = P.module_morphism(function=left_mult_by_self, codomain=P)
-        return FiniteDimensionalEuclideanJordanAlgebraOperator(
-                 P,
-                 P,
-                 L.matrix() )
+        return FiniteDimensionalEJAOperator(P, P, L.matrix() )
 
 
     def quadratic_representation(self, other=None):
@@ -1303,13 +1358,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             [(0, f2), (1, f0)]
 
         """
-        A = self.subalgebra_generated_by(orthonormalize_basis=True)
+        A = self.subalgebra_generated_by(orthonormalize=True)
         result = []
         for (evalue, proj) in A(self).operator().spectral_decomposition():
             result.append( (evalue, proj(A.one()).superalgebra_element()) )
         return result
 
-    def subalgebra_generated_by(self, orthonormalize_basis=False):
+    def subalgebra_generated_by(self, **kwargs):
         """
         Return the associative subalgebra of the parent EJA generated
         by this element.
@@ -1356,8 +1411,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         """
-        from mjo.eja.eja_element_subalgebra import FiniteDimensionalEuclideanJordanElementSubalgebra
-        return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
+        from mjo.eja.eja_element_subalgebra import FiniteDimensionalEJAElementSubalgebra
+        return FiniteDimensionalEJAElementSubalgebra(self, **kwargs)
 
 
     def subalgebra_idempotent(self):