- from sage.structure.element import is_Matrix
- basis_is_matrices = False
-
- degree = 0
- if n > 0:
- if is_Matrix(basis[0]):
- basis_is_matrices = True
- vector_basis = tuple( map(_mat2vec,basis) )
- degree = basis[0].nrows()**2
- else:
- degree = basis[0].degree()
-
- V = VectorSpace(field, degree)
-
- # If we were asked to orthonormalize, and if the orthonormal
- # basis is different from the given one, then we also want to
- # compute multiplication and inner-product tables for the
- # deorthonormalized basis. These can be used later to
- # construct a deorthonormalized copy of this algebra over QQ
- # in which several operations are much faster.
- self._deortho_multiplication_table = None
- self._deortho_inner_product_table = None
-
- if orthonormalize:
- from mjo.eja.eja_utils import gram_schmidt
- vector_basis = gram_schmidt(vector_basis, inner_product)
- W = V.span_of_basis( vector_basis )
-
- # Normalize the "matrix" basis, too!
- basis = vector_basis
-
- if basis_is_matrices:
- from mjo.eja.eja_utils import _vec2mat
- basis = tuple( map(_vec2mat,basis) )
-
- W = V.span_of_basis( vector_basis )
-
- mult_table = [ [0 for i in range(n)] for j in range(n) ]
- ip_table = [ [0 for i in range(n)] for j in range(n) ]
-
- # Note: the Jordan and inner- products are defined in terms
- # of the ambient basis. It's important that their arguments
- # are in ambient coordinates as well.
- for i in range(n):
- for j in range(i+1):
- # ortho basis w.r.t. ambient coords
- q_i = vector_basis[i]
- q_j = vector_basis[j]
-
- if basis_is_matrices:
- q_i = _vec2mat(q_i)
- q_j = _vec2mat(q_j)
-
- elt = jordan_product(q_i, q_j)
- ip = inner_product(q_i, q_j)
-
- if basis_is_matrices:
- # do another mat2vec because the multiplication
- # table is in terms of vectors
- elt = _mat2vec(elt)
-
- elt = W.coordinate_vector(elt)
- mult_table[i][j] = elt
- mult_table[j][i] = elt
- ip_table[i][j] = ip
- ip_table[j][i] = ip
-
- if basis_is_matrices:
- for m in basis:
- m.set_immutable()
- else:
- basis = tuple( x.column() for x in basis )
-
- super().__init__(field,
- mult_table,
- ip_table,
- prefix,
- category,
- basis, # matrix basis
- check_field,
- check_axioms)
-
-class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
- r"""
- Algebras whose basis consists of vectors with rational
- entries. Equivalently, algebras whose multiplication tables
- contain only rational coefficients.
-
- When an EJA has a basis that can be made rational, we can speed up
- the computation of its characteristic polynomial by doing it over
- ``QQ``. All of the named EJA constructors that we provide fall
- into this category.
- """
- @cached_method
- def _charpoly_coefficients(self):
- r"""
- Override the parent method with something that tries to compute
- over a faster (non-extension) field.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import 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.base_ring() is QQ:
- # There's no need to construct *another* algebra over the
- # rationals if this one is already over the rationals.
- superclass = super(RationalBasisEuclideanJordanAlgebra, self)
- return superclass._charpoly_coefficients()
-
- mult_table = tuple(
- tuple(map(lambda x: x.to_vector(), ls))
- for ls in self._multiplication_table
- )
-
- # Do the computation over the rationals. The answer will be
- # the same, because our basis coordinates are (essentially)
- # rational.
- J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
- mult_table,
- check_field=False,
- check_axioms=False)
- a = J._charpoly_coefficients()
- return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
-
-
-class ConcreteEuclideanJordanAlgebra:
- 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 ConcreteEuclideanJordanAlgebra
-
- TESTS:
-
- Our basis is normalized with respect to the algebra's inner
- product, unless we specify otherwise::
-
- sage: set_random_seed()
- sage: J = ConcreteEuclideanJordanAlgebra.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 = ConcreteEuclideanJordanAlgebra.random_instance()
- sage: x = J.random_element()
- sage: x.operator().is_self_adjoint()
- True
- """
-
- @staticmethod
- def _max_random_instance_size():
- """
- Return an integer "size" that is an upper bound on the size of
- this algebra when it is used in a random test
- case. Unfortunately, the term "size" is ambiguous -- when
- dealing with `R^n` under either the Hadamard or Jordan spin
- product, the "size" refers to the dimension `n`. When dealing
- with a matrix algebra (real symmetric or complex/quaternion
- Hermitian), it refers to the size of the matrix, which is far
- less than the dimension of the underlying vector space.
-
- This method must be implemented in each subclass.
- """
- raise NotImplementedError
-
- @classmethod
- def random_instance(cls, field=AA, **kwargs):
- """
- Return a random instance of this type of algebra.
-
- This method should be implemented in each subclass.
- """
- from sage.misc.prandom import choice
- eja_class = choice(cls.__subclasses__())
- return eja_class.random_instance(field)
-
-
-class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
-
- def __init__(self, field, basis, normalize_basis=True, **kwargs):
- """
- Compared to the superclass constructor, we take a basis instead of
- a multiplication table because the latter can be computed in terms
- of the former when the product is known (like it is here).
- """
- # Used in this class's fast _charpoly_coefficients() override.
- self._basis_normalizers = None
-
- # We're going to loop through this a few times, so now's a good
- # time to ensure that it isn't a generator expression.
- basis = tuple(basis)
-
- algebra_dim = len(basis)
- degree = 0 # size of the matrices
- if algebra_dim > 0:
- degree = basis[0].nrows()
-
- if algebra_dim > 1 and normalize_basis:
- # We'll need sqrt(2) to normalize the basis, and this
- # winds up in the multiplication table, so the whole
- # algebra needs to be over the field extension.
- R = PolynomialRing(field, 'z')
- z = R.gen()
- p = z**2 - 2
- if p.is_irreducible():
- field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
- basis = tuple( s.change_ring(field) for s in basis )
- self._basis_normalizers = tuple(
- ~(self.matrix_inner_product(s,s).sqrt()) for s in basis )
- basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
-
- # Now compute the multiplication and inner product tables.
- # We have to do this *after* normalizing the basis, because
- # scaling affects the answers.
- V = VectorSpace(field, degree**2)
- W = V.span_of_basis( _mat2vec(s) for s in basis )
- mult_table = [[W.zero() for j in range(algebra_dim)]
- for i in range(algebra_dim)]
- ip_table = [[field.zero() for j in range(algebra_dim)]
- for i in range(algebra_dim)]
- for i in range(algebra_dim):
- for j in range(algebra_dim):
- mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
- mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
-
- try:
- # HACK: ignore the error here if we don't need the
- # inner product (as is the case when we construct
- # a dummy QQ-algebra for fast charpoly coefficients.
- ip_table[i][j] = self.matrix_inner_product(basis[i],
- basis[j])
- except:
- pass
-
- super(MatrixEuclideanJordanAlgebra, self).__init__(field,
- mult_table,
- ip_table,
- matrix_basis=basis,
- **kwargs)
-
- if algebra_dim == 0:
- self.one.set_cache(self.zero())