]> 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 3124132b61e5f196268cb027ce1bd26028aaa942..32481621975ab3aabc217eae5ef72d9abf191980 100644 (file)
@@ -21,7 +21,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                               assume_associative=False,
                               category=None,
                               rank=None,
-                              natural_basis=None):
+                              natural_basis=None,
+                              inner_product=None):
         n = len(mult_table)
         mult_table = [b.base_extend(field) for b in mult_table]
         for b in mult_table:
@@ -45,16 +46,19 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                                  names=names,
                                  category=cat,
                                  rank=rank,
-                                 natural_basis=natural_basis)
+                                 natural_basis=natural_basis,
+                                 inner_product=inner_product)
 
 
-    def __init__(self, field,
+    def __init__(self,
+                 field,
                  mult_table,
                  names='e',
                  assume_associative=False,
                  category=None,
                  rank=None,
-                 natural_basis=None):
+                 natural_basis=None,
+                 inner_product=None):
         """
         EXAMPLES:
 
@@ -70,6 +74,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         """
         self._rank = rank
         self._natural_basis = natural_basis
+        self._inner_product = inner_product
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
                      mult_table,
@@ -85,6 +90,34 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         return fmt.format(self.degree(), self.base_ring())
 
 
+    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.
+
+        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 self._inner_product is None:
+            return x.trace_inner_product(y)
+        else:
+            return self._inner_product(x,y)
+
+
     def natural_basis(self):
         """
         Return a more-natural representation of this algebra's basis.
@@ -100,7 +133,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         EXAMPLES::
 
-            sage: J = RealSymmetricSimpleEJA(2)
+            sage: J = RealSymmetricEJA(2)
             sage: J.basis()
             Family (e0, e1, e2)
             sage: J.natural_basis()
@@ -111,7 +144,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         ::
 
-            sage: J = JordanSpinSimpleEJA(2)
+            sage: J = JordanSpinEJA(2)
             sage: J.basis()
             Family (e0, e1)
             sage: J.natural_basis()
@@ -142,6 +175,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``.
@@ -208,6 +286,68 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 raise NotImplementedError('irregular element')
 
 
+        def inner_product(self, other):
+            """
+            Return the parent algebra's inner product of myself and ``other``.
+
+            EXAMPLES:
+
+            The inner product in the Jordan spin algebra is the usual
+            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 = JordanSpinEJA(3)
+                sage: x = vector(QQ,[1,2,3])
+                sage: y = vector(QQ,[4,5,6])
+                sage: x.inner_product(y)
+                32
+                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
+            it gives us back a real number::
+
+                sage: set_random_seed()
+                sage: J = random_eja()
+                sage: x = J.random_element()
+                sage: y = J.random_element()
+                sage: x.inner_product(y) in RR
+                True
+
+            """
+            P = self.parent()
+            if not other in P:
+                raise TypeError("'other' must live in the same algebra")
+
+            return P.inner_product(self, other)
+
+
         def operator_commutes_with(self, other):
             """
             Return whether or not this element operator-commutes
@@ -238,7 +378,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()
@@ -251,12 +391,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()
@@ -285,7 +425,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()
@@ -328,7 +468,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
@@ -415,7 +555,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
@@ -440,7 +580,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()
@@ -452,7 +592,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
@@ -461,6 +601,111 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             return self.span_of_powers().dimension()
 
 
+        def minimal_polynomial(self):
+            """
+            EXAMPLES::
+
+                sage: set_random_seed()
+                sage: x = random_eja().random_element()
+                sage: x.degree() == x.minimal_polynomial().degree()
+                True
+
+            ::
+
+                sage: set_random_seed()
+                sage: x = random_eja().random_element()
+                sage: x.degree() == x.minimal_polynomial().degree()
+                True
+
+            The minimal polynomial and the characteristic polynomial coincide
+            and are known (see Alizadeh, Example 11.11) for all elements of
+            the spin factor algebra that aren't scalar multiples of the
+            identity::
+
+                sage: set_random_seed()
+                sage: n = ZZ.random_element(2,10)
+                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: 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()))
+
+            # Recursive call, but should work since elt lives in an
+            # associative algebra.
+            return elt.minimal_polynomial()
+
+
+        def natural_representation(self):
+            """
+            Return a more-natural representation of this element.
+
+            Every finite-dimensional Euclidean Jordan Algebra is a
+            direct sum of five simple algebras, four of which comprise
+            Hermitian matrices. This method returns the original
+            "natural" representation of this element as a Hermitian
+            matrix, if it has one. If not, you get the usual representation.
+
+            EXAMPLES::
+
+                sage: J = ComplexHermitianEJA(3)
+                sage: J.one()
+                e0 + e5 + e8
+                sage: J.one().natural_representation()
+                [1 0 0 0 0 0]
+                [0 1 0 0 0 0]
+                [0 0 1 0 0 0]
+                [0 0 0 1 0 0]
+                [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()
+            return W.linear_combination(zip(self.vector(), B))
+
+
         def operator_matrix(self):
             """
             Return the matrix that represents left- (or right-)
@@ -528,64 +773,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             return fda_elt.matrix().transpose()
 
 
-
-        def minimal_polynomial(self):
-            """
-            EXAMPLES::
-
-                sage: set_random_seed()
-                sage: x = random_eja().random_element()
-                sage: x.degree() == x.minimal_polynomial().degree()
-                True
-
-            ::
-
-                sage: set_random_seed()
-                sage: x = random_eja().random_element()
-                sage: x.degree() == x.minimal_polynomial().degree()
-                True
-
-            The minimal polynomial and the characteristic polynomial coincide
-            and are known (see Alizadeh, Example 11.11) for all elements of
-            the spin factor algebra that aren't scalar multiples of the
-            identity::
-
-                sage: set_random_seed()
-                sage: n = ZZ.random_element(2,10)
-                sage: J = JordanSpinSimpleEJA(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: 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()))
-
-            # Recursive call, but should work since elt lives in an
-            # associative algebra.
-            return elt.minimal_polynomial()
-
-
         def quadratic_representation(self, other=None):
             """
             Return the quadratic representation of this element.
@@ -597,7 +784,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]
@@ -656,7 +843,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()
@@ -744,7 +931,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 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
@@ -800,7 +987,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()
@@ -819,7 +1006,7 @@ 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()
 
@@ -857,7 +1044,10 @@ def eja_rn(dimension, field=QQ):
     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
            for i in xrange(dimension) ]
 
-    return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
+    return FiniteDimensionalEuclideanJordanAlgebra(field,
+                                                   Qs,
+                                                   rank=dimension,
+                                                   inner_product=_usual_ip)
 
 
 
@@ -878,6 +1068,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.
 
@@ -889,9 +1085,10 @@ def random_eja():
     """
     n = ZZ.random_element(1,5)
     constructor = choice([eja_rn,
-                          JordanSpinSimpleEJA,
-                          RealSymmetricSimpleEJA,
-                          ComplexHermitianSimpleEJA])
+                          JordanSpinEJA,
+                          RealSymmetricEJA,
+                          ComplexHermitianEJA,
+                          QuaternionHermitianEJA])
     return constructor(n, field=QQ)
 
 
@@ -952,6 +1149,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
@@ -972,19 +1217,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:
@@ -997,7 +1236,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)
@@ -1018,24 +1257,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)
@@ -1052,14 +1305,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()
@@ -1071,16 +1335,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
@@ -1088,7 +1473,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
@@ -1103,21 +1488,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,
@@ -1130,36 +1538,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
 
-    """
-    S = _complex_hermitian_basis(n)
-    (Qs, T) = _multiplication_table_from_matrix_basis(S)
-    return FiniteDimensionalEuclideanJordanAlgebra(field,
-                                                   Qs,
-                                                   rank=n,
-                                                   natural_basis=T)
+    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
 
-def QuaternionHermitianSimpleEJA(n):
+    """
+    @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 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.
-    """
-    n = 3
-    pass
+    TESTS:
+
+    The degree of this algebra is `n^2`::
 
-def JordanSpinSimpleEJA(n, field=QQ):
+        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
+
+    """
+    @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 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 =
@@ -1170,7 +1652,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
@@ -1187,28 +1669,34 @@ def JordanSpinSimpleEJA(n, field=QQ):
         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))
+        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)