]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_element.py
eja: speed up minimal_polynomial(), in theory.
[sage.d.git] / mjo / eja / eja_element.py
index 66138b2089cfcc2190c4282393c221e2c7c188fc..52933e2decdcf3d9dc1218a568bf29bcff6cf5f1 100644 (file)
@@ -4,7 +4,7 @@ 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):
     """
@@ -664,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()
@@ -910,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.
         #
@@ -1047,19 +1047,30 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         """
         if self.is_zero():
-            # We would generate a zero-dimensional subalgebra
-            # where the minimal polynomial would be constant.
-            # That might be correct, but only if *this* algebra
-            # is trivial too.
-            if not self.parent().is_trivial():
-                # Pretty sure we know what the minimal polynomial of
-                # the zero operator is going to be. This ensures
-                # consistency of e.g. the polynomial variable returned
-                # in the "normal" case without us having to think about it.
-                return self.operator().minimal_polynomial()
-
+            # Pretty sure we know what the minimal polynomial of
+            # the zero operator is going to be. This ensures
+            # consistency of e.g. the polynomial variable returned
+            # in the "normal" case without us having to think about it.
+            return self.operator().minimal_polynomial()
+
+        # If we don't orthonormalize the subalgebra's basis, then the
+        # first two monomials in the subalgebra will be self^0 and
+        # self^1... assuming that self^1 is not a scalar multiple of
+        # self^0 (the unit element). We special case these to avoid
+        # having to solve a system to coerce self into the subalgebra.
         A = self.subalgebra_generated_by(orthonormalize=False)
-        return A(self).operator().minimal_polynomial()
+
+        if A.dimension() == 1:
+            # Does a solve to find the scalar multiple alpha such that
+            # alpha*unit = self. We have to do this because the basis
+            # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
+            unit = self.parent().one()
+            alpha = self.to_vector() / unit.to_vector()
+            return (unit.operator()*alpha).minimal_polynomial()
+        else:
+            # If the dimension of the subalgebra is >= 2, then we just
+            # use the second basis element.
+            return A.monomial(1).operator().minimal_polynomial()
 
 
 
@@ -1077,7 +1088,9 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
-            ....:                                  QuaternionHermitianEJA)
+            ....:                                  HadamardEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  RealSymmetricEJA)
 
         EXAMPLES::
 
@@ -1107,14 +1120,35 @@ 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 hasattr(W, 'cartesian_factors'):
+            # Aaaaand linear combinations don't work in Cartesian
+            # product spaces, even though they provide a method with
+            # that name. This is hidden behind an "if" because the
+            # _scale() function is slow.
+            pairs = zip(B, self.to_vector())
+            return W.sum( _scale(b, alpha) for (b,alpha) in pairs )
+        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):
@@ -1379,7 +1413,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:
 
@@ -1412,12 +1459,12 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             True
 
         """
-        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 = self.parent().subalgebra(powers,
+                                     associative=True,
+                                     check_field=False,
+                                     check_axioms=False,
+                                     **kwargs)
         A.one.set_cache(A(self.parent().one()))
         return A
 
@@ -1526,6 +1573,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()