eja: work towards fixing element subalgebras.

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

diff --git a/mjo/eja/eja_element_subalgebra.py b/mjo/eja/eja_element_subalgebra.py
index 7bc4b3a9a79c5cae8dbd60db3f63d581d215fc23..ee8e0c67db2773537aefc04d82ee58b7a0b4cca0 100644 (file)
--- a/mjo/eja/eja_element_subalgebra.py
+++ b/mjo/eja/eja_element_subalgebra.py
@@ -2,66 +2,24 @@ from sage.matrix.constructor import matrix
 from sage.misc.cachefunc import cached_method
 from sage.rings.all import QQ
 
-from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
-
-
-class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclideanJordanSubalgebra):
-    def __init__(self, elt, orthonormalize_basis):
-        self._superalgebra = elt.parent()
-        category = self._superalgebra.category().Associative()
-        V = self._superalgebra.vector_space()
-        field = self._superalgebra.base_ring()
-
-        # This list is guaranteed to contain all independent powers,
-        # because it's the maximal set of powers that could possibly
-        # be independent (by a dimension argument).
-        powers = [ elt**k for k in range(V.dimension()) ]
-        power_vectors = [ p.to_vector() for p in powers ]
-        P = matrix(field, power_vectors)
-
-        if orthonormalize_basis == False:
-            # Echelonize the matrix ourselves, because otherwise the
-            # call to P.pivot_rows() below can choose a non-optimal
-            # row-reduction algorithm. In particular, scaling can
-            # 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
-            if field is not QQ:
-                algo = "scaled_partial_pivoting"
-            P.echelonize(algorithm=algo)
-
-            # In this case, we just need to figure out which elements
-            # of the "powers" list are redundant... First compute the
-            # vector subspace spanned by the powers of the given
-            # element.
-
-            # Figure out which powers form a linearly-independent set.
-            ind_rows = P.pivot_rows()
-
-            # Pick those out of the list of all powers.
-            superalgebra_basis = tuple(map(powers.__getitem__, ind_rows))
-        else:
-            # If we're going to orthonormalize the basis anyway, we
-            # might as well just do Gram-Schmidt on the whole list of
-            # powers. The redundant ones will get zero'd out. If this
-            # looks like a roundabout way to orthonormalize, it is.
-            # But converting everything from algebra elements to vectors
-            # to matrices and then back again turns out to be about
-            # as fast as reimplementing our own Gram-Schmidt that
-            # works in an EJA.
-            G,_ = P.gram_schmidt(orthonormal=True)
-            basis_vectors = [ g for g in G.rows() if not g.is_zero() ]
-            superalgebra_basis = [ self._superalgebra.from_vector(b)
-                                   for b in basis_vectors ]
-
-        fdeja = super(FiniteDimensionalEuclideanJordanElementSubalgebra, self)
-        fdeja.__init__(self._superalgebra,
-                       superalgebra_basis,
-                       category=category,
-                       check_axioms=False)
+from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+
+
+class FiniteDimensionalEJAElementSubalgebra(FiniteDimensionalEJASubalgebra):
+    def __init__(self, elt, orthonormalize=True, **kwargs):
+        superalgebra = elt.parent()
+
+        # TODO: going up to the superalgebra dimension here is
+        # overkill.  We should append p vectors as rows to a matrix
+        # and continually rref() it until the rank stops going
+        # up. When n=10 but the dimension of the algebra is 1, that
+        # can save a shitload of time (especially over AA).
+        powers = tuple( elt**k for k in range(elt.degree()) )
+
+        super().__init__(superalgebra,
+                         powers,
+                         associative=True,
+                         **kwargs)
 
         # The rank is the highest possible degree of a minimal
         # polynomial, and is bounded above by the dimension. We know
@@ -105,7 +63,7 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         The identity element acts like the identity over the rationals::
 
             sage: set_random_seed()
-            sage: x = random_eja(field=QQ).random_element()
+            sage: x = random_eja(field=QQ,orthonormalize=False).random_element()
             sage: A = x.subalgebra_generated_by()
             sage: x = A.random_element()
             sage: A.one()*x == x and x*A.one() == x
@@ -125,7 +83,7 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         the rationals::
 
             sage: set_random_seed()
-            sage: x = random_eja(field=QQ).random_element()
+            sage: x = random_eja(field=QQ,orthonormalize=False).random_element()
             sage: A = x.subalgebra_generated_by()
             sage: actual = A.one().operator().matrix()
             sage: expected = matrix.identity(A.base_ring(), A.dimension())