eja: fix random_instance() for JordanSpinEJA.

[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 a26381b12dbe649d5dad9bc0cd056adbb1606f62..7bc4b3a9a79c5cae8dbd60db3f63d581d215fc23 100644 (file)
--- a/mjo/eja/eja_element_subalgebra.py
+++ b/mjo/eja/eja_element_subalgebra.py
@@ -1,4 +1,6 @@
 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
 
@@ -18,6 +20,19 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         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
@@ -28,7 +43,6 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
 
             # Pick those out of the list of all powers.
             superalgebra_basis = tuple(map(powers.__getitem__, ind_rows))
-            basis_vectors = map(power_vectors.__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
@@ -58,6 +72,7 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         self.rank.set_cache(self.dimension())
 
 
+    @cached_method
     def one(self):
         """
         Return the multiplicative identity element of this algebra.
@@ -65,7 +80,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::
@@ -131,12 +146,6 @@ class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclide
         """
         if self.dimension() == 0:
             return self.zero()
-        else:
-            sa_one = self.superalgebra().one().to_vector()
-            # The extra hackery is because foo.to_vector() might not
-            # live in foo.parent().vector_space()!
-            coords = sum( a*b for (a,b)
-                          in zip(sa_one,
-                                 self.superalgebra().vector_space().basis()) )
-            return self.from_vector(self.vector_space().coordinate_vector(coords))
+
+        return self(self.superalgebra().one())