]> 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 79efd3876b6fadd0d2a7494a45856bf46bf3a621..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,9 +167,10 @@ 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_utils import _all2list, _mat2vec
+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):
     r"""
 
 def EuclideanJordanAlgebras(field):
     r"""
@@ -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()
@@ -1281,7 +1266,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         #
         # Of course, matrices aren't vectors in sage, so we have to
         # appeal to the "long vectors" isometry.
         #
         # Of course, matrices aren't vectors in sage, so we have to
         # appeal to the "long vectors" isometry.
-        oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
+
+        V = VectorSpace(self.base_ring(), self.dimension()**2)
+        oper_vecs = [ V(g.operator().matrix().list()) for g in self.gens() ]
 
         # Now we use basic linear algebra to find the coefficients,
         # of the matrices-as-vectors-linear-combination, which should
 
         # Now we use basic linear algebra to find the coefficients,
         # of the matrices-as-vectors-linear-combination, which should
@@ -1291,7 +1278,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # We used the isometry on the left-hand side already, but we
         # still need to do it for the right-hand side. Recall that we
         # wanted something that summed to the identity matrix.
         # We used the isometry on the left-hand side already, but we
         # still need to do it for the right-hand side. Recall that we
         # wanted something that summed to the identity matrix.
-        b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
+        b = V( matrix.identity(self.base_ring(), self.dimension()).list() )
 
         # Now if there's an identity element in the algebra, this
         # should work. We solve on the left to avoid having to
 
         # Now if there's an identity element in the algebra, this
         # should work. We solve on the left to avoid having to
@@ -1376,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
@@ -1389,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():
@@ -1463,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):
@@ -1515,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"""
@@ -1530,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
@@ -1538,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():
@@ -1628,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
@@ -1639,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
@@ -1654,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):
@@ -1677,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.
 
@@ -1732,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,
@@ -1743,6 +1763,15 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        check_field=False,
                                        check_axioms=False)
 
                                        check_field=False,
                                        check_axioms=False)
 
+    def rational_algebra(self):
+        # Using None as a flag here (rather than just assigning "self"
+        # to self._rational_algebra by default) feels a little bit
+        # more sane to me in a garbage-collected environment.
+        if self._rational_algebra is None:
+            return self
+        else:
+            return self._rational_algebra
+
     @cached_method
     def _charpoly_coefficients(self):
         r"""
     @cached_method
     def _charpoly_coefficients(self):
         r"""
@@ -1759,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
@@ -1767,21 +1796,17 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             Algebraic Real Field
 
         """
             Algebraic Real Field
 
         """
-        if self._rational_algebra is None:
-            # There's no need to construct *another* algebra over the
-            # rationals if this one is already over the
-            # rationals. Likewise, if we never orthonormalized our
-            # basis, we might as well just use the given one.
+        if self.rational_algebra() is self:
+            # Bypass the hijinks if they won't benefit us.
             return super()._charpoly_coefficients()
 
             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())
         a = ( a_i.change_ring(self.base_ring())
-              for a_i in self._rational_algebra._charpoly_coefficients() )
+              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())
@@ -1790,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.
 
@@ -1808,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
@@ -1819,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()
@@ -1893,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`
@@ -1910,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
 
@@ -2039,7 +2058,6 @@ class MatrixEJA(FiniteDimensionalEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-
         super().__init__(self._denormalized_basis(matrix_space),
                          self.jordan_product,
                          self.trace_inner_product,
         super().__init__(self._denormalized_basis(matrix_space),
                          self.jordan_product,
                          self.trace_inner_product,
@@ -2050,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
@@ -2083,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)
@@ -2092,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()
@@ -2134,24 +2150,17 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         return cls(n, **kwargs)
 
     def __init__(self, n, field=AA, **kwargs):
         return cls(n, **kwargs)
 
     def __init__(self, n, field=AA, **kwargs):
-        # We know this is a valid EJA, but will double-check
-        # if the user passes check_axioms=True.
-        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
-
         A = MatrixSpace(field, n)
         super().__init__(A, **kwargs)
 
         from mjo.eja.eja_cache import real_symmetric_eja_coeffs
         a = real_symmetric_eja_coeffs(self)
         if a is not None:
         A = MatrixSpace(field, n)
         super().__init__(A, **kwargs)
 
         from mjo.eja.eja_cache import real_symmetric_eja_coeffs
         a = real_symmetric_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 
 
 
 
-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()
@@ -2224,10 +2222,6 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     """
     def __init__(self, n, field=AA, **kwargs):
 
     """
     def __init__(self, n, field=AA, **kwargs):
-        # We know this is a valid EJA, but will double-check
-        # if the user passes check_axioms=True.
-        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
-
         from mjo.hurwitz import ComplexMatrixAlgebra
         A = ComplexMatrixAlgebra(n, scalars=field)
         super().__init__(A, **kwargs)
         from mjo.hurwitz import ComplexMatrixAlgebra
         A = ComplexMatrixAlgebra(n, scalars=field)
         super().__init__(A, **kwargs)
@@ -2235,10 +2229,7 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
         a = complex_hermitian_eja_coeffs(self)
         if a is not None:
         from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
         a = complex_hermitian_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
     @staticmethod
     def _max_random_instance_size(max_dimension):
 
     @staticmethod
     def _max_random_instance_size(max_dimension):
@@ -2259,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
@@ -2284,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)
@@ -2293,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()
@@ -2317,10 +2306,6 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     """
     def __init__(self, n, field=AA, **kwargs):
 
     """
     def __init__(self, n, field=AA, **kwargs):
-        # We know this is a valid EJA, but will double-check
-        # if the user passes check_axioms=True.
-        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
-
         from mjo.hurwitz import QuaternionMatrixAlgebra
         A = QuaternionMatrixAlgebra(n, scalars=field)
         super().__init__(A, **kwargs)
         from mjo.hurwitz import QuaternionMatrixAlgebra
         A = QuaternionMatrixAlgebra(n, scalars=field)
         super().__init__(A, **kwargs)
@@ -2328,10 +2313,7 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs
         a = quaternion_hermitian_eja_coeffs(self)
         if a is not None:
         from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs
         a = quaternion_hermitian_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 
 
 
 
@@ -2356,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
 
@@ -2382,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,
@@ -2446,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.
@@ -2487,10 +2469,7 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs
         a = octonion_hermitian_eja_coeffs(self)
         if a is not None:
         from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs
         a = octonion_hermitian_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 class AlbertEJA(OctonionHermitianEJA):
 
 
 class AlbertEJA(OctonionHermitianEJA):
@@ -2676,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)
@@ -2838,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)
@@ -2937,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
@@ -2949,6 +2926,7 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
         sage: from mjo.eja.eja_algebra import (random_eja,
         ....:                                  CartesianProductEJA,
 
         sage: from mjo.eja.eja_algebra import (random_eja,
         ....:                                  CartesianProductEJA,
+        ....:                                  ComplexHermitianEJA,
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
@@ -2958,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])
@@ -3060,6 +3037,28 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         | b2 || 0  | 0  | b2 |
         +----++----+----+----+
 
         | b2 || 0  | 0  | b2 |
         +----++----+----+----+
 
+    The "matrix space" of a Cartesian product always consists of
+    ordered pairs (or triples, or...) whose components are the
+    matrix spaces of its factors::
+
+            sage: J1 = HadamardEJA(2)
+            sage: J2 = ComplexHermitianEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.matrix_space()
+            The Cartesian product of (Full MatrixSpace of 2 by 1 dense
+            matrices over Algebraic Real Field, Module of 2 by 2 matrices
+            with entries in Algebraic Field over the scalar ring Algebraic
+            Real Field)
+            sage: J.one().to_matrix()[0]
+            [1]
+            [1]
+            sage: J.one().to_matrix()[1]
+            +---+---+
+            | 1 | 0 |
+            +---+---+
+            | 0 | 1 |
+            +---+---+
+
     TESTS:
 
     All factors must share the same base field::
     TESTS:
 
     All factors must share the same base field::
@@ -3073,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
@@ -3082,8 +3080,8 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         sage: expected = J.one()             # long time
         sage: actual == expected             # long time
         True
         sage: expected = J.one()             # long time
         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:
@@ -3101,8 +3099,13 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         if associative: category = category.Associative()
         category = category.join([category, category.CartesianProducts()])
 
         if associative: category = category.Associative()
         category = category.join([category, category.CartesianProducts()])
 
-        # Compute my matrix space. This category isn't perfect, but
-        # is good enough for what we need to do.
+        # Compute my matrix space.  We don't simply use the
+        # ``cartesian_product()`` functor here because it acts
+        # differently on SageMath MatrixSpaces and our custom
+        # MatrixAlgebras, which are CombinatorialFreeModules. We
+        # always want the result to be represented (and indexed) as an
+        # ordered tuple. This category isn't perfect, but is good
+        # enough for what we need to do.
         MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
         MS_cat = MS_cat.Unital().CartesianProducts()
         MS_factors = tuple( J.matrix_space() for J in factors )
         MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
         MS_cat = MS_cat.Unital().CartesianProducts()
         MS_factors = tuple( J.matrix_space() for J in factors )
@@ -3153,6 +3156,8 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         self._inner_product_matrix = matrix.block_diagonal(
             [J._inner_product_matrix for J in factors]
         )
         self._inner_product_matrix = matrix.block_diagonal(
             [J._inner_product_matrix for J in factors]
         )
+        self._inner_product_matrix._cache = {'hermitian': True}
+        self._inner_product_matrix.set_immutable()
 
         # Building the multiplication table is a bit more tricky
         # because we have to embed the entries of the factors'
 
         # Building the multiplication table is a bit more tricky
         # because we have to embed the entries of the factors'
@@ -3180,6 +3185,34 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         ones = tuple(J.one().to_matrix() for J in factors)
         self.one.set_cache(self(ones))
 
         ones = tuple(J.one().to_matrix() for J in factors)
         self.one.set_cache(self(ones))
 
+    def _sets_keys(self):
+        r"""
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  RealSymmetricEJA)
+
+        TESTS:
+
+        The superclass uses ``_sets_keys()`` to implement its
+        ``cartesian_factors()`` method::
+
+            sage: J1 = RealSymmetricEJA(2,
+            ....:                       field=QQ,
+            ....:                       orthonormalize=False,
+            ....:                       prefix="a")
+            sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J = cartesian_product([J1,J2])
+            sage: x = sum(i*J.gens()[i] for i in range(len(J.gens())))
+            sage: x.cartesian_factors()
+            (a1 + 2*a2, 3*b0 + 4*b1 + 5*b2 + 6*b3)
+
+        """
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct,
+        # but returning a tuple instead of a list.
+        return tuple(range(len(self.cartesian_factors())))
+
     def cartesian_factors(self):
         # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
         return self._sets
     def cartesian_factors(self):
         # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
         return self._sets
@@ -3196,65 +3229,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         return cartesian_product.symbol.join("%s" % factor
                                              for factor in self._sets)
 
         return cartesian_product.symbol.join("%s" % factor
                                              for factor in self._sets)
 
-    def matrix_space(self):
-        r"""
-        Return the space that our matrix basis lives in as a Cartesian
-        product.
-
-        We don't simply use the ``cartesian_product()`` functor here
-        because it acts differently on SageMath MatrixSpaces and our
-        custom MatrixAlgebras, which are CombinatorialFreeModules. We
-        always want the result to be represented (and indexed) as
-        an ordered tuple.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
-            ....:                                  HadamardEJA,
-            ....:                                  OctonionHermitianEJA,
-            ....:                                  RealSymmetricEJA)
-
-        EXAMPLES::
-
-            sage: J1 = HadamardEJA(1)
-            sage: J2 = RealSymmetricEJA(2)
-            sage: J = cartesian_product([J1,J2])
-            sage: J.matrix_space()
-            The Cartesian product of (Full MatrixSpace of 1 by 1 dense
-            matrices over Algebraic Real Field, Full MatrixSpace of 2
-            by 2 dense matrices over Algebraic Real Field)
-
-        ::
-
-            sage: J1 = ComplexHermitianEJA(1)
-            sage: J2 = ComplexHermitianEJA(1)
-            sage: J = cartesian_product([J1,J2])
-            sage: J.one().to_matrix()[0]
-            +---+
-            | 1 |
-            +---+
-            sage: J.one().to_matrix()[1]
-            +---+
-            | 1 |
-            +---+
-
-        ::
-
-            sage: J1 = OctonionHermitianEJA(1)
-            sage: J2 = OctonionHermitianEJA(1)
-            sage: J = cartesian_product([J1,J2])
-            sage: J.one().to_matrix()[0]
-            +----+
-            | e0 |
-            +----+
-            sage: J.one().to_matrix()[1]
-            +----+
-            | e0 |
-            +----+
-
-        """
-        return super().matrix_space()
-
 
     @cached_method
     def cartesian_projection(self, i):
 
     @cached_method
     def cartesian_projection(self, i):
@@ -3316,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])
@@ -3332,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):
@@ -3406,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])
@@ -3419,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])
@@ -3442,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::
 
 
+            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)
 
 
-FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
+EJA.CartesianProduct = CartesianProductEJA
 
 class RationalBasisCartesianProductEJA(CartesianProductEJA,
                                        RationalBasisEJA):
 
 class RationalBasisCartesianProductEJA(CartesianProductEJA,
                                        RationalBasisEJA):
@@ -3456,9 +3455,9 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     SETUP::
 
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        sage: from mjo.eja.eja_algebra import (EJA,
+        ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  JordanSpinEJA,
-        ....:                                  OctonionHermitianEJA,
         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:
         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:
@@ -3479,28 +3478,38 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     The ``cartesian_product()`` function only uses the first factor to
     decide where the result will live; thus we have to be careful to
 
     The ``cartesian_product()`` function only uses the first factor to
     decide where the result will live; thus we have to be careful to
-    check that all factors do indeed have a `_rational_algebra` member
-    before we try to access it::
-
-        sage: J1 = OctonionHermitianEJA(1) # no rational basis
-        sage: J2 = HadamardEJA(2)
-        sage: cartesian_product([J1,J2])
-        Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
-        (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
-        sage: cartesian_product([J2,J1])
-        Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
-        (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+    check that all factors do indeed have a ``rational_algebra()`` method
+    before we construct an algebra that claims to have a rational basis::
+
+        sage: J1 = HadamardEJA(2)
+        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 = 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
+        Real Field
+        sage: cartesian_product([J1,J2]) # factor one is RationalBasisEJA
+        Traceback (most recent call last):
+        ...
+        ValueError: factor not a RationalBasisEJA
 
     """
     def __init__(self, algebras, **kwargs):
 
     """
     def __init__(self, algebras, **kwargs):
+        if not all( hasattr(r, "rational_algebra") for r in algebras ):
+            raise ValueError("factor not a RationalBasisEJA")
+
         CartesianProductEJA.__init__(self, algebras, **kwargs)
 
         CartesianProductEJA.__init__(self, algebras, **kwargs)
 
-        self._rational_algebra = None
-        if self.vector_space().base_field() is not QQ:
-            if all( hasattr(r, "_rational_algebra") for r in algebras ):
-                self._rational_algebra = cartesian_product([
-                    r._rational_algebra for r in algebras
-                ])
+    @cached_method
+    def rational_algebra(self):
+        if self.base_ring() is QQ:
+            return self
+
+        return cartesian_product([
+            r.rational_algebra() for r in self.cartesian_factors()
+        ])
 
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
 
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
@@ -3514,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
@@ -3540,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) )