]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/euclidean_jordan_algebra.py
eja: use 't' for the minimal polynomial variable name.
[sage.d.git] / mjo / eja / euclidean_jordan_algebra.py
index 7c7c2d6ba3f0d54ab34c89d37e47e6bbd5844e40..8f508bc305304d249440294267cc56f6099c631b 100644 (file)
@@ -21,8 +21,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                               assume_associative=False,
                               category=None,
                               rank=None,
-                              natural_basis=None,
-                              inner_product=None):
+                              natural_basis=None):
         n = len(mult_table)
         mult_table = [b.base_extend(field) for b in mult_table]
         for b in mult_table:
@@ -46,18 +45,17 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                                  names=names,
                                  category=cat,
                                  rank=rank,
-                                 natural_basis=natural_basis,
-                                 inner_product=inner_product)
+                                 natural_basis=natural_basis)
 
 
-    def __init__(self, field,
+    def __init__(self,
+                 field,
                  mult_table,
                  names='e',
                  assume_associative=False,
                  category=None,
                  rank=None,
-                 natural_basis=None,
-                 inner_product=None):
+                 natural_basis=None):
         """
         EXAMPLES:
 
@@ -73,7 +71,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         """
         self._rank = rank
         self._natural_basis = natural_basis
-        self._inner_product = inner_product
+        self._multiplication_table = mult_table
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
                      mult_table,
@@ -89,11 +87,61 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         return fmt.format(self.degree(), self.base_ring())
 
 
+    def characteristic_polynomial(self):
+        r = self.rank()
+        n = self.dimension()
+
+        names = ['X' + str(i) for i in range(1,n+1)]
+        R = PolynomialRing(self.base_ring(), names)
+        J = FiniteDimensionalEuclideanJordanAlgebra(R,
+                                                    self._multiplication_table,
+                                                    rank=r)
+        x0 = J.zero()
+        c = 1
+        for g in J.gens():
+            x0 += c*g
+            c +=1
+        if not x0.is_regular():
+            raise ValueError("don't know a regular element")
+
+        # Get the vector space (as opposed to module) so that
+        # span_of_basis() works.
+        V = x0.vector().parent().ambient_vector_space()
+        V1 = V.span_of_basis( (x0**k).vector() for k in range(r) )
+        B = V1.basis() + V1.complement().basis()
+        W = V.span_of_basis(B)
+
+        def e(k):
+            # The coordinates of e_k with respect to the basis B.
+            # But, the e_k are elements of B...
+            return identity_matrix(J.base_ring(), n).column(k-1).column()
+
+        # A matrix implementation 1
+        x = J(vector(R, R.gens()))
+        l1 = [column_matrix(W.coordinates((x**k).vector())) for k in range(r)]
+        l2 = [e(k) for k in range(r+1, n+1)]
+        A_of_x = block_matrix(1, n, (l1 + l2))
+        xr = W.coordinates((x**r).vector())
+        a = []
+        for i in range(n):
+            A_cols = A.columns()
+            A_cols[i] = xr
+            numerator = column_matrix(A.base_ring(), A_cols).det()
+            denominator = A.det()
+            ai = numerator/denominator
+            a.append(ai)
+
+        # Note: all entries past the rth should be zero.
+        return a
+
+
     def inner_product(self, x, y):
         """
         The inner product associated with this Euclidean Jordan algebra.
 
-        Will default to the trace inner product if nothing else.
+        Defaults to the trace inner product, but can be overridden by
+        subclasses if they are sure that the necessary properties are
+        satisfied.
 
         EXAMPLES:
 
@@ -111,10 +159,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         """
         if (not x in self) or (not y in self):
             raise TypeError("arguments must live in this algebra")
-        if self._inner_product is None:
-            return x.trace_inner_product(y)
-        else:
-            return self._inner_product(x,y)
+        return x.trace_inner_product(y)
 
 
     def natural_basis(self):
@@ -132,7 +177,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         EXAMPLES::
 
-            sage: J = RealSymmetricSimpleEJA(2)
+            sage: J = RealSymmetricEJA(2)
             sage: J.basis()
             Family (e0, e1, e2)
             sage: J.natural_basis()
@@ -143,7 +188,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         ::
 
-            sage: J = JordanSpinSimpleEJA(2)
+            sage: J = JordanSpinEJA(2)
             sage: J.basis()
             Family (e0, e1)
             sage: J.natural_basis()
@@ -180,14 +225,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             The identity in `S^n` is converted to the identity in the EJA::
 
-                sage: J = RealSymmetricSimpleEJA(3)
+                sage: J = RealSymmetricEJA(3)
                 sage: I = identity_matrix(QQ,3)
                 sage: J(I) == J.one()
                 True
 
             This skew-symmetric matrix can't be represented in the EJA::
 
-                sage: J = RealSymmetricSimpleEJA(3)
+                sage: J = RealSymmetricEJA(3)
                 sage: A = matrix(QQ,3, lambda i,j: i-j)
                 sage: J(A)
                 Traceback (most recent call last):
@@ -295,7 +340,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             inner product on `R^n` (this example only works because the
             basis for the Jordan algebra is the standard basis in `R^n`)::
 
-                sage: J = JordanSpinSimpleEJA(3)
+                sage: J = JordanSpinEJA(3)
                 sage: x = vector(QQ,[1,2,3])
                 sage: y = vector(QQ,[4,5,6])
                 sage: x.inner_product(y)
@@ -308,7 +353,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             so the inner product of the identity matrix with itself
             should be the `n`::
 
-                sage: J = RealSymmetricSimpleEJA(3)
+                sage: J = RealSymmetricEJA(3)
                 sage: J.one().inner_product(J.one())
                 3
 
@@ -317,7 +362,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             part because the product of Hermitian matrices may not be
             Hermitian::
 
-                sage: J = ComplexHermitianSimpleEJA(3)
+                sage: J = ComplexHermitianEJA(3)
+                sage: J.one().inner_product(J.one())
+                3
+
+            Ditto for the quaternions::
+
+                sage: J = QuaternionHermitianEJA(3)
                 sage: J.one().inner_product(J.one())
                 3
 
@@ -384,12 +435,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             EXAMPLES::
 
-                sage: J = JordanSpinSimpleEJA(2)
+                sage: J = JordanSpinEJA(2)
                 sage: e0,e1 = J.gens()
                 sage: x = e0 + e1
                 sage: x.det()
                 0
-                sage: J = JordanSpinSimpleEJA(3)
+                sage: J = JordanSpinEJA(3)
                 sage: e0,e1,e2 = J.gens()
                 sage: x = e0 + e1 + e2
                 sage: x.det()
@@ -418,7 +469,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
                 sage: set_random_seed()
                 sage: n = ZZ.random_element(1,10)
-                sage: J = JordanSpinSimpleEJA(n)
+                sage: J = JordanSpinEJA(n)
                 sage: x = J.random_element()
                 sage: while x.is_zero():
                 ....:     x = J.random_element()
@@ -490,8 +541,36 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             We can't use the superclass method because it relies on
             the algebra being associative.
+
+            ALGORITHM:
+
+            The usual way to do this is to check if the determinant is
+            zero, but we need the characteristic polynomial for the
+            determinant. The minimal polynomial is a lot easier to get,
+            so we use Corollary 2 in Chapter V of Koecher to check
+            whether or not the paren't algebra's zero element is a root
+            of this element's minimal polynomial.
+
+            TESTS:
+
+            The identity element is always invertible::
+
+                sage: set_random_seed()
+                sage: J = random_eja()
+                sage: J.one().is_invertible()
+                True
+
+            The zero element is never invertible::
+
+                sage: set_random_seed()
+                sage: J = random_eja()
+                sage: J.zero().is_invertible()
+                False
+
             """
-            return not self.det().is_zero()
+            zero = self.parent().zero()
+            p = self.minimal_polynomial()
+            return not (p(zero) == zero)
 
 
         def is_nilpotent(self):
@@ -548,7 +627,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             The identity element always has degree one, but any element
             linearly-independent from it is regular::
 
-                sage: J = JordanSpinSimpleEJA(5)
+                sage: J = JordanSpinEJA(5)
                 sage: J.one().is_regular()
                 False
                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
@@ -573,7 +652,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             EXAMPLES::
 
-                sage: J = JordanSpinSimpleEJA(4)
+                sage: J = JordanSpinEJA(4)
                 sage: J.one().degree()
                 1
                 sage: e0,e1,e2,e3 = J.gens()
@@ -585,7 +664,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
                 sage: set_random_seed()
                 sage: n = ZZ.random_element(1,10)
-                sage: J = JordanSpinSimpleEJA(n)
+                sage: J = JordanSpinEJA(n)
                 sage: x = J.random_element()
                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
                 True
@@ -596,14 +675,30 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         def minimal_polynomial(self):
             """
-            EXAMPLES::
+            Return the minimal polynomial of this element,
+            as a function of the variable `t`.
+
+            ALGORITHM:
+
+            We restrict ourselves to the associative subalgebra
+            generated by this element, and then return the minimal
+            polynomial of this element's operator matrix (in that
+            subalgebra). This works by Baes Proposition 2.3.16.
+
+            TESTS:
+
+            The minimal polynomial of the identity and zero elements are
+            always the same::
 
                 sage: set_random_seed()
-                sage: x = random_eja().random_element()
-                sage: x.degree() == x.minimal_polynomial().degree()
-                True
+                sage: J = random_eja()
+                sage: J.one().minimal_polynomial()
+                t - 1
+                sage: J.zero().minimal_polynomial()
+                t
 
-            ::
+            The degree of an element is (by one definition) the degree
+            of its minimal polynomial::
 
                 sage: set_random_seed()
                 sage: x = random_eja().random_element()
@@ -617,38 +712,30 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
                 sage: set_random_seed()
                 sage: n = ZZ.random_element(2,10)
-                sage: J = JordanSpinSimpleEJA(n)
+                sage: J = JordanSpinEJA(n)
                 sage: y = J.random_element()
                 sage: while y == y.coefficient(0)*J.one():
                 ....:     y = J.random_element()
                 sage: y0 = y.vector()[0]
                 sage: y_bar = y.vector()[1:]
                 sage: actual = y.minimal_polynomial()
-                sage: x = SR.symbol('x', domain='real')
-                sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
+                sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
+                sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
                 sage: bool(actual == expected)
                 True
 
             """
-            # The element we're going to call "minimal_polynomial()" on.
-            # Either myself, interpreted as an element of a finite-
-            # dimensional algebra, or an element of an associative
-            # subalgebra.
-            elt = None
-
-            if self.parent().is_associative():
-                elt = FiniteDimensionalAlgebraElement(self.parent(), self)
-            else:
-                V = self.span_of_powers()
-                assoc_subalg = self.subalgebra_generated_by()
-                # Mis-design warning: the basis used for span_of_powers()
-                # and subalgebra_generated_by() must be the same, and in
-                # the same order!
-                elt = assoc_subalg(V.coordinates(self.vector()))
+            V = self.span_of_powers()
+            assoc_subalg = self.subalgebra_generated_by()
+            # Mis-design warning: the basis used for span_of_powers()
+            # and subalgebra_generated_by() must be the same, and in
+            # the same order!
+            elt = assoc_subalg(V.coordinates(self.vector()))
 
-            # Recursive call, but should work since elt lives in an
-            # associative algebra.
-            return elt.minimal_polynomial()
+            # We get back a symbolic polynomial in 'x' but want a real
+            # polynomial in 't'.
+            p_of_x = elt.operator_matrix().minimal_polynomial()
+            return p_of_x.change_variable_name('t')
 
 
         def natural_representation(self):
@@ -663,7 +750,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             EXAMPLES::
 
-                sage: J = ComplexHermitianSimpleEJA(3)
+                sage: J = ComplexHermitianEJA(3)
                 sage: J.one()
                 e0 + e5 + e8
                 sage: J.one().natural_representation()
@@ -674,6 +761,25 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 [0 0 0 0 1 0]
                 [0 0 0 0 0 1]
 
+            ::
+
+                sage: J = QuaternionHermitianEJA(3)
+                sage: J.one()
+                e0 + e9 + e14
+                sage: J.one().natural_representation()
+                [1 0 0 0 0 0 0 0 0 0 0 0]
+                [0 1 0 0 0 0 0 0 0 0 0 0]
+                [0 0 1 0 0 0 0 0 0 0 0 0]
+                [0 0 0 1 0 0 0 0 0 0 0 0]
+                [0 0 0 0 1 0 0 0 0 0 0 0]
+                [0 0 0 0 0 1 0 0 0 0 0 0]
+                [0 0 0 0 0 0 1 0 0 0 0 0]
+                [0 0 0 0 0 0 0 1 0 0 0 0]
+                [0 0 0 0 0 0 0 0 1 0 0 0]
+                [0 0 0 0 0 0 0 0 0 1 0 0]
+                [0 0 0 0 0 0 0 0 0 0 1 0]
+                [0 0 0 0 0 0 0 0 0 0 0 1]
+
             """
             B = self.parent().natural_basis()
             W = B[0].matrix_space()
@@ -758,7 +864,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
                 sage: set_random_seed()
                 sage: n = ZZ.random_element(1,10)
-                sage: J = JordanSpinSimpleEJA(n)
+                sage: J = JordanSpinEJA(n)
                 sage: x = J.random_element()
                 sage: x_vec = x.vector()
                 sage: x0 = x_vec[0]
@@ -832,7 +938,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             # The dimension of the subalgebra can't be greater than
             # the big algebra, so just put everything into a list
             # and let span() get rid of the excess.
-            V = self.vector().parent()
+            #
+            # We do the extra ambient_vector_space() in case we're messing
+            # with polynomials and the direct parent is a module.
+            V = self.vector().parent().ambient_vector_space()
             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 
 
@@ -901,11 +1010,11 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             TESTS::
 
                 sage: set_random_seed()
-                sage: J = eja_rn(5)
+                sage: J = RealCartesianProductEJA(5)
                 sage: c = J.random_element().subalgebra_idempotent()
                 sage: c^2 == c
                 True
-                sage: J = JordanSpinSimpleEJA(5)
+                sage: J = JordanSpinEJA(5)
                 sage: c = J.random_element().subalgebra_idempotent()
                 sage: c^2 == c
                 True
@@ -961,7 +1070,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             EXAMPLES::
 
-                sage: J = JordanSpinSimpleEJA(3)
+                sage: J = JordanSpinEJA(3)
                 sage: e0,e1,e2 = J.gens()
                 sage: x = e0 + e1 + e2
                 sage: x.trace()
@@ -985,16 +1094,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             return (self*other).trace()
 
 
-def eja_rn(dimension, field=QQ):
+class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
 
+    Note: this is nothing more than the Cartesian product of ``n``
+    copies of the spin algebra. Once Cartesian product algebras
+    are implemented, this can go.
+
     EXAMPLES:
 
     This multiplication table can be verified by hand::
 
-        sage: J = eja_rn(3)
+        sage: J = RealCartesianProductEJA(3)
         sage: e0,e1,e2 = J.gens()
         sage: e0*e0
         e0
@@ -1010,19 +1123,21 @@ def eja_rn(dimension, field=QQ):
         e2
 
     """
-    # The FiniteDimensionalAlgebra constructor takes a list of
-    # matrices, the ith representing right multiplication by the ith
-    # basis element in the vector space. So if e_1 = (1,0,0), then
-    # right (Hadamard) multiplication of x by e_1 picks out the first
-    # component of x; and likewise for the ith basis element e_i.
-    Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
-           for i in xrange(dimension) ]
-
-    return FiniteDimensionalEuclideanJordanAlgebra(field,
-                                                   Qs,
-                                                   rank=dimension,
-                                                   inner_product=_usual_ip)
+    @staticmethod
+    def __classcall_private__(cls, n, field=QQ):
+        # The FiniteDimensionalAlgebra constructor takes a list of
+        # matrices, the ith representing right multiplication by the ith
+        # basis element in the vector space. So if e_1 = (1,0,0), then
+        # right (Hadamard) multiplication of x by e_1 picks out the first
+        # component of x; and likewise for the ith basis element e_i.
+        Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
+               for i in xrange(n) ]
+
+        fdeja = super(RealCartesianProductEJA, cls)
+        return fdeja.__classcall_private__(cls, field, Qs, rank=n)
 
+    def inner_product(self, x, y):
+        return _usual_ip(x,y)
 
 
 def random_eja():
@@ -1042,6 +1157,12 @@ def random_eja():
       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
         product.
 
+      * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
+        in the space of ``2n``-by-``2n`` real symmetric matrices.
+
+      * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
+        in the space of ``4n``-by-``4n`` real symmetric matrices.
+
     Later this might be extended to return Cartesian products of the
     EJAs above.
 
@@ -1052,10 +1173,11 @@ def random_eja():
 
     """
     n = ZZ.random_element(1,5)
-    constructor = choice([eja_rn,
-                          JordanSpinSimpleEJA,
-                          RealSymmetricSimpleEJA,
-                          ComplexHermitianSimpleEJA])
+    constructor = choice([RealCartesianProductEJA,
+                          JordanSpinEJA,
+                          RealSymmetricEJA,
+                          ComplexHermitianEJA,
+                          QuaternionHermitianEJA])
     return constructor(n, field=QQ)
 
 
@@ -1116,6 +1238,48 @@ def _complex_hermitian_basis(n, field=QQ):
     return tuple(S)
 
 
+def _quaternion_hermitian_basis(n, field=QQ):
+    """
+    Returns a basis for the space of quaternion Hermitian n-by-n matrices.
+
+    TESTS::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
+        True
+
+    """
+    Q = QuaternionAlgebra(QQ,-1,-1)
+    I,J,K = Q.gens()
+
+    # This is like the symmetric case, but we need to be careful:
+    #
+    #   * We want conjugate-symmetry, not just symmetry.
+    #   * The diagonal will (as a result) be real.
+    #
+    S = []
+    for i in xrange(n):
+        for j in xrange(i+1):
+            Eij = matrix(Q, n, lambda k,l: k==i and l==j)
+            if i == j:
+                Sij = _embed_quaternion_matrix(Eij)
+                S.append(Sij)
+            else:
+                # Beware, orthogonal but not normalized! The second,
+                # third, and fourth ones have a minus because they're
+                # conjugated.
+                Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
+                S.append(Sij_real)
+                Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
+                S.append(Sij_I)
+                Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
+                S.append(Sij_J)
+                Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
+                S.append(Sij_K)
+    return tuple(S)
+
+
 def _mat2vec(m):
         return vector(m.base_ring(), m.list())
 
@@ -1190,6 +1354,20 @@ def _embed_complex_matrix(M):
         [ 0 -1| 6  0]
         [ 1  0| 0  6]
 
+    TESTS:
+
+    Embedding is a homomorphism (isomorphism, in fact)::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(5)
+        sage: F = QuadraticField(-1, 'i')
+        sage: X = random_matrix(F, n)
+        sage: Y = random_matrix(F, n)
+        sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
+        sage: expected = _embed_complex_matrix(X*Y)
+        sage: actual == expected
+        True
+
     """
     n = M.nrows()
     if M.ncols() != n:
@@ -1219,7 +1397,9 @@ def _unembed_complex_matrix(M):
         [  2*i + 1   4*i + 3]
         [ 10*i + 9 12*i + 11]
 
-    TESTS::
+    TESTS:
+
+    Unembedding is the inverse of embedding::
 
         sage: set_random_seed()
         sage: F = QuadraticField(-1, 'i')
@@ -1273,6 +1453,18 @@ def _embed_quaternion_matrix(M):
         [-3  4  1 -2]
         [-4 -3  2  1]
 
+    Embedding is a homomorphism (isomorphism, in fact)::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(5)
+        sage: Q = QuaternionAlgebra(QQ,-1,-1)
+        sage: X = random_matrix(Q, n)
+        sage: Y = random_matrix(Q, n)
+        sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
+        sage: expected = _embed_quaternion_matrix(X*Y)
+        sage: actual == expected
+        True
+
     """
     quaternions = M.base_ring()
     n = M.nrows()
@@ -1311,7 +1503,9 @@ def _unembed_quaternion_matrix(M):
         sage: _unembed_quaternion_matrix(M)
         [1 + 2*i + 3*j + 4*k]
 
-    TESTS::
+    TESTS:
+
+    Unembedding is the inverse of embedding::
 
         sage: set_random_seed()
         sage: Q = QuaternionAlgebra(QQ, -1, -1)
@@ -1360,7 +1554,7 @@ def _matrix_ip(X,Y):
     return (X_mat*Y_mat).trace()
 
 
-def RealSymmetricSimpleEJA(n, field=QQ):
+class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -1368,7 +1562,7 @@ def RealSymmetricSimpleEJA(n, field=QQ):
 
     EXAMPLES::
 
-        sage: J = RealSymmetricSimpleEJA(2)
+        sage: J = RealSymmetricEJA(2)
         sage: e0, e1, e2 = J.gens()
         sage: e0*e0
         e0
@@ -1383,22 +1577,44 @@ def RealSymmetricSimpleEJA(n, field=QQ):
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
-        sage: J = RealSymmetricSimpleEJA(n)
+        sage: J = RealSymmetricEJA(n)
         sage: J.degree() == (n^2 + n)/2
         True
 
+    The Jordan multiplication is what we think it is::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = RealSymmetricEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: actual = (x*y).natural_representation()
+        sage: X = x.natural_representation()
+        sage: Y = y.natural_representation()
+        sage: expected = (X*Y + Y*X)/2
+        sage: actual == expected
+        True
+        sage: J(expected) == x*y
+        True
+
     """
-    S = _real_symmetric_basis(n, field=field)
-    (Qs, T) = _multiplication_table_from_matrix_basis(S)
+    @staticmethod
+    def __classcall_private__(cls, n, field=QQ):
+        S = _real_symmetric_basis(n, field=field)
+        (Qs, T) = _multiplication_table_from_matrix_basis(S)
 
-    return FiniteDimensionalEuclideanJordanAlgebra(field,
-                                                   Qs,
-                                                   rank=n,
-                                                   natural_basis=T,
-                                                   inner_product=_matrix_ip)
+        fdeja = super(RealSymmetricEJA, cls)
+        return fdeja.__classcall_private__(cls,
+                                           field,
+                                           Qs,
+                                           rank=n,
+                                           natural_basis=T)
+
+    def inner_product(self, x, y):
+        return _matrix_ip(x,y)
 
 
-def ComplexHermitianSimpleEJA(n, field=QQ):
+class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -1411,47 +1627,110 @@ def ComplexHermitianSimpleEJA(n, field=QQ):
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
-        sage: J = ComplexHermitianSimpleEJA(n)
+        sage: J = ComplexHermitianEJA(n)
         sage: J.degree() == n^2
         True
 
+    The Jordan multiplication is what we think it is::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = ComplexHermitianEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: actual = (x*y).natural_representation()
+        sage: X = x.natural_representation()
+        sage: Y = y.natural_representation()
+        sage: expected = (X*Y + Y*X)/2
+        sage: actual == expected
+        True
+        sage: J(expected) == x*y
+        True
+
     """
-    S = _complex_hermitian_basis(n)
-    (Qs, T) = _multiplication_table_from_matrix_basis(S)
+    @staticmethod
+    def __classcall_private__(cls, n, field=QQ):
+        S = _complex_hermitian_basis(n)
+        (Qs, T) = _multiplication_table_from_matrix_basis(S)
 
-    # Since a+bi on the diagonal is represented as
-    #
-    #   a + bi  = [  a  b  ]
-    #             [ -b  a  ],
-    #
-    # we'll double-count the "a" entries if we take the trace of
-    # the embedding.
-    ip = lambda X,Y: _matrix_ip(X,Y)/2
+        fdeja = super(ComplexHermitianEJA, cls)
+        return fdeja.__classcall_private__(cls,
+                                           field,
+                                           Qs,
+                                           rank=n,
+                                           natural_basis=T)
 
-    return FiniteDimensionalEuclideanJordanAlgebra(field,
-                                                   Qs,
-                                                   rank=n,
-                                                   natural_basis=T,
-                                                   inner_product=ip)
+    def inner_product(self, x, y):
+        # Since a+bi on the diagonal is represented as
+        #
+        #   a + bi  = [  a  b  ]
+        #             [ -b  a  ],
+        #
+        # we'll double-count the "a" entries if we take the trace of
+        # the embedding.
+        return _matrix_ip(x,y)/2
 
 
-def QuaternionHermitianSimpleEJA(n):
+class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
     real-part-of-trace inner product. It has dimension `2n^2 - n` over
     the reals.
-    """
-    pass
 
-def OctonionHermitianSimpleEJA(n):
-    """
-    This shit be crazy. It has dimension 27 over the reals.
+    TESTS:
+
+    The degree of this algebra is `n^2`::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = QuaternionHermitianEJA(n)
+        sage: J.degree() == 2*(n^2) - n
+        True
+
+    The Jordan multiplication is what we think it is::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = QuaternionHermitianEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: actual = (x*y).natural_representation()
+        sage: X = x.natural_representation()
+        sage: Y = y.natural_representation()
+        sage: expected = (X*Y + Y*X)/2
+        sage: actual == expected
+        True
+        sage: J(expected) == x*y
+        True
+
     """
-    n = 3
-    pass
+    @staticmethod
+    def __classcall_private__(cls, n, field=QQ):
+        S = _quaternion_hermitian_basis(n)
+        (Qs, T) = _multiplication_table_from_matrix_basis(S)
+
+        fdeja = super(QuaternionHermitianEJA, cls)
+        return fdeja.__classcall_private__(cls,
+                                           field,
+                                           Qs,
+                                           rank=n,
+                                           natural_basis=T)
 
-def JordanSpinSimpleEJA(n, field=QQ):
+    def inner_product(self, x, y):
+        # Since a+bi+cj+dk on the diagonal is represented as
+        #
+        #   a + bi +cj + dk = [  a  b  c  d]
+        #                     [ -b  a -d  c]
+        #                     [ -c  d  a -b]
+        #                     [ -d -c  b  a],
+        #
+        # we'll quadruple-count the "a" entries if we take the trace of
+        # the embedding.
+        return _matrix_ip(x,y)/4
+
+
+class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the usual inner product and jordan product ``x*y =
@@ -1462,7 +1741,7 @@ def JordanSpinSimpleEJA(n, field=QQ):
 
     This multiplication table can be verified by hand::
 
-        sage: J = JordanSpinSimpleEJA(4)
+        sage: J = JordanSpinEJA(4)
         sage: e0,e1,e2,e3 = J.gens()
         sage: e0*e0
         e0
@@ -1479,31 +1758,27 @@ def JordanSpinSimpleEJA(n, field=QQ):
         sage: e2*e3
         0
 
-    In one dimension, this is the reals under multiplication::
-
-        sage: J1 = JordanSpinSimpleEJA(1)
-        sage: J2 = eja_rn(1)
-        sage: J1 == J2
-        True
-
     """
-    Qs = []
-    id_matrix = identity_matrix(field, n)
-    for i in xrange(n):
-        ei = id_matrix.column(i)
-        Qi = zero_matrix(field, n)
-        Qi.set_row(0, ei)
-        Qi.set_column(0, ei)
-        Qi += diagonal_matrix(n, [ei[0]]*n)
-        # The addition of the diagonal matrix adds an extra ei[0] in the
-        # upper-left corner of the matrix.
-        Qi[0,0] = Qi[0,0] * ~field(2)
-        Qs.append(Qi)
-
-    # The rank of the spin factor algebra is two, UNLESS we're in a
-    # one-dimensional ambient space (the rank is bounded by the
-    # ambient dimension).
-    return FiniteDimensionalEuclideanJordanAlgebra(field,
-                                                   Qs,
-                                                   rank=min(n,2),
-                                                   inner_product=_usual_ip)
+    @staticmethod
+    def __classcall_private__(cls, n, field=QQ):
+        Qs = []
+        id_matrix = identity_matrix(field, n)
+        for i in xrange(n):
+            ei = id_matrix.column(i)
+            Qi = zero_matrix(field, n)
+            Qi.set_row(0, ei)
+            Qi.set_column(0, ei)
+            Qi += diagonal_matrix(n, [ei[0]]*n)
+            # The addition of the diagonal matrix adds an extra ei[0] in the
+            # upper-left corner of the matrix.
+            Qi[0,0] = Qi[0,0] * ~field(2)
+            Qs.append(Qi)
+
+        # The rank of the spin algebra is two, unless we're in a
+        # one-dimensional ambient space (because the rank is bounded by
+        # the ambient dimension).
+        fdeja = super(JordanSpinEJA, cls)
+        return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
+
+    def inner_product(self, x, y):
+        return _usual_ip(x,y)