]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/euclidean_jordan_algebra.py
eja: turn the other simple EJA constructors into classes.
[sage.d.git] / mjo / eja / euclidean_jordan_algebra.py
index e30c34128b564499d0485c3d0f9fcb68a826e6c7..32481621975ab3aabc217eae5ef72d9abf191980 100644 (file)
@@ -50,7 +50,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                                  inner_product=inner_product)
 
 
                                  inner_product=inner_product)
 
 
-    def __init__(self, field,
+    def __init__(self,
+                 field,
                  mult_table,
                  names='e',
                  assume_associative=False,
                  mult_table,
                  names='e',
                  assume_associative=False,
@@ -94,6 +95,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         The inner product associated with this Euclidean Jordan algebra.
 
         Will default to the trace inner product if nothing else.
         The inner product associated with this Euclidean Jordan algebra.
 
         Will default to the trace inner product if nothing else.
+
+        EXAMPLES:
+
+        The inner product must satisfy its axiom for this algebra to truly
+        be a Euclidean Jordan Algebra::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: z = J.random_element()
+            sage: (x*y).inner_product(z) == y.inner_product(x*z)
+            True
+
         """
         if (not x in self) or (not y in self):
             raise TypeError("arguments must live in this algebra")
         """
         if (not x in self) or (not y in self):
             raise TypeError("arguments must live in this algebra")
@@ -118,7 +133,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         EXAMPLES::
 
 
         EXAMPLES::
 
-            sage: J = RealSymmetricSimpleEJA(2)
+            sage: J = RealSymmetricEJA(2)
             sage: J.basis()
             Family (e0, e1, e2)
             sage: J.natural_basis()
             sage: J.basis()
             Family (e0, e1, e2)
             sage: J.natural_basis()
@@ -129,7 +144,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         ::
 
 
         ::
 
-            sage: J = JordanSpinSimpleEJA(2)
+            sage: J = JordanSpinEJA(2)
             sage: J.basis()
             Family (e0, e1)
             sage: J.natural_basis()
             sage: J.basis()
             Family (e0, e1)
             sage: J.natural_basis()
@@ -160,6 +175,51 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         An element of a Euclidean Jordan algebra.
         """
 
         An element of a Euclidean Jordan algebra.
         """
 
+        def __init__(self, A, elt=None):
+            """
+            EXAMPLES:
+
+            The identity in `S^n` is converted to the identity in the EJA::
+
+                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 = RealSymmetricEJA(3)
+                sage: A = matrix(QQ,3, lambda i,j: i-j)
+                sage: J(A)
+                Traceback (most recent call last):
+                ...
+                ArithmeticError: vector is not in free module
+
+            """
+            # Goal: if we're given a matrix, and if it lives in our
+            # parent algebra's "natural ambient space," convert it
+            # into an algebra element.
+            #
+            # The catch is, we make a recursive call after converting
+            # the given matrix into a vector that lives in the algebra.
+            # This we need to try the parent class initializer first,
+            # to avoid recursing forever if we're given something that
+            # already fits into the algebra, but also happens to live
+            # in the parent's "natural ambient space" (this happens with
+            # vectors in R^n).
+            try:
+                FiniteDimensionalAlgebraElement.__init__(self, A, elt)
+            except ValueError:
+                natural_basis = A.natural_basis()
+                if elt in natural_basis[0].matrix_space():
+                    # Thanks for nothing! Matrix spaces aren't vector
+                    # spaces in Sage, so we have to figure out its
+                    # natural-basis coordinates ourselves.
+                    V = VectorSpace(elt.base_ring(), elt.nrows()**2)
+                    W = V.span( _mat2vec(s) for s in natural_basis )
+                    coords =  W.coordinates(_mat2vec(elt))
+                    FiniteDimensionalAlgebraElement.__init__(self, A, coords)
+
         def __pow__(self, n):
             """
             Return ``self`` raised to the power ``n``.
         def __pow__(self, n):
             """
             Return ``self`` raised to the power ``n``.
@@ -236,7 +296,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`)::
 
             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)
                 sage: x = vector(QQ,[1,2,3])
                 sage: y = vector(QQ,[4,5,6])
                 sage: x.inner_product(y)
@@ -249,7 +309,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             so the inner product of the identity matrix with itself
             should be the `n`::
 
             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
 
                 sage: J.one().inner_product(J.one())
                 3
 
@@ -258,7 +318,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             part because the product of Hermitian matrices may not be
             Hermitian::
 
             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
 
                 sage: J.one().inner_product(J.one())
                 3
 
@@ -325,12 +391,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             EXAMPLES::
 
 
             EXAMPLES::
 
-                sage: J = JordanSpinSimpleEJA(2)
+                sage: J = JordanSpinEJA(2)
                 sage: e0,e1 = J.gens()
                 sage: x = e0 + e1
                 sage: x.det()
                 0
                 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()
                 sage: e0,e1,e2 = J.gens()
                 sage: x = e0 + e1 + e2
                 sage: x.det()
@@ -359,7 +425,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
                 sage: set_random_seed()
                 sage: n = ZZ.random_element(1,10)
 
                 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()
                 sage: x = J.random_element()
                 sage: while x.is_zero():
                 ....:     x = J.random_element()
@@ -489,7 +555,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             The identity element always has degree one, but any element
             linearly-independent from it is regular::
 
             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
                 sage: J.one().is_regular()
                 False
                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
@@ -514,7 +580,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             EXAMPLES::
 
 
             EXAMPLES::
 
-                sage: J = JordanSpinSimpleEJA(4)
+                sage: J = JordanSpinEJA(4)
                 sage: J.one().degree()
                 1
                 sage: e0,e1,e2,e3 = J.gens()
                 sage: J.one().degree()
                 1
                 sage: e0,e1,e2,e3 = J.gens()
@@ -526,7 +592,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
                 sage: set_random_seed()
                 sage: n = ZZ.random_element(1,10)
 
                 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
                 sage: x = J.random_element()
                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
                 True
@@ -558,7 +624,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
                 sage: set_random_seed()
                 sage: n = ZZ.random_element(2,10)
 
                 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: y = J.random_element()
                 sage: while y == y.coefficient(0)*J.one():
                 ....:     y = J.random_element()
@@ -604,7 +670,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             EXAMPLES::
 
 
             EXAMPLES::
 
-                sage: J = ComplexHermitianSimpleEJA(3)
+                sage: J = ComplexHermitianEJA(3)
                 sage: J.one()
                 e0 + e5 + e8
                 sage: J.one().natural_representation()
                 sage: J.one()
                 e0 + e5 + e8
                 sage: J.one().natural_representation()
@@ -615,6 +681,25 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 [0 0 0 0 1 0]
                 [0 0 0 0 0 1]
 
                 [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()
             """
             B = self.parent().natural_basis()
             W = B[0].matrix_space()
@@ -699,7 +784,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
                 sage: set_random_seed()
                 sage: n = ZZ.random_element(1,10)
 
                 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]
                 sage: x = J.random_element()
                 sage: x_vec = x.vector()
                 sage: x0 = x_vec[0]
@@ -846,7 +931,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: c = J.random_element().subalgebra_idempotent()
                 sage: c^2 == c
                 True
                 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
                 sage: c = J.random_element().subalgebra_idempotent()
                 sage: c^2 == c
                 True
@@ -902,7 +987,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             EXAMPLES::
 
 
             EXAMPLES::
 
-                sage: J = JordanSpinSimpleEJA(3)
+                sage: J = JordanSpinEJA(3)
                 sage: e0,e1,e2 = J.gens()
                 sage: x = e0 + e1 + e2
                 sage: x.trace()
                 sage: e0,e1,e2 = J.gens()
                 sage: x = e0 + e1 + e2
                 sage: x.trace()
@@ -983,6 +1068,12 @@ def random_eja():
       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
         product.
 
       * 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.
 
     Later this might be extended to return Cartesian products of the
     EJAs above.
 
@@ -994,9 +1085,10 @@ def random_eja():
     """
     n = ZZ.random_element(1,5)
     constructor = choice([eja_rn,
     """
     n = ZZ.random_element(1,5)
     constructor = choice([eja_rn,
-                          JordanSpinSimpleEJA,
-                          RealSymmetricSimpleEJA,
-                          ComplexHermitianSimpleEJA])
+                          JordanSpinEJA,
+                          RealSymmetricEJA,
+                          ComplexHermitianEJA,
+                          QuaternionHermitianEJA])
     return constructor(n, field=QQ)
 
 
     return constructor(n, field=QQ)
 
 
@@ -1057,6 +1149,48 @@ def _complex_hermitian_basis(n, field=QQ):
     return tuple(S)
 
 
     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())
 
 def _mat2vec(m):
         return vector(m.base_ring(), m.list())
 
@@ -1123,13 +1257,27 @@ def _embed_complex_matrix(M):
         sage: x2 = F(1 + 2*i)
         sage: x3 = F(-i)
         sage: x4 = F(6)
         sage: x2 = F(1 + 2*i)
         sage: x3 = F(-i)
         sage: x4 = F(6)
-        sage: M = matrix(F,2,[x1,x2,x3,x4])
+        sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
         sage: _embed_complex_matrix(M)
         sage: _embed_complex_matrix(M)
-        [ 4  2| 1 -2]
-        [-2  4| 2  1]
+        [ 4 -2| 1  2]
+        [ 2  4|-2  1]
         [-----+-----]
         [-----+-----]
-        [ 0  1| 6  0]
-        [-1  0| 0  6]
+        [ 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()
 
     """
     n = M.nrows()
@@ -1140,7 +1288,7 @@ def _embed_complex_matrix(M):
     for z in M.list():
         a = z.real()
         b = z.imag()
     for z in M.list():
         a = z.real()
         b = z.imag()
-        blocks.append(matrix(field, 2, [[a,-b],[b,a]]))
+        blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
 
     # We can drop the imaginaries here.
     return block_matrix(field.base_ring(), n, blocks)
 
     # We can drop the imaginaries here.
     return block_matrix(field.base_ring(), n, blocks)
@@ -1157,8 +1305,19 @@ def _unembed_complex_matrix(M):
         ....:                 [ 9,  10, 11, 12],
         ....:                 [-10, 9, -12, 11] ])
         sage: _unembed_complex_matrix(A)
         ....:                 [ 9,  10, 11, 12],
         ....:                 [-10, 9, -12, 11] ])
         sage: _unembed_complex_matrix(A)
-        [  -2*i + 1   -4*i + 3]
-        [ -10*i + 9 -12*i + 11]
+        [  2*i + 1   4*i + 3]
+        [ 10*i + 9 12*i + 11]
+
+    TESTS:
+
+    Unembedding is the inverse of embedding::
+
+        sage: set_random_seed()
+        sage: F = QuadraticField(-1, 'i')
+        sage: M = random_matrix(F, 3)
+        sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
+        True
+
     """
     n = ZZ(M.nrows())
     if M.ncols() != n:
     """
     n = ZZ(M.nrows())
     if M.ncols() != n:
@@ -1176,14 +1335,123 @@ def _unembed_complex_matrix(M):
         for j in xrange(n/2):
             submat = M[2*k:2*k+2,2*j:2*j+2]
             if submat[0,0] != submat[1,1]:
         for j in xrange(n/2):
             submat = M[2*k:2*k+2,2*j:2*j+2]
             if submat[0,0] != submat[1,1]:
-                raise ValueError('bad real submatrix')
+                raise ValueError('bad on-diagonal submatrix')
             if submat[0,1] != -submat[1,0]:
             if submat[0,1] != -submat[1,0]:
-                raise ValueError('bad imag submatrix')
-            z = submat[0,0] + submat[1,0]*i
+                raise ValueError('bad off-diagonal submatrix')
+            z = submat[0,0] + submat[0,1]*i
             elements.append(z)
 
     return matrix(F, n/2, elements)
 
             elements.append(z)
 
     return matrix(F, n/2, elements)
 
+
+def _embed_quaternion_matrix(M):
+    """
+    Embed the n-by-n quaternion matrix ``M`` into the space of real
+    matrices of size 4n-by-4n by first sending each quaternion entry
+    `z = a + bi + cj + dk` to the block-complex matrix
+    ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
+    a real matrix.
+
+    EXAMPLES::
+
+        sage: Q = QuaternionAlgebra(QQ,-1,-1)
+        sage: i,j,k = Q.gens()
+        sage: x = 1 + 2*i + 3*j + 4*k
+        sage: M = matrix(Q, 1, [[x]])
+        sage: _embed_quaternion_matrix(M)
+        [ 1  2  3  4]
+        [-2  1 -4  3]
+        [-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()
+    if M.ncols() != n:
+        raise ValueError("the matrix 'M' must be square")
+
+    F = QuadraticField(-1, 'i')
+    i = F.gen()
+
+    blocks = []
+    for z in M.list():
+        t = z.coefficient_tuple()
+        a = t[0]
+        b = t[1]
+        c = t[2]
+        d = t[3]
+        cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
+                                    [-c + d*i, a - b*i]])
+        blocks.append(_embed_complex_matrix(cplx_matrix))
+
+    # We should have real entries by now, so use the realest field
+    # we've got for the return value.
+    return block_matrix(quaternions.base_ring(), n, blocks)
+
+
+def _unembed_quaternion_matrix(M):
+    """
+    The inverse of _embed_quaternion_matrix().
+
+    EXAMPLES::
+
+        sage: M = matrix(QQ, [[ 1,  2,  3,  4],
+        ....:                 [-2,  1, -4,  3],
+        ....:                 [-3,  4,  1, -2],
+        ....:                 [-4, -3,  2,  1]])
+        sage: _unembed_quaternion_matrix(M)
+        [1 + 2*i + 3*j + 4*k]
+
+    TESTS:
+
+    Unembedding is the inverse of embedding::
+
+        sage: set_random_seed()
+        sage: Q = QuaternionAlgebra(QQ, -1, -1)
+        sage: M = random_matrix(Q, 3)
+        sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
+        True
+
+    """
+    n = ZZ(M.nrows())
+    if M.ncols() != n:
+        raise ValueError("the matrix 'M' must be square")
+    if not n.mod(4).is_zero():
+        raise ValueError("the matrix 'M' must be a complex embedding")
+
+    Q = QuaternionAlgebra(QQ,-1,-1)
+    i,j,k = Q.gens()
+
+    # Go top-left to bottom-right (reading order), converting every
+    # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
+    # quaternion block.
+    elements = []
+    for l in xrange(n/4):
+        for m in xrange(n/4):
+            submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
+            if submat[0,0] != submat[1,1].conjugate():
+                raise ValueError('bad on-diagonal submatrix')
+            if submat[0,1] != -submat[1,0].conjugate():
+                raise ValueError('bad off-diagonal submatrix')
+            z  = submat[0,0].real() + submat[0,0].imag()*i
+            z += submat[0,1].real()*j + submat[0,1].imag()*k
+            elements.append(z)
+
+    return matrix(Q, n/4, elements)
+
+
 # The usual inner product on R^n.
 def _usual_ip(x,y):
     return x.vector().inner_product(y.vector())
 # The usual inner product on R^n.
 def _usual_ip(x,y):
     return x.vector().inner_product(y.vector())
@@ -1197,7 +1465,7 @@ def _matrix_ip(X,Y):
     return (X_mat*Y_mat).trace()
 
 
     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
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -1205,7 +1473,7 @@ def RealSymmetricSimpleEJA(n, field=QQ):
 
     EXAMPLES::
 
 
     EXAMPLES::
 
-        sage: J = RealSymmetricSimpleEJA(2)
+        sage: J = RealSymmetricEJA(2)
         sage: e0, e1, e2 = J.gens()
         sage: e0*e0
         e0
         sage: e0, e1, e2 = J.gens()
         sage: e0*e0
         e0
@@ -1220,22 +1488,44 @@ def RealSymmetricSimpleEJA(n, field=QQ):
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
 
         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
 
         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,
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -1248,47 +1538,110 @@ def ComplexHermitianSimpleEJA(n, field=QQ):
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
 
         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
 
         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.
     """
     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 =
     """
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the usual inner product and jordan product ``x*y =
@@ -1299,7 +1652,7 @@ def JordanSpinSimpleEJA(n, field=QQ):
 
     This multiplication table can be verified by hand::
 
 
     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
         sage: e0,e1,e2,e3 = J.gens()
         sage: e0*e0
         e0
@@ -1316,31 +1669,34 @@ def JordanSpinSimpleEJA(n, field=QQ):
         sage: e2*e3
         0
 
         sage: e2*e3
         0
 
-    In one dimension, this is the reals under multiplication::
+    """
+    @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)
+
+        fdeja = super(JordanSpinEJA, cls)
+        return fdeja.__classcall_private__(cls, field, Qs)
 
 
-        sage: J1 = JordanSpinSimpleEJA(1)
-        sage: J2 = eja_rn(1)
-        sage: J1 == J2
-        True
+    def rank(self):
+        """
+        Return the rank of this Jordan Spin Algebra.
 
 
-    """
-    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)
+        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).
+        """
+        return min(self.dimension(),2)
+
+    def inner_product(self, x, y):
+        return _usual_ip(x,y)