eja: refactor the class hierarchy to separate the matrix EJAs.

diff --git a/mjo/eja/TODO b/mjo/eja/TODO
index f6d7743782613d241c50e05eee665d59278bd4b4..89b4d329b76d848c0caf1528c1c78b7faec1ff9a 100644 (file)
--- a/mjo/eja/TODO
+++ b/mjo/eja/TODO
@@ -12,10 +12,6 @@
 6. Refactor the current ungodly fast charpoly hack (relies on the
    theory to ensure that the charpolys are equal.)
 
 6. Refactor the current ungodly fast charpoly hack (relies on the
    theory to ensure that the charpolys are equal.)
 

-7. If we factor out a "matrix algebra" class, then it would make sense
-   to replace the custom embedding/unembedding functions with static
-   _real_embedding() and _real_unembedding() methods.
+7. Implement random_instance() for the main EJA class.

 
 

-8. Implement random_instance() for the main EJA class.
-
-9. Implement random_instance() for the subalgebra class.
+8. Implement random_instance() for the subalgebra class.

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index f687e46215a1f1aefb5a5955acffd4fc4c6be301..a207250b4092f97e07d93a62c0ed7e23d9f9536d 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -890,414 +890,123 @@ def random_eja():
 
 
 
 
 
 

-def _real_symmetric_basis(n, field):
-    """
-    Return a basis for the space of real symmetric n-by-n matrices.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import _real_symmetric_basis
-
-    TESTS::
-
-        sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: B = _real_symmetric_basis(n, QQ)
-        sage: all( M.is_symmetric() for M in  B)
-        True
-
-    """
-    # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
-    # coordinates.
-    S = []
-    for i in xrange(n):
-        for j in xrange(i+1):
-            Eij = matrix(field, n, lambda k,l: k==i and l==j)
-            if i == j:
-                Sij = Eij
-            else:
-                Sij = Eij + Eij.transpose()
-            S.append(Sij)
-    return tuple(S)
-
-
-def _complex_hermitian_basis(n, field):
-    """
-    Returns a basis for the space of complex Hermitian n-by-n matrices.
-
-    Why do we embed these? Basically, because all of numerical linear
-    algebra assumes that you're working with vectors consisting of `n`
-    entries from a field and scalars from the same field. There's no way
-    to tell SageMath that (for example) the vectors contain complex
-    numbers, while the scalar field is real.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
-
-    TESTS::
-
-        sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: field = QuadraticField(2, 'sqrt2')
-        sage: B = _complex_hermitian_basis(n, field)
-        sage: all( M.is_symmetric() for M in  B)
-        True
-
-    """
-    R = PolynomialRing(field, 'z')
-    z = R.gen()
-    F = NumberField(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
-    I = F.gen()
-
-    # 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(F, n, lambda k,l: k==i and l==j)
-            if i == j:
-                Sij = _embed_complex_matrix(Eij)
-                S.append(Sij)
-            else:
-                # The second one has a minus because it's conjugated.
-                Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
-                S.append(Sij_real)
-                Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
-                S.append(Sij_imag)

 
 

-    # Since we embedded these, we can drop back to the "field" that we
-    # started with instead of the complex extension "F".
-    return tuple( s.change_ring(field) for s in S )

 
 
 
 

+class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+    @staticmethod
+    def _max_test_case_size():
+        # Play it safe, since this will be squared and the underlying
+        # field can have dimension 4 (quaternions) too.
+        return 3

 
 

-def _quaternion_hermitian_basis(n, field):
-    """
-    Returns a basis for the space of quaternion Hermitian n-by-n matrices.
-
-    Why do we embed these? Basically, because all of numerical linear
-    algebra assumes that you're working with vectors consisting of `n`
-    entries from a field and scalars from the same field. There's no way
-    to tell SageMath that (for example) the vectors contain complex
-    numbers, while the scalar field is real.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
-
-    TESTS::
-
-        sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: B = _quaternion_hermitian_basis(n, QQ)
-        sage: all( M.is_symmetric() for M in B )
-        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 _multiplication_table_from_matrix_basis(basis):
-    """
-    At least three of the five simple Euclidean Jordan algebras have the
-    symmetric multiplication (A,B) |-> (AB + BA)/2, where the
-    multiplication on the right is matrix multiplication. Given a basis
-    for the underlying matrix space, this function returns a
-    multiplication table (obtained by looping through the basis
-    elements) for an algebra of those matrices.
-    """
-    # In S^2, for example, we nominally have four coordinates even
-    # though the space is of dimension three only. The vector space V
-    # is supposed to hold the entire long vector, and the subspace W
-    # of V will be spanned by the vectors that arise from symmetric
-    # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
-    field = basis[0].base_ring()
-    dimension = basis[0].nrows()
-
-    V = VectorSpace(field, dimension**2)
-    W = V.span_of_basis( _mat2vec(s) for s in basis )
-    n = len(basis)
-    mult_table = [[W.zero() for j in range(n)] for i in range(n)]
-    for i in range(n):
-        for j in range(n):
-            mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
-            mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
-
-    return mult_table
-
-
-def _embed_complex_matrix(M):
-    """
-    Embed the n-by-n complex matrix ``M`` into the space of real
-    matrices of size 2n-by-2n via the map the sends each entry `z = a +
-    bi` to the block matrix ``[[a,b],[-b,a]]``.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
-        ....:                                  ComplexHermitianEJA)
-
-    EXAMPLES::
-
-        sage: F = QuadraticField(-1, 'i')
-        sage: x1 = F(4 - 2*i)
-        sage: x2 = F(1 + 2*i)
-        sage: x3 = F(-i)
-        sage: x4 = F(6)
-        sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
-        sage: _embed_complex_matrix(M)
-        [ 4 -2| 1  2]
-        [ 2  4|-2  1]
-        [-----+-----]
-        [ 0 -1| 6  0]
-        [ 1  0| 0  6]
-
-    TESTS:
-
-    Embedding is a homomorphism (isomorphism, in fact)::
-
-        sage: set_random_seed()
-        sage: n_max = ComplexHermitianEJA._max_test_case_size()
-        sage: n = ZZ.random_element(n_max)
-        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 ValueError("the matrix 'M' must be square")
-    field = M.base_ring()
-    blocks = []
-    for z in M.list():
-        a = z.vector()[0] # real part, I guess
-        b = z.vector()[1] # imag part, I guess
-        blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
-
-    # We can drop the imaginaries here.
-    return matrix.block(field.base_ring(), n, blocks)
-
-
-def _unembed_complex_matrix(M):
-    """
-    The inverse of _embed_complex_matrix().
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
-        ....:                                  _unembed_complex_matrix)
-
-    EXAMPLES::
-
-        sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
-        ....:                 [-2,  1,  -4,  3],
-        ....:                 [ 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]
-
-    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
+    @classmethod
+    def _denormalized_basis(cls, n, field):
+        raise NotImplementedError

 
 

-    """
-    n = ZZ(M.nrows())
-    if M.ncols() != n:
-        raise ValueError("the matrix 'M' must be square")
-    if not n.mod(2).is_zero():
-        raise ValueError("the matrix 'M' must be a complex embedding")
-
-    field = M.base_ring() # This should already have sqrt2
-    R = PolynomialRing(field, 'z')
-    z = R.gen()
-    F = NumberField(z**2 + 1,'i', embedding=CLF(-1).sqrt())
-    i = F.gen()
-
-    # Go top-left to bottom-right (reading order), converting every
-    # 2-by-2 block we see to a single complex element.
-    elements = []
-    for k in xrange(n/2):
-        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 on-diagonal submatrix')
-            if submat[0,1] != -submat[1,0]:
-                raise ValueError('bad off-diagonal submatrix')
-            z = submat[0,0] + submat[0,1]*i
-            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.
+    def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
+        S = self._denormalized_basis(n, field)

 
 

-    SETUP::
+        if n > 1 and normalize_basis:
+            # We'll need sqrt(2) to normalize the basis, and this
+            # winds up in the multiplication table, so the whole
+            # algebra needs to be over the field extension.
+            R = PolynomialRing(field, 'z')
+            z = R.gen()
+            p = z**2 - 2
+            if p.is_irreducible():
+                field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
+                S = [ s.change_ring(field) for s in S ]
+            self._basis_normalizers = tuple(
+                ~(self.natural_inner_product(s,s).sqrt()) for s in S )
+            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )

 
 

-        sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
-        ....:                                  QuaternionHermitianEJA)
+        Qs = self.multiplication_table_from_matrix_basis(S)

 
 

-    EXAMPLES::
+        fdeja = super(MatrixEuclideanJordanAlgebra, self)
+        return fdeja.__init__(field,
+                              Qs,
+                              rank=n,
+                              natural_basis=S,
+                              **kwargs)

 
 

-        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)::
+    @staticmethod
+    def multiplication_table_from_matrix_basis(basis):
+        """
+        At least three of the five simple Euclidean Jordan algebras have the
+        symmetric multiplication (A,B) |-> (AB + BA)/2, where the
+        multiplication on the right is matrix multiplication. Given a basis
+        for the underlying matrix space, this function returns a
+        multiplication table (obtained by looping through the basis
+        elements) for an algebra of those matrices.
+        """
+        # In S^2, for example, we nominally have four coordinates even
+        # though the space is of dimension three only. The vector space V
+        # is supposed to hold the entire long vector, and the subspace W
+        # of V will be spanned by the vectors that arise from symmetric
+        # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
+        field = basis[0].base_ring()
+        dimension = basis[0].nrows()
+
+        V = VectorSpace(field, dimension**2)
+        W = V.span_of_basis( _mat2vec(s) for s in basis )
+        n = len(basis)
+        mult_table = [[W.zero() for j in range(n)] for i in range(n)]
+        for i in range(n):
+            for j in range(n):
+                mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
+                mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))

 
 

-        sage: set_random_seed()
-        sage: n_max = QuaternionHermitianEJA._max_test_case_size()
-        sage: n = ZZ.random_element(n_max)
-        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
+        return mult_table

 
 

-    """
-    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 matrix.block(quaternions.base_ring(), n, blocks)
-
-
-def _unembed_quaternion_matrix(M):
-    """
-    The inverse of _embed_quaternion_matrix().

 
 

-    SETUP::
+    @staticmethod
+    def real_embed(M):
+        """
+        Embed the matrix ``M`` into a space of real matrices.

 
 

-        sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
-        ....:                                  _unembed_quaternion_matrix)
+        The matrix ``M`` can have entries in any field at the moment:
+        the real numbers, complex numbers, or quaternions. And although
+        they are not a field, we can probably support octonions at some
+        point, too. This function returns a real matrix that "acts like"
+        the original with respect to matrix multiplication; i.e.

 
 

-    EXAMPLES::
+          real_embed(M*N) = real_embed(M)*real_embed(N)

 
 

-        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]
+        """
+        raise NotImplementedError

 
 

-    TESTS:

 
 

-    Unembedding is the inverse of embedding::
+    @staticmethod
+    def real_unembed(M):
+        """
+        The inverse of :meth:`real_embed`.
+        """
+        raise NotImplementedError

 
 

-        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")
-
-    # Use the base ring of the matrix to ensure that its entries can be
-    # multiplied by elements of the quaternion algebra.
-    field = M.base_ring()
-    Q = QuaternionAlgebra(field,-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].vector()[0]   # real part
-            z += submat[0,0].vector()[1]*i # imag part
-            z += submat[0,1].vector()[0]*j # real part
-            z += submat[0,1].vector()[1]*k # imag part
-            elements.append(z)
-
-    return matrix(Q, n/4, elements)
-
-
-
-class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
+    @classmethod
+    def natural_inner_product(cls,X,Y):
+        Xu = cls.real_unembed(X)
+        Yu = cls.real_unembed(Y)
+        tr = (Xu*Yu).trace()
+        if tr in RLF:
+            # It's real already.
+            return tr
+
+        # Otherwise, try the thing that works for complex numbers; and
+        # if that doesn't work, the thing that works for quaternions.
+        try:
+            return tr.vector()[0] # real part, imag part is index 1
+        except AttributeError:
+            # A quaternions doesn't have a vector() method, but does
+            # have coefficient_tuple() method that returns the
+            # coefficients of 1, i, j, and k -- in that order.
+            return tr.coefficient_tuple()[0]
+
+
+class RealSymmetricEJA(MatrixEuclideanJordanAlgebra):

     """
     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


@@ -1379,38 +1088,188 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
         True
 
     """
         True
 
     """

-    def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
-        S = _real_symmetric_basis(n, field)
+    @classmethod
+    def _denormalized_basis(cls, n, field):
+        """
+        Return a basis for the space of real symmetric n-by-n matrices.

 
 

-        if n > 1 and normalize_basis:
-            # We'll need sqrt(2) to normalize the basis, and this
-            # winds up in the multiplication table, so the whole
-            # algebra needs to be over the field extension.
-            R = PolynomialRing(field, 'z')
-            z = R.gen()
-            p = z**2 - 2
-            if p.is_irreducible():
-                field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
-                S = [ s.change_ring(field) for s in S ]
-            self._basis_normalizers = tuple(
-                ~(self.natural_inner_product(s,s).sqrt()) for s in S )
-            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import RealSymmetricEJA

 
 

-        Qs = _multiplication_table_from_matrix_basis(S)
+        TESTS::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
+            sage: all( M.is_symmetric() for M in  B)
+            True
+
+        """
+        # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
+        # coordinates.
+        S = []
+        for i in xrange(n):
+            for j in xrange(i+1):
+                Eij = matrix(field, n, lambda k,l: k==i and l==j)
+                if i == j:
+                    Sij = Eij
+                else:
+                    Sij = Eij + Eij.transpose()
+                S.append(Sij)
+        return tuple(S)

 
 

-        fdeja = super(RealSymmetricEJA, self)
-        return fdeja.__init__(field,
-                              Qs,
-                              rank=n,
-                              natural_basis=S,
-                              **kwargs)

 
     @staticmethod
     def _max_test_case_size():
 
     @staticmethod
     def _max_test_case_size():

-        return 5
+        return 5 # Dimension 10
+
+    @staticmethod
+    def real_embed(M):
+        """
+        Embed the matrix ``M`` into a space of real matrices.
+
+        The matrix ``M`` can have entries in any field at the moment:
+        the real numbers, complex numbers, or quaternions. And although
+        they are not a field, we can probably support octonions at some
+        point, too. This function returns a real matrix that "acts like"
+        the original with respect to matrix multiplication; i.e.
+
+          real_embed(M*N) = real_embed(M)*real_embed(N)
+
+        """
+        return M
+
+
+    @staticmethod
+    def real_unembed(M):
+        """
+        The inverse of :meth:`real_embed`.
+        """
+        return M
+
+

 
 

+class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
+    @staticmethod
+    def real_embed(M):
+        """
+        Embed the n-by-n complex matrix ``M`` into the space of real
+        matrices of size 2n-by-2n via the map the sends each entry `z = a +
+        bi` to the block matrix ``[[a,b],[-b,a]]``.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import \
+            ....:   ComplexMatrixEuclideanJordanAlgebra
+
+        EXAMPLES::
+
+            sage: F = QuadraticField(-1, 'i')
+            sage: x1 = F(4 - 2*i)
+            sage: x2 = F(1 + 2*i)
+            sage: x3 = F(-i)
+            sage: x4 = F(6)
+            sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
+            sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
+            [ 4 -2| 1  2]
+            [ 2  4|-2  1]
+            [-----+-----]
+            [ 0 -1| 6  0]
+            [ 1  0| 0  6]
+
+        TESTS:

 
 

-class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
+        Embedding is a homomorphism (isomorphism, in fact)::
+
+            sage: set_random_seed()
+            sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
+            sage: n = ZZ.random_element(n_max)
+            sage: F = QuadraticField(-1, 'i')
+            sage: X = random_matrix(F, n)
+            sage: Y = random_matrix(F, n)
+            sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
+            sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
+            sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
+            sage: Xe*Ye == XYe
+            True
+
+        """
+        n = M.nrows()
+        if M.ncols() != n:
+            raise ValueError("the matrix 'M' must be square")
+        field = M.base_ring()
+        blocks = []
+        for z in M.list():
+            a = z.vector()[0] # real part, I guess
+            b = z.vector()[1] # imag part, I guess
+            blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
+
+        # We can drop the imaginaries here.
+        return matrix.block(field.base_ring(), n, blocks)
+
+
+    @staticmethod
+    def real_unembed(M):
+        """
+        The inverse of _embed_complex_matrix().
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import \
+            ....:   ComplexMatrixEuclideanJordanAlgebra
+
+        EXAMPLES::
+
+            sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
+            ....:                 [-2,  1,  -4,  3],
+            ....:                 [ 9,  10, 11, 12],
+            ....:                 [-10, 9, -12, 11] ])
+            sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
+            [  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: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
+            sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
+            True
+
+        """
+        n = ZZ(M.nrows())
+        if M.ncols() != n:
+            raise ValueError("the matrix 'M' must be square")
+        if not n.mod(2).is_zero():
+            raise ValueError("the matrix 'M' must be a complex embedding")
+
+        field = M.base_ring() # This should already have sqrt2
+        R = PolynomialRing(field, 'z')
+        z = R.gen()
+        F = NumberField(z**2 + 1,'i', embedding=CLF(-1).sqrt())
+        i = F.gen()
+
+        # Go top-left to bottom-right (reading order), converting every
+        # 2-by-2 block we see to a single complex element.
+        elements = []
+        for k in xrange(n/2):
+            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 on-diagonal submatrix')
+                if submat[0,1] != -submat[1,0]:
+                    raise ValueError('bad off-diagonal submatrix')
+                z = submat[0,0] + submat[0,1]*i
+                elements.append(z)
+
+        return matrix(F, n/2, elements)
+
+
+class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):

     """
     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,


@@ -1482,47 +1341,195 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         True
 
     """
         True
 
     """

-    def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
-        S = _complex_hermitian_basis(n, field)
+    @classmethod
+    def _denormalized_basis(cls, n, field):
+        """
+        Returns a basis for the space of complex Hermitian n-by-n matrices.

 
 

-        if n > 1 and normalize_basis:
-            # We'll need sqrt(2) to normalize the basis, and this
-            # winds up in the multiplication table, so the whole
-            # algebra needs to be over the field extension.
-            R = PolynomialRing(field, 'z')
-            z = R.gen()
-            p = z**2 - 2
-            if p.is_irreducible():
-                field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
-                S = [ s.change_ring(field) for s in S ]
-            self._basis_normalizers = tuple(
-                ~(self.natural_inner_product(s,s).sqrt()) for s in S )
-            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )
+        Why do we embed these? Basically, because all of numerical linear
+        algebra assumes that you're working with vectors consisting of `n`
+        entries from a field and scalars from the same field. There's no way
+        to tell SageMath that (for example) the vectors contain complex
+        numbers, while the scalar field is real.

 
 

-        Qs = _multiplication_table_from_matrix_basis(S)
+        SETUP::

 
 

-        fdeja = super(ComplexHermitianEJA, self)
-        return fdeja.__init__(field,
-                              Qs,
-                              rank=n,
-                              natural_basis=S,
-                              **kwargs)
+            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: field = QuadraticField(2, 'sqrt2')
+            sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
+            sage: all( M.is_symmetric() for M in  B)
+            True

 
 

+        """
+        R = PolynomialRing(field, 'z')
+        z = R.gen()
+        F = NumberField(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+        I = F.gen()

 
 

+        # 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(F, n, lambda k,l: k==i and l==j)
+                if i == j:
+                    Sij = cls.real_embed(Eij)
+                    S.append(Sij)
+                else:
+                    # The second one has a minus because it's conjugated.
+                    Sij_real = cls.real_embed(Eij + Eij.transpose())
+                    S.append(Sij_real)
+                    Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
+                    S.append(Sij_imag)
+
+        # Since we embedded these, we can drop back to the "field" that we
+        # started with instead of the complex extension "F".
+        return tuple( s.change_ring(field) for s in S )
+
+
+
+class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):

     @staticmethod
     @staticmethod

-    def _max_test_case_size():
-        return 4
+    def real_embed(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.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import \
+            ....:   QuaternionMatrixEuclideanJordanAlgebra
+
+        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: QuaternionMatrixEuclideanJordanAlgebra.real_embed(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_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
+            sage: n = ZZ.random_element(n_max)
+            sage: Q = QuaternionAlgebra(QQ,-1,-1)
+            sage: X = random_matrix(Q, n)
+            sage: Y = random_matrix(Q, n)
+            sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
+            sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
+            sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
+            sage: Xe*Ye == XYe
+            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]
+            cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
+                                 [-c + d*i, a - b*i]])
+            realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM)
+            blocks.append(realM)
+
+        # We should have real entries by now, so use the realest field
+        # we've got for the return value.
+        return matrix.block(quaternions.base_ring(), n, blocks)
+
+

 
     @staticmethod
 
     @staticmethod

-    def natural_inner_product(X,Y):
-        Xu = _unembed_complex_matrix(X)
-        Yu = _unembed_complex_matrix(Y)
-        # The trace need not be real; consider Xu = (i*I) and Yu = I.
-        return ((Xu*Yu).trace()).vector()[0] # real part, I guess
+    def real_unembed(M):
+        """
+        The inverse of _embed_quaternion_matrix().
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import \
+            ....:   QuaternionMatrixEuclideanJordanAlgebra
+
+        EXAMPLES::
+
+            sage: M = matrix(QQ, [[ 1,  2,  3,  4],
+            ....:                 [-2,  1, -4,  3],
+            ....:                 [-3,  4,  1, -2],
+            ....:                 [-4, -3,  2,  1]])
+            sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(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: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
+            sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
+            True

 
 

-class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
+        """
+        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")
+
+        # Use the base ring of the matrix to ensure that its entries can be
+        # multiplied by elements of the quaternion algebra.
+        field = M.base_ring()
+        Q = QuaternionAlgebra(field,-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 = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
+                    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].vector()[0]   # real part
+                z += submat[0,0].vector()[1]*i # imag part
+                z += submat[0,1].vector()[0]*j # real part
+                z += submat[0,1].vector()[1]*k # imag part
+                elements.append(z)
+
+        return matrix(Q, n/4, elements)
+
+
+
+class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):

     """
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
     """
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the


@@ -1594,45 +1601,57 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         True
 
     """
         True
 
     """

-    def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
-        S = _quaternion_hermitian_basis(n, field)
+    @classmethod
+    def _denormalized_basis(cls, n, field):
+        """
+        Returns a basis for the space of quaternion Hermitian n-by-n matrices.

 
 

-        if n > 1 and normalize_basis:
-            # We'll need sqrt(2) to normalize the basis, and this
-            # winds up in the multiplication table, so the whole
-            # algebra needs to be over the field extension.
-            R = PolynomialRing(field, 'z')
-            z = R.gen()
-            p = z**2 - 2
-            if p.is_irreducible():
-                field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
-                S = [ s.change_ring(field) for s in S ]
-            self._basis_normalizers = tuple(
-                ~(self.natural_inner_product(s,s).sqrt()) for s in S )
-            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )
+        Why do we embed these? Basically, because all of numerical
+        linear algebra assumes that you're working with vectors consisting
+        of `n` entries from a field and scalars from the same field. There's
+        no way to tell SageMath that (for example) the vectors contain
+        complex numbers, while the scalar field is real.

 
 

-        Qs = _multiplication_table_from_matrix_basis(S)
+        SETUP::

 
 

-        fdeja = super(QuaternionHermitianEJA, self)
-        return fdeja.__init__(field,
-                              Qs,
-                              rank=n,
-                              natural_basis=S,
-                              **kwargs)
+            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA

 
 

-    @staticmethod
-    def _max_test_case_size():
-        return 3
+        TESTS::

 
 

-    @staticmethod
-    def natural_inner_product(X,Y):
-        Xu = _unembed_quaternion_matrix(X)
-        Yu = _unembed_quaternion_matrix(Y)
-        # The trace need not be real; consider Xu = (i*I) and Yu = I.
-        # The result will be a quaternion algebra element, which doesn't
-        # have a "vector" method, but does have coefficient_tuple() method
-        # that returns the coefficients of 1, i, j, and k -- in that order.
-        return ((Xu*Yu).trace()).coefficient_tuple()[0]
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
+            sage: all( M.is_symmetric() for M in B )
+            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 = cls.real_embed(Eij)
+                    S.append(Sij)
+                else:
+                    # The second, third, and fourth ones have a minus
+                    # because they're conjugated.
+                    Sij_real = cls.real_embed(Eij + Eij.transpose())
+                    S.append(Sij_real)
+                    Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
+                    S.append(Sij_I)
+                    Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
+                    S.append(Sij_J)
+                    Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
+                    S.append(Sij_K)
+        return tuple(S)