]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: rename operator_inner_product -> operator_trace inner_product.
[sage.d.git] / mjo / eja / eja_algebra.py
index 2bad32c2f500193e4126b7c5e209c0acb3116ede..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
@@ -34,12 +34,13 @@ for these simple algebras:
   * :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:`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
@@ -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::
 
-    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: 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: J2 = FiniteDimensionalEJA((b1,), jp, ip2)
+    sage: J2 = EJA((b1,), jp, ip2)
     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::
 
-    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
 
@@ -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: 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()
@@ -167,8 +168,8 @@ from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
 from mjo.eja.eja_element import (CartesianProductEJAElement,
-                                 FiniteDimensionalEJAElement)
-from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
+                                 EJAElement)
+from mjo.eja.eja_operator import EJAOperator
 from mjo.eja.eja_utils import _all2list
 
 def EuclideanJordanAlgebras(field):
@@ -182,7 +183,7 @@ def EuclideanJordanAlgebras(field):
     category = category.WithBasis().Unital().Commutative()
     return category
 
-class FiniteDimensionalEJA(CombinatorialFreeModule):
+class EJA(CombinatorialFreeModule):
     r"""
     A finite-dimensional Euclidean Jordan algebra.
 
@@ -230,7 +231,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
     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)
@@ -238,7 +238,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         sage: J.subalgebra(basis, orthonormalize=False).is_associative()
         True
     """
-    Element = FiniteDimensionalEJAElement
+    Element = EJAElement
 
     @staticmethod
     def _check_input_field(field):
@@ -367,7 +367,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         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
@@ -432,7 +432,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         TESTS::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: J(1)
             Traceback (most recent call last):
@@ -457,7 +456,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         TESTS::
 
-            sage: set_random_seed()
             sage: J = random_eja()
             sage: n = J.dimension()
             sage: bi = J.zero()
@@ -499,7 +497,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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)
@@ -510,7 +507,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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)
@@ -523,7 +519,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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()
@@ -636,7 +631,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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
@@ -758,7 +752,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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
@@ -871,7 +864,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             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
@@ -920,13 +913,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             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()
-            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):
@@ -948,7 +941,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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)
@@ -959,7 +951,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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)
@@ -972,7 +963,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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()
@@ -1200,19 +1190,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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: 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: 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
@@ -1226,14 +1214,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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: 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
@@ -1241,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())
@@ -1257,7 +1243,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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()
@@ -1266,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()
@@ -1379,7 +1363,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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
@@ -1392,7 +1375,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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():
@@ -1466,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()
-        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):
@@ -1518,6 +1487,64 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                       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"""
@@ -1533,7 +1560,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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
@@ -1541,16 +1567,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         """
         n = self.dimension()
         R = self.coordinate_polynomial_ring()
-        vars = R.gens()
         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():
@@ -1631,7 +1650,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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
@@ -1642,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::
 
-            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
@@ -1657,8 +1674,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         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):
@@ -1680,7 +1697,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 
 
-class RationalBasisEJA(FiniteDimensionalEJA):
+class RationalBasisEJA(EJA):
     r"""
     Algebras whose supplied basis elements have all rational entries.
 
@@ -1735,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.
-            self._rational_algebra = FiniteDimensionalEJA(
+            self._rational_algebra = EJA(
                                        basis,
                                        jordan_product,
                                        inner_product,
@@ -1771,7 +1788,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
 
             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
@@ -1783,14 +1800,13 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             # 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() )
 
-        # 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())
@@ -1799,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 )
 
-class ConcreteEJA(FiniteDimensionalEJA):
+class ConcreteEJA(EJA):
     r"""
     A class for the Euclidean Jordan algebras that we know by name.
 
@@ -1817,7 +1833,6 @@ class ConcreteEJA(FiniteDimensionalEJA):
     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
@@ -1828,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::
 
-        sage: set_random_seed()
         sage: J = ConcreteEJA.random_instance()
         sage: x = J.random_element()
         sage: x.operator().is_self_adjoint()
@@ -1902,11 +1916,11 @@ class ConcreteEJA(FiniteDimensionalEJA):
         return eja_class.random_instance(max_dimension, *args, **kwargs)
 
 
-class MatrixEJA(FiniteDimensionalEJA):
+class HermitianMatrixEJA(EJA):
     @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`
@@ -1919,41 +1933,37 @@ class MatrixEJA(FiniteDimensionalEJA):
             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::
 
-            sage: set_random_seed()
             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: set_random_seed()
             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: set_random_seed()
             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: set_random_seed()
             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
 
@@ -2058,7 +2068,7 @@ class MatrixEJA(FiniteDimensionalEJA):
         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
@@ -2091,7 +2101,6 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     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)
@@ -2100,7 +2109,6 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     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()
@@ -2152,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,
@@ -2179,20 +2187,10 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         ...
         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`::
 
-        sage: set_random_seed()
         sage: d = ComplexHermitianEJA._max_random_instance_dimension()
         sage: n = ComplexHermitianEJA._max_random_instance_size(d)
         sage: J = ComplexHermitianEJA(n)
@@ -2201,7 +2199,6 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     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()
@@ -2253,7 +2250,7 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         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
@@ -2278,7 +2275,6 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     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)
@@ -2287,7 +2283,6 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     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()
@@ -2343,11 +2338,11 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
-class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
+class OctonionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA):
     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
 
@@ -2369,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: J = FiniteDimensionalEJA(basis,
+        sage: J = EJA(basis,
         ....:                          jp,
         ....:                          ip,
         ....:                          field=QQ,
@@ -2433,7 +2428,7 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
     @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.
@@ -2660,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::
 
-        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)
@@ -2822,7 +2816,6 @@ class JordanSpinEJA(BilinearFormEJA):
 
         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)
@@ -2921,7 +2914,7 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA):
         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
@@ -2943,7 +2936,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
     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])
@@ -3080,7 +3072,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
     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
@@ -3299,7 +3290,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
         The answer never changes::
 
-            sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
@@ -3315,7 +3305,7 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         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):
@@ -3389,7 +3379,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
         The answer never changes::
 
-            sage: set_random_seed()
             sage: J1 = random_eja()
             sage: J2 = random_eja()
             sage: J = cartesian_product([J1,J2])
@@ -3402,7 +3391,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         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])
@@ -3425,11 +3413,39 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         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):
@@ -3439,7 +3455,7 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+        sage: from mjo.eja.eja_algebra import (EJA,
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
@@ -3469,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: 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
@@ -3507,7 +3523,6 @@ def random_eja(max_dimension=None, *args, **kwargs):
 
     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
@@ -3533,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])
+
+
+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) )