eja: speed up det() for Cartesian product elements.

[sage.d.git] / mjo / eja / eja_element.py

from sage.matrix.constructor import matrix
from sage.misc.cachefunc import cached_method
from sage.modules.free_module import VectorSpace
from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement

from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
from mjo.eja.eja_utils import _scale


class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
    """
    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 __pow__(self, n):
        """
        Return ``self`` raised to the power ``n``.

        Jordan algebras are always power-associative; see for
        example Faraut and Korányi, Proposition II.1.2 (ii).

        We have to override this because our superclass uses row
        vectors instead of column vectors! We, on the other hand,
        assume column vectors everywhere.

        SETUP::

            sage: from mjo.eja.eja_algebra import random_eja

        TESTS:

        The definition of `x^2` is the unambiguous `x*x`::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: x*x == (x^2)
            True

        A few examples of power-associativity::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: x*(x*x)*(x*x) == x^5
            True
            sage: (x*x)*(x*x*x) == x^5
            True

        We also know that powers operator-commute (Koecher, Chapter
        III, Corollary 1)::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: m = ZZ.random_element(0,10)
            sage: n = ZZ.random_element(0,10)
            sage: Lxm = (x^m).operator()
            sage: Lxn = (x^n).operator()
            sage: Lxm*Lxn == Lxn*Lxm
            True

        """
        if n == 0:
            return self.parent().one()
        elif n == 1:
            return self
        else:
            return (self**(n-1))*self


    def apply_univariate_polynomial(self, p):
        """
        Apply the univariate polynomial ``p`` to this element.

        A priori, SageMath won't allow us to apply a univariate
        polynomial to an element of an EJA, because we don't know
        that EJAs are rings (they are usually not associative). Of
        course, we know that EJAs are power-associative, so the
        operation is ultimately kosher. This function sidesteps
        the CAS to get the answer we want and expect.

        SETUP::

            sage: from mjo.eja.eja_algebra import (HadamardEJA,
            ....:                                  random_eja)

        EXAMPLES::

            sage: R = PolynomialRing(QQ, 't')
            sage: t = R.gen(0)
            sage: p = t^4 - t^3 + 5*t - 2
            sage: J = HadamardEJA(5)
            sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
            True

        TESTS:

        We should always get back an element of the algebra::

            sage: set_random_seed()
            sage: p = PolynomialRing(AA, 't').random_element()
            sage: J = random_eja()
            sage: x = J.random_element()
            sage: x.apply_univariate_polynomial(p) in J
            True

        """
        if len(p.variables()) > 1:
            raise ValueError("not a univariate polynomial")
        P = self.parent()
        R = P.base_ring()
        # Convert the coeficcients to the parent's base ring,
        # because a priori they might live in an (unnecessarily)
        # larger ring for which P.sum() would fail below.
        cs = [ R(c) for c in p.coefficients(sparse=False) ]
        return P.sum( cs[k]*(self**k) for k in range(len(cs)) )


    def characteristic_polynomial(self):
        """
        Return the characteristic polynomial of this element.

        SETUP::

            sage: from mjo.eja.eja_algebra import (random_eja,
            ....:                                  HadamardEJA)

        EXAMPLES:

        The rank of `R^3` is three, and the minimal polynomial of
        the identity element is `(t-1)` from which it follows that
        the characteristic polynomial should be `(t-1)^3`::

            sage: J = HadamardEJA(3)
            sage: J.one().characteristic_polynomial()
            t^3 - 3*t^2 + 3*t - 1

        Likewise, the characteristic of the zero element in the
        rank-three algebra `R^{n}` should be `t^{3}`::

            sage: J = HadamardEJA(3)
            sage: J.zero().characteristic_polynomial()
            t^3

        TESTS:

        The characteristic polynomial of an element should evaluate
        to zero on that element::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: p = x.characteristic_polynomial()
            sage: x.apply_univariate_polynomial(p).is_zero()
            True

        The characteristic polynomials of the zero and unit elements
        should be what we think they are in a subalgebra, too::

            sage: J = HadamardEJA(3)
            sage: p1 = J.one().characteristic_polynomial()
            sage: q1 = J.zero().characteristic_polynomial()
            sage: b0,b1,b2 = J.gens()
            sage: A = (b0 + 2*b1 + 3*b2).subalgebra_generated_by() # dim 3
            sage: p2 = A.one().characteristic_polynomial()
            sage: q2 = A.zero().characteristic_polynomial()
            sage: p1 == p2
            True
            sage: q1 == q2
            True

        """
        p = self.parent().characteristic_polynomial_of()
        return p(*self.to_vector())


    def inner_product(self, other):
        """
        Return the parent algebra's inner product of myself and ``other``.

        SETUP::

            sage: from mjo.eja.eja_algebra import (
            ....:   ComplexHermitianEJA,
            ....:   JordanSpinEJA,
            ....:   QuaternionHermitianEJA,
            ....:   RealSymmetricEJA,
            ....:   random_eja)

        EXAMPLES:

        The inner product in the Jordan spin algebra is the usual
        inner product on `R^n` (this example only works because the
        basis for the Jordan algebra is the standard basis in `R^n`)::

            sage: J = JordanSpinEJA(3)
            sage: x = vector(QQ,[1,2,3])
            sage: y = vector(QQ,[4,5,6])
            sage: x.inner_product(y)
            32
            sage: J.from_vector(x).inner_product(J.from_vector(y))
            32

        The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
        multiplication is the usual matrix multiplication in `S^n`,
        so the inner product of the identity matrix with itself
        should be the `n`::

            sage: J = RealSymmetricEJA(3)
            sage: J.one().inner_product(J.one())
            3

        Likewise, the inner product on `C^n` is `<X,Y> =
        Re(trace(X*Y))`, where we must necessarily take the real
        part because the product of Hermitian matrices may not be
        Hermitian::

            sage: J = ComplexHermitianEJA(3)
            sage: J.one().inner_product(J.one())
            3

        Ditto for the quaternions::

            sage: J = QuaternionHermitianEJA(2)
            sage: J.one().inner_product(J.one())
            2

        TESTS:

        Ensure that we can always compute an inner product, and that
        it gives us back a real number::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x,y = J.random_elements(2)
            sage: x.inner_product(y) in RLF
            True

        """
        P = self.parent()
        if not other in P:
            raise TypeError("'other' must live in the same algebra")

        return P.inner_product(self, other)


    def operator_commutes_with(self, other):
        """
        Return whether or not this element operator-commutes
        with ``other``.

        SETUP::

            sage: from mjo.eja.eja_algebra import random_eja

        EXAMPLES:

        The definition of a Jordan algebra says that any element
        operator-commutes with its square::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: x.operator_commutes_with(x^2)
            True

        TESTS:

        Test Lemma 1 from Chapter III of Koecher::

            sage: set_random_seed()
            sage: u,v = random_eja().random_elements(2)
            sage: lhs = u.operator_commutes_with(u*v)
            sage: rhs = v.operator_commutes_with(u^2)
            sage: lhs == rhs
            True

        Test the first polarization identity from my notes, Koecher
        Chapter III, or from Baes (2.3)::

            sage: set_random_seed()
            sage: x,y = random_eja().random_elements(2)
            sage: Lx = x.operator()
            sage: Ly = y.operator()
            sage: Lxx = (x*x).operator()
            sage: Lxy = (x*y).operator()
            sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
            True

        Test the second polarization identity from my notes or from
        Baes (2.4)::

            sage: set_random_seed()
            sage: x,y,z = random_eja().random_elements(3)
            sage: Lx = x.operator()
            sage: Ly = y.operator()
            sage: Lz = z.operator()
            sage: Lzy = (z*y).operator()
            sage: Lxy = (x*y).operator()
            sage: Lxz = (x*z).operator()
            sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
            True

        Test the third polarization identity from my notes or from
        Baes (2.5)::

            sage: set_random_seed()
            sage: u,y,z = random_eja().random_elements(3)
            sage: Lu = u.operator()
            sage: Ly = y.operator()
            sage: Lz = z.operator()
            sage: Lzy = (z*y).operator()
            sage: Luy = (u*y).operator()
            sage: Luz = (u*z).operator()
            sage: Luyz = (u*(y*z)).operator()
            sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
            sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
            sage: bool(lhs == rhs)
            True

        """
        if not other in self.parent():
            raise TypeError("'other' must live in the same algebra")

        A = self.operator()
        B = other.operator()
        return (A*B == B*A)


    def det(self):
        """
        Return my determinant, the product of my eigenvalues.

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  TrivialEJA,
            ....:                                  RealSymmetricEJA,
            ....:                                  ComplexHermitianEJA,
            ....:                                  random_eja)

        EXAMPLES::

            sage: J = JordanSpinEJA(2)
            sage: x = sum( J.gens() )
            sage: x.det()
            0

        ::

            sage: J = JordanSpinEJA(3)
            sage: x = sum( J.gens() )
            sage: x.det()
            -1

        The determinant of the sole element in the rank-zero trivial
        algebra is ``1``, by three paths of reasoning. First, its
        characteristic polynomial is a constant ``1``, so the constant
        term in that polynomial is ``1``. Second, the characteristic
        polynomial evaluated at zero is again ``1``. And finally, the
        (empty) product of its eigenvalues is likewise just unity::

            sage: J = TrivialEJA()
            sage: J.zero().det()
            1

        TESTS:

        An element is invertible if and only if its determinant is
        non-zero::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: x.is_invertible() == (x.det() != 0)
            True

        Ensure that the determinant is multiplicative on an associative
        subalgebra as in Faraut and Korányi's Proposition II.2.2::

            sage: set_random_seed()
            sage: J = random_eja().random_element().subalgebra_generated_by()
            sage: x,y = J.random_elements(2)
            sage: (x*y).det() == x.det()*y.det()
            True

        The determinant in real matrix algebras is the usual determinant::

            sage: set_random_seed()
            sage: X = matrix.random(QQ,3)
            sage: X = X + X.T
            sage: J1 = RealSymmetricEJA(3)
            sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
            sage: expected = X.det()
            sage: actual1 = J1(X).det()
            sage: actual2 = J2(X).det()
            sage: actual1 == expected
            True
            sage: actual2 == expected
            True

        """
        P = self.parent()
        r = P.rank()

        if r == 0:
            # Special case, since we don't get the a0=1
            # coefficient when the rank of the algebra
            # is zero.
            return P.base_ring().one()

        p = P._charpoly_coefficients()[0]
        # The _charpoly_coeff function already adds the factor of -1
        # to ensure that _charpoly_coefficients()[0] is really what
        # appears in front of t^{0} in the charpoly. However, we want
        # (-1)^r times THAT for the determinant.
        return ((-1)**r)*p(*self.to_vector())


    @cached_method
    def inverse(self):
        """
        Return the Jordan-multiplicative inverse of this element.

        ALGORITHM:

        In general we appeal to the quadratic representation as in
        Koecher's Theorem 12 in Chapter III, Section 5. But if the
        parent algebra's "characteristic polynomial of" coefficients
        happen to be cached, then we use Proposition II.2.4 in Faraut
        and Korányi which gives a formula for the inverse based on the
        characteristic polynomial and the Cayley-Hamilton theorem for
        Euclidean Jordan algebras::

        SETUP::

            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
            ....:                                  JordanSpinEJA,
            ....:                                  random_eja)

        EXAMPLES:

        The inverse in the spin factor algebra is given in Alizadeh's
        Example 11.11::

            sage: set_random_seed()
            sage: J = JordanSpinEJA.random_instance()
            sage: x = J.random_element()
            sage: while not x.is_invertible():
            ....:     x = J.random_element()
            sage: x_vec = x.to_vector()
            sage: x0 = x_vec[:1]
            sage: x_bar = x_vec[1:]
            sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
            sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
            sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
            sage: x.inverse() == J.from_vector(x_inverse)
            True

        Trying to invert a non-invertible element throws an error:

            sage: JordanSpinEJA(3).zero().inverse()
            Traceback (most recent call last):
            ...
            ZeroDivisionError: element is not invertible

        TESTS:

        The identity element is its own inverse::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: J.one().inverse() == J.one()
            True

        If an element has an inverse, it acts like one::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x = J.random_element()
            sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
            True

        The inverse of the inverse is what we started with::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x = J.random_element()
            sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
            True

        Proposition II.2.3 in Faraut and Korányi says that the inverse
        of an element is the inverse of its left-multiplication operator
        applied to the algebra's identity, when that inverse exists::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x = J.random_element()
            sage: (not x.operator().is_invertible()) or (
            ....:    x.operator().inverse()(J.one()) == x.inverse() )
            True

        Check that the fast (cached) and slow algorithms give the same
        answer::

            sage: set_random_seed()                              # long time
            sage: J = random_eja(field=QQ, orthonormalize=False) # long time
            sage: x = J.random_element()                         # long time
            sage: while not x.is_invertible():                   # long time
            ....:     x = J.random_element()                     # long time
            sage: slow = x.inverse()                             # long time
            sage: _ = J._charpoly_coefficients()                 # long time
            sage: fast = x.inverse()                             # long time
            sage: slow == fast                                   # long time
            True
        """
        not_invertible_msg = "element is not invertible"
        if self.parent()._charpoly_coefficients.is_in_cache():
            # We can invert using our charpoly if it will be fast to
            # compute. If the coefficients are cached, our rank had
            # better be too!
            if self.det().is_zero():
                raise ZeroDivisionError(not_invertible_msg)
            r = self.parent().rank()
            a = self.characteristic_polynomial().coefficients(sparse=False)
            return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det()

        try:
            inv = (~self.quadratic_representation())(self)
            self.is_invertible.set_cache(True)
            return inv
        except ZeroDivisionError:
            self.is_invertible.set_cache(False)
            raise ZeroDivisionError(not_invertible_msg)


    @cached_method
    def is_invertible(self):
        """
        Return whether or not this element is invertible.

        ALGORITHM:

        If computing my determinant will be fast, we do so and compare
        with zero (Proposition II.2.4 in Faraut and
        Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
        reduces the problem to the invertibility of my quadratic
        representation.

        SETUP::

            sage: from mjo.eja.eja_algebra import random_eja

        TESTS:

        The identity element is always invertible::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: J.one().is_invertible()
            True

        The zero element is never invertible in a non-trivial algebra::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: (not J.is_trivial()) and J.zero().is_invertible()
            False

        Test that the fast (cached) and slow algorithms give the same
        answer::

            sage: set_random_seed()                              # long time
            sage: J = random_eja(field=QQ, orthonormalize=False) # long time
            sage: x = J.random_element()                         # long time
            sage: slow = x.is_invertible()                       # long time
            sage: _ = J._charpoly_coefficients()                 # long time
            sage: fast = x.is_invertible()                       # long time
            sage: slow == fast                                   # long time
            True
        """
        if self.is_zero():
            if self.parent().is_trivial():
                return True
            else:
                return False

        if self.parent()._charpoly_coefficients.is_in_cache():
            # The determinant will be quicker than inverting the
            # quadratic representation, most likely.
            return (not self.det().is_zero())

        # The easiest way to determine if I'm invertible is to try.
        try:
            inv = (~self.quadratic_representation())(self)
            self.inverse.set_cache(inv)
            return True
        except ZeroDivisionError:
            return False


    def is_primitive_idempotent(self):
        """
        Return whether or not this element is a primitive (or minimal)
        idempotent.

        A primitive idempotent is a non-zero idempotent that is not
        the sum of two other non-zero idempotents. Remark 2.7.15 in
        Baes shows that this is what he refers to as a "minimal
        idempotent."

        An element of a Euclidean Jordan algebra is a minimal idempotent
        if it :meth:`is_idempotent` and if its Peirce subalgebra
        corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
        Proposition 2.7.17).

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  RealSymmetricEJA,
            ....:                                  TrivialEJA,
            ....:                                  random_eja)

        WARNING::

        This method is sloooooow.

        EXAMPLES:

        The spectral decomposition of a non-regular element should always
        contain at least one non-minimal idempotent::

            sage: J = RealSymmetricEJA(3)
            sage: x = sum(J.gens())
            sage: x.is_regular()
            False
            sage: [ c.is_primitive_idempotent()
            ....:   for (l,c) in x.spectral_decomposition() ]
            [False, True]

        On the other hand, the spectral decomposition of a regular
        element should always be in terms of minimal idempotents::

            sage: J = JordanSpinEJA(4)
            sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
            sage: x.is_regular()
            True
            sage: [ c.is_primitive_idempotent()
            ....:   for (l,c) in x.spectral_decomposition() ]
            [True, True]

        TESTS:

        The identity element is minimal only in an EJA of rank one::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
            True

        A non-idempotent cannot be a minimal idempotent::

            sage: set_random_seed()
            sage: J = JordanSpinEJA(4)
            sage: x = J.random_element()
            sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
            False

        Proposition 2.7.19 in Baes says that an element is a minimal
        idempotent if and only if it's idempotent with trace equal to
        unity::

            sage: set_random_seed()
            sage: J = JordanSpinEJA(4)
            sage: x = J.random_element()
            sage: expected = (x.is_idempotent() and x.trace() == 1)
            sage: actual = x.is_primitive_idempotent()
            sage: actual == expected
            True

        Primitive idempotents must be non-zero::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: J.zero().is_idempotent()
            True
            sage: J.zero().is_primitive_idempotent()
            False

        As a consequence of the fact that primitive idempotents must
        be non-zero, there are no primitive idempotents in a trivial
        Euclidean Jordan algebra::

            sage: J = TrivialEJA()
            sage: J.one().is_idempotent()
            True
            sage: J.one().is_primitive_idempotent()
            False

        """
        if not self.is_idempotent():
            return False

        if self.is_zero():
            return False

        (_,_,J1) = self.parent().peirce_decomposition(self)
        return (J1.dimension() == 1)


    def is_nilpotent(self):
        """
        Return whether or not some power of this element is zero.

        ALGORITHM:

        We use Theorem 5 in Chapter III of Koecher, which says that
        an element ``x`` is nilpotent if and only if ``x.operator()``
        is nilpotent. And it is a basic fact of linear algebra that
        an operator on an `n`-dimensional space is nilpotent if and
        only if, when raised to the `n`th power, it equals the zero
        operator (for example, see Axler Corollary 8.8).

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  random_eja)

        EXAMPLES::

            sage: J = JordanSpinEJA(3)
            sage: x = sum(J.gens())
            sage: x.is_nilpotent()
            False

        TESTS:

        The identity element is never nilpotent, except in a trivial EJA::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: J.one().is_nilpotent() and not J.is_trivial()
            False

        The additive identity is always nilpotent::

            sage: set_random_seed()
            sage: random_eja().zero().is_nilpotent()
            True

        """
        P = self.parent()
        zero_operator = P.zero().operator()
        return self.operator()**P.dimension() == zero_operator


    def is_regular(self):
        """
        Return whether or not this is a regular element.

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  random_eja)

        EXAMPLES:

        The identity element always has degree one, but any element
        linearly-independent from it is regular::

            sage: J = JordanSpinEJA(5)
            sage: J.one().is_regular()
            False
            sage: b0, b1, b2, b3, b4 = J.gens()
            sage: b0 == J.one()
            True
            sage: for x in J.gens():
            ....:     (J.one() + x).is_regular()
            False
            True
            True
            True
            True

        TESTS:

        The zero element should never be regular, unless the parent
        algebra has dimension less than or equal to one::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: J.dimension() <= 1 or not J.zero().is_regular()
            True

        The unit element isn't regular unless the algebra happens to
        consist of only its scalar multiples::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: J.dimension() <= 1 or not J.one().is_regular()
            True

        """
        return self.degree() == self.parent().rank()


    def degree(self):
        """
        Return the degree of this element, which is defined to be
        the degree of its minimal polynomial.

        ALGORITHM:

        .........

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  random_eja)

        EXAMPLES::

            sage: J = JordanSpinEJA(4)
            sage: J.one().degree()
            1
            sage: b0,b1,b2,b3 = J.gens()
            sage: (b0 - b1).degree()
            2

        In the spin factor algebra (of rank two), all elements that
        aren't multiples of the identity are regular::

            sage: set_random_seed()
            sage: J = JordanSpinEJA.random_instance()
            sage: n = J.dimension()
            sage: x = J.random_element()
            sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
            True

        TESTS:

        The zero and unit elements are both of degree one in nontrivial
        algebras::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: d = J.zero().degree()
            sage: (J.is_trivial() and d == 0) or d == 1
            True
            sage: d = J.one().degree()
            sage: (J.is_trivial() and d == 0) or d == 1
            True

        Our implementation agrees with the definition::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: x.degree() == x.minimal_polynomial().degree()
            True

        """
        n = self.parent().dimension()

        if n == 0:
            # The minimal polynomial is an empty product, i.e. the
            # constant polynomial "1" having degree zero.
            return 0
        elif self.is_zero():
            # The minimal polynomial of zero in a nontrivial algebra
            # is "t", and is of degree one.
            return 1
        elif n == 1:
            # If this is a nonzero element of a nontrivial algebra, it
            # has degree at least one. It follows that, in an algebra
            # of dimension one, the degree must be actually one.
            return 1

        # BEWARE: The subalgebra_generated_by() method uses the result
        # of this method to construct a basis for the subalgebra. That
        # means, in particular, that we cannot implement this method
        # as ``self.subalgebra_generated_by().dimension()``.

        # Algorithm: keep appending (vector representations of) powers
        # self as rows to a matrix and echelonizing it. When its rank
        # stops increasing, we've reached a redundancy.

        # Given the special cases above, we can assume that "self" is
        # nonzero, the algebra is nontrivial, and that its dimension
        # is at least two.
        M = matrix([(self.parent().one()).to_vector()])
        old_rank = 1

        # Specifying the row-reduction algorithm can e.g. help over
        # AA because it avoids the RecursionError that gets thrown
        # when we have to look too hard for a root.
        #
        # Beware: QQ supports an entirely different set of "algorithm"
        # keywords than do AA and RR.
        algo = None
        from sage.rings.all import QQ
        if self.parent().base_ring() is not QQ:
            algo = "scaled_partial_pivoting"

        for d in range(1,n):
            M = matrix(M.rows() + [(self**d).to_vector()])
            M.echelonize(algo)
            new_rank = M.rank()
            if new_rank == old_rank:
                return new_rank
            else:
                old_rank = new_rank

        return n



    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,
        as a function of the variable `t`.

        ALGORITHM:

        We restrict ourselves to the associative subalgebra
        generated by this element, and then return the minimal
        polynomial of this element's operator matrix (in that
        subalgebra). This works by Baes Proposition 2.3.16.

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  RealSymmetricEJA,
            ....:                                  TrivialEJA,
            ....:                                  random_eja)

        EXAMPLES:

        Keeping in mind that the polynomial ``1`` evaluates the identity
        element (also the zero element) of the trivial algebra, it is clear
        that the polynomial ``1`` is the minimal polynomial of the only
        element in a trivial algebra::

            sage: J = TrivialEJA()
            sage: J.one().minimal_polynomial()
            1
            sage: J.zero().minimal_polynomial()
            1

        TESTS:

        The minimal polynomial of the identity and zero elements are
        always the same, except in trivial algebras where the minimal
        polynomial of the unit/zero element is ``1``::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: mu = J.one().minimal_polynomial()
            sage: t = mu.parent().gen()
            sage: mu + int(J.is_trivial())*(t-2)
            t - 1
            sage: mu = J.zero().minimal_polynomial()
            sage: t = mu.parent().gen()
            sage: mu + int(J.is_trivial())*(t-1)
            t

        The degree of an element is (by one definition) the degree
        of its minimal polynomial::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: x.degree() == x.minimal_polynomial().degree()
            True

        The minimal polynomial and the characteristic polynomial coincide
        and are known (see Alizadeh, Example 11.11) for all elements of
        the spin factor algebra that aren't scalar multiples of the
        identity. We require the dimension of the algebra to be at least
        two here so that said elements actually exist::

            sage: set_random_seed()
            sage: d_max = JordanSpinEJA._max_random_instance_dimension()
            sage: n = ZZ.random_element(2, max(2,d_max))
            sage: J = JordanSpinEJA(n)
            sage: y = J.random_element()
            sage: while y == y.coefficient(0)*J.one():
            ....:     y = J.random_element()
            sage: y0 = y.to_vector()[0]
            sage: y_bar = y.to_vector()[1:]
            sage: actual = y.minimal_polynomial()
            sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
            sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
            sage: bool(actual == expected)
            True

        The minimal polynomial should always kill its element::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: p = x.minimal_polynomial()
            sage: x.apply_univariate_polynomial(p)
            0

        The minimal polynomial is invariant under a change of basis,
        and in particular, a re-scaling of the basis::

            sage: set_random_seed()
            sage: d_max = RealSymmetricEJA._max_random_instance_dimension()
            sage: d = ZZ.random_element(1, d_max)
            sage: n = RealSymmetricEJA._max_random_instance_size(d)
            sage: J1 = RealSymmetricEJA(n)
            sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
            sage: X = random_matrix(AA,n)
            sage: X = X*X.transpose()
            sage: x1 = J1(X)
            sage: x2 = J2(X)
            sage: x1.minimal_polynomial() == x2.minimal_polynomial()
            True

        """
        if self.is_zero():
            # Pretty sure we know what the minimal polynomial of
            # the zero operator is going to be. This ensures
            # consistency of e.g. the polynomial variable returned
            # in the "normal" case without us having to think about it.
            return self.operator().minimal_polynomial()

        # If we don't orthonormalize the subalgebra's basis, then the
        # first two monomials in the subalgebra will be self^0 and
        # self^1... assuming that self^1 is not a scalar multiple of
        # self^0 (the unit element). We special case these to avoid
        # having to solve a system to coerce self into the subalgebra.
        A = self.subalgebra_generated_by(orthonormalize=False)

        if A.dimension() == 1:
            # Does a solve to find the scalar multiple alpha such that
            # alpha*unit = self. We have to do this because the basis
            # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
            unit = self.parent().one()
            alpha = self.to_vector() / unit.to_vector()
            return (unit.operator()*alpha).minimal_polynomial()
        else:
            # If the dimension of the subalgebra is >= 2, then we just
            # use the second basis element.
            return A.monomial(1).operator().minimal_polynomial()



    def to_matrix(self):
        """
        Return an (often more natural) representation of this element as a
        matrix.

        Every finite-dimensional Euclidean Jordan Algebra is a direct
        sum of five simple algebras, four of which comprise Hermitian
        matrices. This method returns a "natural" matrix
        representation of this element as either a Hermitian matrix or
        column vector.

        SETUP::

            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
            ....:                                  HadamardEJA,
            ....:                                  QuaternionHermitianEJA,
            ....:                                  RealSymmetricEJA)

        EXAMPLES::

            sage: J = ComplexHermitianEJA(3)
            sage: J.one()
            b0 + b3 + b8
            sage: J.one().to_matrix()
            +---+---+---+
            | 1 | 0 | 0 |
            +---+---+---+
            | 0 | 1 | 0 |
            +---+---+---+
            | 0 | 0 | 1 |
            +---+---+---+

        ::

            sage: J = QuaternionHermitianEJA(2)
            sage: J.one()
            b0 + b5
            sage: J.one().to_matrix()
            +---+---+
            | 1 | 0 |
            +---+---+
            | 0 | 1 |
            +---+---+

        This also works in Cartesian product algebras::

            sage: J1 = HadamardEJA(1)
            sage: J2 = RealSymmetricEJA(2)
            sage: J = cartesian_product([J1,J2])
            sage: x = sum(J.gens())
            sage: x.to_matrix()[0]
            [1]
            sage: x.to_matrix()[1]
            [                  1 0.7071067811865475?]
            [0.7071067811865475?                   1]

        """
        B = self.parent().matrix_basis()
        W = self.parent().matrix_space()

        # This is just a manual "from_vector()", but of course
        # matrix spaces aren't vector spaces in sage, so they
        # don't have a from_vector() method.
        return W.linear_combination( zip(B, self.to_vector()) )



    def norm(self):
        """
        The norm of this element with respect to :meth:`inner_product`.

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  HadamardEJA)

        EXAMPLES::

            sage: J = HadamardEJA(2)
            sage: x = sum(J.gens())
            sage: x.norm()
            1.414213562373095?

        ::

            sage: J = JordanSpinEJA(4)
            sage: x = sum(J.gens())
            sage: x.norm()
            2

        """
        return self.inner_product(self).sqrt()


    def operator(self):
        """
        Return the left-multiplication-by-this-element
        operator on the ambient algebra.

        SETUP::

            sage: from mjo.eja.eja_algebra import random_eja

        TESTS::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x,y = J.random_elements(2)
            sage: x.operator()(y) == x*y
            True
            sage: y.operator()(x) == x*y
            True

        """
        P = self.parent()
        left_mult_by_self = lambda y: self*y
        L = P.module_morphism(function=left_mult_by_self, codomain=P)
        return FiniteDimensionalEJAOperator(P, P, L.matrix() )


    def quadratic_representation(self, other=None):
        """
        Return the quadratic representation of this element.

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  random_eja)

        EXAMPLES:

        The explicit form in the spin factor algebra is given by
        Alizadeh's Example 11.12::

            sage: set_random_seed()
            sage: x = JordanSpinEJA.random_instance().random_element()
            sage: x_vec = x.to_vector()
            sage: Q = matrix.identity(x.base_ring(), 0)
            sage: n = x_vec.degree()
            sage: if n > 0:
            ....:     x0 = x_vec[0]
            ....:     x_bar = x_vec[1:]
            ....:     A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
            ....:     B = 2*x0*x_bar.row()
            ....:     C = 2*x0*x_bar.column()
            ....:     D = matrix.identity(x.base_ring(), n-1)
            ....:     D = (x0^2 - x_bar.inner_product(x_bar))*D
            ....:     D = D + 2*x_bar.tensor_product(x_bar)
            ....:     Q = matrix.block(2,2,[A,B,C,D])
            sage: Q == x.quadratic_representation().matrix()
            True

        Test all of the properties from Theorem 11.2 in Alizadeh::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x,y = J.random_elements(2)
            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: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
            True

        Property 2 (multiply on the right for :trac:`28272`):

            sage: alpha = J.base_ring().random_element()
            sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
            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*Lx*Qex - Qx )
            True

            sage: 2*Lx*Qex - Qx == Lxx
            True

        Property 5:

            sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
            True

        Property 6:

            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

        """
        if other is None:
            other=self
        elif not other in self.parent():
            raise TypeError("'other' must live in the same algebra")

        L = self.operator()
        M = other.operator()
        return ( L*M + M*L - (self*other).operator() )



    def spectral_decomposition(self):
        """
        Return the unique spectral decomposition of this element.

        ALGORITHM:

        Following Faraut and Korányi's Theorem III.1.1, we restrict this
        element's left-multiplication-by operator to the subalgebra it
        generates. We then compute the spectral decomposition of that
        operator, and the spectral projectors we get back must be the
        left-multiplication-by operators for the idempotents we
        seek. Thus applying them to the identity element gives us those
        idempotents.

        Since the eigenvalues are required to be distinct, we take
        the spectral decomposition of the zero element to be zero
        times the identity element of the algebra (which is idempotent,
        obviously).

        SETUP::

            sage: from mjo.eja.eja_algebra import RealSymmetricEJA

        EXAMPLES:

        The spectral decomposition of the identity is ``1`` times itself,
        and the spectral decomposition of zero is ``0`` times the identity::

            sage: J = RealSymmetricEJA(3)
            sage: J.one()
            b0 + b2 + b5
            sage: J.one().spectral_decomposition()
            [(1, b0 + b2 + b5)]
            sage: J.zero().spectral_decomposition()
            [(0, b0 + b2 + b5)]

        TESTS::

            sage: J = RealSymmetricEJA(4)
            sage: x = sum(J.gens())
            sage: sd = x.spectral_decomposition()
            sage: l0 = sd[0][0]
            sage: l1 = sd[1][0]
            sage: c0 = sd[0][1]
            sage: c1 = sd[1][1]
            sage: c0.inner_product(c1) == 0
            True
            sage: c0.is_idempotent()
            True
            sage: c1.is_idempotent()
            True
            sage: c0 + c1 == J.one()
            True
            sage: l0*c0 + l1*c1 == x
            True

        The spectral decomposition should work in subalgebras, too::

            sage: J = RealSymmetricEJA(4)
            sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
            sage: A = 2*b5 - 2*b8
            sage: (lambda1, c1) = A.spectral_decomposition()[1]
            sage: (J0, J5, J1) = J.peirce_decomposition(c1)
            sage: (f0, f1, f2) = J1.gens()
            sage: f0.spectral_decomposition()
            [(0, 1.000000000000000?*c2), (1, 1.000000000000000?*c0)]

        """
        A = self.subalgebra_generated_by(orthonormalize=True)
        result = []
        for (evalue, proj) in A(self).operator().spectral_decomposition():
            result.append( (evalue, proj(A.one()).superalgebra_element()) )
        return result

    def subalgebra_generated_by(self, **kwargs):
        """
        Return the associative subalgebra of the parent EJA generated
        by this element.

        Since our parent algebra is unital, we want "subalgebra" to mean
        "unital subalgebra" as well; thus the subalgebra that an element
        generates will itself be a Euclidean Jordan algebra after
        restricting the algebra operations appropriately. This is the
        subalgebra that Faraut and Korányi work with in section II.2, for
        example.

        SETUP::

            sage: from mjo.eja.eja_algebra import (random_eja,
            ....:                                  HadamardEJA,
            ....:                                  RealSymmetricEJA)

        EXAMPLES:

        We can create subalgebras of Cartesian product EJAs that are not
        themselves Cartesian product EJAs (they're just "regular" EJAs)::

            sage: J1 = HadamardEJA(3)
            sage: J2 = RealSymmetricEJA(2)
            sage: J = cartesian_product([J1,J2])
            sage: J.one().subalgebra_generated_by()
            Euclidean Jordan algebra of dimension 1 over Algebraic Real Field

        TESTS:

        This subalgebra, being composed of only powers, is associative::

            sage: set_random_seed()
            sage: x0 = random_eja().random_element()
            sage: A = x0.subalgebra_generated_by()
            sage: x,y,z = A.random_elements(3)
            sage: (x*y)*z == x*(y*z)
            True

        Squaring in the subalgebra should work the same as in
        the superalgebra::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: A = x.subalgebra_generated_by()
            sage: A(x^2) == A(x)*A(x)
            True

        By definition, the subalgebra generated by the zero element is
        the one-dimensional algebra generated by the identity
        element... unless the original algebra was trivial, in which
        case the subalgebra is trivial too::

            sage: set_random_seed()
            sage: A = random_eja().zero().subalgebra_generated_by()
            sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
            True

        """
        powers = tuple( self**k for k in range(self.degree()) )
        A = self.parent().subalgebra(powers,
                                     associative=True,
                                     check_field=False,
                                     check_axioms=False,
                                     **kwargs)
        A.one.set_cache(A(self.parent().one()))
        return A


    def subalgebra_idempotent(self):
        """
        Find an idempotent in the associative subalgebra I generate
        using Proposition 2.3.5 in Baes.

        SETUP::

            sage: from mjo.eja.eja_algebra import random_eja

        TESTS:

        Ensure that we can find an idempotent in a non-trivial algebra
        where there are non-nilpotent elements, or that we get the dumb
        solution in the trivial algebra::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x = J.random_element()
            sage: while x.is_nilpotent() and not J.is_trivial():
            ....:     x = J.random_element()
            sage: c = x.subalgebra_idempotent()
            sage: c^2 == c
            True

        """
        if self.parent().is_trivial():
            return self

        if self.is_nilpotent():
            raise ValueError("this only works with non-nilpotent elements!")

        J = self.subalgebra_generated_by()
        u = J(self)

        # The image of the matrix of left-u^m-multiplication
        # will be minimal for some natural number s...
        s = 0
        minimal_dim = J.dimension()
        for i in range(1, minimal_dim):
            this_dim = (u**i).operator().matrix().image().dimension()
            if this_dim < minimal_dim:
                minimal_dim = this_dim
                s = i

        # Now minimal_matrix should correspond to the smallest
        # non-zero subspace in Baes's (or really, Koecher's)
        # proposition.
        #
        # However, we need to restrict the matrix to work on the
        # subspace... or do we? Can't we just solve, knowing that
        # A(c) = u^(s+1) should have a solution in the big space,
        # too?
        #
        # Beware, solve_right() means that we're using COLUMN vectors.
        # Our FiniteDimensionalAlgebraElement superclass uses rows.
        u_next = u**(s+1)
        A = u_next.operator().matrix()
        c = J.from_vector(A.solve_right(u_next.to_vector()))

        # Now c is the idempotent we want, but it still lives in the subalgebra.
        return c.superalgebra_element()


    def trace(self):
        """
        Return my trace, the sum of my eigenvalues.

        In a trivial algebra, however you want to look at it, the trace is
        an empty sum for which we declare the result to be zero.

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  HadamardEJA,
            ....:                                  TrivialEJA,
            ....:                                  random_eja)

        EXAMPLES::

            sage: J = TrivialEJA()
            sage: J.zero().trace()
            0

        ::
            sage: J = JordanSpinEJA(3)
            sage: x = sum(J.gens())
            sage: x.trace()
            2

        ::

            sage: J = HadamardEJA(5)
            sage: J.one().trace()
            5

        TESTS:

        The trace of an element is a real number::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: J.random_element().trace() in RLF
            True

        The trace is linear::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x,y = J.random_elements(2)
            sage: alpha = J.base_ring().random_element()
            sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
            True

        """
        P = self.parent()
        r = P.rank()

        if r == 0:
            # Special case for the trivial algebra where
            # the trace is an empty sum.
            return P.base_ring().zero()

        p = P._charpoly_coefficients()[r-1]
        # The _charpoly_coeff function already adds the factor of
        # -1 to ensure that _charpoly_coeff(r-1) is really what
        # appears in front of t^{r-1} in the charpoly. However,
        # we want the negative of THAT for the trace.
        return -p(*self.to_vector())


    def trace_inner_product(self, other):
        """
        Return the trace inner product of myself and ``other``.

        SETUP::

            sage: from mjo.eja.eja_algebra import random_eja

        TESTS:

        The trace inner product is commutative, bilinear, and associative::

            sage: set_random_seed()
            sage: J = random_eja()
            sage: x,y,z = J.random_elements(3)
            sage: # commutative
            sage: x.trace_inner_product(y) == y.trace_inner_product(x)
            True
            sage: # bilinear
            sage: a = J.base_ring().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
            sage: # associative
            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")

        return (self*other).trace()


    def trace_norm(self):
        """
        The norm of this element with respect to :meth:`trace_inner_product`.

        SETUP::

            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
            ....:                                  HadamardEJA)

        EXAMPLES::

            sage: J = HadamardEJA(2)
            sage: x = sum(J.gens())
            sage: x.trace_norm()
            1.414213562373095?

        ::

            sage: J = JordanSpinEJA(4)
            sage: x = sum(J.gens())
            sage: x.trace_norm()
            2.828427124746190?

        """
        return self.trace_inner_product(self).sqrt()


class CartesianProductEJAElement(FiniteDimensionalEJAElement):
    def det(self):
        r"""
        Compute the determinant of this product-element using the
        determianants of its factors.

        This result Follows from the spectral decomposition of (say)
        the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
        0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
        """
        from sage.misc.misc_c import prod
        return prod( f.det() for f in self.cartesian_factors() )

    def to_matrix(self):
        # An override is necessary to call our custom _scale().
        B = self.parent().matrix_basis()
        W = self.parent().matrix_space()

        # Aaaaand linear combinations don't work in Cartesian
        # product spaces, even though they provide a method with
        # that name. This is hidden behind an "if" because the
        # _scale() function is slow.
        pairs = zip(B, self.to_vector())
        return W.sum( _scale(b, alpha) for (b,alpha) in pairs )