]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: cache the charpoly coefficients for the AlbertEJA.
[sage.d.git] / mjo / eja / eja_algebra.py
index 3027075374a177ba5f72bbd690a2c07dad02d136..587d8e339463adaafb0a390ec164e76a6e5ca0c4 100644 (file)
@@ -104,6 +104,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         product. This will be applied to ``basis`` to compute an
         inner-product table (basically a matrix) for this algebra.
 
+      - ``matrix_space`` -- the space that your matrix basis lives in,
+        or ``None`` (the default). So long as your basis does not have
+        length zero you can omit this. But in trivial algebras, it is
+        required.
+
       - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
         field for the algebra.
 
@@ -128,7 +133,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         sage: basis = tuple(b.superalgebra_element() for b in A.basis())
         sage: J.subalgebra(basis, orthonormalize=False).is_associative()
         True
-
     """
     Element = FiniteDimensionalEJAElement
 
@@ -137,6 +141,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  jordan_product,
                  inner_product,
                  field=AA,
+                 matrix_space=None,
                  orthonormalize=True,
                  associative=None,
                  cartesian_product=False,
@@ -172,6 +177,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         category = MagmaticAlgebras(field).FiniteDimensional()
         category = category.WithBasis().Unital().Commutative()
 
+        if n <= 1:
+            # All zero- and one-dimensional algebras are just the real
+            # numbers with (some positive multiples of) the usual
+            # multiplication as its Jordan and inner-product.
+            associative = True
         if associative is None:
             # We should figure it out. As with check_axioms, we have to do
             # this without the help of the _jordan_product_is_associative()
@@ -231,8 +241,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             basis = tuple(gram_schmidt(basis, inner_product))
 
         # Save the (possibly orthonormalized) matrix basis for
-        # later...
+        # later, as well as the space that its elements live in.
+        # In most cases we can deduce the matrix space, but when
+        # n == 0 (that is, there are no basis elements) we cannot.
         self._matrix_basis = basis
+        if matrix_space is None:
+            self._matrix_space = self._matrix_basis[0].parent()
+        else:
+            self._matrix_space = matrix_space
 
         # Now create the vector space for the algebra, which will have
         # its own set of non-ambient coordinates (in terms of the
@@ -1016,7 +1032,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
             sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
             sage: J.matrix_space()
-            Full MatrixSpace of 4 by 4 dense matrices over Rational Field
+            Module of 2 by 2 matrices with entries in Algebraic Field over
+            the scalar ring Rational Field
             sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
             sage: J.matrix_space()
             Module of 1 by 1 matrices with entries in Quaternion
@@ -1024,10 +1041,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             the scalar ring Rational Field
 
         """
-        if self.is_trivial():
-            return MatrixSpace(self.base_ring(), 0)
-        else:
-            return self.matrix_basis()[0].parent()
+        return self._matrix_space
 
 
     @cached_method
@@ -1604,6 +1618,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        jordan_product,
                                        inner_product,
                                        field=QQ,
+                                       matrix_space=self.matrix_space(),
                                        associative=self.is_associative(),
                                        orthonormalize=False,
                                        check_field=False,
@@ -1731,7 +1746,87 @@ class ConcreteEJA(FiniteDimensionalEJA):
         return eja_class.random_instance(*args, **kwargs)
 
 
-class MatrixEJA:
+class MatrixEJA(FiniteDimensionalEJA):
+    @staticmethod
+    def _denormalized_basis(A):
+        """
+        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.hurwitz import (ComplexMatrixAlgebra,
+            ....:                          QuaternionMatrixAlgebra,
+            ....:                          OctonionMatrixAlgebra)
+            sage: from mjo.eja.eja_algebra import MatrixEJA
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = MatrixSpace(QQ, n)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in  B)
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in  B)
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in B )
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in B )
+            True
+
+        """
+        # These work for real MatrixSpace, whose monomials only have
+        # two coordinates (because the last one would always be "1").
+        es = A.base_ring().gens()
+        gen = lambda A,m: A.monomial(m[:2])
+
+        if hasattr(A, 'entry_algebra_gens'):
+            # We've got a MatrixAlgebra, and its monomials will have
+            # three coordinates.
+            es = A.entry_algebra_gens()
+            gen = lambda A,m: A.monomial(m)
+
+        basis = []
+        for i in range(A.nrows()):
+            for j in range(i+1):
+                if i == j:
+                    E_ii = gen(A, (i,j,es[0]))
+                    basis.append(E_ii)
+                else:
+                    for e in es:
+                        E_ij  = gen(A, (i,j,e))
+                        E_ij += E_ij.conjugate_transpose()
+                        basis.append(E_ij)
+
+        return tuple( basis )
+
     @staticmethod
     def jordan_product(X,Y):
         return (X*Y + Y*X)/2
@@ -1741,11 +1836,74 @@ class MatrixEJA:
         r"""
         A trace inner-product for matrices that aren't embedded in the
         reals. It takes MATRICES as arguments, not EJA elements.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  OctonionHermitianEJA)
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        ::
+
+            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        ::
+
+            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
         """
-        return (X*Y).trace().real()
+        tr = (X*Y).trace()
+        if hasattr(tr, 'coefficient'):
+            # Works for octonions, and has to come first because they
+            # also have a "real()" method that doesn't return an
+            # element of the scalar ring.
+            return tr.coefficient(0)
+        elif hasattr(tr, 'coefficient_tuple'):
+            # Works for quaternions.
+            return tr.coefficient_tuple()[0]
 
+        # Works for real and complex numbers.
+        return tr.real()
 
-class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
+
+    def __init__(self, matrix_space, **kwargs):
+        # We know this is a valid EJA, but will double-check
+        # if the user passes check_axioms=True.
+        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+
+        super().__init__(self._denormalized_basis(matrix_space),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         field=matrix_space.base_ring(),
+                         matrix_space=matrix_space,
+                         **kwargs)
+
+        self.rank.set_cache(matrix_space.nrows())
+        self.one.set_cache( self(matrix_space.one()) )
+
+class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -1810,38 +1968,6 @@ class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
-    @classmethod
-    def _denormalized_basis(cls, n, field):
-        """
-        Return a basis for the space of real symmetric n-by-n matrices.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import RealSymmetricEJA
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
-            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 range(n):
-            for j in range(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)
-
-
     @staticmethod
     def _max_random_instance_size():
         return 4 # Dimension 10
@@ -1859,27 +1985,12 @@ class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        associative = False
-        if n <= 1:
-            associative = True
-
-        super().__init__(self._denormalized_basis(n,field),
-                         self.jordan_product,
-                         self.trace_inner_product,
-                         field=field,
-                         associative=associative,
-                         **kwargs)
-
-        # TODO: this could be factored out somehow, but is left here
-        # because the MatrixEJA is not presently a subclass of the
-        # FDEJA class that defines rank() and one().
-        self.rank.set_cache(n)
-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
+        A = MatrixSpace(field, n)
+        super().__init__(A, **kwargs)
 
 
 
-class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
+class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -1892,13 +2003,28 @@ class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
 
     EXAMPLES:
 
-    In theory, our "field" can be any subfield of the reals::
+    In theory, our "field" can be any subfield of the reals, but we
+    can't use inexact real fields at the moment because SageMath
+    doesn't know how to convert their elements into complex numbers,
+    or even into algebraic reals::
 
-        sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
-        Euclidean Jordan algebra of dimension 4 over Real Double Field
-        sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
-        Euclidean Jordan algebra of dimension 4 over Real Field with
-        53 bits of precision
+        sage: QQbar(RDF(1))
+        Traceback (most recent call last):
+        ...
+        TypeError: Illegal initializer for algebraic number
+        sage: AA(RR(1))
+        Traceback (most recent call last):
+        ...
+        TypeError: Illegal initializer for algebraic number
+
+    This causes the following error when we try to scale a matrix of
+    complex numbers by an inexact real number::
+
+        sage: ComplexHermitianEJA(2,field=RR)
+        Traceback (most recent call last):
+        ...
+        TypeError: Unable to coerce entries (=(1.00000000000000,
+        -0.000000000000000)) to coefficients in Algebraic Real Field
 
     TESTS:
 
@@ -1934,95 +2060,16 @@ class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
 
         sage: ComplexHermitianEJA(0)
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
-
     """
-
-    @classmethod
-    def _denormalized_basis(cls, 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 ComplexHermitianEJA
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = ComplexHermitianEJA._denormalized_basis(n,QQ)
-            sage: all( M.is_hermitian() for M in  B)
-            True
-
-        """
-        from mjo.hurwitz import ComplexMatrixAlgebra
-        A = ComplexMatrixAlgebra(n, scalars=field)
-        es = A.entry_algebra_gens()
-
-        basis = []
-        for i in range(n):
-            for j in range(i+1):
-                if i == j:
-                    E_ii = A.monomial( (i,j,es[0]) )
-                    basis.append(E_ii)
-                else:
-                    for e in es:
-                        E_ij  = A.monomial( (i,j,e)             )
-                        ec = e.conjugate()
-                        # If the conjugate has a negative sign in front
-                        # of it, (j,i,ec) won't be a monomial!
-                        if (j,i,ec) in A.indices():
-                            E_ij += A.monomial( (j,i,ec) )
-                        else:
-                            E_ij -= A.monomial( (j,i,-ec) )
-                        basis.append(E_ij)
-
-        return tuple( basis )
-
-    @staticmethod
-    def trace_inner_product(X,Y):
-        r"""
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
-
-        TESTS::
-
-            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
-            sage: I = J.one().to_matrix()
-            sage: J.trace_inner_product(I, -I)
-            -2
-
-        """
-        return (X*Y).trace().real()
-
     def __init__(self, n, field=AA, **kwargs):
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        associative = False
-        if n <= 1:
-            associative = True
+        from mjo.hurwitz import ComplexMatrixAlgebra
+        A = ComplexMatrixAlgebra(n, scalars=field)
+        super().__init__(A, **kwargs)
 
-        super().__init__(self._denormalized_basis(n,field),
-                         self.jordan_product,
-                         self.trace_inner_product,
-                         field=field,
-                         associative=associative,
-                         **kwargs)
-        # TODO: this could be factored out somehow, but is left here
-        # because the MatrixEJA is not presently a subclass of the
-        # FDEJA class that defines rank() and one().
-        self.rank.set_cache(n)
-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
 
     @staticmethod
     def _max_random_instance_size():
@@ -2037,7 +2084,7 @@ class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         return cls(n, **kwargs)
 
 
-class QuaternionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
+class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
     r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
@@ -2094,97 +2141,14 @@ class QuaternionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
-    @classmethod
-    def _denormalized_basis(cls, 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 QuaternionHermitianEJA
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
-            sage: all( M.is_hermitian() for M in B )
-            True
-
-        """
-        from mjo.hurwitz import QuaternionMatrixAlgebra
-        A = QuaternionMatrixAlgebra(n, scalars=field)
-        es = A.entry_algebra_gens()
-
-        basis = []
-        for i in range(n):
-            for j in range(i+1):
-                if i == j:
-                    E_ii = A.monomial( (i,j,es[0]) )
-                    basis.append(E_ii)
-                else:
-                    for e in es:
-                        E_ij  = A.monomial( (i,j,e)             )
-                        ec = e.conjugate()
-                        # If the conjugate has a negative sign in front
-                        # of it, (j,i,ec) won't be a monomial!
-                        if (j,i,ec) in A.indices():
-                            E_ij += A.monomial( (j,i,ec) )
-                        else:
-                            E_ij -= A.monomial( (j,i,-ec) )
-                        basis.append(E_ij)
-
-        return tuple( basis )
-
-
-    @staticmethod
-    def trace_inner_product(X,Y):
-        r"""
-        Overload the superclass method because the quaternions are weird
-        and we need to use ``coefficient_tuple()`` to get the realpart.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
-
-        TESTS::
-
-            sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
-            sage: I = J.one().to_matrix()
-            sage: J.trace_inner_product(I, -I)
-            -2
-
-        """
-        return (X*Y).trace().coefficient_tuple()[0]
-
     def __init__(self, n, field=AA, **kwargs):
         # We know this is a valid EJA, but will double-check
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        associative = False
-        if n <= 1:
-            associative = True
-
-        super().__init__(self._denormalized_basis(n,field),
-                         self.jordan_product,
-                         self.trace_inner_product,
-                         field=field,
-                         associative=associative,
-                         **kwargs)
-
-        # TODO: this could be factored out somehow, but is left here
-        # because the MatrixEJA is not presently a subclass of the
-        # FDEJA class that defines rank() and one().
-        self.rank.set_cache(n)
-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
+        from mjo.hurwitz import QuaternionMatrixAlgebra
+        A = QuaternionMatrixAlgebra(n, scalars=field)
+        super().__init__(A, **kwargs)
 
 
     @staticmethod
@@ -2202,12 +2166,13 @@ class QuaternionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, **kwargs)
 
-class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
+class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
     r"""
     SETUP::
 
         sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
         ....:                                  OctonionHermitianEJA)
+        sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
 
     EXAMPLES:
 
@@ -2222,7 +2187,8 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
     After a change-of-basis, the 2-by-2 algebra has the same
     multiplication table as the ten-dimensional Jordan spin algebra::
 
-        sage: b = OctonionHermitianEJA._denormalized_basis(2,QQ)
+        sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
+        sage: b = OctonionHermitianEJA._denormalized_basis(A)
         sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
         sage: jp = OctonionHermitianEJA.jordan_product
         sage: ip = OctonionHermitianEJA.trace_inner_product
@@ -2311,82 +2277,17 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        super().__init__(self._denormalized_basis(n,field),
-                         self.jordan_product,
-                         self.trace_inner_product,
-                         field=field,
-                         **kwargs)
-
-        # TODO: this could be factored out somehow, but is left here
-        # because the MatrixEJA is not presently a subclass of the
-        # FDEJA class that defines rank() and one().
-        self.rank.set_cache(n)
-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
-
-
-    @classmethod
-    def _denormalized_basis(cls, n, field):
-        """
-        Returns a basis for the space of octonion Hermitian n-by-n
-        matrices.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
-
-        EXAMPLES::
-
-            sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
-            sage: all( M.is_hermitian() for M in B )
-            True
-            sage: len(B)
-            27
-
-        """
         from mjo.hurwitz import OctonionMatrixAlgebra
         A = OctonionMatrixAlgebra(n, scalars=field)
-        es = A.entry_algebra_gens()
-
-        basis = []
-        for i in range(n):
-            for j in range(i+1):
-                if i == j:
-                    E_ii = A.monomial( (i,j,es[0]) )
-                    basis.append(E_ii)
-                else:
-                    for e in es:
-                        E_ij  = A.monomial( (i,j,e)             )
-                        ec = e.conjugate()
-                        # If the conjugate has a negative sign in front
-                        # of it, (j,i,ec) won't be a monomial!
-                        if (j,i,ec) in A.indices():
-                            E_ij += A.monomial( (j,i,ec) )
-                        else:
-                            E_ij -= A.monomial( (j,i,-ec) )
-                        basis.append(E_ij)
-
-        return tuple( basis )
+        super().__init__(A, **kwargs)
 
-    @staticmethod
-    def trace_inner_product(X,Y):
-        r"""
-        The octonions don't know that the reals are embedded in them,
-        so we have to take the e0 component ourselves.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
-
-        TESTS::
-
-            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
-            sage: I = J.one().to_matrix()
-            sage: J.trace_inner_product(I, -I)
-            -2
-
-        """
-        return (X*Y).trace().coefficient(0)
+        if n == 3:
+            from mjo.eja.eja_cache import albert_eja_coeffs
+            a = albert_eja_coeffs(self.coordinate_polynomial_ring())
+            if self._rational_algebra is None:
+                self._charpoly_coefficients.set_cache(a)
+            else:
+                self._rational_algebra._charpoly_coefficients.set_cache(a)
 
 
 class AlbertEJA(OctonionHermitianEJA):
@@ -2451,13 +2352,14 @@ class HadamardEJA(RationalBasisEJA, ConcreteEJA):
         (r0, r1, r2)
     """
     def __init__(self, n, field=AA, **kwargs):
+        MS = MatrixSpace(field, n, 1)
+
         if n == 0:
             jordan_product = lambda x,y: x
             inner_product = lambda x,y: x
         else:
             def jordan_product(x,y):
-                P = x.parent()
-                return P( xi*yi for (xi,yi) in zip(x,y) )
+                return MS( xi*yi for (xi,yi) in zip(x,y) )
 
             def inner_product(x,y):
                 return (x.T*y)[0,0]
@@ -2471,20 +2373,17 @@ class HadamardEJA(RationalBasisEJA, ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        column_basis = tuple( b.column()
-                              for b in FreeModule(field, n).basis() )
+        column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )
         super().__init__(column_basis,
                          jordan_product,
                          inner_product,
                          field=field,
+                         matrix_space=MS,
                          associative=True,
                          **kwargs)
         self.rank.set_cache(n)
 
-        if n == 0:
-            self.one.set_cache( self.zero() )
-        else:
-            self.one.set_cache( sum(self.gens()) )
+        self.one.set_cache( self.sum(self.gens()) )
 
     @staticmethod
     def _max_random_instance_size():
@@ -2598,22 +2497,22 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
         # verify things, we'll skip the rest of the checks.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
+        n = B.nrows()
+        MS = MatrixSpace(field, n, 1)
+
         def inner_product(x,y):
             return (y.T*B*x)[0,0]
 
         def jordan_product(x,y):
-            P = x.parent()
             x0 = x[0,0]
             xbar = x[1:,0]
             y0 = y[0,0]
             ybar = y[1:,0]
             z0 = inner_product(y,x)
             zbar = y0*xbar + x0*ybar
-            return P([z0] + zbar.list())
+            return MS([z0] + zbar.list())
 
-        n = B.nrows()
-        column_basis = tuple( b.column()
-                              for b in FreeModule(field, n).basis() )
+        column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )
 
         # TODO: I haven't actually checked this, but it seems legit.
         associative = False
@@ -2624,6 +2523,7 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
                          jordan_product,
                          inner_product,
                          field=field,
+                         matrix_space=MS,
                          associative=associative,
                          **kwargs)
 
@@ -2631,7 +2531,6 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
         # one-dimensional ambient space (because the rank is bounded
         # by the ambient dimension).
         self.rank.set_cache(min(n,2))
-
         if n == 0:
             self.one.set_cache( self.zero() )
         else:
@@ -2779,10 +2678,11 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA):
         0
 
     """
-    def __init__(self, **kwargs):
+    def __init__(self, field=AA, **kwargs):
         jordan_product = lambda x,y: x
-        inner_product = lambda x,y: 0
+        inner_product = lambda x,y: field.zero()
         basis = ()
+        MS = MatrixSpace(field,0)
 
         # New defaults for keyword arguments
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
@@ -2792,6 +2692,8 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA):
                          jordan_product,
                          inner_product,
                          associative=True,
+                         field=field,
+                         matrix_space=MS,
                          **kwargs)
 
         # The rank is zero using my definition, namely the dimension of the
@@ -2969,7 +2871,14 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
         associative = all( f.is_associative() for f in factors )
 
-        MS = self.matrix_space()
+        # Compute my matrix space. This category isn't perfect, but
+        # is good enough for what we need to do.
+        MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
+        MS_cat = MS_cat.Unital().CartesianProducts()
+        MS_factors = tuple( J.matrix_space() for J in factors )
+        from sage.sets.cartesian_product import CartesianProduct
+        MS = CartesianProduct(MS_factors, MS_cat)
+
         basis = []
         zero = MS.zero()
         for i in range(m):
@@ -3004,15 +2913,16 @@ class CartesianProductEJA(FiniteDimensionalEJA):
                                       jordan_product,
                                       inner_product,
                                       field=field,
+                                      matrix_space=MS,
                                       orthonormalize=False,
                                       associative=associative,
                                       cartesian_product=True,
                                       check_field=False,
                                       check_axioms=False)
 
+        self.rank.set_cache(sum(J.rank() for J in factors))
         ones = tuple(J.one().to_matrix() for J in factors)
         self.one.set_cache(self(ones))
-        self.rank.set_cache(sum(J.rank() for J in factors))
 
     def cartesian_factors(self):
         # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
@@ -3064,11 +2974,13 @@ class CartesianProductEJA(FiniteDimensionalEJA):
             sage: J2 = ComplexHermitianEJA(1)
             sage: J = cartesian_product([J1,J2])
             sage: J.one().to_matrix()[0]
-            [1 0]
-            [0 1]
+            +---+
+            | 1 |
+            +---+
             sage: J.one().to_matrix()[1]
-            [1 0]
-            [0 1]
+            +---+
+            | 1 |
+            +---+
 
         ::
 
@@ -3085,16 +2997,7 @@ class CartesianProductEJA(FiniteDimensionalEJA):
             +----+
 
         """
-        scalars = self.cartesian_factor(0).base_ring()
-
-        # This category isn't perfect, but is good enough for what we
-        # need to do.
-        cat = MagmaticAlgebras(scalars).FiniteDimensional().WithBasis()
-        cat = cat.Unital().CartesianProducts()
-        factors = tuple( J.matrix_space() for J in self.cartesian_factors() )
-
-        from sage.sets.cartesian_product import CartesianProduct
-        return CartesianProduct(factors, cat)
+        return super().matrix_space()
 
 
     @cached_method