]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
COPYING,LICENSE: add (AGPL-3.0+)
[sage.d.git] / mjo / eja / eja_algebra.py
index 907f40d72135df42d211cdd9bb8923df6bf462e6..adcc3436b1302e09cd20007d0525aee08e32a48f 100644 (file)
@@ -1,4 +1,4 @@
-"""
+r"""
 Representations and constructions for Euclidean Jordan algebras.
 
 A Euclidean Jordan algebra is a Jordan algebra that has some
 Representations and constructions for Euclidean Jordan algebras.
 
 A Euclidean Jordan algebra is a Jordan algebra that has some
@@ -34,12 +34,13 @@ for these simple algebras:
   * :class:`QuaternionHermitianEJA`
   * :class:`OctonionHermitianEJA`
 
   * :class:`QuaternionHermitianEJA`
   * :class:`OctonionHermitianEJA`
 
-In addition to these, we provide two other example constructions,
+In addition to these, we provide a few other example constructions,
 
   * :class:`JordanSpinEJA`
   * :class:`HadamardEJA`
   * :class:`AlbertEJA`
   * :class:`TrivialEJA`
 
   * :class:`JordanSpinEJA`
   * :class:`HadamardEJA`
   * :class:`AlbertEJA`
   * :class:`TrivialEJA`
+  * :class:`ComplexSkewSymmetricEJA`
 
 The Jordan spin algebra is a bilinear form algebra where the bilinear
 form is the identity. The Hadamard EJA is simply a Cartesian product
 
 The Jordan spin algebra is a bilinear form algebra where the bilinear
 form is the identity. The Hadamard EJA is simply a Cartesian product
@@ -71,18 +72,18 @@ matrix, whereas the inner product must return a scalar. Our basis for
 the one-by-one matrices is of course the set consisting of a single
 matrix with its sole entry non-zero::
 
 the one-by-one matrices is of course the set consisting of a single
 matrix with its sole entry non-zero::
 
-    sage: from mjo.eja.eja_algebra import FiniteDimensionalEJA
+    sage: from mjo.eja.eja_algebra import EJA
     sage: jp = lambda X,Y: X*Y
     sage: ip = lambda X,Y: X[0,0]*Y[0,0]
     sage: b1 = matrix(AA, [[1]])
     sage: jp = lambda X,Y: X*Y
     sage: ip = lambda X,Y: X[0,0]*Y[0,0]
     sage: b1 = matrix(AA, [[1]])
-    sage: J1 = FiniteDimensionalEJA((b1,), jp, ip)
+    sage: J1 = EJA((b1,), jp, ip)
     sage: J1
     Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
 
 In fact, any positive scalar multiple of that inner-product would work::
 
     sage: ip2 = lambda X,Y: 16*ip(X,Y)
     sage: J1
     Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
 
 In fact, any positive scalar multiple of that inner-product would work::
 
     sage: ip2 = lambda X,Y: 16*ip(X,Y)
-    sage: J2 = FiniteDimensionalEJA((b1,), jp, ip2)
+    sage: J2 = EJA((b1,), jp, ip2)
     sage: J2
     Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
 
     sage: J2
     Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
 
@@ -90,7 +91,7 @@ But beware that your basis will be orthonormalized _with respect to the
 given inner-product_ unless you pass ``orthonormalize=False`` to the
 constructor. For example::
 
 given inner-product_ unless you pass ``orthonormalize=False`` to the
 constructor. For example::
 
-    sage: J3 = FiniteDimensionalEJA((b1,), jp, ip2, orthonormalize=False)
+    sage: J3 = EJA((b1,), jp, ip2, orthonormalize=False)
     sage: J3
     Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
 
     sage: J3
     Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
 
@@ -117,7 +118,7 @@ Another option for your basis is to use elemebts of a
 
     sage: from mjo.matrix_algebra import MatrixAlgebra
     sage: A = MatrixAlgebra(1,AA,AA)
 
     sage: from mjo.matrix_algebra import MatrixAlgebra
     sage: A = MatrixAlgebra(1,AA,AA)
-    sage: J4 = FiniteDimensionalEJA(A.gens(), jp, ip)
+    sage: J4 = EJA(A.gens(), jp, ip)
     sage: J4
     Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
     sage: J4.basis()[0].to_matrix()
     sage: J4
     Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
     sage: J4.basis()[0].to_matrix()
@@ -166,8 +167,9 @@ from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
-from mjo.eja.eja_element import FiniteDimensionalEJAElement
-from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
+from mjo.eja.eja_element import (CartesianProductEJAElement,
+                                 EJAElement)
+from mjo.eja.eja_operator import EJAOperator
 from mjo.eja.eja_utils import _all2list
 
 def EuclideanJordanAlgebras(field):
 from mjo.eja.eja_utils import _all2list
 
 def EuclideanJordanAlgebras(field):
@@ -181,7 +183,7 @@ def EuclideanJordanAlgebras(field):
     category = category.WithBasis().Unital().Commutative()
     return category
 
     category = category.WithBasis().Unital().Commutative()
     return category
 
-class FiniteDimensionalEJA(CombinatorialFreeModule):
+class EJA(CombinatorialFreeModule):
     r"""
     A finite-dimensional Euclidean Jordan algebra.
 
     r"""
     A finite-dimensional Euclidean Jordan algebra.
 
@@ -229,7 +231,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
     We should compute that an element subalgebra is associative even
     if we circumvent the element method::
 
     We should compute that an element subalgebra is associative even
     if we circumvent the element method::
 
-        sage: set_random_seed()
         sage: J = random_eja(field=QQ,orthonormalize=False)
         sage: x = J.random_element()
         sage: A = x.subalgebra_generated_by(orthonormalize=False)
         sage: J = random_eja(field=QQ,orthonormalize=False)
         sage: x = J.random_element()
         sage: A = x.subalgebra_generated_by(orthonormalize=False)
@@ -237,7 +238,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         sage: J.subalgebra(basis, orthonormalize=False).is_associative()
         True
     """
         sage: J.subalgebra(basis, orthonormalize=False).is_associative()
         True
     """
-    Element = FiniteDimensionalEJAElement
+    Element = EJAElement
 
     @staticmethod
     def _check_input_field(field):
 
     @staticmethod
     def _check_input_field(field):
@@ -366,7 +367,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         if orthonormalize:
             # Now "self._matrix_span" is the vector space of our
 
         if orthonormalize:
             # Now "self._matrix_span" is the vector space of our
-            # algebra coordinates. The variables "X1", "X2",...  refer
+            # algebra coordinates. The variables "X0", "X1",...  refer
             # to the entries of vectors in self._matrix_span. Thus to
             # convert back and forth between the orthonormal
             # coordinates and the given ones, we need to stick the
             # to the entries of vectors in self._matrix_span. Thus to
             # convert back and forth between the orthonormal
             # coordinates and the given ones, we need to stick the
@@ -431,7 +432,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         TESTS::
 
 
         TESTS::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J(1)
             Traceback (most recent call last):
             sage: J = random_eja()
             sage: J(1)
             Traceback (most recent call last):
@@ -456,7 +456,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         TESTS::
 
 
         TESTS::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: n = J.dimension()
             sage: bi = J.zero()
             sage: J = random_eja()
             sage: n = J.dimension()
             sage: bi = J.zero()
@@ -498,7 +497,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Our inner product is "associative," which means the following for
         a symmetric bilinear form::
 
         Our inner product is "associative," which means the following for
         a symmetric bilinear form::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
@@ -509,7 +507,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Ensure that this is the usual inner product for the algebras
         over `R^n`::
 
         Ensure that this is the usual inner product for the algebras
         over `R^n`::
 
-            sage: set_random_seed()
             sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
             sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
@@ -522,7 +519,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         one). This is in Faraut and Koranyi, and also my "On the
         symmetry..." paper::
 
         one). This is in Faraut and Koranyi, and also my "On the
         symmetry..." paper::
 
-            sage: set_random_seed()
             sage: J = BilinearFormEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
             sage: J = BilinearFormEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
@@ -635,7 +631,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         The values we've presupplied to the constructors agree with
         the computation::
 
         The values we've presupplied to the constructors agree with
         the computation::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J.is_associative() == J._jordan_product_is_associative()
             True
             sage: J = random_eja()
             sage: J.is_associative() == J._jordan_product_is_associative()
             True
@@ -757,7 +752,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Ensure that we can convert any element back and forth
         faithfully between its matrix and algebra representations::
 
         Ensure that we can convert any element back and forth
         faithfully between its matrix and algebra representations::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J(x.to_matrix()) == x
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J(x.to_matrix()) == x
@@ -870,7 +864,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = JordanSpinEJA(3)
             sage: p = J.characteristic_polynomial_of(); p
 
             sage: J = JordanSpinEJA(3)
             sage: p = J.characteristic_polynomial_of(); p
-            X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
+            X0^2 - X1^2 - X2^2 + (-2*t)*X0 + t^2
             sage: xvec = J.one().to_vector()
             sage: p(*xvec)
             t^2 - 2*t + 1
             sage: xvec = J.one().to_vector()
             sage: p(*xvec)
             t^2 - 2*t + 1
@@ -919,13 +913,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = HadamardEJA(2)
             sage: J.coordinate_polynomial_ring()
 
             sage: J = HadamardEJA(2)
             sage: J.coordinate_polynomial_ring()
-            Multivariate Polynomial Ring in X1, X2...
+            Multivariate Polynomial Ring in X0, X1...
             sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
             sage: J.coordinate_polynomial_ring()
             sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
             sage: J.coordinate_polynomial_ring()
-            Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
+            Multivariate Polynomial Ring in X0, X1, X2, X3, X4, X5...
 
         """
 
         """
-        var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
+        var_names = tuple( "X%d" % z for z in range(self.dimension()) )
         return PolynomialRing(self.base_ring(), var_names)
 
     def inner_product(self, x, y):
         return PolynomialRing(self.base_ring(), var_names)
 
     def inner_product(self, x, y):
@@ -947,7 +941,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Our inner product is "associative," which means the following for
         a symmetric bilinear form::
 
         Our inner product is "associative," which means the following for
         a symmetric bilinear form::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
@@ -958,7 +951,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Ensure that this is the usual inner product for the algebras
         over `R^n`::
 
         Ensure that this is the usual inner product for the algebras
         over `R^n`::
 
-            sage: set_random_seed()
             sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
             sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
@@ -971,7 +963,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         one). This is in Faraut and Koranyi, and also my "On the
         symmetry..." paper::
 
         one). This is in Faraut and Koranyi, and also my "On the
         symmetry..." paper::
 
-            sage: set_random_seed()
             sage: J = BilinearFormEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
             sage: J = BilinearFormEJA.random_instance()
             sage: n = J.dimension()
             sage: x = J.random_element()
@@ -1199,19 +1190,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         The identity element acts like the identity, regardless of
         whether or not we orthonormalize::
 
         The identity element acts like the identity, regardless of
         whether or not we orthonormalize::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             True
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             True
-            sage: A = x.subalgebra_generated_by()
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: y = A.random_element()
             sage: A.one()*y == y and y*A.one() == y
             True
 
         ::
 
             sage: y = A.random_element()
             sage: A.one()*y == y and y*A.one() == y
             True
 
         ::
 
-            sage: set_random_seed()
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: x = J.random_element()
             sage: J.one()*x == x and x*J.one() == x
@@ -1225,14 +1214,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         regardless of the base field and whether or not we
         orthonormalize::
 
         regardless of the base field and whether or not we
         orthonormalize::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: actual = J.one().operator().matrix()
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: actual == expected
             True
             sage: x = J.random_element()
             sage: J = random_eja()
             sage: actual = J.one().operator().matrix()
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: actual == expected
             True
             sage: x = J.random_element()
-            sage: A = x.subalgebra_generated_by()
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
             sage: actual = A.one().operator().matrix()
             sage: expected = matrix.identity(A.base_ring(), A.dimension())
             sage: actual == expected
             sage: actual = A.one().operator().matrix()
             sage: expected = matrix.identity(A.base_ring(), A.dimension())
             sage: actual == expected
@@ -1240,7 +1228,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         ::
 
 
         ::
 
-            sage: set_random_seed()
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: actual = J.one().operator().matrix()
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: actual = J.one().operator().matrix()
             sage: expected = matrix.identity(J.base_ring(), J.dimension())
@@ -1256,7 +1243,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Ensure that the cached unit element (often precomputed by
         hand) agrees with the computed one::
 
         Ensure that the cached unit element (often precomputed by
         hand) agrees with the computed one::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: cached = J.one()
             sage: J.one.clear_cache()
             sage: J = random_eja()
             sage: cached = J.one()
             sage: J.one.clear_cache()
@@ -1265,7 +1251,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         ::
 
 
         ::
 
-            sage: set_random_seed()
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: cached = J.one()
             sage: J.one.clear_cache()
             sage: J = random_eja(field=QQ, orthonormalize=False)
             sage: cached = J.one()
             sage: J.one.clear_cache()
@@ -1378,7 +1363,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Every algebra decomposes trivially with respect to its identity
         element::
 
         Every algebra decomposes trivially with respect to its identity
         element::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J0,J5,J1 = J.peirce_decomposition(J.one())
             sage: J0.dimension() == 0 and J5.dimension() == 0
             sage: J = random_eja()
             sage: J0,J5,J1 = J.peirce_decomposition(J.one())
             sage: J0.dimension() == 0 and J5.dimension() == 0
@@ -1391,7 +1375,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         elements in the two subalgebras are the projections onto their
         respective subspaces of the superalgebra's identity element::
 
         elements in the two subalgebras are the projections onto their
         respective subspaces of the superalgebra's identity element::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: if not J.is_trivial():
             sage: J = random_eja()
             sage: x = J.random_element()
             sage: if not J.is_trivial():
@@ -1465,26 +1448,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # For a general base ring... maybe we can trust this to do the
         # right thing? Unlikely, but.
         V = self.vector_space()
         # For a general base ring... maybe we can trust this to do the
         # right thing? Unlikely, but.
         V = self.vector_space()
-        v = V.random_element()
-
-        if self.base_ring() is AA:
-            # The "random element" method of the algebraic reals is
-            # stupid at the moment, and only returns integers between
-            # -2 and 2, inclusive:
-            #
-            #   https://trac.sagemath.org/ticket/30875
-            #
-            # Instead, we implement our own "random vector" method,
-            # and then coerce that into the algebra. We use the vector
-            # space degree here instead of the dimension because a
-            # subalgebra could (for example) be spanned by only two
-            # vectors, each with five coordinates.  We need to
-            # generate all five coordinates.
-            if thorough:
-                v *= QQbar.random_element().real()
-            else:
-                v *= QQ.random_element()
+        if self.base_ring() is AA and not thorough:
+            # Now that AA generates actually random random elements
+            # (post Trac 30875), we only need to de-thorough the
+            # randomness when asked to.
+            V = V.change_ring(QQ)
 
 
+        v = V.random_element()
         return self.from_vector(V.coordinate_vector(v))
 
     def random_elements(self, count, thorough=False):
         return self.from_vector(V.coordinate_vector(v))
 
     def random_elements(self, count, thorough=False):
@@ -1517,6 +1487,64 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                       for idx in range(count) )
 
 
                       for idx in range(count) )
 
 
+    def operator_polynomial_matrix(self):
+        r"""
+        Return the matrix of polynomials (over this algebra's
+        :meth:`coordinate_polynomial_ring`) that, when evaluated at
+        the basis coordinates of an element `x`, produces the basis
+        representation of `L_{x}`.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  JordanSpinEJA)
+
+        EXAMPLES::
+
+            sage: J = HadamardEJA(4)
+            sage: L_x = J.operator_polynomial_matrix()
+            sage: L_x
+            [X0  0  0  0]
+            [ 0 X1  0  0]
+            [ 0  0 X2  0]
+            [ 0  0  0 X3]
+            sage: x = J.one()
+            sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
+            sage: L_x.subs(dict(d))
+            [1 0 0 0]
+            [0 1 0 0]
+            [0 0 1 0]
+            [0 0 0 1]
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: L_x = J.operator_polynomial_matrix()
+            sage: L_x
+            [X0 X1 X2 X3]
+            [X1 X0  0  0]
+            [X2  0 X0  0]
+            [X3  0  0 X0]
+            sage: x = J.one()
+            sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
+            sage: L_x.subs(dict(d))
+            [1 0 0 0]
+            [0 1 0 0]
+            [0 0 1 0]
+            [0 0 0 1]
+
+        """
+        R = self.coordinate_polynomial_ring()
+
+        def L_x_i_j(i,j):
+            # From a result in my book, these are the entries of the
+            # basis representation of L_x.
+            return sum( v*self.monomial(k).operator().matrix()[i,j]
+                        for (k,v) in enumerate(R.gens()) )
+
+        n = self.dimension()
+        return matrix(R, n, n, L_x_i_j)
+
     @cached_method
     def _charpoly_coefficients(self):
         r"""
     @cached_method
     def _charpoly_coefficients(self):
         r"""
@@ -1532,7 +1560,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         The theory shows that these are all homogeneous polynomials of
         a known degree::
 
         The theory shows that these are all homogeneous polynomials of
         a known degree::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
             True
             sage: J = random_eja()
             sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
             True
@@ -1540,16 +1567,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()
         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()
-        vars = R.gens()
         F = R.fraction_field()
 
         F = R.fraction_field()
 
-        def L_x_i_j(i,j):
-            # From a result in my book, these are the entries of the
-            # basis representation of L_x.
-            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
-                        for k in range(n) )
-
-        L_x = matrix(F, n, n, L_x_i_j)
+        L_x = self.operator_polynomial_matrix()
 
         r = None
         if self.rank.is_in_cache():
 
         r = None
         if self.rank.is_in_cache():
@@ -1630,7 +1650,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         positive integer rank, unless the algebra is trivial in
         which case its rank will be zero::
 
         positive integer rank, unless the algebra is trivial in
         which case its rank will be zero::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: r = J.rank()
             sage: r in ZZ
             sage: J = random_eja()
             sage: r = J.rank()
             sage: r in ZZ
@@ -1641,7 +1660,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         Ensure that computing the rank actually works, since the ranks
         of all simple algebras are known and will be cached by default::
 
         Ensure that computing the rank actually works, since the ranks
         of all simple algebras are known and will be cached by default::
 
-            sage: set_random_seed()    # long time
             sage: J = random_eja()     # long time
             sage: cached = J.rank()    # long time
             sage: J.rank.clear_cache() # long time
             sage: J = random_eja()     # long time
             sage: cached = J.rank()    # long time
             sage: J.rank.clear_cache() # long time
@@ -1656,8 +1674,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         r"""
         Create a subalgebra of this algebra from the given basis.
         """
         r"""
         Create a subalgebra of this algebra from the given basis.
         """
-        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
-        return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+        from mjo.eja.eja_subalgebra import EJASubalgebra
+        return EJASubalgebra(self, basis, **kwargs)
 
 
     def vector_space(self):
 
 
     def vector_space(self):
@@ -1679,7 +1697,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 
 
 
 
 
-class RationalBasisEJA(FiniteDimensionalEJA):
+class RationalBasisEJA(EJA):
     r"""
     Algebras whose supplied basis elements have all rational entries.
 
     r"""
     Algebras whose supplied basis elements have all rational entries.
 
@@ -1734,7 +1752,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             # Note: the same Jordan and inner-products work here,
             # because they are necessarily defined with respect to
             # ambient coordinates and not any particular basis.
             # Note: the same Jordan and inner-products work here,
             # because they are necessarily defined with respect to
             # ambient coordinates and not any particular basis.
-            self._rational_algebra = FiniteDimensionalEJA(
+            self._rational_algebra = EJA(
                                        basis,
                                        jordan_product,
                                        inner_product,
                                        basis,
                                        jordan_product,
                                        inner_product,
@@ -1770,7 +1788,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
 
             sage: J = JordanSpinEJA(3)
             sage: J._charpoly_coefficients()
 
             sage: J = JordanSpinEJA(3)
             sage: J._charpoly_coefficients()
-            (X1^2 - X2^2 - X3^2, -2*X1)
+            (X0^2 - X1^2 - X2^2, -2*X0)
             sage: a0 = J._charpoly_coefficients()[0]
             sage: J.base_ring()
             Algebraic Real Field
             sage: a0 = J._charpoly_coefficients()[0]
             sage: J.base_ring()
             Algebraic Real Field
@@ -1782,14 +1800,13 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             # Bypass the hijinks if they won't benefit us.
             return super()._charpoly_coefficients()
 
             # Bypass the hijinks if they won't benefit us.
             return super()._charpoly_coefficients()
 
-        # Do the computation over the rationals. The answer will be
-        # the same, because all we've done is a change of basis.
-        # Then, change back from QQ to our real base ring
+        # Do the computation over the rationals.
         a = ( a_i.change_ring(self.base_ring())
               for a_i in self.rational_algebra()._charpoly_coefficients() )
 
         a = ( a_i.change_ring(self.base_ring())
               for a_i in self.rational_algebra()._charpoly_coefficients() )
 
-        # Otherwise, convert the coordinate variables back to the
-        # deorthonormalized ones.
+        # Convert our coordinate variables into deorthonormalized ones
+        # and substitute them into the deorthonormalized charpoly
+        # coefficients.
         R = self.coordinate_polynomial_ring()
         from sage.modules.free_module_element import vector
         X = vector(R, R.gens())
         R = self.coordinate_polynomial_ring()
         from sage.modules.free_module_element import vector
         X = vector(R, R.gens())
@@ -1798,7 +1815,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         subs_dict = { X[i]: BX[i] for i in range(len(X)) }
         return tuple( a_i.subs(subs_dict) for a_i in a )
 
         subs_dict = { X[i]: BX[i] for i in range(len(X)) }
         return tuple( a_i.subs(subs_dict) for a_i in a )
 
-class ConcreteEJA(FiniteDimensionalEJA):
+class ConcreteEJA(EJA):
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
@@ -1816,7 +1833,6 @@ class ConcreteEJA(FiniteDimensionalEJA):
     Our basis is normalized with respect to the algebra's inner
     product, unless we specify otherwise::
 
     Our basis is normalized with respect to the algebra's inner
     product, unless we specify otherwise::
 
-        sage: set_random_seed()
         sage: J = ConcreteEJA.random_instance()
         sage: all( b.norm() == 1 for b in J.gens() )
         True
         sage: J = ConcreteEJA.random_instance()
         sage: all( b.norm() == 1 for b in J.gens() )
         True
@@ -1827,7 +1843,6 @@ class ConcreteEJA(FiniteDimensionalEJA):
     natural->EJA basis representation is an isometry and within the
     EJA the operator is self-adjoint by the Jordan axiom::
 
     natural->EJA basis representation is an isometry and within the
     EJA the operator is self-adjoint by the Jordan axiom::
 
-        sage: set_random_seed()
         sage: J = ConcreteEJA.random_instance()
         sage: x = J.random_element()
         sage: x.operator().is_self_adjoint()
         sage: J = ConcreteEJA.random_instance()
         sage: x = J.random_element()
         sage: x.operator().is_self_adjoint()
@@ -1901,11 +1916,11 @@ class ConcreteEJA(FiniteDimensionalEJA):
         return eja_class.random_instance(max_dimension, *args, **kwargs)
 
 
         return eja_class.random_instance(max_dimension, *args, **kwargs)
 
 
-class MatrixEJA(FiniteDimensionalEJA):
+class HermitianMatrixEJA(EJA):
     @staticmethod
     def _denormalized_basis(A):
         """
     @staticmethod
     def _denormalized_basis(A):
         """
-        Returns a basis for the space of complex Hermitian n-by-n matrices.
+        Returns a basis for the given Hermitian matrix space.
 
         Why do we embed these? Basically, because all of numerical linear
         algebra assumes that you're working with vectors consisting of `n`
 
         Why do we embed these? Basically, because all of numerical linear
         algebra assumes that you're working with vectors consisting of `n`
@@ -1918,41 +1933,37 @@ class MatrixEJA(FiniteDimensionalEJA):
             sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
             ....:                          QuaternionMatrixAlgebra,
             ....:                          OctonionMatrixAlgebra)
             sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
             ....:                          QuaternionMatrixAlgebra,
             ....:                          OctonionMatrixAlgebra)
-            sage: from mjo.eja.eja_algebra import MatrixEJA
+            sage: from mjo.eja.eja_algebra import HermitianMatrixEJA
 
         TESTS::
 
 
         TESTS::
 
-            sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
             sage: A = MatrixSpace(QQ, n)
             sage: n = ZZ.random_element(1,5)
             sage: A = MatrixSpace(QQ, n)
-            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: B = HermitianMatrixEJA._denormalized_basis(A)
             sage: all( M.is_hermitian() for M in  B)
             True
 
         ::
 
             sage: all( M.is_hermitian() for M in  B)
             True
 
         ::
 
-            sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
             sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
             sage: n = ZZ.random_element(1,5)
             sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
-            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: B = HermitianMatrixEJA._denormalized_basis(A)
             sage: all( M.is_hermitian() for M in  B)
             True
 
         ::
 
             sage: all( M.is_hermitian() for M in  B)
             True
 
         ::
 
-            sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
             sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
             sage: n = ZZ.random_element(1,5)
             sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
-            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: B = HermitianMatrixEJA._denormalized_basis(A)
             sage: all( M.is_hermitian() for M in B )
             True
 
         ::
 
             sage: all( M.is_hermitian() for M in B )
             True
 
         ::
 
-            sage: set_random_seed()
             sage: n = ZZ.random_element(1,5)
             sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
             sage: n = ZZ.random_element(1,5)
             sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
-            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: B = HermitianMatrixEJA._denormalized_basis(A)
             sage: all( M.is_hermitian() for M in B )
             True
 
             sage: all( M.is_hermitian() for M in B )
             True
 
@@ -2057,7 +2068,7 @@ class MatrixEJA(FiniteDimensionalEJA):
         self.rank.set_cache(matrix_space.nrows())
         self.one.set_cache( self(matrix_space.one()) )
 
         self.rank.set_cache(matrix_space.nrows())
         self.one.set_cache( self(matrix_space.one()) )
 
-class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
+class RealSymmetricEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -2090,7 +2101,6 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     The dimension of this algebra is `(n^2 + n) / 2`::
 
 
     The dimension of this algebra is `(n^2 + n) / 2`::
 
-        sage: set_random_seed()
         sage: d = RealSymmetricEJA._max_random_instance_dimension()
         sage: n = RealSymmetricEJA._max_random_instance_size(d)
         sage: J = RealSymmetricEJA(n)
         sage: d = RealSymmetricEJA._max_random_instance_dimension()
         sage: n = RealSymmetricEJA._max_random_instance_size(d)
         sage: J = RealSymmetricEJA(n)
@@ -2099,7 +2109,6 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     The Jordan multiplication is what we think it is::
 
 
     The Jordan multiplication is what we think it is::
 
-        sage: set_random_seed()
         sage: J = RealSymmetricEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
         sage: J = RealSymmetricEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
@@ -2151,7 +2160,7 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
 
 
 
 
 
-class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
+class ComplexHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -2178,20 +2187,10 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         ...
         TypeError: Illegal initializer for algebraic number
 
         ...
         TypeError: Illegal initializer for algebraic number
 
-    This causes the following error when we try to scale a matrix of
-    complex numbers by an inexact real number::
-
-        sage: ComplexHermitianEJA(2,field=RR)
-        Traceback (most recent call last):
-        ...
-        TypeError: Unable to coerce entries (=(1.00000000000000,
-        -0.000000000000000)) to coefficients in Algebraic Real Field
-
     TESTS:
 
     The dimension of this algebra is `n^2`::
 
     TESTS:
 
     The dimension of this algebra is `n^2`::
 
-        sage: set_random_seed()
         sage: d = ComplexHermitianEJA._max_random_instance_dimension()
         sage: n = ComplexHermitianEJA._max_random_instance_size(d)
         sage: J = ComplexHermitianEJA(n)
         sage: d = ComplexHermitianEJA._max_random_instance_dimension()
         sage: n = ComplexHermitianEJA._max_random_instance_size(d)
         sage: J = ComplexHermitianEJA(n)
@@ -2200,7 +2199,6 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     The Jordan multiplication is what we think it is::
 
 
     The Jordan multiplication is what we think it is::
 
-        sage: set_random_seed()
         sage: J = ComplexHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
         sage: J = ComplexHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
@@ -2252,7 +2250,7 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         return cls(n, **kwargs)
 
 
         return cls(n, **kwargs)
 
 
-class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
+class QuaternionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
@@ -2277,7 +2275,6 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     The dimension of this algebra is `2*n^2 - n`::
 
 
     The dimension of this algebra is `2*n^2 - n`::
 
-        sage: set_random_seed()
         sage: d = QuaternionHermitianEJA._max_random_instance_dimension()
         sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
         sage: J = QuaternionHermitianEJA(n)
         sage: d = QuaternionHermitianEJA._max_random_instance_dimension()
         sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
         sage: J = QuaternionHermitianEJA(n)
@@ -2286,7 +2283,6 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     The Jordan multiplication is what we think it is::
 
 
     The Jordan multiplication is what we think it is::
 
-        sage: set_random_seed()
         sage: J = QuaternionHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
         sage: J = QuaternionHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
         sage: actual = (x*y).to_matrix()
@@ -2342,11 +2338,11 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
         n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
-class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
+class OctonionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
     r"""
     SETUP::
 
     r"""
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+        sage: from mjo.eja.eja_algebra import (EJA,
         ....:                                  OctonionHermitianEJA)
         sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
 
         ....:                                  OctonionHermitianEJA)
         sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
 
@@ -2368,7 +2364,7 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
         sage: jp = OctonionHermitianEJA.jordan_product
         sage: ip = OctonionHermitianEJA.trace_inner_product
         sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
         sage: jp = OctonionHermitianEJA.jordan_product
         sage: ip = OctonionHermitianEJA.trace_inner_product
-        sage: J = FiniteDimensionalEJA(basis,
+        sage: J = EJA(basis,
         ....:                          jp,
         ....:                          ip,
         ....:                          field=QQ,
         ....:                          jp,
         ....:                          ip,
         ....:                          field=QQ,
@@ -2432,7 +2428,7 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
     @staticmethod
     def _max_random_instance_size(max_dimension):
         r"""
     @staticmethod
     def _max_random_instance_size(max_dimension):
         r"""
-        The maximum rank of a random QuaternionHermitianEJA.
+        The maximum rank of a random OctonionHermitianEJA.
         """
         # There's certainly a formula for this, but with only four
         # cases to worry about, I'm not that motivated to derive it.
         """
         # There's certainly a formula for this, but with only four
         # cases to worry about, I'm not that motivated to derive it.
@@ -2659,7 +2655,6 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
     matrix.  We opt not to orthonormalize the basis, because if we
     did, we would have to normalize the `s_{i}` in a similar manner::
 
     matrix.  We opt not to orthonormalize the basis, because if we
     did, we would have to normalize the `s_{i}` in a similar manner::
 
-        sage: set_random_seed()
         sage: n = ZZ.random_element(5)
         sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
         sage: B11 = matrix.identity(QQ,1)
         sage: n = ZZ.random_element(5)
         sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
         sage: B11 = matrix.identity(QQ,1)
@@ -2821,7 +2816,6 @@ class JordanSpinEJA(BilinearFormEJA):
 
         Ensure that we have the usual inner product on `R^n`::
 
 
         Ensure that we have the usual inner product on `R^n`::
 
-            sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
             sage: J = JordanSpinEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: actual = x.inner_product(y)
@@ -2920,7 +2914,7 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA):
         return cls(**kwargs)
 
 
         return cls(**kwargs)
 
 
-class CartesianProductEJA(FiniteDimensionalEJA):
+class CartesianProductEJA(EJA):
     r"""
     The external (orthogonal) direct sum of two or more Euclidean
     Jordan algebras. Every Euclidean Jordan algebra decomposes into an
     r"""
     The external (orthogonal) direct sum of two or more Euclidean
     Jordan algebras. Every Euclidean Jordan algebra decomposes into an
@@ -2942,7 +2936,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
     The Jordan product is inherited from our factors and implemented by
     our CombinatorialFreeModule Cartesian product superclass::
 
     The Jordan product is inherited from our factors and implemented by
     our CombinatorialFreeModule Cartesian product superclass::
 
-        sage: set_random_seed()
         sage: J1 = HadamardEJA(2)
         sage: J2 = RealSymmetricEJA(2)
         sage: J = cartesian_product([J1,J2])
         sage: J1 = HadamardEJA(2)
         sage: J2 = RealSymmetricEJA(2)
         sage: J = cartesian_product([J1,J2])
@@ -3079,7 +3072,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
     The cached unit element is the same one that would be computed::
 
 
     The cached unit element is the same one that would be computed::
 
-        sage: set_random_seed()              # long time
         sage: J1 = random_eja()              # long time
         sage: J2 = random_eja()              # long time
         sage: J = cartesian_product([J1,J2]) # long time
         sage: J1 = random_eja()              # long time
         sage: J2 = random_eja()              # long time
         sage: J = cartesian_product([J1,J2]) # long time
@@ -3089,6 +3081,7 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         sage: actual == expected             # long time
         True
     """
         sage: actual == expected             # long time
         True
     """
+    Element = CartesianProductEJAElement
     def __init__(self, factors, **kwargs):
         m = len(factors)
         if m == 0:
     def __init__(self, factors, **kwargs):
         m = len(factors)
         if m == 0:
@@ -3297,7 +3290,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
         The answer never changes::
 
 
         The answer never changes::
 
-            sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
@@ -3313,7 +3305,7 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
                                    codomain=Ji)
 
         Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
                                    codomain=Ji)
 
-        return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
+        return EJAOperator(self,Ji,Pi.matrix())
 
     @cached_method
     def cartesian_embedding(self, i):
 
     @cached_method
     def cartesian_embedding(self, i):
@@ -3387,7 +3379,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
         The answer never changes::
 
 
         The answer never changes::
 
-            sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
@@ -3400,7 +3391,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         produce the identity map, and mismatching them should produce
         the zero map::
 
         produce the identity map, and mismatching them should produce
         the zero map::
 
-            sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
@@ -3423,11 +3413,39 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         Ji = self.cartesian_factor(i)
         Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
                                  codomain=self)
         Ji = self.cartesian_factor(i)
         Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
                                  codomain=self)
-        return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
+        return EJAOperator(Ji,self,Ei.matrix())
+
+
+    def subalgebra(self, basis, **kwargs):
+        r"""
+        Create a subalgebra of this algebra from the given basis.
+
+        Only overridden to allow us to use a special Cartesian product
+        subalgebra class.
 
 
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  QuaternionHermitianEJA)
+
+        EXAMPLES:
 
 
+        Subalgebras of Cartesian product EJAs have a different class
+        than those of non-Cartesian-product EJAs::
 
 
-FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
+            sage: J1 = HadamardEJA(2,field=QQ,orthonormalize=False)
+            sage: J2 = QuaternionHermitianEJA(0,field=QQ,orthonormalize=False)
+            sage: J = cartesian_product([J1,J2])
+            sage: K1 = J1.subalgebra((J1.one(),), orthonormalize=False)
+            sage: K = J.subalgebra((J.one(),), orthonormalize=False)
+            sage: K1.__class__ is K.__class__
+            False
+
+        """
+        from mjo.eja.eja_subalgebra import CartesianProductEJASubalgebra
+        return CartesianProductEJASubalgebra(self, basis, **kwargs)
+
+EJA.CartesianProduct = CartesianProductEJA
 
 class RationalBasisCartesianProductEJA(CartesianProductEJA,
                                        RationalBasisEJA):
 
 class RationalBasisCartesianProductEJA(CartesianProductEJA,
                                        RationalBasisEJA):
@@ -3437,7 +3455,7 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     SETUP::
 
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+        sage: from mjo.eja.eja_algebra import (EJA,
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
@@ -3467,7 +3485,7 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
         sage: jp = lambda X,Y: X*Y
         sage: ip = lambda X,Y: X[0,0]*Y[0,0]
         sage: b1 = matrix(QQ, [[1]])
         sage: jp = lambda X,Y: X*Y
         sage: ip = lambda X,Y: X[0,0]*Y[0,0]
         sage: b1 = matrix(QQ, [[1]])
-        sage: J2 = FiniteDimensionalEJA((b1,), jp, ip)
+        sage: J2 = EJA((b1,), jp, ip)
         sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA
         Euclidean Jordan algebra of dimension 1 over Algebraic Real
         Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic
         sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA
         Euclidean Jordan algebra of dimension 1 over Algebraic Real
         Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic
@@ -3505,7 +3523,6 @@ def random_eja(max_dimension=None, *args, **kwargs):
 
     TESTS::
 
 
     TESTS::
 
-        sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
         sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
         sage: J.dimension() <= n
         sage: n = ZZ.random_element(1,5)
         sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
         sage: J.dimension() <= n
@@ -3531,3 +3548,234 @@ def random_eja(max_dimension=None, *args, **kwargs):
         # if the sub-call also Decides on a cartesian product.
         J2 = random_eja(new_max_dimension, *args, **kwargs)
         return cartesian_product([J1,J2])
         # if the sub-call also Decides on a cartesian product.
         J2 = random_eja(new_max_dimension, *args, **kwargs)
         return cartesian_product([J1,J2])
+
+
+class ComplexSkewSymmetricEJA(RationalBasisEJA, ConcreteEJA):
+    r"""
+    The skew-symmetric EJA of order `2n` described in Faraut and
+    Koranyi's Exercise III.1.b. It has dimension `2n^2 - n`.
+
+    It is (not obviously) isomorphic to the QuaternionHermitianEJA of
+    order `n`, as can be inferred by comparing rank/dimension or
+    explicitly from their "characteristic polynomial of" functions,
+    which just so happen to align nicely.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (ComplexSkewSymmetricEJA,
+        ....:                                  QuaternionHermitianEJA)
+        sage: from mjo.eja.eja_operator import EJAOperator
+
+    EXAMPLES:
+
+    This EJA is isomorphic to the quaternions::
+
+        sage: J = ComplexSkewSymmetricEJA(2, field=QQ, orthonormalize=False)
+        sage: K = QuaternionHermitianEJA(2, field=QQ, orthonormalize=False)
+        sage: jordan_isom_matrix = matrix.diagonal(QQ,[-1,1,1,1,1,-1])
+        sage: phi = EJAOperator(J,K,jordan_isom_matrix)
+        sage: all( phi(x*y) == phi(x)*phi(y)
+        ....:      for x in J.gens()
+        ....:      for y in J.gens() )
+        True
+        sage: x,y = J.random_elements(2)
+        sage: phi(x*y) == phi(x)*phi(y)
+        True
+
+    TESTS:
+
+    Random elements should satisfy the same conditions that the basis
+    elements do::
+
+        sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ,
+        ....:                                             orthonormalize=False)
+        sage: x,y = K.random_elements(2)
+        sage: z = x*y
+        sage: x = x.to_matrix()
+        sage: y = y.to_matrix()
+        sage: z = z.to_matrix()
+        sage: all( e.is_skew_symmetric() for e in (x,y,z) )
+        True
+        sage: J = -K.one().to_matrix()
+        sage: all( e*J == J*e.conjugate() for e in (x,y,z) )
+        True
+
+    The power law in Faraut & Koranyi's II.7.a is satisfied.
+    We're in a subalgebra of theirs, but powers are still
+    defined the same::
+
+        sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ,
+        ....:                                             orthonormalize=False)
+        sage: x = K.random_element()
+        sage: k = ZZ.random_element(5)
+        sage: actual = x^k
+        sage: J = -K.one().to_matrix()
+        sage: expected = K(-J*(J*x.to_matrix())^k)
+        sage: actual == expected
+        True
+
+    """
+    @staticmethod
+    def _max_random_instance_size(max_dimension):
+        # Obtained by solving d = 2n^2 - n, which comes from noticing
+        # that, in 2x2 block form, any element of this algebra has a
+        # free skew-symmetric top-left block, a Hermitian top-right
+        # block, and two bottom blocks that are determined by the top.
+        # The ZZ-int-ZZ thing is just "floor."
+        return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))
+
+    @classmethod
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
+        return cls(n, **kwargs)
+
+    @staticmethod
+    def _denormalized_basis(A):
+        """
+        SETUP::
+
+            sage: from mjo.hurwitz import ComplexMatrixAlgebra
+            sage: from mjo.eja.eja_algebra import ComplexSkewSymmetricEJA
+
+        TESTS:
+
+        The basis elements are all skew-Hermitian::
+
+            sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension()
+            sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max)
+            sage: n = ZZ.random_element(n_max + 1)
+            sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
+            sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A)
+            sage: all( M.is_skew_symmetric() for M in  B)
+            True
+
+        The basis elements ``b`` all satisfy ``b*J == J*b.conjugate()``,
+        as in the definition of the algebra::
+
+            sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension()
+            sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max)
+            sage: n = ZZ.random_element(n_max + 1)
+            sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
+            sage: I_n = matrix.identity(ZZ, n)
+            sage: J = matrix.block(ZZ, 2, 2, (0, I_n, -I_n, 0), subdivide=False)
+            sage: J = A.from_list(J.rows())
+            sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A)
+            sage: all( b*J == J*b.conjugate()  for b in B )
+            True
+
+        """
+        es = A.entry_algebra_gens()
+        gen = lambda A,m: A.monomial(m)
+
+        basis = []
+
+        # The size of the blocks. We're going to treat these thing as
+        # 2x2 block matrices,
+        #
+        #   [  x1        x2      ]
+        #   [ -x2-conj   x1-conj ]
+        #
+        # where x1 is skew-symmetric and x2 is Hermitian.
+        #
+        m = A.nrows()/2
+
+        # We only loop through the top half of the matrix, because the
+        # bottom can be constructed from the top.
+        for i in range(m):
+            # First do the top-left block, which is skew-symmetric.
+            # We can compute the bottom-right block in the process.
+            for j in range(i+1):
+                if i != j:
+                    # Skew-symmetry implies zeros for (i == j).
+                    for e in es:
+                        # Top-left block's entry.
+                        E_ij  = gen(A, (i,j,e))
+                        E_ij -= gen(A, (j,i,e))
+
+                        # Bottom-right block's entry.
+                        F_ij  = gen(A, (i+m,j+m,e)).conjugate()
+                        F_ij -= gen(A, (j+m,i+m,e)).conjugate()
+
+                        basis.append(E_ij + F_ij)
+
+            # Now do the top-right block, which is Hermitian, and compute
+            # the bottom-left block along the way.
+            for j in range(m,i+m+1):
+                if (i+m) == j:
+                    # Hermitian matrices have real diagonal entries.
+                    # Top-right block's entry.
+                    E_ii = gen(A, (i,j,es[0]))
+
+                    # Bottom-left block's entry. Don't conjugate
+                    # 'cause it's real.
+                    E_ii -= gen(A, (i+m,j-m,es[0]))
+                    basis.append(E_ii)
+                else:
+                    for e in es:
+                        # Top-right block's entry. BEWARE! We're not
+                        # reflecting across the main diagonal as in
+                        # (i,j)~(j,i). We're only reflecting across
+                        # the diagonal for the top-right block.
+                        E_ij  = gen(A, (i,j,e))
+
+                        # Shift it back to non-offset coords, transpose,
+                        # conjugate, and put it back:
+                        #
+                        # (i,j) -> (i,j-m) -> (j-m, i) -> (j-m, i+m)
+                        E_ij += gen(A, (j-m,i+m,e)).conjugate()
+
+                        # Bottom-left's block's below-diagonal entry.
+                        # Just shift the top-right coords down m and
+                        # left m.
+                        F_ij  = -gen(A, (i+m,j-m,e)).conjugate()
+                        F_ij += -gen(A, (j,i,e)) # double-conjugate cancels
+
+                        basis.append(E_ij + F_ij)
+
+        return tuple( basis )
+
+    @staticmethod
+    @cached_method
+    def _J_matrix(matrix_space):
+        n = matrix_space.nrows() // 2
+        F = matrix_space.base_ring()
+        I_n = matrix.identity(F, n)
+        J = matrix.block(F, 2, 2, (0, I_n, -I_n, 0), subdivide=False)
+        return matrix_space.from_list(J.rows())
+
+    def J_matrix(self):
+        return ComplexSkewSymmetricEJA._J_matrix(self.matrix_space())
+
+    def __init__(self, n, field=AA, **kwargs):
+        # New code; always check the axioms.
+        #if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+        from mjo.hurwitz import ComplexMatrixAlgebra
+        A = ComplexMatrixAlgebra(2*n, scalars=field)
+        J = ComplexSkewSymmetricEJA._J_matrix(A)
+
+        def jordan_product(X,Y):
+            return (X*J*Y + Y*J*X)/2
+
+        def inner_product(X,Y):
+            return (X.conjugate_transpose()*Y).trace().real()
+
+        super().__init__(self._denormalized_basis(A),
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         matrix_space=A,
+                         **kwargs)
+
+        # This algebra is conjectured (by me) to be isomorphic to
+        # the quaternion Hermitian EJA of size n, and the rank
+        # would follow from that.
+        #self.rank.set_cache(n)
+        self.one.set_cache( self(-J) )