]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_element.py
eja: complete the CFM_CartesianProduct purge.
[sage.d.git] / mjo / eja / eja_element.py
index f0f6da756e841e65b7d9272f5feb09404c9f1a09..e8c183b51d522efd14aba57b779b4b5ebfa4b896 100644 (file)
@@ -1,9 +1,10 @@
 from sage.matrix.constructor import matrix
+from sage.misc.cachefunc import cached_method
 from sage.modules.free_module import VectorSpace
 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
 
 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
-from mjo.eja.eja_utils import _mat2vec
+from mjo.eja.eja_utils import _mat2vec, _scale
 
 class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
     """
@@ -438,6 +439,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         return ((-1)**r)*p(*self.to_vector())
 
 
+    @cached_method
     def inverse(self):
         """
         Return the Jordan-multiplicative inverse of this element.
@@ -482,7 +484,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: JordanSpinEJA(3).zero().inverse()
             Traceback (most recent call last):
             ...
-            ValueError: element is not invertible
+            ZeroDivisionError: element is not invertible
 
         TESTS:
 
@@ -534,40 +536,38 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: slow == fast                                   # long time
             True
         """
-        if not self.is_invertible():
-            raise ValueError("element is not invertible")
-
+        not_invertible_msg = "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!
+            if self.det().is_zero():
+                raise ZeroDivisionError(not_invertible_msg)
             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)
+        try:
+            inv = (~self.quadratic_representation())(self)
+            self.is_invertible.set_cache(True)
+            return inv
+        except ZeroDivisionError:
+            self.is_invertible.set_cache(False)
+            raise ZeroDivisionError(not_invertible_msg)
 
 
+    @cached_method
     def is_invertible(self):
         """
         Return whether or not this element is invertible.
 
         ALGORITHM:
 
-        The usual way to do this is to check if the determinant is
-        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 parent algebra's zero element is a root
-        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.
+        If computing my determinant will be fast, we do so and compare
+        with zero (Proposition II.2.4 in Faraut and
+        Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
+        reduces the problem to the invertibility of my quadratic
+        representation.
 
         SETUP::
 
@@ -600,7 +600,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: fast = x.is_invertible()                       # long time
             sage: slow == fast                                   # long time
             True
-
         """
         if self.is_zero():
             if self.parent().is_trivial():
@@ -609,15 +608,17 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
                 return False
 
         if self.parent()._charpoly_coefficients.is_in_cache():
-            # The determinant will be quicker than computing the minimal
-            # polynomial from scratch, most likely.
+            # The determinant will be quicker than inverting the
+            # quadratic representation, 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()
-        p = self.minimal_polynomial()
-        return not (p(zero) == zero)
+        # The easiest way to determine if I'm invertible is to try.
+        try:
+            inv = (~self.quadratic_representation())(self)
+            self.inverse.set_cache(inv)
+            return True
+        except ZeroDivisionError:
+            return False
 
 
     def is_primitive_idempotent(self):
@@ -663,7 +664,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         element should always be in terms of minimal idempotents::
 
             sage: J = JordanSpinEJA(4)
-            sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
+            sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
             sage: x.is_regular()
             True
             sage: [ c.is_primitive_idempotent()
@@ -909,7 +910,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         M = matrix([(self.parent().one()).to_vector()])
         old_rank = 1
 
-        # Specifying the row-reduction algorithm can e.g.  help over
+        # 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.
         #
@@ -1057,7 +1058,7 @@ class FiniteDimensionalEJAElement(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()
 
 
@@ -1076,7 +1077,9 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
-            ....:                                  QuaternionHermitianEJA)
+            ....:                                  HadamardEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  RealSymmetricEJA)
 
         EXAMPLES::
 
@@ -1106,14 +1109,36 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             [0 0 0 0 0 0 1 0]
             [0 0 0 0 0 0 0 1]
 
+        This also works in Cartesian product algebras::
+
+            sage: J1 = HadamardEJA(1)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: x = sum(J.gens())
+            sage: x.to_matrix()[0]
+            [1]
+            sage: x.to_matrix()[1]
+            [                  1 0.7071067811865475?]
+            [0.7071067811865475?                   1]
+
         """
         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()) )
+        if self.parent()._matrix_basis_is_cartesian:
+            # Aaaaand linear combinations don't work in Cartesian
+            # product spaces, even though they provide a method
+            # with that name. This is special-cased because the
+            # _scale() function is slow.
+            pairs = zip(B, self.to_vector())
+            return sum( ( _scale(b, alpha) for (b,alpha) in pairs ),
+                        W.zero())
+        else:
+            # 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):
@@ -1358,7 +1383,7 @@ class FiniteDimensionalEJAElement(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()) )
@@ -1378,7 +1403,20 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import random_eja
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  HadamardEJA,
+            ....:                                  RealSymmetricEJA)
+
+        EXAMPLES:
+
+        We can create subalgebras of Cartesian product EJAs that are not
+        themselves Cartesian product EJAs (they're just "regular" EJAs)::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.one().subalgebra_generated_by()
+            Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
 
         TESTS:
 
@@ -1411,8 +1449,14 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             True
 
         """
-        from mjo.eja.eja_element_subalgebra import FiniteDimensionalEJAElementSubalgebra
-        return FiniteDimensionalEJAElementSubalgebra(self, **kwargs)
+        powers = tuple( self**k for k in range(self.degree()) )
+        A = self.parent().subalgebra(powers,
+                                     associative=True,
+                                     check_field=False,
+                                     check_axioms=False,
+                                     **kwargs)
+        A.one.set_cache(A(self.parent().one()))
+        return A
 
 
     def subalgebra_idempotent(self):
@@ -1519,6 +1563,15 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: J.random_element().trace() in RLF
             True
 
+        The trace is linear::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x,y = J.random_elements(2)
+            sage: alpha = J.base_ring().random_element()
+            sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
+            True
+
         """
         P = self.parent()
         r = P.rank()