]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/euclidean_jordan_algebra.py
eja: begin to stub out a characteristic_polynomial() for EJAs.
[sage.d.git] / mjo / eja / euclidean_jordan_algebra.py
index 7c1f249b4632766e4b1808acc0275b2e55b8b506..29248fcfaabfad33567070eaa65efffef7df87ff 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,18 +87,79 @@ 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:
+
+        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 ArgumentError("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)
+            raise TypeError("arguments must live in this algebra")
+        return x.trace_inner_product(y)
 
 
     def natural_basis(self):
@@ -118,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()
@@ -129,7 +188,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         ::
 
-            sage: J = JordanSpinSimpleEJA(2)
+            sage: J = JordanSpinEJA(2)
             sage: J.basis()
             Family (e0, e1)
             sage: J.natural_basis()
@@ -160,6 +219,51 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         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``.
@@ -236,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)
@@ -244,6 +348,30 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: J(x).inner_product(J(y))
                 32
 
+            The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
+            multiplication is the usual matrix multiplication in `S^n`,
+            so the inner product of the identity matrix with itself
+            should be the `n`::
+
+                sage: J = RealSymmetricEJA(3)
+                sage: J.one().inner_product(J.one())
+                3
+
+            Likewise, the inner product on `C^n` is `<X,Y> =
+            Re(trace(X*Y))`, where we must necessarily take the real
+            part because the product of Hermitian matrices may not be
+            Hermitian::
+
+                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
+
             TESTS:
 
             Ensure that we can always compute an inner product, and that
@@ -259,7 +387,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             """
             P = self.parent()
             if not other in P:
-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")
 
             return P.inner_product(self, other)
 
@@ -294,7 +422,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             """
             if not other in self.parent():
-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")
 
             A = self.operator_matrix()
             B = other.operator_matrix()
@@ -307,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()
@@ -341,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()
@@ -384,7 +512,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             # TODO: we can do better once the call to is_invertible()
             # doesn't crash on irregular elements.
             #if not self.is_invertible():
-            #    raise ArgumentError('element is not invertible')
+            #    raise ValueError('element is not invertible')
 
             # We do this a little different than the usual recursive
             # call to a finite-dimensional algebra element, because we
@@ -413,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):
@@ -471,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
@@ -496,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()
@@ -508,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
@@ -519,6 +675,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         def minimal_polynomial(self):
             """
+            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.
+
             EXAMPLES::
 
                 sage: set_random_seed()
@@ -540,7 +703,7 @@ 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()
@@ -553,25 +716,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 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()))
-
-            # Recursive call, but should work since elt lives in an
-            # associative algebra.
-            return elt.minimal_polynomial()
+            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()))
+            return elt.operator_matrix().minimal_polynomial()
 
 
         def natural_representation(self):
@@ -586,7 +737,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             EXAMPLES::
 
-                sage: J = ComplexHermitianSimpleEJA(3)
+                sage: J = ComplexHermitianEJA(3)
                 sage: J.one()
                 e0 + e5 + e8
                 sage: J.one().natural_representation()
@@ -597,6 +748,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()
@@ -681,7 +851,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]
@@ -740,7 +910,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             if other is None:
                 other=self
             elif not other in self.parent():
-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")
 
             L = self.operator_matrix()
             M = other.operator_matrix()
@@ -755,7 +925,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()) )
 
 
@@ -824,11 +997,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
@@ -884,7 +1057,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()
@@ -903,21 +1076,25 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             Return the trace inner product of myself and ``other``.
             """
             if not other in self.parent():
-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")
 
             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
@@ -933,22 +1110,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) ]
-
-    # The usual inner product on R^n.
-    ip = lambda x, y: x.vector().inner_product(y.vector())
-
-    return FiniteDimensionalEuclideanJordanAlgebra(field,
-                                                   Qs,
-                                                   rank=dimension,
-                                                   inner_product=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():
@@ -968,6 +1144,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.
 
@@ -978,10 +1160,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)
 
 
@@ -1042,6 +1225,54 @@ 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())
+
+def _vec2mat(v):
+        return matrix(v.base_ring(), sqrt(v.degree()), v.list())
+
 def _multiplication_table_from_matrix_basis(basis):
     """
     At least three of the five simple Euclidean Jordan algebras have the
@@ -1062,19 +1293,13 @@ def _multiplication_table_from_matrix_basis(basis):
     field = basis[0].base_ring()
     dimension = basis[0].nrows()
 
-    def mat2vec(m):
-        return vector(field, m.list())
-
-    def vec2mat(v):
-        return matrix(field, dimension, v.list())
-
     V = VectorSpace(field, dimension**2)
-    W = V.span( mat2vec(s) for s in basis )
+    W = V.span( _mat2vec(s) for s in basis )
 
     # Taking the span above reorders our basis (thanks, jerk!) so we
     # need to put our "matrix basis" in the same order as the
     # (reordered) vector basis.
-    S = tuple( vec2mat(b) for b in W.basis() )
+    S = tuple( _vec2mat(b) for b in W.basis() )
 
     Qs = []
     for s in S:
@@ -1087,7 +1312,7 @@ def _multiplication_table_from_matrix_basis(basis):
         # why we're computing rows here and not columns.
         Q_rows = []
         for t in S:
-            this_row = mat2vec((s*t + t*s)/2)
+            this_row = _mat2vec((s*t + t*s)/2)
             Q_rows.append(W.coordinates(this_row))
         Q = matrix(field, W.dimension(), Q_rows)
         Qs.append(Q)
@@ -1108,24 +1333,38 @@ def _embed_complex_matrix(M):
         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)
-        [ 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()
     if M.ncols() != n:
-        raise ArgumentError("the matrix 'M' must be square")
+        raise ValueError("the matrix 'M' must be square")
     field = M.base_ring()
     blocks = []
     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)
@@ -1142,14 +1381,25 @@ def _unembed_complex_matrix(M):
         ....:                 [ 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:
-        raise ArgumentError("the matrix 'M' must be square")
+        raise ValueError("the matrix 'M' must be square")
     if not n.mod(2).is_zero():
-        raise ArgumentError("the matrix 'M' must be a complex embedding")
+        raise ValueError("the matrix 'M' must be a complex embedding")
 
     F = QuadraticField(-1, 'i')
     i = F.gen()
@@ -1161,16 +1411,137 @@ 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]:
-                raise ArgumentError('bad real submatrix')
+                raise ValueError('bad on-diagonal submatrix')
             if submat[0,1] != -submat[1,0]:
-                raise ArgumentError('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)
 
 
-def RealSymmetricSimpleEJA(n, field=QQ):
+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 inner product used for the real symmetric simple EJA.
+# We keep it as a separate function because e.g. the complex
+# algebra uses the same inner product, except divided by 2.
+def _matrix_ip(X,Y):
+    X_mat = X.natural_representation()
+    Y_mat = Y.natural_representation()
+    return (X_mat*Y_mat).trace()
+
+
+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
@@ -1178,7 +1549,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
@@ -1193,21 +1564,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)
+        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,
@@ -1220,36 +1614,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)
-    return FiniteDimensionalEuclideanJordanAlgebra(field,
-                                                   Qs,
-                                                   rank=n,
-                                                   natural_basis=T)
+    @staticmethod
+    def __classcall_private__(cls, n, field=QQ):
+        S = _complex_hermitian_basis(n)
+        (Qs, T) = _multiplication_table_from_matrix_basis(S)
 
+        fdeja = super(ComplexHermitianEJA, cls)
+        return fdeja.__classcall_private__(cls,
+                                           field,
+                                           Qs,
+                                           rank=n,
+                                           natural_basis=T)
 
-def QuaternionHermitianSimpleEJA(n):
+    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
+
+
+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 =
@@ -1260,7 +1728,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
@@ -1278,26 +1746,26 @@ def JordanSpinSimpleEJA(n, field=QQ):
         0
 
     """
-    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 usual inner product on R^n.
-    ip = lambda x, y: x.vector().inner_product(y.vector())
-
-    # 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=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)