eja: eliminate the special element subalgebra class.

[sage.d.git] / mjo / eja / eja_element.py

diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py
index a6230a40c2fa0ce88830366b7690eb57ef69963b..e30dbb13f39d06b6d49c6bfd34c829ab30d6112d 100644 (file)
--- a/mjo/eja/eja_element.py
+++ b/mjo/eja/eja_element.py
@@ -2,15 +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 sage.modules.free_module import VectorSpace
 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
 

-# TODO: make this unnecessary somehow.
-from sage.misc.lazy_import import lazy_import
-lazy_import('mjo.eja.eja_algebra', 'FiniteDimensionalEuclideanJordanAlgebra')
-lazy_import('mjo.eja.eja_element_subalgebra',
-            'FiniteDimensionalEuclideanJordanElementSubalgebra')
-from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
+from mjo.eja.eja_operator import FiniteDimensionalEJAOperator

 from mjo.eja.eja_utils import _mat2vec
 
 from mjo.eja.eja_utils import _mat2vec
 

-class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
+class FiniteDimensionalEJAElement(IndexedFreeModuleElement):

     """
     An element of a Euclidean Jordan algebra.
     """
     """
     An element of a Euclidean Jordan algebra.
     """


@@ -181,7 +176,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         """
             True
 
         """

-        p = self.parent().characteristic_polynomial()
+        p = self.parent().characteristic_polynomial_of()

         return p(*self.to_vector())
 
 
         return p(*self.to_vector())
 
 


@@ -232,9 +227,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         Ditto for the quaternions::
 
 
         Ditto for the quaternions::
 

-            sage: J = QuaternionHermitianEJA(3)
+            sage: J = QuaternionHermitianEJA(2)

             sage: J.one().inner_product(J.one())
             sage: J.one().inner_product(J.one())

-            3
+            2

 
         TESTS:
 
 
         TESTS:
 


@@ -345,6 +340,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
             ....:                                  TrivialEJA,
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
             ....:                                  TrivialEJA,

+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,

             ....:                                  random_eja)
 
         EXAMPLES::
             ....:                                  random_eja)
 
         EXAMPLES::


@@ -392,6 +389,37 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: x,y = J.random_elements(2)
             sage: (x*y).det() == x.det()*y.det()
             True
             sage: x,y = J.random_elements(2)
             sage: (x*y).det() == x.det()*y.det()
             True

+
+        The determinant in matrix algebras is just the usual determinant::
+
+            sage: set_random_seed()
+            sage: X = matrix.random(QQ,3)
+            sage: X = X + X.T
+            sage: J1 = RealSymmetricEJA(3)
+            sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
+            sage: expected = X.det()
+            sage: actual1 = J1(X).det()
+            sage: actual2 = J2(X).det()
+            sage: actual1 == expected
+            True
+            sage: actual2 == expected
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J1 = ComplexHermitianEJA(2)
+            sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: X = matrix.random(GaussianIntegers(), 2)
+            sage: X = X + X.H
+            sage: expected = AA(X.det())
+            sage: actual1 = J1(J1.real_embed(X)).det()
+            sage: actual2 = J2(J2.real_embed(X)).det()
+            sage: expected == actual1
+            True
+            sage: expected == actual2
+            True
+

         """
         P = self.parent()
         r = P.rank()
         """
         P = self.parent()
         r = P.rank()


@@ -416,8 +444,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ALGORITHM:
 
 
         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::
 
 
         SETUP::
 


@@ -487,26 +520,31 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             ....:    x.operator().inverse()(J.one()) == x.inverse() )
             True
 
             ....:    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
             True

-

         """
         if not self.is_invertible():
             raise ValueError("element is not invertible")
 
         """
         if not self.is_invertible():
             raise ValueError("element is not invertible")
 

+        if self.parent()._charpoly_coefficients.is_in_cache():
+            # We can invert using our charpoly if it will be fast to
+            # compute. If the coefficients are cached, our rank had
+            # better be too!
+            r = self.parent().rank()
+            a = self.characteristic_polynomial().coefficients(sparse=False)
+            return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det()
+

         return (~self.quadratic_representation())(self)
 
 
         return (~self.quadratic_representation())(self)
 
 


@@ -520,9 +558,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         zero, but we need the characteristic polynomial for the
         determinant. The minimal polynomial is a lot easier to get,
         so we use Corollary 2 in Chapter V of Koecher to check
         zero, but we need the characteristic polynomial for the
         determinant. The minimal polynomial is a lot easier to get,
         so we use Corollary 2 in Chapter V of Koecher to check

-        whether or not the paren't algebra's zero element is a root
+        whether or not the parent algebra's zero element is a root

         of this element's minimal polynomial.
 
         of this element's minimal polynomial.
 

+        That is... unless the coefficients of our algebra's
+        "characteristic polynomial of" function are already cached!
+        In that case, we just use the determinant (which will be fast
+        as a result).
+

         Beware that we can't use the superclass method, because it
         relies on the algebra being associative.
 
         Beware that we can't use the superclass method, because it
         relies on the algebra being associative.
 


@@ -546,6 +589,18 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: (not J.is_trivial()) and J.zero().is_invertible()
             False
 
             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():
         """
         if self.is_zero():
             if self.parent().is_trivial():


@@ -553,6 +608,11 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             else:
                 return False
 
             else:
                 return False
 

+        if self.parent()._charpoly_coefficients.is_in_cache():
+            # The determinant will be quicker than computing the minimal
+            # polynomial from scratch, most likely.
+            return (not self.det().is_zero())
+

         # In fact, we only need to know if the constant term is non-zero,
         # so we can pass in the field's zero element instead.
         zero = self.base_ring().zero()
         # In fact, we only need to know if the constant term is non-zero,
         # so we can pass in the field's zero element instead.
         zero = self.base_ring().zero()


@@ -770,10 +830,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ALGORITHM:
 
 
         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::
 
 
         SETUP::
 


@@ -821,12 +878,59 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         """
             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
             # 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
+        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
             return 1

-        return self.subalgebra_generated_by().dimension()
+
+        # 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):
 
 
     def left_matrix(self):


@@ -903,7 +1007,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         two here so that said elements actually exist::
 
             sage: set_random_seed()
         two here so that said elements actually exist::
 
             sage: set_random_seed()

-            sage: n_max = max(2, JordanSpinEJA._max_test_case_size())
+            sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())

             sage: n = ZZ.random_element(2, n_max)
             sage: J = JordanSpinEJA(n)
             sage: y = J.random_element()
             sage: n = ZZ.random_element(2, n_max)
             sage: J = JordanSpinEJA(n)
             sage: y = J.random_element()


@@ -929,10 +1033,10 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         and in particular, a re-scaling of the basis::
 
             sage: set_random_seed()
         and in particular, a re-scaling of the basis::
 
             sage: set_random_seed()

-            sage: n_max = RealSymmetricEJA._max_test_case_size()
+            sage: n_max = RealSymmetricEJA._max_random_instance_size()

             sage: n = ZZ.random_element(1, n_max)
             sage: J1 = RealSymmetricEJA(n)
             sage: n = ZZ.random_element(1, n_max)
             sage: J1 = RealSymmetricEJA(n)

-            sage: J2 = RealSymmetricEJA(n,normalize_basis=False)
+            sage: J2 = RealSymmetricEJA(n,orthonormalize=False)

             sage: X = random_matrix(AA,n)
             sage: X = X*X.transpose()
             sage: x1 = J1(X)
             sage: X = random_matrix(AA,n)
             sage: X = X*X.transpose()
             sage: x1 = J1(X)


@@ -953,20 +1057,21 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
                 # in the "normal" case without us having to think about it.
                 return self.operator().minimal_polynomial()
 
                 # in the "normal" case without us having to think about it.
                 return self.operator().minimal_polynomial()
 

-        A = self.subalgebra_generated_by()
+        A = self.subalgebra_generated_by(orthonormalize=False)

         return A(self).operator().minimal_polynomial()
 
 
 
         return A(self).operator().minimal_polynomial()
 
 
 

-    def natural_representation(self):
+    def to_matrix(self):

         """
         """

-        Return a more-natural representation of this element.
+        Return an (often more natural) representation of this element as a
+        matrix.

 
 

-        Every finite-dimensional Euclidean Jordan Algebra is a
-        direct sum of five simple algebras, four of which comprise
-        Hermitian matrices. This method returns the original
-        "natural" representation of this element as a Hermitian
-        matrix, if it has one. If not, you get the usual representation.
+        Every finite-dimensional Euclidean Jordan Algebra is a direct
+        sum of five simple algebras, four of which comprise Hermitian
+        matrices. This method returns a "natural" matrix
+        representation of this element as either a Hermitian matrix or
+        column vector.

 
         SETUP::
 
 
         SETUP::
 


@@ -978,7 +1083,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: J = ComplexHermitianEJA(3)
             sage: J.one()
             e0 + e3 + e8
             sage: J = ComplexHermitianEJA(3)
             sage: J.one()
             e0 + e3 + e8

-            sage: J.one().natural_representation()
+            sage: J.one().to_matrix()

             [1 0 0 0 0 0]
             [0 1 0 0 0 0]
             [0 0 1 0 0 0]
             [1 0 0 0 0 0]
             [0 1 0 0 0 0]
             [0 0 1 0 0 0]


@@ -988,27 +1093,27 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ::
 
 
         ::
 

-            sage: J = QuaternionHermitianEJA(3)
+            sage: J = QuaternionHermitianEJA(2)

             sage: J.one()
             sage: J.one()

-            e0 + e5 + e14
-            sage: J.one().natural_representation()
-            [1 0 0 0 0 0 0 0 0 0 0 0]
-            [0 1 0 0 0 0 0 0 0 0 0 0]
-            [0 0 1 0 0 0 0 0 0 0 0 0]
-            [0 0 0 1 0 0 0 0 0 0 0 0]
-            [0 0 0 0 1 0 0 0 0 0 0 0]
-            [0 0 0 0 0 1 0 0 0 0 0 0]
-            [0 0 0 0 0 0 1 0 0 0 0 0]
-            [0 0 0 0 0 0 0 1 0 0 0 0]
-            [0 0 0 0 0 0 0 0 1 0 0 0]
-            [0 0 0 0 0 0 0 0 0 1 0 0]
-            [0 0 0 0 0 0 0 0 0 0 1 0]
-            [0 0 0 0 0 0 0 0 0 0 0 1]
+            e0 + e5
+            sage: J.one().to_matrix()
+            [1 0 0 0 0 0 0 0]
+            [0 1 0 0 0 0 0 0]
+            [0 0 1 0 0 0 0 0]
+            [0 0 0 1 0 0 0 0]
+            [0 0 0 0 1 0 0 0]
+            [0 0 0 0 0 1 0 0]
+            [0 0 0 0 0 0 1 0]
+            [0 0 0 0 0 0 0 1]

 
         """
 
         """

-        B = self.parent().natural_basis()
-        W = self.parent().natural_basis_space()
-        return W.linear_combination(zip(B,self.to_vector()))
+        B = self.parent().matrix_basis()
+        W = self.parent().matrix_space()
+
+        # This is just a manual "from_vector()", but of course
+        # matrix spaces aren't vector spaces in sage, so they
+        # don't have a from_vector() method.
+        return W.linear_combination( zip(B, self.to_vector()) )

 
 
     def norm(self):
 
 
     def norm(self):


@@ -1061,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)
         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):
 
 
     def quadratic_representation(self, other=None):


@@ -1253,16 +1355,16 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: (J0, J5, J1) = J.peirce_decomposition(c1)
             sage: (f0, f1, f2) = J1.gens()
             sage: f0.spectral_decomposition()
             sage: (J0, J5, J1) = J.peirce_decomposition(c1)
             sage: (f0, f1, f2) = J1.gens()
             sage: f0.spectral_decomposition()

-            [(0, 1.000000000000000?*f2), (1, 1.000000000000000?*f0)]
+            [(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
 
         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.
         """
         Return the associative subalgebra of the parent EJA generated
         by this element.


@@ -1309,7 +1411,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         """
             True
 
         """

-        return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
+        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+        powers = tuple( self**k for k in range(self.degree()) )
+        A = FiniteDimensionalEJASubalgebra(self.parent(),
+                                           powers,
+                                           associative=True,
+                                           **kwargs)
+        A.one.set_cache(A(self.parent().one()))
+        return A

 
 
     def subalgebra_idempotent(self):
 
 
     def subalgebra_idempotent(self):