+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: actual = J.one().operator().matrix()
+ sage: expected = matrix.identity(J.base_ring(), J.dimension())
+ sage: actual == expected
+ True
+ sage: x = J.random_element()
+ sage: A = x.subalgebra_generated_by()
+ sage: actual = A.one().operator().matrix()
+ sage: expected = matrix.identity(A.base_ring(), A.dimension())
+ sage: actual == expected
+ True
+
+ ::
+
+ sage: set_random_seed()
+ sage: J = random_eja(field=QQ, orthonormalize=False)
+ sage: actual = J.one().operator().matrix()
+ sage: expected = matrix.identity(J.base_ring(), J.dimension())
+ sage: actual == expected
+ True
+ sage: x = J.random_element()
+ sage: A = x.subalgebra_generated_by(orthonormalize=False)
+ sage: actual = A.one().operator().matrix()
+ sage: expected = matrix.identity(A.base_ring(), A.dimension())
+ sage: actual == expected
+ True
+
+ Ensure that the cached unit element (often precomputed by
+ hand) agrees with the computed one::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: cached = J.one()
+ sage: J.one.clear_cache()
+ sage: J.one() == cached
+ True
+
+ ::
+
+ sage: set_random_seed()
+ sage: J = random_eja(field=QQ, orthonormalize=False)
+ sage: cached = J.one()
+ sage: J.one.clear_cache()
+ sage: J.one() == cached
+ True
+
+ """
+ # We can brute-force compute the matrices of the operators
+ # that correspond to the basis elements of this algebra.
+ # If some linear combination of those basis elements is the
+ # algebra identity, then the same linear combination of
+ # their matrices has to be the identity matrix.
+ #
+ # Of course, matrices aren't vectors in sage, so we have to
+ # appeal to the "long vectors" isometry.
+ oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
+
+ # Now we use basic linear algebra to find the coefficients,
+ # of the matrices-as-vectors-linear-combination, which should
+ # work for the original algebra basis too.
+ A = matrix(self.base_ring(), oper_vecs)
+
+ # We used the isometry on the left-hand side already, but we
+ # still need to do it for the right-hand side. Recall that we
+ # wanted something that summed to the identity matrix.
+ b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
+
+ # Now if there's an identity element in the algebra, this
+ # should work. We solve on the left to avoid having to
+ # transpose the matrix "A".
+ return self.from_vector(A.solve_left(b))
+
+
+ def peirce_decomposition(self, c):
+ """
+ The Peirce decomposition of this algebra relative to the
+ idempotent ``c``.
+
+ In the future, this can be extended to a complete system of
+ orthogonal idempotents.
+
+ INPUT:
+
+ - ``c`` -- an idempotent of this algebra.
+
+ OUTPUT:
+
+ A triple (J0, J5, J1) containing two subalgebras and one subspace
+ of this algebra,
+
+ - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
+ corresponding to the eigenvalue zero.
+
+ - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
+ corresponding to the eigenvalue one-half.
+
+ - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
+ corresponding to the eigenvalue one.
+
+ These are the only possible eigenspaces for that operator, and this
+ algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
+ orthogonal, and are subalgebras of this algebra with the appropriate
+ restrictions.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
+
+ EXAMPLES:
+
+ The canonical example comes from the symmetric matrices, which
+ decompose into diagonal and off-diagonal parts::
+
+ sage: J = RealSymmetricEJA(3)
+ sage: C = matrix(QQ, [ [1,0,0],
+ ....: [0,1,0],
+ ....: [0,0,0] ])
+ sage: c = J(C)
+ sage: J0,J5,J1 = J.peirce_decomposition(c)
+ sage: J0
+ Euclidean Jordan algebra of dimension 1...
+ sage: J5
+ Vector space of degree 6 and dimension 2...
+ sage: J1
+ Euclidean Jordan algebra of dimension 3...
+ sage: J0.one().to_matrix()
+ [0 0 0]
+ [0 0 0]
+ [0 0 1]
+ sage: orig_df = AA.options.display_format
+ sage: AA.options.display_format = 'radical'
+ sage: J.from_vector(J5.basis()[0]).to_matrix()
+ [ 0 0 1/2*sqrt(2)]
+ [ 0 0 0]
+ [1/2*sqrt(2) 0 0]
+ sage: J.from_vector(J5.basis()[1]).to_matrix()
+ [ 0 0 0]
+ [ 0 0 1/2*sqrt(2)]
+ [ 0 1/2*sqrt(2) 0]
+ sage: AA.options.display_format = orig_df
+ sage: J1.one().to_matrix()
+ [1 0 0]
+ [0 1 0]
+ [0 0 0]
+
+ TESTS:
+
+ Every algebra decomposes trivially with respect to its identity
+ element::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: J0,J5,J1 = J.peirce_decomposition(J.one())
+ sage: J0.dimension() == 0 and J5.dimension() == 0
+ True
+ sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
+ True
+
+ The decomposition is into eigenspaces, and its components are
+ therefore necessarily orthogonal. Moreover, the identity
+ elements in the two subalgebras are the projections onto their
+ respective subspaces of the superalgebra's identity element::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: x = J.random_element()
+ sage: if not J.is_trivial():
+ ....: while x.is_nilpotent():
+ ....: x = J.random_element()
+ sage: c = x.subalgebra_idempotent()
+ sage: J0,J5,J1 = J.peirce_decomposition(c)
+ sage: ipsum = 0
+ sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
+ ....: w = w.superalgebra_element()
+ ....: y = J.from_vector(y)
+ ....: z = z.superalgebra_element()
+ ....: ipsum += w.inner_product(y).abs()
+ ....: ipsum += w.inner_product(z).abs()
+ ....: ipsum += y.inner_product(z).abs()
+ sage: ipsum
+ 0
+ sage: J1(c) == J1.one()
+ True
+ sage: J0(J.one() - c) == J0.one()
+ True
+
+ """
+ if not c.is_idempotent():
+ raise ValueError("element is not idempotent: %s" % c)
+
+ # Default these to what they should be if they turn out to be
+ # trivial, because eigenspaces_left() won't return eigenvalues
+ # corresponding to trivial spaces (e.g. it returns only the
+ # eigenspace corresponding to lambda=1 if you take the
+ # decomposition relative to the identity element).
+ trivial = self.subalgebra((), check_axioms=False)
+ J0 = trivial # eigenvalue zero
+ J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
+ J1 = trivial # eigenvalue one
+
+ for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
+ if eigval == ~(self.base_ring()(2)):
+ J5 = eigspace
+ else:
+ gens = tuple( self.from_vector(b) for b in eigspace.basis() )
+ subalg = self.subalgebra(gens, check_axioms=False)
+ if eigval == 0:
+ J0 = subalg
+ elif eigval == 1:
+ J1 = subalg
+ else:
+ raise ValueError("unexpected eigenvalue: %s" % eigval)
+
+ return (J0, J5, J1)
+
+
+ def random_element(self, thorough=False):
+ r"""
+ Return a random element of this algebra.
+
+ Our algebra superclass method only returns a linear
+ combination of at most two basis elements. We instead
+ want the vector space "random element" method that
+ returns a more diverse selection.
+
+ INPUT:
+
+ - ``thorough`` -- (boolean; default False) whether or not we
+ should generate irrational coefficients for the random
+ element when our base ring is irrational; this slows the
+ algebra operations to a crawl, but any truly random method
+ should include them
+
+ """
+ # For a general base ring... maybe we can trust this to do the
+ # right thing? Unlikely, but.
+ V = self.vector_space()
+ v = V.random_element()
+
+ if self.base_ring() is AA:
+ # The "random element" method of the algebraic reals is
+ # stupid at the moment, and only returns integers between
+ # -2 and 2, inclusive:
+ #
+ # https://trac.sagemath.org/ticket/30875
+ #
+ # Instead, we implement our own "random vector" method,
+ # and then coerce that into the algebra. We use the vector
+ # space degree here instead of the dimension because a
+ # subalgebra could (for example) be spanned by only two
+ # vectors, each with five coordinates. We need to
+ # generate all five coordinates.
+ if thorough:
+ v *= QQbar.random_element().real()
+ else:
+ v *= QQ.random_element()
+
+ return self.from_vector(V.coordinate_vector(v))
+
+ def random_elements(self, count, thorough=False):
+ """
+ Return ``count`` random elements as a tuple.
+
+ INPUT:
+
+ - ``thorough`` -- (boolean; default False) whether or not we
+ should generate irrational coefficients for the random
+ elements when our base ring is irrational; this slows the
+ algebra operations to a crawl, but any truly random method
+ should include them
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+ EXAMPLES::
+
+ sage: J = JordanSpinEJA(3)
+ sage: x,y,z = J.random_elements(3)
+ sage: all( [ x in J, y in J, z in J ])
+ True
+ sage: len( J.random_elements(10) ) == 10
+ True
+
+ """
+ return tuple( self.random_element(thorough)
+ for idx in range(count) )
+
+
+ @cached_method
+ def _charpoly_coefficients(self):
+ r"""
+ The `r` polynomial coefficients of the "characteristic polynomial
+ of" function.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import random_eja
+
+ TESTS:
+
+ The theory shows that these are all homogeneous polynomials of
+ a known degree::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
+ True
+
+ """
+ n = self.dimension()
+ R = self.coordinate_polynomial_ring()
+ vars = R.gens()
+ F = R.fraction_field()
+
+ def L_x_i_j(i,j):
+ # From a result in my book, these are the entries of the
+ # basis representation of L_x.
+ return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
+ for k in range(n) )
+
+ L_x = matrix(F, n, n, L_x_i_j)
+
+ r = None
+ if self.rank.is_in_cache():
+ r = self.rank()
+ # There's no need to pad the system with redundant
+ # columns if we *know* they'll be redundant.
+ n = r
+
+ # Compute an extra power in case the rank is equal to
+ # the dimension (otherwise, we would stop at x^(r-1)).
+ x_powers = [ (L_x**k)*self.one().to_vector()
+ for k in range(n+1) ]
+ A = matrix.column(F, x_powers[:n])
+ AE = A.extended_echelon_form()
+ E = AE[:,n:]
+ A_rref = AE[:,:n]
+ if r is None:
+ r = A_rref.rank()
+ b = x_powers[r]
+
+ # The theory says that only the first "r" coefficients are
+ # nonzero, and they actually live in the original polynomial
+ # ring and not the fraction field. We negate them because in
+ # the actual characteristic polynomial, they get moved to the
+ # other side where x^r lives. We don't bother to trim A_rref
+ # down to a square matrix and solve the resulting system,
+ # because the upper-left r-by-r portion of A_rref is
+ # guaranteed to be the identity matrix, so e.g.
+ #
+ # A_rref.solve_right(Y)
+ #
+ # would just be returning Y.
+ return (-E*b)[:r].change_ring(R)
+
+ @cached_method
+ def rank(self):
+ r"""
+ Return the rank of this EJA.
+
+ This is a cached method because we know the rank a priori for
+ all of the algebras we can construct. Thus we can avoid the
+ expensive ``_charpoly_coefficients()`` call unless we truly
+ need to compute the whole characteristic polynomial.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (HadamardEJA,
+ ....: JordanSpinEJA,
+ ....: RealSymmetricEJA,
+ ....: ComplexHermitianEJA,
+ ....: QuaternionHermitianEJA,
+ ....: random_eja)
+
+ EXAMPLES:
+
+ The rank of the Jordan spin algebra is always two::
+
+ sage: JordanSpinEJA(2).rank()
+ 2
+ sage: JordanSpinEJA(3).rank()
+ 2
+ sage: JordanSpinEJA(4).rank()
+ 2
+
+ The rank of the `n`-by-`n` Hermitian real, complex, or
+ quaternion matrices is `n`::
+
+ sage: RealSymmetricEJA(4).rank()
+ 4
+ sage: ComplexHermitianEJA(3).rank()
+ 3
+ sage: QuaternionHermitianEJA(2).rank()
+ 2
+
+ TESTS:
+
+ Ensure that every EJA that we know how to construct has a
+ positive integer rank, unless the algebra is trivial in
+ which case its rank will be zero::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: r = J.rank()
+ sage: r in ZZ
+ True
+ sage: r > 0 or (r == 0 and J.is_trivial())
+ True
+
+ Ensure that computing the rank actually works, since the ranks
+ of all simple algebras are known and will be cached by default::
+
+ sage: set_random_seed() # long time
+ sage: J = random_eja() # long time
+ sage: cached = J.rank() # long time
+ sage: J.rank.clear_cache() # long time
+ sage: J.rank() == cached # long time
+ True
+
+ """
+ return len(self._charpoly_coefficients())
+
+
+ def subalgebra(self, basis, **kwargs):
+ r"""
+ Create a subalgebra of this algebra from the given basis.
+ """
+ from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+ return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+
+
+ def vector_space(self):
+ """
+ Return the vector space that underlies this algebra.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import RealSymmetricEJA
+
+ EXAMPLES::
+
+ sage: J = RealSymmetricEJA(2)
+ sage: J.vector_space()
+ Vector space of dimension 3 over...
+
+ """
+ return self.zero().to_vector().parent().ambient_vector_space()
+
+
+
+class RationalBasisEJA(FiniteDimensionalEJA):
+ r"""
+ Algebras whose supplied basis elements have all rational entries.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import BilinearFormEJA
+
+ EXAMPLES:
+
+ The supplied basis is orthonormalized by default::
+
+ sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
+ sage: J = BilinearFormEJA(B)
+ sage: J.matrix_basis()
+ (
+ [1] [ 0] [ 0]
+ [0] [1/5] [32/5]
+ [0], [ 0], [ 5]
+ )
+
+ """
+ def __init__(self,
+ basis,
+ jordan_product,
+ inner_product,
+ field=AA,
+ check_field=True,
+ **kwargs):
+
+ if check_field:
+ # Abuse the check_field parameter to check that the entries of
+ # out basis (in ambient coordinates) are in the field QQ.
+ # Use _all2list to get the vector coordinates of octonion
+ # entries and not the octonions themselves (which are not
+ # rational).
+ if not all( all(b_i in QQ for b_i in _all2list(b))
+ for b in basis ):
+ raise TypeError("basis not rational")
+
+ super().__init__(basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ check_field=check_field,
+ **kwargs)
+
+ self._rational_algebra = None
+ if field is not QQ:
+ # There's no point in constructing the extra algebra if this
+ # one is already rational.
+ #
+ # Note: the same Jordan and inner-products work here,
+ # because they are necessarily defined with respect to
+ # ambient coordinates and not any particular basis.
+ self._rational_algebra = FiniteDimensionalEJA(
+ basis,
+ jordan_product,
+ inner_product,
+ field=QQ,
+ matrix_space=self.matrix_space(),
+ associative=self.is_associative(),
+ orthonormalize=False,
+ check_field=False,
+ check_axioms=False)
+
+ def rational_algebra(self):
+ # Using None as a flag here (rather than just assigning "self"
+ # to self._rational_algebra by default) feels a little bit
+ # more sane to me in a garbage-collected environment.
+ if self._rational_algebra is None:
+ return self
+ else:
+ return self._rational_algebra
+
+ @cached_method
+ def _charpoly_coefficients(self):
+ r"""
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+ ....: JordanSpinEJA)
+
+ EXAMPLES:
+
+ The base ring of the resulting polynomial coefficients is what
+ it should be, and not the rationals (unless the algebra was
+ already over the rationals)::
+
+ sage: J = JordanSpinEJA(3)
+ sage: J._charpoly_coefficients()
+ (X1^2 - X2^2 - X3^2, -2*X1)
+ sage: a0 = J._charpoly_coefficients()[0]
+ sage: J.base_ring()
+ Algebraic Real Field
+ sage: a0.base_ring()
+ Algebraic Real Field
+
+ """
+ if self.rational_algebra() is self:
+ # Bypass the hijinks if they won't benefit us.
+ return super()._charpoly_coefficients()
+
+ # Do the computation over the rationals. The answer will be
+ # the same, because all we've done is a change of basis.
+ # Then, change back from QQ to our real base ring
+ a = ( a_i.change_ring(self.base_ring())
+ for a_i in self.rational_algebra()._charpoly_coefficients() )
+
+ # Otherwise, convert the coordinate variables back to the
+ # deorthonormalized ones.
+ R = self.coordinate_polynomial_ring()
+ from sage.modules.free_module_element import vector
+ X = vector(R, R.gens())
+ BX = self._deortho_matrix*X
+
+ subs_dict = { X[i]: BX[i] for i in range(len(X)) }
+ return tuple( a_i.subs(subs_dict) for a_i in a )
+
+class ConcreteEJA(FiniteDimensionalEJA):
+ r"""
+ A class for the Euclidean Jordan algebras that we know by name.
+
+ These are the Jordan algebras whose basis, multiplication table,
+ rank, and so on are known a priori. More to the point, they are
+ the Euclidean Jordan algebras for which we are able to conjure up
+ a "random instance."
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import ConcreteEJA
+
+ TESTS:
+
+ Our basis is normalized with respect to the algebra's inner
+ product, unless we specify otherwise::
+
+ sage: set_random_seed()
+ sage: J = ConcreteEJA.random_instance()
+ sage: all( b.norm() == 1 for b in J.gens() )
+ True
+
+ Since our basis is orthonormal with respect to the algebra's inner
+ product, and since we know that this algebra is an EJA, any
+ left-multiplication operator's matrix will be symmetric because
+ natural->EJA basis representation is an isometry and within the
+ EJA the operator is self-adjoint by the Jordan axiom::
+
+ sage: set_random_seed()
+ sage: J = ConcreteEJA.random_instance()
+ sage: x = J.random_element()
+ sage: x.operator().is_self_adjoint()
+ True
+ """
+
+ @staticmethod
+ def _max_random_instance_dimension():
+ r"""
+ The maximum dimension of any random instance. Ten dimensions seems
+ to be about the point where everything takes a turn for the
+ worse. And dimension ten (but not nine) allows the 4-by-4 real
+ Hermitian matrices, the 2-by-2 quaternion Hermitian matrices,
+ and the 2-by-2 octonion Hermitian matrices.
+ """
+ return 10
+
+ @staticmethod
+ def _max_random_instance_size(max_dimension):
+ """
+ Return an integer "size" that is an upper bound on the size of
+ this algebra when it is used in a random test case. This size
+ (which can be passed to the algebra's constructor) is itself
+ based on the ``max_dimension`` parameter.
+
+ This method must be implemented in each subclass.
+ """
+ raise NotImplementedError
+
+ @classmethod
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
+ """
+ Return a random instance of this type of algebra whose dimension
+ is less than or equal to the lesser of ``max_dimension`` and
+ the value returned by ``_max_random_instance_dimension()``. If
+ the dimension bound is omitted, then only the
+ ``_max_random_instance_dimension()`` is used as a bound.
+
+ This method should be implemented in each subclass.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import ConcreteEJA
+
+ TESTS:
+
+ Both the class bound and the ``max_dimension`` argument are upper
+ bounds on the dimension of the algebra returned::
+
+ sage: from sage.misc.prandom import choice
+ sage: eja_class = choice(ConcreteEJA.__subclasses__())
+ sage: class_max_d = eja_class._max_random_instance_dimension()
+ sage: J = eja_class.random_instance(max_dimension=20,
+ ....: field=QQ,
+ ....: orthonormalize=False)
+ sage: J.dimension() <= class_max_d
+ True
+ sage: J = eja_class.random_instance(max_dimension=2,
+ ....: field=QQ,
+ ....: orthonormalize=False)
+ sage: J.dimension() <= 2
+ True
+
+ """
+ from sage.misc.prandom import choice
+ eja_class = choice(cls.__subclasses__())
+
+ # These all bubble up to the RationalBasisEJA superclass
+ # constructor, so any (kw)args valid there are also valid
+ # here.
+ return eja_class.random_instance(max_dimension, *args, **kwargs)
+
+
+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
+
+ @staticmethod
+ def trace_inner_product(X,Y):
+ 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
+
+ """
+ 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()
+
+
+ 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
+ product. It has dimension `(n^2 + n)/2` over the reals.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import RealSymmetricEJA
+
+ EXAMPLES::
+
+ sage: J = RealSymmetricEJA(2)
+ sage: b0, b1, b2 = J.gens()
+ sage: b0*b0
+ b0
+ sage: b1*b1
+ 1/2*b0 + 1/2*b2
+ sage: b2*b2
+ b2
+
+ In theory, our "field" can be any subfield of the reals::
+
+ sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
+ Euclidean Jordan algebra of dimension 3 over Real Double Field
+ sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
+ Euclidean Jordan algebra of dimension 3 over Real Field with
+ 53 bits of precision
+
+ TESTS:
+
+ The dimension of this algebra is `(n^2 + n) / 2`::
+
+ sage: set_random_seed()
+ sage: d = RealSymmetricEJA._max_random_instance_dimension()
+ sage: n = RealSymmetricEJA._max_random_instance_size(d)
+ sage: J = RealSymmetricEJA(n)
+ sage: J.dimension() == (n^2 + n)/2
+ True
+
+ The Jordan multiplication is what we think it is::
+
+ sage: set_random_seed()
+ sage: J = RealSymmetricEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: actual = (x*y).to_matrix()
+ sage: X = x.to_matrix()
+ sage: Y = y.to_matrix()
+ sage: expected = (X*Y + Y*X)/2
+ sage: actual == expected
+ True
+ sage: J(expected) == x*y
+ True
+
+ We can change the generator prefix::
+
+ sage: RealSymmetricEJA(3, prefix='q').gens()
+ (q0, q1, q2, q3, q4, q5)
+
+ We can construct the (trivial) algebra of rank zero::
+
+ sage: RealSymmetricEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
+ """
+ @staticmethod
+ def _max_random_instance_size(max_dimension):
+ # Obtained by solving d = (n^2 + n)/2.
+ # The ZZ-int-ZZ thing is just "floor."
+ return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/2))
+
+ @classmethod
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ """
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
+ return cls(n, **kwargs)
+
+ 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
+
+ A = MatrixSpace(field, n)
+ super().__init__(A, **kwargs)
+
+ from mjo.eja.eja_cache import real_symmetric_eja_coeffs
+ a = real_symmetric_eja_coeffs(self)
+ if a is not None:
+ self.rational_algebra()._charpoly_coefficients.set_cache(a)
+
+
+
+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,
+ and the real-part-of-trace inner product. It has dimension `n^2` over
+ the reals.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+
+ EXAMPLES:
+
+ 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: 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:
+
+ The dimension of this algebra is `n^2`::
+
+ sage: set_random_seed()
+ sage: d = ComplexHermitianEJA._max_random_instance_dimension()
+ sage: n = ComplexHermitianEJA._max_random_instance_size(d)
+ sage: J = ComplexHermitianEJA(n)
+ sage: J.dimension() == n^2
+ True
+
+ The Jordan multiplication is what we think it is::
+
+ sage: set_random_seed()
+ sage: J = ComplexHermitianEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: actual = (x*y).to_matrix()
+ sage: X = x.to_matrix()
+ sage: Y = y.to_matrix()
+ sage: expected = (X*Y + Y*X)/2
+ sage: actual == expected
+ True
+ sage: J(expected) == x*y
+ True
+
+ We can change the generator prefix::
+
+ sage: ComplexHermitianEJA(2, prefix='z').gens()
+ (z0, z1, z2, z3)
+
+ We can construct the (trivial) algebra of rank zero::
+
+ sage: ComplexHermitianEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
+ """
+ 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
+
+ from mjo.hurwitz import ComplexMatrixAlgebra
+ A = ComplexMatrixAlgebra(n, scalars=field)
+ super().__init__(A, **kwargs)
+
+ from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
+ a = complex_hermitian_eja_coeffs(self)
+ if a is not None:
+ self.rational_algebra()._charpoly_coefficients.set_cache(a)
+
+ @staticmethod
+ def _max_random_instance_size(max_dimension):
+ # Obtained by solving d = n^2.
+ # The ZZ-int-ZZ thing is just "floor."
+ return ZZ(int(ZZ(max_dimension).sqrt()))
+
+ @classmethod
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ """
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
+ return cls(n, **kwargs)
+
+
+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
+ real-part-of-trace inner product. It has dimension `2n^2 - n` over
+ the reals.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
+
+ EXAMPLES:
+
+ In theory, our "field" can be any subfield of the reals::
+
+ sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
+ Euclidean Jordan algebra of dimension 6 over Real Double Field
+ sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
+ Euclidean Jordan algebra of dimension 6 over Real Field with
+ 53 bits of precision
+
+ TESTS:
+
+ The dimension of this algebra is `2*n^2 - n`::
+
+ sage: set_random_seed()
+ sage: d = QuaternionHermitianEJA._max_random_instance_dimension()
+ sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
+ sage: J = QuaternionHermitianEJA(n)
+ sage: J.dimension() == 2*(n^2) - n
+ True
+
+ The Jordan multiplication is what we think it is::
+
+ sage: set_random_seed()
+ sage: J = QuaternionHermitianEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: actual = (x*y).to_matrix()
+ sage: X = x.to_matrix()
+ sage: Y = y.to_matrix()
+ sage: expected = (X*Y + Y*X)/2
+ sage: actual == expected
+ True
+ sage: J(expected) == x*y
+ True
+
+ We can change the generator prefix::
+
+ sage: QuaternionHermitianEJA(2, prefix='a').gens()
+ (a0, a1, a2, a3, a4, a5)
+
+ We can construct the (trivial) algebra of rank zero::
+
+ sage: QuaternionHermitianEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
+ """
+ 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
+
+ from mjo.hurwitz import QuaternionMatrixAlgebra
+ A = QuaternionMatrixAlgebra(n, scalars=field)
+ super().__init__(A, **kwargs)
+
+ from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs
+ a = quaternion_hermitian_eja_coeffs(self)
+ if a is not None:
+ self.rational_algebra()._charpoly_coefficients.set_cache(a)
+
+
+
+ @staticmethod
+ def _max_random_instance_size(max_dimension):
+ r"""
+ The maximum rank of a random QuaternionHermitianEJA.
+ """
+ # Obtained by solving d = 2n^2 - n.
+ # The ZZ-int-ZZ thing is just "floor."
+ return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))
+
+ @classmethod
+ def random_instance(cls, max_dimension=None, *args, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ """
+ class_max_d = cls._max_random_instance_dimension()
+ if (max_dimension is None or max_dimension > class_max_d):
+ max_dimension = class_max_d
+ max_size = cls._max_random_instance_size(max_dimension)
+ n = ZZ.random_element(max_size + 1)
+ return cls(n, **kwargs)
+
+class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
+ r"""
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+ ....: OctonionHermitianEJA)
+ sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra