eja: speed up minimal_polynomial(), in theory.

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

diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py
index 9a770ae5f68b3e19f3946ca7716f299a3ff82685..52933e2decdcf3d9dc1218a568bf29bcff6cf5f1 100644 (file)
--- a/mjo/eja/eja_element.py
+++ b/mjo/eja/eja_element.py
@@ -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()
@@ -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()
 
 
 
@@ -1125,14 +1136,13 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         B = self.parent().matrix_basis()
         W = self.parent().matrix_space()
 
-        if self.parent()._matrix_basis_is_cartesian:
+        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 special-cased because the
+            # 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 sum( ( _scale(b, alpha) for (b,alpha) in pairs ),
-                        W.zero())
+            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