]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_element_subalgebra.py
eja: begin major overhaul of class hierarchy and naming.
[sage.d.git] / mjo / eja / eja_element_subalgebra.py
index 67342c6be3df71a63eb4cf82337f52ff7edf5fea..dceb3b405a4c5a663c61a966993ce890ed516b49 100644 (file)
@@ -1,53 +1,50 @@
 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
+from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
 
 
-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()
+class FiniteDimensionalEJAElementSubalgebra(FiniteDimensionalEJASubalgebra):
+    def __init__(self, elt, orthonormalize=True, **kwargs):
+        superalgebra = elt.parent()
 
-        # 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)
+        powers = tuple( elt**k for k in range(superalgebra.dimension()) )
+        power_vectors = ( p.to_vector() for p in powers )
+        P = matrix(superalgebra.base_ring(), power_vectors)
 
-        if orthonormalize_basis == False:
-            # 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))
-
-            # If our superalgebra is a subalgebra of something else, then
-            # these vectors won't have the right coordinates for
-            # V.span_of_basis() unless we use V.from_vector() on them.
-            basis_vectors = map(power_vectors.__getitem__, ind_rows)
+        if orthonormalize:
+            basis = powers # let god sort 'em out
         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 ]
-
-        W = V.span_of_basis( V.from_vector(v) for v in basis_vectors )
+            # 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 superalgebra.base_ring() 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.
+                basis = tuple(map(powers.__getitem__, ind_rows))
+
+
+        super().__init__(superalgebra,
+                         basis,
+                         associative=True,
+                         **kwargs)
 
         # The rank is the highest possible degree of a minimal
         # polynomial, and is bounded above by the dimension. We know
@@ -55,39 +52,10 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         # polynomial has the same degree as the space's dimension
         # (remember how we constructed the space?), so that must be
         # its rank too.
-        rank = W.dimension()
-
-        fdeja = super(FiniteDimensionalEuclideanJordanElementSubalgebra, self)
-        return fdeja.__init__(self._superalgebra,
-                              superalgebra_basis,
-                              rank=rank,
-                              category=category)
-
-
-    def _a_regular_element(self):
-        """
-        Override the superalgebra method to return the one
-        regular element that is sure to exist in this
-        subalgebra, namely the element that generated it.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import random_eja
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: J = random_eja().random_element().subalgebra_generated_by()
-            sage: J._a_regular_element().is_regular()
-            True
-
-        """
-        if self.dimension() == 0:
-            return self.zero()
-        else:
-            return self.monomial(1)
+        self.rank.set_cache(self.dimension())
 
 
+    @cached_method
     def one(self):
         """
         Return the multiplicative identity element of this algebra.
@@ -95,7 +63,7 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         The superclass method computes the identity element, which is
         beyond overkill in this case: the superalgebra identity
         restricted to this algebra is its identity. Note that we can't
-        count on the first basis element being the identity -- it migth
+        count on the first basis element being the identity -- it might
         have been scaled if we orthonormalized the basis.
 
         SETUP::
@@ -120,7 +88,7 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         The identity element acts like the identity over the rationals::
 
             sage: set_random_seed()
-            sage: x = random_eja().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
@@ -140,7 +108,7 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         the rationals::
 
             sage: set_random_seed()
-            sage: x = random_eja().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())
@@ -161,7 +129,6 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         """
         if self.dimension() == 0:
             return self.zero()
-        else:
-            sa_one = self.superalgebra().one().to_vector()
-            sa_coords = self.vector_space().coordinate_vector(sa_one)
-            return self.from_vector(sa_coords)
+
+        return self(self.superalgebra().one())
+