]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_element.py
README: rewrite it, it was rather out-of-date
[sage.d.git] / mjo / eja / eja_element.py
index c98d9a2b64651562928a6489ef88e03fdd2b0dc9..f4d5995ce9302d2fdc518dac2b0ca20f0fadf6d1 100644 (file)
@@ -3,10 +3,11 @@ from sage.misc.cachefunc import cached_method
 from sage.modules.free_module import VectorSpace
 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
 
 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, _scale
+from mjo.eja.eja_operator import EJAOperator
+from mjo.eja.eja_utils import _scale
 
 
-class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
+
+class EJAElement(IndexedFreeModuleElement):
     """
     An element of a Euclidean Jordan algebra.
     """
     """
     An element of a Euclidean Jordan algebra.
     """
@@ -42,14 +43,12 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         The definition of `x^2` is the unambiguous `x*x`::
 
 
         The definition of `x^2` is the unambiguous `x*x`::
 
-            sage: set_random_seed()
             sage: x = random_eja().random_element()
             sage: x*x == (x^2)
             True
 
         A few examples of power-associativity::
 
             sage: x = random_eja().random_element()
             sage: x*x == (x^2)
             True
 
         A few examples of power-associativity::
 
-            sage: set_random_seed()
             sage: x = random_eja().random_element()
             sage: x*(x*x)*(x*x) == x^5
             True
             sage: x = random_eja().random_element()
             sage: x*(x*x)*(x*x) == x^5
             True
@@ -59,7 +58,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         We also know that powers operator-commute (Koecher, Chapter
         III, Corollary 1)::
 
         We also know that powers operator-commute (Koecher, Chapter
         III, Corollary 1)::
 
-            sage: set_random_seed()
             sage: x = random_eja().random_element()
             sage: m = ZZ.random_element(0,10)
             sage: n = ZZ.random_element(0,10)
             sage: x = random_eja().random_element()
             sage: m = ZZ.random_element(0,10)
             sage: n = ZZ.random_element(0,10)
@@ -106,7 +104,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         We should always get back an element of the algebra::
 
 
         We should always get back an element of the algebra::
 
-            sage: set_random_seed()
             sage: p = PolynomialRing(AA, 't').random_element()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: p = PolynomialRing(AA, 't').random_element()
             sage: J = random_eja()
             sage: x = J.random_element()
@@ -131,7 +128,8 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         SETUP::
 
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import HadamardEJA
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  HadamardEJA)
 
         EXAMPLES:
 
 
         EXAMPLES:
 
@@ -155,11 +153,10 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         The characteristic polynomial of an element should evaluate
         to zero on that element::
 
         The characteristic polynomial of an element should evaluate
         to zero on that element::
 
-            sage: set_random_seed()
-            sage: x = HadamardEJA(3).random_element()
+            sage: x = random_eja().random_element()
             sage: p = x.characteristic_polynomial()
             sage: p = x.characteristic_polynomial()
-            sage: x.apply_univariate_polynomial(p)
-            0
+            sage: x.apply_univariate_polynomial(p).is_zero()
+            True
 
         The characteristic polynomials of the zero and unit elements
         should be what we think they are in a subalgebra, too::
 
         The characteristic polynomials of the zero and unit elements
         should be what we think they are in a subalgebra, too::
@@ -237,7 +234,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         Ensure that we can always compute an inner product, and that
         it gives us back a real number::
 
         Ensure that we can always compute an inner product, and that
         it gives us back a real number::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x,y = J.random_elements(2)
             sage: x.inner_product(y) in RLF
             sage: J = random_eja()
             sage: x,y = J.random_elements(2)
             sage: x.inner_product(y) in RLF
@@ -265,7 +261,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         The definition of a Jordan algebra says that any element
         operator-commutes with its square::
 
         The definition of a Jordan algebra says that any element
         operator-commutes with its square::
 
-            sage: set_random_seed()
             sage: x = random_eja().random_element()
             sage: x.operator_commutes_with(x^2)
             True
             sage: x = random_eja().random_element()
             sage: x.operator_commutes_with(x^2)
             True
@@ -274,7 +269,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         Test Lemma 1 from Chapter III of Koecher::
 
 
         Test Lemma 1 from Chapter III of Koecher::
 
-            sage: set_random_seed()
             sage: u,v = random_eja().random_elements(2)
             sage: lhs = u.operator_commutes_with(u*v)
             sage: rhs = v.operator_commutes_with(u^2)
             sage: u,v = random_eja().random_elements(2)
             sage: lhs = u.operator_commutes_with(u*v)
             sage: rhs = v.operator_commutes_with(u^2)
@@ -284,7 +278,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         Test the first polarization identity from my notes, Koecher
         Chapter III, or from Baes (2.3)::
 
         Test the first polarization identity from my notes, Koecher
         Chapter III, or from Baes (2.3)::
 
-            sage: set_random_seed()
             sage: x,y = random_eja().random_elements(2)
             sage: Lx = x.operator()
             sage: Ly = y.operator()
             sage: x,y = random_eja().random_elements(2)
             sage: Lx = x.operator()
             sage: Ly = y.operator()
@@ -296,32 +289,32 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         Test the second polarization identity from my notes or from
         Baes (2.4)::
 
         Test the second polarization identity from my notes or from
         Baes (2.4)::
 
-            sage: set_random_seed()
-            sage: x,y,z = random_eja().random_elements(3)
-            sage: Lx = x.operator()
-            sage: Ly = y.operator()
-            sage: Lz = z.operator()
-            sage: Lzy = (z*y).operator()
-            sage: Lxy = (x*y).operator()
-            sage: Lxz = (x*z).operator()
-            sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
+            sage: x,y,z = random_eja().random_elements(3)  # long time
+            sage: Lx = x.operator()                        # long time
+            sage: Ly = y.operator()                        # long time
+            sage: Lz = z.operator()                        # long time
+            sage: Lzy = (z*y).operator()                   # long time
+            sage: Lxy = (x*y).operator()                   # long time
+            sage: Lxz = (x*z).operator()                   # long time
+            sage: lhs = Lx*Lzy + Lz*Lxy + Ly*Lxz           # long time
+            sage: rhs = Lzy*Lx + Lxy*Lz + Lxz*Ly           # long time
+            sage: bool(lhs == rhs)                         # long time
             True
 
         Test the third polarization identity from my notes or from
         Baes (2.5)::
 
             True
 
         Test the third polarization identity from my notes or from
         Baes (2.5)::
 
-            sage: set_random_seed()
-            sage: u,y,z = random_eja().random_elements(3)
-            sage: Lu = u.operator()
-            sage: Ly = y.operator()
-            sage: Lz = z.operator()
-            sage: Lzy = (z*y).operator()
-            sage: Luy = (u*y).operator()
-            sage: Luz = (u*z).operator()
-            sage: Luyz = (u*(y*z)).operator()
-            sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
-            sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
-            sage: bool(lhs == rhs)
+            sage: u,y,z = random_eja().random_elements(3)  # long time
+            sage: Lu = u.operator()                        # long time
+            sage: Ly = y.operator()                        # long time
+            sage: Lz = z.operator()                        # long time
+            sage: Lzy = (z*y).operator()                   # long time
+            sage: Luy = (u*y).operator()                   # long time
+            sage: Luz = (u*z).operator()                   # long time
+            sage: Luyz = (u*(y*z)).operator()              # long time
+            sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz           # long time
+            sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly         # long time
+            sage: bool(lhs == rhs)                         # long time
             True
 
         """
             True
 
         """
@@ -339,7 +332,8 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         SETUP::
 
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (AlbertEJA,
+            ....:                                  JordanSpinEJA,
             ....:                                  TrivialEJA,
             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
             ....:                                  TrivialEJA,
             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
@@ -375,7 +369,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         An element is invertible if and only if its determinant is
         non-zero::
 
         An element is invertible if and only if its determinant is
         non-zero::
 
-            sage: set_random_seed()
             sage: x = random_eja().random_element()
             sage: x.is_invertible() == (x.det() != 0)
             True
             sage: x = random_eja().random_element()
             sage: x.is_invertible() == (x.det() != 0)
             True
@@ -383,15 +376,14 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         Ensure that the determinant is multiplicative on an associative
         subalgebra as in Faraut and Korányi's Proposition II.2.2::
 
         Ensure that the determinant is multiplicative on an associative
         subalgebra as in Faraut and Korányi's Proposition II.2.2::
 
-            sage: set_random_seed()
-            sage: J = random_eja().random_element().subalgebra_generated_by()
+            sage: x0 = random_eja().random_element()
+            sage: J = x0.subalgebra_generated_by(orthonormalize=False)
             sage: x,y = J.random_elements(2)
             sage: (x*y).det() == x.det()*y.det()
             True
 
         The determinant in real matrix algebras is the usual determinant::
 
             sage: x,y = J.random_elements(2)
             sage: (x*y).det() == x.det()*y.det()
             True
 
         The determinant in real matrix algebras is the usual determinant::
 
-            sage: set_random_seed()
             sage: X = matrix.random(QQ,3)
             sage: X = X + X.T
             sage: J1 = RealSymmetricEJA(3)
             sage: X = matrix.random(QQ,3)
             sage: X = X + X.T
             sage: J1 = RealSymmetricEJA(3)
@@ -404,6 +396,22 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: actual2 == expected
             True
 
             sage: actual2 == expected
             True
 
+        There's a formula for the determinant of the Albert algebra
+        (Yokota, Section 2.1)::
+
+            sage: def albert_det(x):
+            ....:     X = x.to_matrix()
+            ....:     res  = X[0,0]*X[1,1]*X[2,2]
+            ....:     res += 2*(X[1,2]*X[2,0]*X[0,1]).real()
+            ....:     res -= X[0,0]*X[1,2]*X[2,1]
+            ....:     res -= X[1,1]*X[2,0]*X[0,2]
+            ....:     res -= X[2,2]*X[0,1]*X[1,0]
+            ....:     return res.leading_coefficient()
+            sage: J = AlbertEJA(field=QQ, orthonormalize=False)
+            sage: xs = J.random_elements(10)
+            sage: all( albert_det(x) == x.det() for x in xs )
+            True
+
         """
         P = self.parent()
         r = P.rank()
         """
         P = self.parent()
         r = P.rank()
@@ -448,7 +456,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         The inverse in the spin factor algebra is given in Alizadeh's
         Example 11.11::
 
         The inverse in the spin factor algebra is given in Alizadeh's
         Example 11.11::
 
-            sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
             sage: x = J.random_element()
             sage: while not x.is_invertible():
             sage: J = JordanSpinEJA.random_instance()
             sage: x = J.random_element()
             sage: while not x.is_invertible():
@@ -473,14 +480,12 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         The identity element is its own inverse::
 
 
         The identity element is its own inverse::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J.one().inverse() == J.one()
             True
 
         If an element has an inverse, it acts like one::
 
             sage: J = random_eja()
             sage: J.one().inverse() == J.one()
             True
 
         If an element has an inverse, it acts like one::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
@@ -488,7 +493,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         The inverse of the inverse is what we started with::
 
 
         The inverse of the inverse is what we started with::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
@@ -498,17 +502,17 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         of an element is the inverse of its left-multiplication operator
         applied to the algebra's identity, when that inverse exists::
 
         of an element is the inverse of its left-multiplication operator
         applied to the algebra's identity, when that inverse exists::
 
-            sage: set_random_seed()
-            sage: J = random_eja()
-            sage: x = J.random_element()
-            sage: (not x.operator().is_invertible()) or (
-            ....:    x.operator().inverse()(J.one()) == x.inverse() )
+            sage: J = random_eja()                         # long time
+            sage: x = J.random_element()                   # long time
+            sage: (not x.operator().is_invertible()) or (  # long time
+            ....:    x.operator().inverse()(J.one())       # long time
+            ....:    ==                                    # long time
+            ....:    x.inverse() )                         # long time
             True
 
         Check that the fast (cached) and slow algorithms give the same
         answer::
 
             True
 
         Check that the fast (cached) and slow algorithms give the same
         answer::
 
-            sage: set_random_seed()                              # long time
             sage: J = random_eja(field=QQ, orthonormalize=False) # long time
             sage: x = J.random_element()                         # long time
             sage: while not x.is_invertible():                   # long time
             sage: J = random_eja(field=QQ, orthonormalize=False) # long time
             sage: x = J.random_element()                         # long time
             sage: while not x.is_invertible():                   # long time
@@ -520,15 +524,18 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             True
         """
         not_invertible_msg = "element is not invertible"
             True
         """
         not_invertible_msg = "element is not invertible"
-        if self.parent()._charpoly_coefficients.is_in_cache():
+
+        algebra = self.parent()
+        if algebra._charpoly_coefficients.is_in_cache():
             # We can invert using our charpoly if it will be fast to
             # compute. If the coefficients are cached, our rank had
             # better be too!
             if self.det().is_zero():
                 raise ZeroDivisionError(not_invertible_msg)
             # We can invert using our charpoly if it will be fast to
             # compute. If the coefficients are cached, our rank had
             # better be too!
             if self.det().is_zero():
                 raise ZeroDivisionError(not_invertible_msg)
-            r = self.parent().rank()
+            r = algebra.rank()
             a = self.characteristic_polynomial().coefficients(sparse=False)
             a = self.characteristic_polynomial().coefficients(sparse=False)
-            return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det()
+            return (-1)**(r+1)*algebra.sum(a[i+1]*self**i
+                                           for i in range(r))/self.det()
 
         try:
             inv = (~self.quadratic_representation())(self)
 
         try:
             inv = (~self.quadratic_representation())(self)
@@ -548,7 +555,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         If computing my determinant will be fast, we do so and compare
         with zero (Proposition II.2.4 in Faraut and
 
         If computing my determinant will be fast, we do so and compare
         with zero (Proposition II.2.4 in Faraut and
-        Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
+        Korányi). Otherwise, Proposition II.3.2 in Faraut and Korányi
         reduces the problem to the invertibility of my quadratic
         representation.
 
         reduces the problem to the invertibility of my quadratic
         representation.
 
@@ -560,14 +567,12 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         The identity element is always invertible::
 
 
         The identity element is always invertible::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J.one().is_invertible()
             True
 
         The zero element is never invertible in a non-trivial algebra::
 
             sage: J = random_eja()
             sage: J.one().is_invertible()
             True
 
         The zero element is never invertible in a non-trivial algebra::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: (not J.is_trivial()) and J.zero().is_invertible()
             False
             sage: J = random_eja()
             sage: (not J.is_trivial()) and J.zero().is_invertible()
             False
@@ -575,7 +580,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         Test that the fast (cached) and slow algorithms give the same
         answer::
 
         Test that the fast (cached) and slow algorithms give the same
         answer::
 
-            sage: set_random_seed()                              # long time
             sage: J = random_eja(field=QQ, orthonormalize=False) # long time
             sage: x = J.random_element()                         # long time
             sage: slow = x.is_invertible()                       # long time
             sage: J = random_eja(field=QQ, orthonormalize=False) # long time
             sage: x = J.random_element()                         # long time
             sage: slow = x.is_invertible()                       # long time
@@ -658,14 +662,12 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         The identity element is minimal only in an EJA of rank one::
 
 
         The identity element is minimal only in an EJA of rank one::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
             True
 
         A non-idempotent cannot be a minimal idempotent::
 
             sage: J = random_eja()
             sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
             True
 
         A non-idempotent cannot be a minimal idempotent::
 
-            sage: set_random_seed()
             sage: J = JordanSpinEJA(4)
             sage: x = J.random_element()
             sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
             sage: J = JordanSpinEJA(4)
             sage: x = J.random_element()
             sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
@@ -675,7 +677,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         idempotent if and only if it's idempotent with trace equal to
         unity::
 
         idempotent if and only if it's idempotent with trace equal to
         unity::
 
-            sage: set_random_seed()
             sage: J = JordanSpinEJA(4)
             sage: x = J.random_element()
             sage: expected = (x.is_idempotent() and x.trace() == 1)
             sage: J = JordanSpinEJA(4)
             sage: x = J.random_element()
             sage: expected = (x.is_idempotent() and x.trace() == 1)
@@ -685,7 +686,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         Primitive idempotents must be non-zero::
 
 
         Primitive idempotents must be non-zero::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J.zero().is_idempotent()
             True
             sage: J = random_eja()
             sage: J.zero().is_idempotent()
             True
@@ -742,14 +742,12 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         The identity element is never nilpotent, except in a trivial EJA::
 
 
         The identity element is never nilpotent, except in a trivial EJA::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J.one().is_nilpotent() and not J.is_trivial()
             False
 
         The additive identity is always nilpotent::
 
             sage: J = random_eja()
             sage: J.one().is_nilpotent() and not J.is_trivial()
             False
 
         The additive identity is always nilpotent::
 
-            sage: set_random_seed()
             sage: random_eja().zero().is_nilpotent()
             True
 
             sage: random_eja().zero().is_nilpotent()
             True
 
@@ -792,7 +790,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         The zero element should never be regular, unless the parent
         algebra has dimension less than or equal to one::
 
         The zero element should never be regular, unless the parent
         algebra has dimension less than or equal to one::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J.dimension() <= 1 or not J.zero().is_regular()
             True
             sage: J = random_eja()
             sage: J.dimension() <= 1 or not J.zero().is_regular()
             True
@@ -800,7 +797,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         The unit element isn't regular unless the algebra happens to
         consist of only its scalar multiples::
 
         The unit element isn't regular unless the algebra happens to
         consist of only its scalar multiples::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J.dimension() <= 1 or not J.one().is_regular()
             True
             sage: J = random_eja()
             sage: J.dimension() <= 1 or not J.one().is_regular()
             True
@@ -816,7 +812,23 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         ALGORITHM:
 
 
         ALGORITHM:
 
-        .........
+        First we handle the special cases where the algebra is
+        trivial, this element is zero, or the dimension of the algebra
+        is one and this element is not zero. With those out of the
+        way, we may assume that ``self`` is nonzero, the algebra is
+        nontrivial, and that the dimension of the algebra is at least
+        two.
+
+        Beginning with the algebra's unit element (power zero), we add
+        successive (basis representations of) powers of this element
+        to a matrix, row-reducing at each step. After row-reducing, we
+        check the rank of the matrix. If adding a row and row-reducing
+        does not increase the rank of the matrix at any point, the row
+        we've just added lives in the span of the previous ones; thus
+        the corresponding power of ``self`` lives in the span of its
+        lesser powers. When that happens, the degree of the minimal
+        polynomial is the rank of the matrix; if it never happens, the
+        degree must be the dimension of the entire space.
 
         SETUP::
 
 
         SETUP::
 
@@ -835,7 +847,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         In the spin factor algebra (of rank two), all elements that
         aren't multiples of the identity are regular::
 
         In the spin factor algebra (of rank two), all elements that
         aren't multiples of the identity are regular::
 
-            sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
             sage: J = JordanSpinEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
@@ -847,7 +858,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         The zero and unit elements are both of degree one in nontrivial
         algebras::
 
         The zero and unit elements are both of degree one in nontrivial
         algebras::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: d = J.zero().degree()
             sage: (J.is_trivial() and d == 0) or d == 1
             sage: J = random_eja()
             sage: d = J.zero().degree()
             sage: (J.is_trivial() and d == 0) or d == 1
@@ -858,11 +868,9 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         Our implementation agrees with the definition::
 
 
         Our implementation agrees with the definition::
 
-            sage: set_random_seed()
             sage: x = random_eja().random_element()
             sage: x.degree() == x.minimal_polynomial().degree()
             True
             sage: x = random_eja().random_element()
             sage: x.degree() == x.minimal_polynomial().degree()
             True
-
         """
         n = self.parent().dimension()
 
         """
         n = self.parent().dimension()
 
@@ -967,7 +975,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         always the same, except in trivial algebras where the minimal
         polynomial of the unit/zero element is ``1``::
 
         always the same, except in trivial algebras where the minimal
         polynomial of the unit/zero element is ``1``::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: mu = J.one().minimal_polynomial()
             sage: t = mu.parent().gen()
             sage: J = random_eja()
             sage: mu = J.one().minimal_polynomial()
             sage: t = mu.parent().gen()
@@ -981,7 +988,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         The degree of an element is (by one definition) the degree
         of its minimal polynomial::
 
         The degree of an element is (by one definition) the degree
         of its minimal polynomial::
 
-            sage: set_random_seed()
             sage: x = random_eja().random_element()
             sage: x.degree() == x.minimal_polynomial().degree()
             True
             sage: x = random_eja().random_element()
             sage: x.degree() == x.minimal_polynomial().degree()
             True
@@ -992,9 +998,8 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         identity. We require the dimension of the algebra to be at least
         two here so that said elements actually exist::
 
         identity. We require the dimension of the algebra to be at least
         two here so that said elements actually exist::
 
-            sage: set_random_seed()
-            sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())
-            sage: n = ZZ.random_element(2, n_max)
+            sage: d_max = JordanSpinEJA._max_random_instance_dimension()
+            sage: n = ZZ.random_element(2, max(2,d_max))
             sage: J = JordanSpinEJA(n)
             sage: y = J.random_element()
             sage: while y == y.coefficient(0)*J.one():
             sage: J = JordanSpinEJA(n)
             sage: y = J.random_element()
             sage: while y == y.coefficient(0)*J.one():
@@ -1009,18 +1014,17 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         The minimal polynomial should always kill its element::
 
 
         The minimal polynomial should always kill its element::
 
-            sage: set_random_seed()
-            sage: x = random_eja().random_element()
-            sage: p = x.minimal_polynomial()
-            sage: x.apply_univariate_polynomial(p)
+            sage: x = random_eja().random_element()  # long time
+            sage: p = x.minimal_polynomial()         # long time
+            sage: x.apply_univariate_polynomial(p)   # long time
             0
 
         The minimal polynomial is invariant under a change of basis,
         and in particular, a re-scaling of the basis::
 
             0
 
         The minimal polynomial is invariant under a change of basis,
         and in particular, a re-scaling of the basis::
 
-            sage: set_random_seed()
-            sage: n_max = RealSymmetricEJA._max_random_instance_size()
-            sage: n = ZZ.random_element(1, n_max)
+            sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
+            sage: d = ZZ.random_element(1, d_max)
+            sage: n = RealSymmetricEJA._max_random_instance_size(d)
             sage: J1 = RealSymmetricEJA(n)
             sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
             sage: X = random_matrix(AA,n)
             sage: J1 = RealSymmetricEJA(n)
             sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
             sage: X = random_matrix(AA,n)
@@ -1119,18 +1123,10 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         B = self.parent().matrix_basis()
         W = self.parent().matrix_space()
 
         B = self.parent().matrix_basis()
         W = self.parent().matrix_space()
 
-        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()) )
+        # 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()) )
 
 
 
 
 
 
@@ -1172,7 +1168,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         TESTS::
 
 
         TESTS::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x,y = J.random_elements(2)
             sage: x.operator()(y) == x*y
             sage: J = random_eja()
             sage: x,y = J.random_elements(2)
             sage: x.operator()(y) == x*y
@@ -1184,7 +1179,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         P = self.parent()
         left_mult_by_self = lambda y: self*y
         L = P.module_morphism(function=left_mult_by_self, codomain=P)
         P = self.parent()
         left_mult_by_self = lambda y: self*y
         L = P.module_morphism(function=left_mult_by_self, codomain=P)
-        return FiniteDimensionalEJAOperator(P, P, L.matrix() )
+        return EJAOperator(P, P, L.matrix() )
 
 
     def quadratic_representation(self, other=None):
 
 
     def quadratic_representation(self, other=None):
@@ -1201,7 +1196,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         The explicit form in the spin factor algebra is given by
         Alizadeh's Example 11.12::
 
         The explicit form in the spin factor algebra is given by
         Alizadeh's Example 11.12::
 
-            sage: set_random_seed()
             sage: x = JordanSpinEJA.random_instance().random_element()
             sage: x_vec = x.to_vector()
             sage: Q = matrix.identity(x.base_ring(), 0)
             sage: x = JordanSpinEJA.random_instance().random_element()
             sage: x_vec = x.to_vector()
             sage: Q = matrix.identity(x.base_ring(), 0)
@@ -1221,7 +1215,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         Test all of the properties from Theorem 11.2 in Alizadeh::
 
 
         Test all of the properties from Theorem 11.2 in Alizadeh::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x,y = J.random_elements(2)
             sage: Lx = x.operator()
             sage: J = random_eja()
             sage: x,y = J.random_elements(2)
             sage: Lx = x.operator()
@@ -1373,7 +1366,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: (J0, J5, J1) = J.peirce_decomposition(c1)
             sage: (f0, f1, f2) = J1.gens()
             sage: f0.spectral_decomposition()
             sage: (J0, J5, J1) = J.peirce_decomposition(c1)
             sage: (f0, f1, f2) = J1.gens()
             sage: f0.spectral_decomposition()
-            [(0, c2), (1, c0)]
+            [(0, 1.000000000000000?*c2), (1, 1.000000000000000?*c0)]
 
         """
         A = self.subalgebra_generated_by(orthonormalize=True)
 
         """
         A = self.subalgebra_generated_by(orthonormalize=True)
@@ -1415,9 +1408,8 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         This subalgebra, being composed of only powers, is associative::
 
 
         This subalgebra, being composed of only powers, is associative::
 
-            sage: set_random_seed()
             sage: x0 = random_eja().random_element()
             sage: x0 = random_eja().random_element()
-            sage: A = x0.subalgebra_generated_by()
+            sage: A = x0.subalgebra_generated_by(orthonormalize=False)
             sage: x,y,z = A.random_elements(3)
             sage: (x*y)*z == x*(y*z)
             True
             sage: x,y,z = A.random_elements(3)
             sage: (x*y)*z == x*(y*z)
             True
@@ -1425,9 +1417,8 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         Squaring in the subalgebra should work the same as in
         the superalgebra::
 
         Squaring in the subalgebra should work the same as in
         the superalgebra::
 
-            sage: set_random_seed()
             sage: x = random_eja().random_element()
             sage: x = random_eja().random_element()
-            sage: A = x.subalgebra_generated_by()
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: A(x^2) == A(x)*A(x)
             True
 
             sage: A(x^2) == A(x)*A(x)
             True
 
@@ -1436,7 +1427,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         element... unless the original algebra was trivial, in which
         case the subalgebra is trivial too::
 
         element... unless the original algebra was trivial, in which
         case the subalgebra is trivial too::
 
-            sage: set_random_seed()
             sage: A = random_eja().zero().subalgebra_generated_by()
             sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
             True
             sage: A = random_eja().zero().subalgebra_generated_by()
             sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
             True
@@ -1467,8 +1457,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         where there are non-nilpotent elements, or that we get the dumb
         solution in the trivial algebra::
 
         where there are non-nilpotent elements, or that we get the dumb
         solution in the trivial algebra::
 
-            sage: set_random_seed()
-            sage: J = random_eja()
+            sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: x = J.random_element()
             sage: while x.is_nilpotent() and not J.is_trivial():
             ....:     x = J.random_element()
             sage: x = J.random_element()
             sage: while x.is_nilpotent() and not J.is_trivial():
             ....:     x = J.random_element()
@@ -1483,7 +1472,10 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         if self.is_nilpotent():
             raise ValueError("this only works with non-nilpotent elements!")
 
         if self.is_nilpotent():
             raise ValueError("this only works with non-nilpotent elements!")
 
-        J = self.subalgebra_generated_by()
+        # The subalgebra is transient (we return an element of the
+        # superalgebra, i.e. this algebra) so why bother
+        # orthonormalizing?
+        J = self.subalgebra_generated_by(orthonormalize=False)
         u = J(self)
 
         # The image of the matrix of left-u^m-multiplication
         u = J(self)
 
         # The image of the matrix of left-u^m-multiplication
@@ -1504,14 +1496,12 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         # subspace... or do we? Can't we just solve, knowing that
         # A(c) = u^(s+1) should have a solution in the big space,
         # too?
         # subspace... or do we? Can't we just solve, knowing that
         # A(c) = u^(s+1) should have a solution in the big space,
         # too?
-        #
-        # Beware, solve_right() means that we're using COLUMN vectors.
-        # Our FiniteDimensionalAlgebraElement superclass uses rows.
         u_next = u**(s+1)
         A = u_next.operator().matrix()
         c = J.from_vector(A.solve_right(u_next.to_vector()))
 
         u_next = u**(s+1)
         A = u_next.operator().matrix()
         c = J.from_vector(A.solve_right(u_next.to_vector()))
 
-        # Now c is the idempotent we want, but it still lives in the subalgebra.
+        # Now c is the idempotent we want, but it still lives in
+        # the subalgebra.
         return c.superalgebra_element()
 
 
         return c.superalgebra_element()
 
 
@@ -1551,20 +1541,24 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         The trace of an element is a real number::
 
 
         The trace of an element is a real number::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J.random_element().trace() in RLF
             True
 
         The trace is linear::
 
             sage: J = random_eja()
             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
 
             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
 
+        The trace of a square is nonnegative::
+
+            sage: x = random_eja().random_element()
+            sage: (x*x).trace() >= 0
+            True
+
         """
         P = self.parent()
         r = P.rank()
         """
         P = self.parent()
         r = P.rank()
@@ -1594,14 +1588,13 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         The trace inner product is commutative, bilinear, and associative::
 
 
         The trace inner product is commutative, bilinear, and associative::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: # commutative
             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
             True
             sage: # bilinear
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: # commutative
             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
             True
             sage: # bilinear
-            sage: a = J.base_ring().random_element();
+            sage: a = J.base_ring().random_element()
             sage: actual = (a*(x+z)).trace_inner_product(y)
             sage: expected = ( a*x.trace_inner_product(y) +
             ....:              a*z.trace_inner_product(y) )
             sage: actual = (a*(x+z)).trace_inner_product(y)
             sage: expected = ( a*x.trace_inner_product(y) +
             ....:              a*z.trace_inner_product(y) )
@@ -1648,3 +1641,187 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         """
         return self.trace_inner_product(self).sqrt()
 
         """
         return self.trace_inner_product(self).sqrt()
+
+
+    def operator_trace_inner_product(self, other):
+        r"""
+        Return the operator inner product of myself and ``other``.
+
+        The "operator inner product," whose name is not standard, is
+        defined be the usual linear-algebraic trace of the
+        ``(x*y).operator()``.
+
+        Proposition III.1.5 in Faraut and Korányi shows that on any
+        Euclidean Jordan algebra, this is another associative inner
+        product under which the cone of squares is symmetric.
+
+        This works even if the basis hasn't been orthonormalized
+        because the eigenvalues of the corresponding matrix don't
+        change when the basis does (they're preserved by any
+        similarity transformation).
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  random_eja)
+
+        EXAMPLES:
+
+        Proposition III.4.2 of Faraut and Korányi shows that on a
+        simple algebra of rank `r` and dimension `n`, this inner
+        product is `n/r` times the canonical
+        :meth:`trace_inner_product`::
+
+            sage: J = JordanSpinEJA(4, field=QQ)
+            sage: x,y = J.random_elements(2)
+            sage: n = J.dimension()
+            sage: r = J.rank()
+            sage: actual = x.operator_trace_inner_product(y)
+            sage: expected = (n/r)*x.trace_inner_product(y)
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: x,y = J.random_elements(2)
+            sage: n = J.dimension()
+            sage: r = J.rank()
+            sage: actual = x.operator_trace_inner_product(y)
+            sage: expected = (n/r)*x.trace_inner_product(y)
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
+            sage: x,y = J.random_elements(2)
+            sage: n = J.dimension()
+            sage: r = J.rank()
+            sage: actual = x.operator_trace_inner_product(y)
+            sage: expected = (n/r)*x.trace_inner_product(y)
+            sage: actual == expected
+            True
+
+        TESTS:
+
+        The operator inner product is commutative, bilinear, and
+        associative::
+
+            sage: J = random_eja()
+            sage: x,y,z = J.random_elements(3)
+            sage: # commutative
+            sage: actual = x.operator_trace_inner_product(y)
+            sage: expected = y.operator_trace_inner_product(x)
+            sage: actual == expected
+            True
+            sage: # bilinear
+            sage: a = J.base_ring().random_element()
+            sage: actual = (a*(x+z)).operator_trace_inner_product(y)
+            sage: expected = ( a*x.operator_trace_inner_product(y) +
+            ....:              a*z.operator_trace_inner_product(y) )
+            sage: actual == expected
+            True
+            sage: actual = x.operator_trace_inner_product(a*(y+z))
+            sage: expected = ( a*x.operator_trace_inner_product(y) +
+            ....:              a*x.operator_trace_inner_product(z) )
+            sage: actual == expected
+            True
+            sage: # associative
+            sage: actual = (x*y).operator_trace_inner_product(z)
+            sage: expected = y.operator_trace_inner_product(x*z)
+            sage: actual == expected
+            True
+
+        Despite the fact that the implementation uses a matrix representation,
+        the answer is independent of the basis used::
+
+            sage: J = RealSymmetricEJA(3, field=QQ, orthonormalize=False)
+            sage: V = RealSymmetricEJA(3)
+            sage: x,y = J.random_elements(2)
+            sage: w = V(x.to_matrix())
+            sage: z = V(y.to_matrix())
+            sage: expected = x.operator_trace_inner_product(y)
+            sage: actual = w.operator_trace_inner_product(z)
+            sage: actual == expected
+            True
+
+        """
+        if not other in self.parent():
+            raise TypeError("'other' must live in the same algebra")
+
+        return (self*other).operator().matrix().trace()
+
+
+    def operator_trace_norm(self):
+        """
+        The norm of this element with respect to
+        :meth:`operator_trace_inner_product`.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  HadamardEJA)
+
+        EXAMPLES:
+
+        On a simple algebra, this will differ from :meth:`trace_norm`
+        by the scalar factor ``(n/r).sqrt()``, where `n` is the
+        dimension of the algebra and `r` its rank. This follows from
+        the corresponding result (Proposition III.4.2 of Faraut and
+        Korányi) for the trace inner product::
+
+            sage: J = HadamardEJA(2)
+            sage: x = sum(J.gens())
+            sage: x.operator_trace_norm()
+            1.414213562373095?
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: x = sum(J.gens())
+            sage: x.operator_trace_norm()
+            4
+
+        """
+        return self.operator_trace_inner_product(self).sqrt()
+
+
+class CartesianProductParentEJAElement(EJAElement):
+    r"""
+    An intermediate class for elements that have a Cartesian
+    product as their parent algebra.
+
+    This is needed because the ``to_matrix`` method (which gives you a
+    representation from the superalgebra) needs to do special stuff
+    for Cartesian products. Specifically, an EJA subalgebra of a
+    Cartesian product EJA will not itself be a Cartesian product (it
+    has its own basis) -- but we want ``to_matrix()`` to be able to
+    give us a Cartesian product representation.
+    """
+    def to_matrix(self):
+        # An override is necessary to call our custom _scale().
+        B = self.parent().matrix_basis()
+        W = self.parent().matrix_space()
+
+        # Aaaaand linear combinations don't work in Cartesian
+        # product spaces, even though they provide a method with
+        # that name. This is hidden in a subclass because the
+        # _scale() function is slow.
+        pairs = zip(B, self.to_vector())
+        return W.sum( _scale(b, alpha) for (b,alpha) in pairs )
+
+class CartesianProductEJAElement(CartesianProductParentEJAElement):
+    def det(self):
+        r"""
+        Compute the determinant of this product-element using the
+        determianants of its factors.
+
+        This result Follows from the spectral decomposition of (say)
+        the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
+        0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
+        """
+        from sage.misc.misc_c import prod
+        return prod( f.det() for f in self.cartesian_factors() )