]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/euclidean_jordan_algebra.py
eja: fix an erroneous test case.
[sage.d.git] / mjo / eja / euclidean_jordan_algebra.py
index 74c0bf1cb11e38460f39aef4920fa70f60d02409..30d04f4a352b7d553470bd456d18d224e2ec5428 100644 (file)
@@ -5,12 +5,405 @@ are used in optimization, and have some additional nice methods beyond
 what can be supported in a general Jordan Algebra.
 """
 
-from sage.categories.magmatic_algebras import MagmaticAlgebras
+from sage.categories.finite_dimensional_algebras_with_basis import FiniteDimensionalAlgebrasWithBasis
+from sage.categories.morphism import SetMorphism
 from sage.structure.element import is_Matrix
 from sage.structure.category_object import normalize_names
 
 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra
 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement
+from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_morphism import FiniteDimensionalAlgebraMorphism, FiniteDimensionalAlgebraHomset
+
+
+class FiniteDimensionalEuclideanJordanAlgebraHomset(FiniteDimensionalAlgebraHomset):
+
+    def has_coerce_map_from(self, S):
+        """
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2)
+            sage: H = J.Hom(J)
+            sage: H.has_coerce_map_from(QQ)
+            True
+
+        """
+        try:
+            # The Homset classes override has_coerce_map_from() with
+            # something that crashes when it's given e.g. QQ.
+            if S.is_subring(self.codomain().base_ring()):
+                return True
+        except:
+            pclass = super(FiniteDimensionalEuclideanJordanAlgebraHomset, self)
+            return pclass.has_coerce_map_from(S)
+
+
+    def _coerce_map_from_(self, S):
+        """
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2)
+            sage: H = J.Hom(J)
+            sage: H.coerce(2)
+            Morphism from Euclidean Jordan algebra of degree 3 over Rational
+            Field to Euclidean Jordan algebra of degree 3 over Rational Field
+            given by matrix
+            [2 0 0]
+            [0 2 0]
+            [0 0 2]
+
+        """
+        C = self.codomain()
+        R = C.base_ring()
+        if S.is_subring(R):
+            h = S.hom(self.codomain())
+            return SetMorphism(Hom(S,C), lambda x: h(x).operator())
+
+
+    def __call__(self, x):
+        """
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2)
+            sage: H = J.Hom(J)
+            sage: H(2)
+            Morphism from Euclidean Jordan algebra of degree 3 over Rational
+            Field to Euclidean Jordan algebra of degree 3 over Rational Field
+            given by matrix
+            [2 0 0]
+            [0 2 0]
+            [0 0 2]
+
+        """
+        if x in self.base_ring():
+            cols = self.domain().dimension()
+            rows = self.codomain().dimension()
+            x = x*identity_matrix(self.codomain().base_ring(), rows, cols)
+        return FiniteDimensionalEuclideanJordanAlgebraMorphism(self, x)
+
+
+    def one(self):
+        """
+        Return the identity morphism, but as a member of the right
+        space (so that we can add it, multiply it, etc.)
+        """
+        cols = self.domain().dimension()
+        rows = self.codomain().dimension()
+        mat = identity_matrix(self.base_ring(), rows, cols)
+        return FiniteDimensionalEuclideanJordanAlgebraMorphism(self, mat)
+
+
+
+class FiniteDimensionalEuclideanJordanAlgebraMorphism(FiniteDimensionalAlgebraMorphism):
+    """
+    A linear map between two finite-dimensional EJAs.
+
+    This is a very thin wrapper around FiniteDimensionalAlgebraMorphism
+    that does only a few things:
+
+      1. Avoids the ``unitary`` and ``check`` arguments to the constructor
+         that will always be ``False``. This is necessary because these
+         are homomorphisms with respect to ADDITION, but the SageMath
+         machinery wants to check that they're homomorphisms with respect
+         to (Jordan) MULTIPLICATION. That obviously doesn't work.
+
+      2. Inputs and outputs the underlying matrix with respect to COLUMN
+         vectors, unlike the parent class.
+
+      3. Allows us to add, subtract, negate, multiply (compose), and
+         invert morphisms in the obvious way.
+
+    If this seems a bit heavyweight, it is. I would have been happy to
+    use a the ring morphism that underlies the finite-dimensional
+    algebra morphism, but they don't seem to be callable on elements of
+    our EJA, and you can't add/multiply/etc. them.
+    """
+    def _add_(self, other):
+        """
+        Add two EJA morphisms in the obvious way.
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: x = J.zero()
+            sage: y = J.one()
+            sage: x.operator() + y.operator()
+            Morphism from Euclidean Jordan algebra of degree 6 over Rational
+            Field to Euclidean Jordan algebra of degree 6 over Rational Field
+            given by matrix
+            [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]
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: (x.operator() + y.operator()) in J.Hom(J)
+            True
+
+        """
+        P = self.parent()
+        if not other in P:
+            raise ValueError("summands must live in the same space")
+
+        return FiniteDimensionalEuclideanJordanAlgebraMorphism(
+                  P,
+                  self.matrix() + other.matrix() )
+
+
+    def __init__(self, parent, f):
+        FiniteDimensionalAlgebraMorphism.__init__(self,
+                                                  parent,
+                                                  f.transpose(),
+                                                  unitary=False,
+                                                  check=False)
+
+
+    def __invert__(self):
+        """
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2)
+            sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
+            sage: x.is_invertible()
+            True
+            sage: ~x.operator()
+            Morphism from Euclidean Jordan algebra of degree 3 over Rational
+            Field to Euclidean Jordan algebra of degree 3 over Rational Field
+            given by matrix
+            [-3/2    2 -1/2]
+            [   1    0    0]
+            [-1/2    0  1/2]
+            sage: x.operator_matrix().inverse()
+            [-3/2    2 -1/2]
+            [   1    0    0]
+            [-1/2    0  1/2]
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: not x.is_invertible() or (
+            ....:   (~x.operator()).matrix() == x.operator_matrix().inverse() )
+            True
+
+        """
+        A = self.matrix()
+        if not A.is_invertible():
+            raise ValueError("morphism is not invertible")
+
+        P = self.parent()
+        return FiniteDimensionalEuclideanJordanAlgebraMorphism(self.parent(),
+                                                                A.inverse())
+
+    def _lmul_(self, right):
+        """
+        Compose two EJA morphisms using multiplicative notation.
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2)
+            sage: x = J.zero()
+            sage: y = J.one()
+            sage: x.operator() * y.operator()
+            Morphism from Euclidean Jordan algebra of degree 3 over Rational
+            Field to Euclidean Jordan algebra of degree 3 over Rational Field
+            given by matrix
+            [0 0 0]
+            [0 0 0]
+            [0 0 0]
+
+        ::
+
+            sage: J = RealSymmetricEJA(2)
+            sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
+            sage: x.operator()
+            Morphism from Euclidean Jordan algebra of degree 3 over Rational
+            Field to Euclidean Jordan algebra of degree 3 over Rational Field
+            given by matrix
+            [  0   1   0]
+            [1/2   1 1/2]
+            [  0   1   2]
+            sage: 2*x.operator()
+            Morphism from Euclidean Jordan algebra of degree 3 over Rational
+            Field to Euclidean Jordan algebra of degree 3 over Rational Field
+            given by matrix
+            [0 2 0]
+            [1 2 1]
+            [0 2 4]
+            sage: x.operator()*2
+            Morphism from Euclidean Jordan algebra of degree 3 over Rational
+            Field to Euclidean Jordan algebra of degree 3 over Rational Field
+            given by matrix
+            [0 2 0]
+            [1 2 1]
+            [0 2 4]
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: (x.operator() * y.operator()) in J.Hom(J)
+            True
+
+        """
+        try:
+            # I think the morphism classes break the coercion framework
+            # somewhere along the way, so we have to do this ourselves.
+            right = self.parent().coerce(right)
+        except:
+            pass
+
+        if not right.codomain() is self.domain():
+            raise ValueError("(co)domains must agree for composition")
+
+        return FiniteDimensionalEuclideanJordanAlgebraMorphism(
+                 self.parent(),
+                 self.matrix()*right.matrix() )
+
+    __mul__ = _lmul_
+
+
+    def __pow__(self, n):
+        """
+
+        TESTS::
+
+            sage: J = JordanSpinEJA(4)
+            sage: e0,e1,e2,e3 = J.gens()
+            sage: x = -5/2*e0 + 1/2*e2 + 20*e3
+            sage: Qx = x.quadratic_representation()
+            sage: Qx^0
+            Morphism from Euclidean Jordan algebra of degree 4 over Rational
+            Field to Euclidean Jordan algebra of degree 4 over Rational Field
+            given by matrix
+            [1 0 0 0]
+            [0 1 0 0]
+            [0 0 1 0]
+            [0 0 0 1]
+            sage: (x^0).quadratic_representation() == Qx^0
+            True
+
+        """
+        if n == 0:
+            # We get back the stupid identity morphism which doesn't
+            # live in the right space.
+            return self.parent().one()
+        elif n == 1:
+            return self
+        else:
+            return FiniteDimensionalAlgebraMorphism.__pow__(self,n)
+
+
+    def _neg_(self):
+        """
+        Negate this morphism.
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2)
+            sage: x = J.one()
+            sage: -x.operator()
+            Morphism from Euclidean Jordan algebra of degree 3 over Rational
+            Field to Euclidean Jordan algebra of degree 3 over Rational Field
+            given by matrix
+            [-1  0  0]
+            [ 0 -1  0]
+            [ 0  0 -1]
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: -x.operator() in J.Hom(J)
+            True
+
+        """
+        return FiniteDimensionalEuclideanJordanAlgebraMorphism(
+                  self.parent(),
+                  -self.matrix() )
+
+
+    def _repr_(self):
+        """
+        We override only the representation that is shown to the user,
+        because we want the matrix to be with respect to COLUMN vectors.
+
+        TESTS:
+
+        Ensure that we see the transpose of the underlying matrix object:
+
+            sage: J = RealSymmetricEJA(3)
+            sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
+            sage: L = x.operator()
+            sage: L
+            Morphism from Euclidean Jordan algebra of degree 6 over Rational
+            Field to Euclidean Jordan algebra of degree 6 over Rational Field
+            given by matrix
+            [  0   1   2   0   0   0]
+            [1/2 3/2   2 1/2   1   0]
+            [  1   2 5/2   0 1/2   1]
+            [  0   1   0   3   4   0]
+            [  0   1 1/2   2   4   2]
+            [  0   0   2   0   4   5]
+            sage: L._matrix
+            [  0 1/2   1   0   0   0]
+            [  1 3/2   2   1   1   0]
+            [  2   2 5/2   0 1/2   2]
+            [  0 1/2   0   3   2   0]
+            [  0   1 1/2   4   4   4]
+            [  0   0   1   0   2   5]
+
+        """
+        return "Morphism from {} to {} given by matrix\n{}".format(
+            self.domain(), self.codomain(), self.matrix())
+
+
+    def __sub__(self, other):
+        """
+        Subtract one morphism from another using addition and negation.
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2)
+            sage: L1 = J.one().operator()
+            sage: L1 - L1
+            Morphism from Euclidean Jordan algebra of degree 3 over Rational
+            Field to Euclidean Jordan algebra of degree 3 over Rational
+            Field given by matrix
+            [0 0 0]
+            [0 0 0]
+            [0 0 0]
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: x.operator() - y.operator() in J.Hom(J)
+            True
+
+        """
+        return self + (-other)
+
+
+    def matrix(self):
+        """
+        Return the matrix of this morphism with respect to a left-action
+        on column vectors.
+        """
+        return FiniteDimensionalAlgebraMorphism.matrix(self).transpose()
+
 
 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
     @staticmethod
@@ -30,7 +423,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 raise ValueError("input is not a multiplication table")
         mult_table = tuple(mult_table)
 
-        cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
+        cat = FiniteDimensionalAlgebrasWithBasis(field)
         cat.or_subcategory(category)
         if assume_associative:
             cat = cat.Associative()
@@ -48,6 +441,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                                  natural_basis=natural_basis)
 
 
+    def _Hom_(self, B, cat):
+        """
+        Construct a homset of ``self`` and ``B``.
+        """
+        return FiniteDimensionalEuclideanJordanAlgebraHomset(self,
+                                                             B,
+                                                             category=cat)
+
+
     def __init__(self,
                  field,
                  mult_table,
@@ -87,6 +489,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         return fmt.format(self.degree(), self.base_ring())
 
 
+    def _a_regular_element(self):
+        """
+        Guess a regular element. Needed to compute the basis for our
+        characteristic polynomial coefficients.
+        """
+        gs = self.gens()
+        z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
+        if not z.is_regular():
+            raise ValueError("don't know a regular element")
+        return z
+
+
+    @cached_method
+    def _charpoly_basis_space(self):
+        """
+        Return the vector space spanned by the basis used in our
+        characteristic polynomial coefficients. This is used not only to
+        compute those coefficients, but also any time we need to
+        evaluate the coefficients (like when we compute the trace or
+        determinant).
+        """
+        z = self._a_regular_element()
+        V = z.vector().parent().ambient_vector_space()
+        V1 = V.span_of_basis( (z**k).vector() for k in range(self.rank()) )
+        b =  (V1.basis() + V1.complement().basis())
+        return V.span_of_basis(b)
+
 
     @cached_method
     def _charpoly_coeff(self, i):
@@ -98,25 +527,38 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         store the trace/determinant (a_{r-1} and a_{0} respectively)
         separate from the entire characteristic polynomial.
         """
-        (A_of_x, x) = self._charpoly_matrix()
+        (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
         R = A_of_x.base_ring()
-        A_cols = A_of_x.columns()
-        A_cols[i] = (x**self.rank()).vector()
-        numerator = column_matrix(A_of_x.base_ring(), A_cols).det()
-        denominator = A_of_x.det()
+        if i >= self.rank():
+            # Guaranteed by theory
+            return R.zero()
+
+        # Danger: the in-place modification is done for performance
+        # reasons (reconstructing a matrix with huge polynomial
+        # entries is slow), but I don't know how cached_method works,
+        # so it's highly possible that we're modifying some global
+        # list variable by reference, here. In other words, you
+        # probably shouldn't call this method twice on the same
+        # algebra, at the same time, in two threads
+        Ai_orig = A_of_x.column(i)
+        A_of_x.set_column(i,xr)
+        numerator = A_of_x.det()
+        A_of_x.set_column(i,Ai_orig)
 
         # We're relying on the theory here to ensure that each a_i is
         # indeed back in R, and the added negative signs are to make
         # the whole charpoly expression sum to zero.
-        return R(-numerator/denominator)
+        return R(-numerator/detA)
 
 
     @cached_method
-    def _charpoly_matrix(self):
+    def _charpoly_matrix_system(self):
         """
         Compute the matrix whose entries A_ij are polynomials in
-        X1,...,XN. This same matrix is used in more than one method and
-        it's not so fast to construct.
+        X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
+        corresponding to `x^r` and the determinent of the matrix A =
+        [A_ij]. In other words, all of the fixed (cachable) data needed
+        to compute the coefficients of the characteristic polynomial.
         """
         r = self.rank()
         n = self.dimension()
@@ -130,16 +572,49 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         idmat = identity_matrix(J.base_ring(), n)
 
+        W = self._charpoly_basis_space()
+        W = W.change_ring(R.fraction_field())
+
+        # Starting with the standard coordinates x = (X1,X2,...,Xn)
+        # and then converting the entries to W-coordinates allows us
+        # to pass in the standard coordinates to the charpoly and get
+        # back the right answer. Specifically, with x = (X1,X2,...,Xn),
+        # we have
+        #
+        #   W.coordinates(x^2) eval'd at (standard z-coords)
+        #     =
+        #   W-coords of (z^2)
+        #     =
+        #   W-coords of (standard coords of x^2 eval'd at std-coords of z)
+        #
+        # We want the middle equivalent thing in our matrix, but use
+        # the first equivalent thing instead so that we can pass in
+        # standard coordinates.
         x = J(vector(R, R.gens()))
-        l1 = [column_matrix((x**k).vector()) for k in range(r)]
+        l1 = [column_matrix(W.coordinates((x**k).vector())) for k in range(r)]
         l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
         A_of_x = block_matrix(R, 1, n, (l1 + l2))
-        return (A_of_x, x)
+        xr = W.coordinates((x**r).vector())
+        return (A_of_x, x, xr, A_of_x.det())
 
 
     @cached_method
     def characteristic_polynomial(self):
         """
+
+        .. WARNING::
+
+            This implementation doesn't guarantee that the polynomial
+            denominator in the coefficients is not identically zero, so
+            theoretically it could crash. The way that this is handled
+            in e.g. Faraut and Koranyi is to use a basis that guarantees
+            the denominator is non-zero. But, doing so requires knowledge
+            of at least one regular element, and we don't even know how
+            to do that. The trade-off is that, if we use the standard basis,
+            the resulting polynomial will accept the "usual" coordinates. In
+            other words, we don't have to do a change of basis before e.g.
+            computing the trace or determinant.
+
         EXAMPLES:
 
         The characteristic polynomial in the spin algebra is given in
@@ -267,6 +742,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         An element of a Euclidean Jordan algebra.
         """
 
+        def __dir__(self):
+            """
+            Oh man, I should not be doing this. This hides the "disabled"
+            methods ``left_matrix`` and ``matrix`` from introspection;
+            in particular it removes them from tab-completion.
+            """
+            return filter(lambda s: s not in ['left_matrix', 'matrix'],
+                          dir(self.__class__) )
+
+
         def __init__(self, A, elt=None):
             """
             EXAMPLES:
@@ -585,8 +1070,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             """
             Return the Jordan-multiplicative inverse of this element.
 
-            We can't use the superclass method because it relies on the
-            algebra being associative.
+            ALGORITHM:
+
+            We appeal to the quadratic representation as in Koecher's
+            Theorem 12 in Chapter III, Section 5.
 
             EXAMPLES:
 
@@ -597,12 +1084,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: n = ZZ.random_element(1,10)
                 sage: J = JordanSpinEJA(n)
                 sage: x = J.random_element()
-                sage: while x.is_zero():
+                sage: while not x.is_invertible():
                 ....:     x = J.random_element()
                 sage: x_vec = x.vector()
                 sage: x0 = x_vec[0]
                 sage: x_bar = x_vec[1:]
-                sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
+                sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
                 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
                 sage: x_inverse = coeff*inv_vec
                 sage: x.inverse() == J(x_inverse)
@@ -625,35 +1112,27 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
                 True
 
-            """
-            if self.parent().is_associative():
-                elt = FiniteDimensionalAlgebraElement(self.parent(), self)
-                return elt.inverse()
+            The inverse of the inverse is what we started with::
 
-            # TODO: we can do better once the call to is_invertible()
-            # doesn't crash on irregular elements.
-            #if not self.is_invertible():
-            #    raise ValueError('element is not invertible')
+                sage: set_random_seed()
+                sage: J = random_eja()
+                sage: x = J.random_element()
+                sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
+                True
 
-            # We do this a little different than the usual recursive
-            # call to a finite-dimensional algebra element, because we
-            # wind up with an inverse that lives in the subalgebra and
-            # we need information about the parent to convert it back.
-            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()))
+            The zero element is never invertible::
 
-            # This will be in the subalgebra's coordinates...
-            fda_elt = FiniteDimensionalAlgebraElement(assoc_subalg, elt)
-            subalg_inverse = fda_elt.inverse()
+                sage: set_random_seed()
+                sage: J = random_eja().zero().inverse()
+                Traceback (most recent call last):
+                ...
+                ValueError: element is not invertible
 
-            # So we have to convert back...
-            basis = [ self.parent(v) for v in V.basis() ]
-            pairs = zip(subalg_inverse.vector(), basis)
-            return self.parent().linear_combination(pairs)
+            """
+            if not self.is_invertible():
+                raise ValueError("element is not invertible")
+
+            return (~self.quadratic_representation())(self)
 
 
         def is_invertible(self):
@@ -794,6 +1273,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             return self.span_of_powers().dimension()
 
 
+        def left_matrix(self):
+            """
+            Our parent class defines ``left_matrix`` and ``matrix``
+            methods whose names are misleading. We don't want them.
+            """
+            raise NotImplementedError("use operator_matrix() instead")
+
+        matrix = left_matrix
+
+
         def minimal_polynomial(self):
             """
             Return the minimal polynomial of this element,
@@ -915,14 +1404,37 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             return W.linear_combination(zip(self.vector(), B))
 
 
+        def operator(self):
+            """
+            Return the left-multiplication-by-this-element
+            operator on the ambient algebra.
+
+            TESTS::
+
+                sage: set_random_seed()
+                sage: J = random_eja()
+                sage: x = J.random_element()
+                sage: y = J.random_element()
+                sage: x.operator()(y) == x*y
+                True
+                sage: y.operator()(x) == x*y
+                True
+
+            """
+            P = self.parent()
+            return FiniteDimensionalEuclideanJordanAlgebraMorphism(
+                     Hom(P,P),
+                     self.operator_matrix() )
+
+
+
         def operator_matrix(self):
             """
             Return the matrix that represents left- (or right-)
             multiplication by this element in the parent algebra.
 
-            We have to override this because the superclass method
-            returns a matrix that acts on row vectors (that is, on
-            the right).
+            We implement this ourselves to work around the fact that
+            our parent class represents everything with row vectors.
 
             EXAMPLES:
 
@@ -1005,7 +1517,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
                 sage: D = D + 2*x_bar.tensor_product(x_bar)
                 sage: Q = block_matrix(2,2,[A,B,C,D])
-                sage: Q == x.quadratic_representation()
+                sage: Q == x.quadratic_representation().matrix()
                 True
 
             Test all of the properties from Theorem 11.2 in Alizadeh::
@@ -1014,38 +1526,76 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: J = random_eja()
                 sage: x = J.random_element()
                 sage: y = J.random_element()
+                sage: Lx = x.operator()
+                sage: Lxx = (x*x).operator()
+                sage: Qx = x.quadratic_representation()
+                sage: Qy = y.quadratic_representation()
+                sage: Qxy = x.quadratic_representation(y)
+                sage: Qex = J.one().quadratic_representation(x)
+                sage: n = ZZ.random_element(10)
+                sage: Qxn = (x^n).quadratic_representation()
 
             Property 1:
 
-                sage: actual = x.quadratic_representation(y)
-                sage: expected = ( (x+y).quadratic_representation()
-                ....:              -x.quadratic_representation()
-                ....:              -y.quadratic_representation() ) / 2
-                sage: actual == expected
+                sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
                 True
 
             Property 2:
 
                 sage: alpha = QQ.random_element()
-                sage: actual = (alpha*x).quadratic_representation()
-                sage: expected = (alpha^2)*x.quadratic_representation()
-                sage: actual == expected
+                sage: (alpha*x).quadratic_representation() == (alpha^2)*Qx
+                True
+
+            Property 3:
+
+                sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
+                True
+
+                sage: not x.is_invertible() or (
+                ....:   ~Qx
+                ....:   ==
+                ....:   x.inverse().quadratic_representation() )
+                True
+
+                sage: Qxy(J.one()) == x*y
+                True
+
+            Property 4:
+
+                sage: not x.is_invertible() or (
+                ....:   x.quadratic_representation(x.inverse())*Qx
+                ....:   == Qx*x.quadratic_representation(x.inverse()) )
+                True
+
+                sage: not x.is_invertible() or (
+                ....:   x.quadratic_representation(x.inverse())*Qx
+                ....:   ==
+                ....:   2*x.operator()*Qex - Qx )
+                True
+
+                sage: 2*x.operator()*Qex - Qx == Lxx
                 True
 
             Property 5:
 
-                sage: Qy = y.quadratic_representation()
-                sage: actual = J(Qy*x.vector()).quadratic_representation()
-                sage: expected = Qy*x.quadratic_representation()*Qy
-                sage: actual == expected
+                sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
                 True
 
             Property 6:
 
-                sage: k = ZZ.random_element(1,10)
-                sage: actual = (x^k).quadratic_representation()
-                sage: expected = (x.quadratic_representation())^k
-                sage: actual == expected
+                sage: Qxn == (Qx)^n
+                True
+
+            Property 7:
+
+                sage: not x.is_invertible() or (
+                ....:   Qx*x.inverse().operator() == Lx )
+                True
+
+            Property 8:
+
+                sage: not x.operator_commutes_with(y) or (
+                ....:   Qx(y)^n == Qxn(y^n) )
                 True
 
             """
@@ -1054,9 +1604,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             elif not other in self.parent():
                 raise TypeError("'other' must live in the same algebra")
 
-            L = self.operator_matrix()
-            M = other.operator_matrix()
-            return ( L*M + M*L - (self*other).operator_matrix() )
+            L = self.operator()
+            M = other.operator()
+            return ( L*M + M*L - (self*other).operator() )
 
 
         def span_of_powers(self):
@@ -1139,12 +1689,11 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             TESTS::
 
                 sage: set_random_seed()
-                sage: J = RealCartesianProductEJA(5)
-                sage: c = J.random_element().subalgebra_idempotent()
-                sage: c^2 == c
-                True
-                sage: J = JordanSpinEJA(5)
-                sage: c = J.random_element().subalgebra_idempotent()
+                sage: J = random_eja()
+                sage: x = J.random_element()
+                sage: while x.is_nilpotent():
+                ....:     x = J.random_element()
+                sage: c = x.subalgebra_idempotent()
                 sage: c^2 == c
                 True
 
@@ -1233,6 +1782,47 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         def trace_inner_product(self, other):
             """
             Return the trace inner product of myself and ``other``.
+
+            TESTS:
+
+            The trace inner product is commutative::
+
+                sage: set_random_seed()
+                sage: J = random_eja()
+                sage: x = J.random_element(); y = J.random_element()
+                sage: x.trace_inner_product(y) == y.trace_inner_product(x)
+                True
+
+            The trace inner product is bilinear::
+
+                sage: set_random_seed()
+                sage: J = random_eja()
+                sage: x = J.random_element()
+                sage: y = J.random_element()
+                sage: z = J.random_element()
+                sage: a = QQ.random_element();
+                sage: actual = (a*(x+z)).trace_inner_product(y)
+                sage: expected = ( a*x.trace_inner_product(y) +
+                ....:              a*z.trace_inner_product(y) )
+                sage: actual == expected
+                True
+                sage: actual = x.trace_inner_product(a*(y+z))
+                sage: expected = ( a*x.trace_inner_product(y) +
+                ....:              a*x.trace_inner_product(z) )
+                sage: actual == expected
+                True
+
+            The trace inner product satisfies the compatibility
+            condition in the definition of 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).trace_inner_product(z) == y.trace_inner_product(x*z)
+                True
+
             """
             if not other in self.parent():
                 raise TypeError("'other' must live in the same algebra")