from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
PolynomialRing,
QuadraticField)
-from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
-from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
+from mjo.eja.eja_element import FiniteDimensionalEJAElement
+from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
from mjo.eja.eja_utils import _mat2vec
-class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
+class FiniteDimensionalEJA(CombinatorialFreeModule):
+ r"""
+ A finite-dimensional Euclidean Jordan algebra.
+
+ INPUT:
+
+ - basis -- a tuple of basis elements in their matrix form.
+
+ - jordan_product -- function of two elements (in matrix form)
+ that returns their jordan product in this algebra; this will
+ be applied to ``basis`` to compute a multiplication table for
+ the algebra.
+
+ - inner_product -- function of two elements (in matrix form) that
+ returns their inner product. This will be applied to ``basis`` to
+ compute an inner-product table (basically a matrix) for this algebra.
+
+ """
+ Element = FiniteDimensionalEJAElement
+
+ def __init__(self,
+ basis,
+ jordan_product,
+ inner_product,
+ field=AA,
+ orthonormalize=True,
+ associative=False,
+ check_field=True,
+ check_axioms=True,
+ prefix='e'):
+
+ if check_field:
+ if not field.is_subring(RR):
+ # Note: this does return true for the real algebraic
+ # field, the rationals, and any quadratic field where
+ # we've specified a real embedding.
+ raise ValueError("scalar field is not real")
+
+ # If the basis given to us wasn't over the field that it's
+ # supposed to be over, fix that. Or, you know, crash.
+ basis = tuple( b.change_ring(field) for b in basis )
+
+ if check_axioms:
+ # Check commutativity of the Jordan and inner-products.
+ # This has to be done before we build the multiplication
+ # and inner-product tables/matrices, because we take
+ # advantage of symmetry in the process.
+ if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
+ for bi in basis
+ for bj in basis ):
+ raise ValueError("Jordan product is not commutative")
+
+ if not all( inner_product(bi,bj) == inner_product(bj,bi)
+ for bi in basis
+ for bj in basis ):
+ raise ValueError("inner-product is not commutative")
+
+
+ category = MagmaticAlgebras(field).FiniteDimensional()
+ category = category.WithBasis().Unital()
+ if associative:
+ # Element subalgebras can take advantage of this.
+ category = category.Associative()
+
+ # Call the superclass constructor so that we can use its from_vector()
+ # method to build our multiplication table.
+ n = len(basis)
+ super().__init__(field,
+ range(n),
+ prefix=prefix,
+ category=category,
+ bracket=False)
+
+ # Now comes all of the hard work. We'll be constructing an
+ # ambient vector space V that our (vectorized) basis lives in,
+ # as well as a subspace W of V spanned by those (vectorized)
+ # basis elements. The W-coordinates are the coefficients that
+ # we see in things like x = 1*e1 + 2*e2.
+ vector_basis = basis
+
+ from sage.structure.element import is_Matrix
+ basis_is_matrices = False
+
+ degree = 0
+ if n > 0:
+ if is_Matrix(basis[0]):
+ if basis[0].is_square():
+ # TODO: this ugly is_square() hack works around the problem
+ # of passing to_matrix()ed vectors in as the basis from a
+ # subalgebra. They aren't REALLY matrices, at least not of
+ # the type that we assume here... Ugh.
+ basis_is_matrices = True
+ from mjo.eja.eja_utils import _vec2mat
+ vector_basis = tuple( map(_mat2vec,basis) )
+ degree = basis[0].nrows()**2
+ else:
+ # convert from column matrices to vectors, yuck
+ basis = tuple( map(_mat2vec,basis) )
+ vector_basis = basis
+ degree = basis[0].degree()
+ else:
+ degree = basis[0].degree()
+
+ # Build an ambient space that fits...
+ V = VectorSpace(field, degree)
+
+ # We overwrite the name "vector_basis" in a second, but never modify it
+ # in place, to this effectively makes a copy of it.
+ deortho_vector_basis = vector_basis
+ self._deortho_matrix = None
+
+ if orthonormalize:
+ from mjo.eja.eja_utils import gram_schmidt
+ if basis_is_matrices:
+ vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y))
+ vector_basis = gram_schmidt(vector_basis, vector_ip)
+ else:
+ vector_basis = gram_schmidt(vector_basis, inner_product)
+
+ # Normalize the "matrix" basis, too!
+ basis = vector_basis
+
+ if basis_is_matrices:
+ basis = tuple( map(_vec2mat,basis) )
+
+ # Save the matrix "basis" for later... this is the last time we'll
+ # reference it in this constructor.
+ if basis_is_matrices:
+ self._matrix_basis = basis
+ else:
+ MS = MatrixSpace(self.base_ring(), degree, 1)
+ self._matrix_basis = tuple( MS(b) for b in basis )
+
+ # Now create the vector space for the algebra...
+ W = V.span_of_basis( vector_basis, check=check_axioms)
+
+ if orthonormalize:
+ # Now "W" is the vector space of our algebra coordinates. The
+ # variables "X1", "X2",... refer to the entries of vectors in
+ # W. Thus to convert back and forth between the orthonormal
+ # coordinates and the given ones, we need to stick the original
+ # basis in W.
+ U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
+ self._deortho_matrix = matrix( U.coordinate_vector(q)
+ for q in vector_basis )
+
+
+ # Now we actually compute the multiplication and inner-product
+ # tables/matrices using the possibly-orthonormalized basis.
+ self._inner_product_matrix = matrix.zero(field, n)
+ self._multiplication_table = [ [0 for j in range(i+1)] for i in range(n) ]
+
+ print("vector_basis:")
+ print(vector_basis)
+ # 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)
+
+ # TODO: the jordan product turns things back into
+ # matrices here even if they're supposed to be
+ # vectors. ugh. Can we get rid of vectors all together
+ # please?
+ elt = W.coordinate_vector(elt)
+ self._multiplication_table[i][j] = self.from_vector(elt)
+ self._inner_product_matrix[i,j] = ip
+ self._inner_product_matrix[j,i] = ip
+
+ self._inner_product_matrix._cache = {'hermitian': True}
+ self._inner_product_matrix.set_immutable()
+
+ if check_axioms:
+ if not self._is_jordanian():
+ raise ValueError("Jordan identity does not hold")
+ if not self._inner_product_is_associative():
+ raise ValueError("inner product is not associative")
+
def _coerce_map_from_base_ring(self):
"""
"""
return None
- def __init__(self,
- field,
- mult_table,
- prefix='e',
- category=None,
- matrix_basis=None,
- check_field=True,
- check_axioms=True):
+
+ def product_on_basis(self, i, j):
+ # We only stored the lower-triangular portion of the
+ # multiplication table.
+ if j <= i:
+ return self._multiplication_table[i][j]
+ else:
+ return self._multiplication_table[j][i]
+
+ def inner_product(self, x, y):
"""
+ The inner product associated with this Euclidean Jordan algebra.
+
+ Defaults to the trace inner product, but can be overridden by
+ subclasses if they are sure that the necessary properties are
+ satisfied.
+
SETUP::
- sage: from mjo.eja.eja_algebra import (
- ....: FiniteDimensionalEuclideanJordanAlgebra,
- ....: JordanSpinEJA,
- ....: random_eja)
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: HadamardEJA,
+ ....: BilinearFormEJA)
EXAMPLES:
- By definition, Jordan multiplication commutes::
+ Our inner product is "associative," which means the following for
+ a symmetric bilinear form::
sage: set_random_seed()
sage: J = random_eja()
- sage: x,y = J.random_elements(2)
- sage: x*y == y*x
+ sage: x,y,z = J.random_elements(3)
+ sage: (x*y).inner_product(z) == y.inner_product(x*z)
True
TESTS:
- The ``field`` we're given must be real with ``check_field=True``::
+ Ensure that this is the usual inner product for the algebras
+ over `R^n`::
- sage: JordanSpinEJA(2,QQbar)
- Traceback (most recent call last):
- ...
- ValueError: scalar field is not real
+ sage: set_random_seed()
+ sage: J = HadamardEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: actual = x.inner_product(y)
+ sage: expected = x.to_vector().inner_product(y.to_vector())
+ sage: actual == expected
+ True
- The multiplication table must be square with ``check_axioms=True``::
+ Ensure that this is one-half of the trace inner-product in a
+ BilinearFormEJA that isn't just the reals (when ``n`` isn't
+ one). This is in Faraut and Koranyi, and also my "On the
+ symmetry..." paper::
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
- Traceback (most recent call last):
- ...
- ValueError: multiplication table is not square
+ sage: set_random_seed()
+ sage: J = BilinearFormEJA.random_instance()
+ sage: n = J.dimension()
+ sage: x = J.random_element()
+ sage: y = J.random_element()
+ sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
+ True
+ """
+ B = self._inner_product_matrix
+ return (B*x.to_vector()).inner_product(y.to_vector())
+
+ def _is_commutative(self):
+ r"""
+ Whether or not this algebra's multiplication table is commutative.
+
+ This method should of course always return ``True``, unless
+ this algebra was constructed with ``check_axioms=False`` and
+ passed an invalid multiplication table.
"""
- if check_field:
- if not field.is_subring(RR):
- # Note: this does return true for the real algebraic
- # field, the rationals, and any quadratic field where
- # we've specified a real embedding.
- raise ValueError("scalar field is not real")
+ return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
+ for i in range(self.dimension())
+ for j in range(self.dimension()) )
- # The multiplication table had better be square
- n = len(mult_table)
- if check_axioms:
- if not all( len(l) == n for l in mult_table ):
- raise ValueError("multiplication table is not square")
-
- self._matrix_basis = matrix_basis
-
- if category is None:
- category = MagmaticAlgebras(field).FiniteDimensional()
- category = category.WithBasis().Unital()
-
- fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
- fda.__init__(field,
- range(n),
- prefix=prefix,
- category=category)
- self.print_options(bracket='')
-
- # The multiplication table we're given is necessarily in terms
- # of vectors, because we don't have an algebra yet for
- # anything to be an element of. However, it's faster in the
- # long run to have the multiplication table be in terms of
- # algebra elements. We do this after calling the superclass
- # constructor so that from_vector() knows what to do.
- self._multiplication_table = [
- list(map(lambda x: self.from_vector(x), ls))
- for ls in mult_table
- ]
+ def _is_jordanian(self):
+ r"""
+ Whether or not this algebra's multiplication table respects the
+ Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
- if check_axioms:
- if not self._is_commutative():
- raise ValueError("algebra is not commutative")
- if not self._is_jordanian():
- raise ValueError("Jordan identity does not hold")
- if not self._inner_product_is_associative():
- raise ValueError("inner product is not associative")
+ We only check one arrangement of `x` and `y`, so for a
+ ``True`` result to be truly true, you should also check
+ :meth:`_is_commutative`. This method should of course always
+ return ``True``, unless this algebra was constructed with
+ ``check_axioms=False`` and passed an invalid multiplication table.
+ """
+ return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
+ ==
+ (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
+ for i in range(self.dimension())
+ for j in range(self.dimension()) )
+
+ def _inner_product_is_associative(self):
+ r"""
+ Return whether or not this algebra's inner product `B` is
+ associative; that is, whether or not `B(xy,z) = B(x,yz)`.
+
+ This method should of course always return ``True``, unless
+ this algebra was constructed with ``check_axioms=False`` and
+ passed an invalid multiplication table.
+ """
+
+ # Used to check whether or not something is zero in an inexact
+ # ring. This number is sufficient to allow the construction of
+ # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
+ epsilon = 1e-16
+
+ for i in range(self.dimension()):
+ for j in range(self.dimension()):
+ for k in range(self.dimension()):
+ x = self.monomial(i)
+ y = self.monomial(j)
+ z = self.monomial(k)
+ diff = (x*y).inner_product(z) - x.inner_product(y*z)
+
+ if self.base_ring().is_exact():
+ if diff != 0:
+ return False
+ else:
+ if diff.abs() > epsilon:
+ return False
+
+ return True
def _element_constructor_(self, elt):
"""
sage: x = J.random_element()
sage: J(x.to_vector().column()) == x
True
+
"""
msg = "not an element of this algebra"
if elt == 0:
# that the integer 3 belongs to the space of 2-by-2 matrices.
raise ValueError(msg)
+ try:
+ elt = elt.column()
+ except (AttributeError, TypeError):
+ # Try to convert a vector into a column-matrix
+ pass
+
if elt not in self.matrix_space():
raise ValueError(msg)
# element's ring because the basis space might be an algebraic
# closure whereas the base ring of the 3-by-3 identity matrix
# could be QQ instead of QQbar.
+ #
+ # We pass check=False because the matrix basis is "guaranteed"
+ # to be linearly independent... right? Ha ha.
V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
- W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() )
+ W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()),
+ check=False)
try:
coords = W.coordinate_vector(_mat2vec(elt))
fmt = "Euclidean Jordan algebra of dimension {} over {}"
return fmt.format(self.dimension(), self.base_ring())
- def product_on_basis(self, i, j):
- return self._multiplication_table[i][j]
-
- def _is_commutative(self):
- r"""
- Whether or not this algebra's multiplication table is commutative.
-
- This method should of course always return ``True``, unless
- this algebra was constructed with ``check_axioms=False`` and
- passed an invalid multiplication table.
- """
- return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
- for i in range(self.dimension())
- for j in range(self.dimension()) )
-
- def _is_jordanian(self):
- r"""
- Whether or not this algebra's multiplication table respects the
- Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
-
- We only check one arrangement of `x` and `y`, so for a
- ``True`` result to be truly true, you should also check
- :meth:`_is_commutative`. This method should of course always
- return ``True``, unless this algebra was constructed with
- ``check_axioms=False`` and passed an invalid multiplication table.
- """
- return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
- ==
- (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
- for i in range(self.dimension())
- for j in range(self.dimension()) )
-
- def _inner_product_is_associative(self):
- r"""
- Return whether or not this algebra's inner product `B` is
- associative; that is, whether or not `B(xy,z) = B(x,yz)`.
-
- This method should of course always return ``True``, unless
- this algebra was constructed with ``check_axioms=False`` and
- passed an invalid multiplication table.
- """
-
- # Used to check whether or not something is zero in an inexact
- # ring. This number is sufficient to allow the construction of
- # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
- epsilon = 1e-16
-
- for i in range(self.dimension()):
- for j in range(self.dimension()):
- for k in range(self.dimension()):
- x = self.monomial(i)
- y = self.monomial(j)
- z = self.monomial(k)
- diff = (x*y).inner_product(z) - x.inner_product(y*z)
-
- if self.base_ring().is_exact():
- if diff != 0:
- return False
- else:
- if diff.abs() > epsilon:
- return False
-
- return True
@cached_method
def characteristic_polynomial_of(self):
sage: J = HadamardEJA(2)
sage: J.coordinate_polynomial_ring()
Multivariate Polynomial Ring in X1, X2...
- sage: J = RealSymmetricEJA(3,QQ)
+ sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
sage: J.coordinate_polynomial_ring()
Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
+----++----+----+----+----+
"""
- M = list(self._multiplication_table) # copy
- for i in range(len(M)):
- # M had better be "square"
- M[i] = [self.monomial(i)] + M[i]
- M = [["*"] + list(self.gens())] + M
+ n = self.dimension()
+ # Prepend the header row.
+ M = [["*"] + list(self.gens())]
+
+ # And to each subsequent row, prepend an entry that belongs to
+ # the left-side "header column."
+ M += [ [self.monomial(i)] + [ self.product_on_basis(i,j)
+ for j in range(n) ]
+ for i in range(n) ]
+
return table(M, header_row=True, header_column=True, frame=True)
[0], [1]
)
"""
- if self._matrix_basis is None:
- M = self.matrix_space()
- return tuple( M(b.to_vector()) for b in self.basis() )
- else:
- return self._matrix_basis
+ return self._matrix_basis
def matrix_space(self):
"""
if self.is_trivial():
return MatrixSpace(self.base_ring(), 0)
- elif self._matrix_basis is None or len(self._matrix_basis) == 0:
- return MatrixSpace(self.base_ring(), self.dimension(), 1)
else:
return self._matrix_basis[0].matrix_space()
if not c.is_idempotent():
raise ValueError("element is not idempotent: %s" % c)
- from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
+ from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
# 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 = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
+ trivial = FiniteDimensionalEJASubalgebra(self, ())
J0 = trivial # eigenvalue zero
J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
J1 = trivial # eigenvalue one
J5 = eigspace
else:
gens = tuple( self.from_vector(b) for b in eigspace.basis() )
- subalg = FiniteDimensionalEuclideanJordanSubalgebra(self,
- gens,
- check_axioms=False)
+ subalg = FiniteDimensionalEJASubalgebra(self,
+ gens,
+ check_axioms=False)
if eigval == 0:
J0 = subalg
elif eigval == 1:
Ensure that computing the rank actually works, since the ranks
of all simple algebras are known and will be cached by default::
- sage: J = HadamardEJA(4)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 4
-
- ::
-
- sage: J = JordanSpinEJA(4)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 2
-
- ::
-
- sage: J = RealSymmetricEJA(3)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 3
-
- ::
-
- sage: J = ComplexHermitianEJA(2)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 2
-
- ::
+ sage: set_random_seed() # long time
+ sage: J = random_eja() # long time
+ sage: caches = J.rank() # long time
+ sage: J.rank.clear_cache() # long time
+ sage: J.rank() == cached # long time
+ True
- sage: J = QuaternionHermitianEJA(2)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 2
"""
return len(self._charpoly_coefficients())
return self.zero().to_vector().parent().ambient_vector_space()
- Element = FiniteDimensionalEuclideanJordanAlgebraElement
+ Element = FiniteDimensionalEJAElement
-
-class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+class RationalBasisEJA(FiniteDimensionalEJA):
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.
+ New class for 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,
+ orthonormalize=True,
+ check_field=True,
+ check_axioms=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.
+ if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
+ raise TypeError("basis not rational")
+
+ 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,
+ orthonormalize=False,
+ check_field=False,
+ check_axioms=False,
+ **kwargs)
+
+ super().__init__(basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ check_field=check_field,
+ check_axioms=check_axioms,
+ **kwargs)
+
@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
+ sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+ ....: JordanSpinEJA)
EXAMPLES:
Algebraic Real Field
"""
- if self.base_ring() is QQ:
+ if self._rational_algebra is None:
# 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(
- map(lambda x: x.to_vector(), ls)
- for ls in self._multiplication_table
- )
+ # rationals if this one is already over the
+ # rationals. Likewise, if we never orthonormalized our
+ # basis, we might as well just use the given one.
+ return super()._charpoly_coefficients()
# 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))
+ # 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() )
+ # Now 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
-class ConcreteEuclideanJordanAlgebra:
+ 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(RationalBasisEJA):
r"""
A class for the Euclidean Jordan algebras that we know by name.
SETUP::
- sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import ConcreteEJA
TESTS:
product, unless we specify otherwise::
sage: set_random_seed()
- sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+ sage: J = ConcreteEJA.random_instance()
sage: all( b.norm() == 1 for b in J.gens() )
True
EJA the operator is self-adjoint by the Jordan axiom::
sage: set_random_seed()
- sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+ sage: J = ConcreteEJA.random_instance()
sage: x = J.random_element()
sage: x.operator().is_self_adjoint()
True
raise NotImplementedError
@classmethod
- def random_instance(cls, field=AA, **kwargs):
+ def random_instance(cls, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
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 = [[W.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
- try:
- # HACK PART DEUX
- self._inner_product_matrix = matrix(field,ip_table)
- except:
- pass
+ # These all bubble up to the RationalBasisEJA superclass
+ # constructor, so any (kw)args valid there are also valid
+ # here.
+ return eja_class.random_instance(*args, **kwargs)
- super(MatrixEuclideanJordanAlgebra, self).__init__(field,
- mult_table,
- matrix_basis=basis,
- **kwargs)
- if algebra_dim == 0:
- self.one.set_cache(self.zero())
- else:
- n = basis[0].nrows()
- # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
- # details of this algebra.
- self.one.set_cache(self(matrix.identity(field,n)))
-
-
- @cached_method
- def _charpoly_coefficients(self):
+class MatrixEJA:
+ @staticmethod
+ def dimension_over_reals():
r"""
- Override the parent method with something that tries to compute
- over a faster (non-extension) field.
- """
- if self._basis_normalizers is None or self.base_ring() is QQ:
- # We didn't normalize, or the basis we started with had
- # entries in a nice field already. Just compute the thing.
- return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
-
- basis = ( (b/n) for (b,n) in zip(self.matrix_basis(),
- self._basis_normalizers) )
-
- # Do this over the rationals and convert back at the end.
- # Only works because we know the entries of the basis are
- # integers. The argument ``check_axioms=False`` is required
- # because the trace inner-product method for this
- # class is a stub and can't actually be checked.
- J = MatrixEuclideanJordanAlgebra(QQ,
- basis,
- normalize_basis=False,
- check_field=False,
- check_axioms=False)
- a = J._charpoly_coefficients()
-
- # Unfortunately, changing the basis does change the
- # coefficients of the characteristic polynomial, but since
- # these are really the coefficients of the "characteristic
- # polynomial of" function, everything is still nice and
- # unevaluated. It's therefore "obvious" how scaling the
- # basis affects the coordinate variables X1, X2, et
- # cetera. Scaling the first basis vector up by "n" adds a
- # factor of 1/n into every "X1" term, for example. So here
- # we simply undo the basis_normalizer scaling that we
- # performed earlier.
- #
- # The a[0] access here is safe because trivial algebras
- # won't have any basis normalizers and therefore won't
- # make it to this "else" branch.
- XS = a[0].parent().gens()
- subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
- for i in range(len(XS)) }
- return tuple( a_i.subs(subs_dict) for a_i in a )
+ The dimension of this matrix's base ring over the reals.
+ The reals are dimension one over themselves, obviously; that's
+ just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
+ have dimension two. Finally, the quaternions have dimension
+ four over the reals.
- @staticmethod
- def real_embed(M):
+ This is used to determine the size of the matrix returned from
+ :meth:`real_embed`, among other things.
+ """
+ raise NotImplementedError
+
+ @classmethod
+ def real_embed(cls,M):
"""
Embed the matrix ``M`` into a space of real matrices.
real_embed(M*N) = real_embed(M)*real_embed(N)
"""
- raise NotImplementedError
+ if M.ncols() != M.nrows():
+ raise ValueError("the matrix 'M' must be square")
+ return M
- @staticmethod
- def real_unembed(M):
+ @classmethod
+ def real_unembed(cls,M):
"""
The inverse of :meth:`real_embed`.
"""
- raise NotImplementedError
+ if M.ncols() != M.nrows():
+ raise ValueError("the matrix 'M' must be square")
+ if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero():
+ raise ValueError("the matrix 'M' must be a real embedding")
+ return M
+
+ @staticmethod
+ def jordan_product(X,Y):
+ return (X*Y + Y*X)/2
@classmethod
- def matrix_inner_product(cls,X,Y):
+ def trace_inner_product(cls,X,Y):
+ r"""
+ Compute the trace inner-product of two real-embeddings.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
+ ....: ComplexHermitianEJA,
+ ....: QuaternionHermitianEJA)
+
+ EXAMPLES::
+
+ This gives the same answer as it would if we computed the trace
+ from the unembedded (original) matrices::
+
+ sage: set_random_seed()
+ sage: J = RealSymmetricEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: Xe = x.to_matrix()
+ sage: Ye = y.to_matrix()
+ sage: X = J.real_unembed(Xe)
+ sage: Y = J.real_unembed(Ye)
+ sage: expected = (X*Y).trace()
+ sage: actual = J.trace_inner_product(Xe,Ye)
+ sage: actual == expected
+ True
+
+ ::
+
+ sage: set_random_seed()
+ sage: J = ComplexHermitianEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: Xe = x.to_matrix()
+ sage: Ye = y.to_matrix()
+ sage: X = J.real_unembed(Xe)
+ sage: Y = J.real_unembed(Ye)
+ sage: expected = (X*Y).trace().real()
+ sage: actual = J.trace_inner_product(Xe,Ye)
+ sage: actual == expected
+ True
+
+ ::
+
+ sage: set_random_seed()
+ sage: J = QuaternionHermitianEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: Xe = x.to_matrix()
+ sage: Ye = y.to_matrix()
+ sage: X = J.real_unembed(Xe)
+ sage: Y = J.real_unembed(Ye)
+ sage: expected = (X*Y).trace().coefficient_tuple()[0]
+ sage: actual = J.trace_inner_product(Xe,Ye)
+ sage: actual == expected
+ True
+
+ """
Xu = cls.real_unembed(X)
Yu = cls.real_unembed(Y)
tr = (Xu*Yu).trace()
return tr.coefficient_tuple()[0]
-class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
+class RealMatrixEJA(MatrixEJA):
@staticmethod
- def real_embed(M):
- """
- The identity function, for embedding real matrices into real
- matrices.
- """
- return M
+ def dimension_over_reals():
+ return 1
- @staticmethod
- def real_unembed(M):
- """
- The identity function, for unembedding real matrices from real
- matrices.
- """
- return M
-
-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
"""
The rank-n simple EJA consisting of real symmetric n-by-n
matrices, the usual symmetric Jordan product, and the trace inner
In theory, our "field" can be any subfield of the reals::
- sage: RealSymmetricEJA(2, RDF)
+ sage: RealSymmetricEJA(2, field=RDF)
Euclidean Jordan algebra of dimension 3 over Real Double Field
- sage: RealSymmetricEJA(2, RR)
+ sage: RealSymmetricEJA(2, field=RR)
Euclidean Jordan algebra of dimension 3 over Real Field with
53 bits of precision
"""
@classmethod
- def _denormalized_basis(cls, n, field):
+ def _denormalized_basis(cls, n):
"""
Return a basis for the space of real symmetric n-by-n matrices.
sage: set_random_seed()
sage: n = ZZ.random_element(1,5)
- sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
+ sage: B = RealSymmetricEJA._denormalized_basis(n)
sage: all( M.is_symmetric() for M in B)
True
S = []
for i in range(n):
for j in range(i+1):
- Eij = matrix(field, n, lambda k,l: k==i and l==j)
+ Eij = matrix(ZZ, n, lambda k,l: k==i and l==j)
if i == j:
Sij = Eij
else:
Sij = Eij + Eij.transpose()
S.append(Sij)
- return S
+ return tuple(S)
@staticmethod
return 4 # Dimension 10
@classmethod
- def random_instance(cls, field=AA, **kwargs):
+ def random_instance(cls, **kwargs):
"""
Return a random instance of this type of algebra.
"""
n = ZZ.random_element(cls._max_random_instance_size() + 1)
- return cls(n, field, **kwargs)
+ return cls(n, **kwargs)
+
+ def __init__(self, n, **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
- def __init__(self, n, field=AA, **kwargs):
- basis = self._denormalized_basis(n, field)
- super(RealSymmetricEJA, self).__init__(field,
- basis,
- check_axioms=False,
+ super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n),
+ self.jordan_product,
+ self.trace_inner_product,
**kwargs)
+
+ # TODO: this could be factored out somehow, but is left here
+ # because the MatrixEJA is not presently a subclass of the
+ # FDEJA class that defines rank() and one().
self.rank.set_cache(n)
+ idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+ self.one.set_cache(self(idV))
-class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
+
+class ComplexMatrixEJA(MatrixEJA):
@staticmethod
- def real_embed(M):
+ def dimension_over_reals():
+ return 2
+
+ @classmethod
+ def real_embed(cls,M):
"""
Embed the n-by-n complex matrix ``M`` into the space of real
matrices of size 2n-by-2n via the map the sends each entry `z = a +
SETUP::
- sage: from mjo.eja.eja_algebra import \
- ....: ComplexMatrixEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
EXAMPLES::
sage: x3 = F(-i)
sage: x4 = F(6)
sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
- sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
+ sage: ComplexMatrixEJA.real_embed(M)
[ 4 -2| 1 2]
[ 2 4|-2 1]
[-----+-----]
sage: F = QuadraticField(-1, 'I')
sage: X = random_matrix(F, n)
sage: Y = random_matrix(F, n)
- sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
- sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
- sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
+ sage: Xe = ComplexMatrixEJA.real_embed(X)
+ sage: Ye = ComplexMatrixEJA.real_embed(Y)
+ sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
sage: Xe*Ye == XYe
True
"""
+ super(ComplexMatrixEJA,cls).real_embed(M)
n = M.nrows()
- if M.ncols() != n:
- raise ValueError("the matrix 'M' must be square")
# We don't need any adjoined elements...
field = M.base_ring().base_ring()
return matrix.block(field, n, blocks)
- @staticmethod
- def real_unembed(M):
+ @classmethod
+ def real_unembed(cls,M):
"""
The inverse of _embed_complex_matrix().
SETUP::
- sage: from mjo.eja.eja_algebra import \
- ....: ComplexMatrixEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
EXAMPLES::
....: [-2, 1, -4, 3],
....: [ 9, 10, 11, 12],
....: [-10, 9, -12, 11] ])
- sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
+ sage: ComplexMatrixEJA.real_unembed(A)
[ 2*I + 1 4*I + 3]
[ 10*I + 9 12*I + 11]
sage: set_random_seed()
sage: F = QuadraticField(-1, 'I')
sage: M = random_matrix(F, 3)
- sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
- sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
+ sage: Me = ComplexMatrixEJA.real_embed(M)
+ sage: ComplexMatrixEJA.real_unembed(Me) == M
True
"""
+ super(ComplexMatrixEJA,cls).real_unembed(M)
n = ZZ(M.nrows())
- if M.ncols() != n:
- raise ValueError("the matrix 'M' must be square")
- if not n.mod(2).is_zero():
- raise ValueError("the matrix 'M' must be a complex embedding")
+ d = cls.dimension_over_reals()
# If "M" was normalized, its base ring might have roots
# adjoined and they can stick around after unembedding.
field = M.base_ring()
R = PolynomialRing(field, 'z')
z = R.gen()
+
+ # Sage doesn't know how to adjoin the complex "i" (the root of
+ # x^2 + 1) to a field in a general way. Here, we just enumerate
+ # all of the cases that I have cared to support so far.
if field is AA:
# Sage doesn't know how to embed AA into QQbar, i.e. how
# to adjoin sqrt(-1) to AA.
F = QQbar
+ elif not field.is_exact():
+ # RDF or RR
+ F = field.complex_field()
else:
+ # Works for QQ and... maybe some other fields.
F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
i = F.gen()
# Go top-left to bottom-right (reading order), converting every
# 2-by-2 block we see to a single complex element.
elements = []
- for k in range(n/2):
- for j in range(n/2):
- submat = M[2*k:2*k+2,2*j:2*j+2]
+ for k in range(n/d):
+ for j in range(n/d):
+ submat = M[d*k:d*k+d,d*j:d*j+d]
if submat[0,0] != submat[1,1]:
raise ValueError('bad on-diagonal submatrix')
if submat[0,1] != -submat[1,0]:
z = submat[0,0] + submat[0,1]*i
elements.append(z)
- return matrix(F, n/2, elements)
-
-
- @classmethod
- def matrix_inner_product(cls,X,Y):
- """
- Compute a matrix inner product in this algebra directly from
- its real embedding.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+ return matrix(F, n/d, elements)
- TESTS:
- This gives the same answer as the slow, default method implemented
- in :class:`MatrixEuclideanJordanAlgebra`::
-
- sage: set_random_seed()
- sage: J = ComplexHermitianEJA.random_instance()
- sage: x,y = J.random_elements(2)
- sage: Xe = x.to_matrix()
- sage: Ye = y.to_matrix()
- sage: X = ComplexHermitianEJA.real_unembed(Xe)
- sage: Y = ComplexHermitianEJA.real_unembed(Ye)
- sage: expected = (X*Y).trace().real()
- sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
- sage: actual == expected
- True
-
- """
- return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2
-
-
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
"""
The rank-n simple EJA consisting of complex Hermitian n-by-n
matrices over the real numbers, the usual symmetric Jordan product,
In theory, our "field" can be any subfield of the reals::
- sage: ComplexHermitianEJA(2, RDF)
+ sage: ComplexHermitianEJA(2, field=RDF)
Euclidean Jordan algebra of dimension 4 over Real Double Field
- sage: ComplexHermitianEJA(2, RR)
+ sage: ComplexHermitianEJA(2, field=RR)
Euclidean Jordan algebra of dimension 4 over Real Field with
53 bits of precision
"""
@classmethod
- def _denormalized_basis(cls, n, field):
+ def _denormalized_basis(cls, n):
"""
Returns a basis for the space of complex Hermitian n-by-n matrices.
sage: set_random_seed()
sage: n = ZZ.random_element(1,5)
sage: field = QuadraticField(2, 'sqrt2')
- sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
+ sage: B = ComplexHermitianEJA._denormalized_basis(n)
sage: all( M.is_symmetric() for M in B)
True
"""
+ field = ZZ
R = PolynomialRing(field, 'z')
z = R.gen()
F = field.extension(z**2 + 1, 'I')
- I = F.gen()
+ I = F.gen(1)
# This is like the symmetric case, but we need to be careful:
#
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the complex extension "F".
- return ( s.change_ring(field) for s in S )
+ return tuple( s.change_ring(field) for s in S )
- def __init__(self, n, field=AA, **kwargs):
- basis = self._denormalized_basis(n,field)
- super(ComplexHermitianEJA,self).__init__(field,
- basis,
- check_axioms=False,
- **kwargs)
+ def __init__(self, n, **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(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n),
+ self.jordan_product,
+ self.trace_inner_product,
+ **kwargs)
+ # TODO: this could be factored out somehow, but is left here
+ # because the MatrixEJA is not presently a subclass of the
+ # FDEJA class that defines rank() and one().
self.rank.set_cache(n)
+ idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+ self.one.set_cache(self(idV))
@staticmethod
def _max_random_instance_size():
return 3 # Dimension 9
@classmethod
- def random_instance(cls, field=AA, **kwargs):
+ def random_instance(cls, **kwargs):
"""
Return a random instance of this type of algebra.
"""
n = ZZ.random_element(cls._max_random_instance_size() + 1)
- return cls(n, field, **kwargs)
+ return cls(n, **kwargs)
-class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
+class QuaternionMatrixEJA(MatrixEJA):
@staticmethod
- def real_embed(M):
+ def dimension_over_reals():
+ return 4
+
+ @classmethod
+ def real_embed(cls,M):
"""
Embed the n-by-n quaternion matrix ``M`` into the space of real
matrices of size 4n-by-4n by first sending each quaternion entry `z
SETUP::
- sage: from mjo.eja.eja_algebra import \
- ....: QuaternionMatrixEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
EXAMPLES::
sage: i,j,k = Q.gens()
sage: x = 1 + 2*i + 3*j + 4*k
sage: M = matrix(Q, 1, [[x]])
- sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
+ sage: QuaternionMatrixEJA.real_embed(M)
[ 1 2 3 4]
[-2 1 -4 3]
[-3 4 1 -2]
sage: Q = QuaternionAlgebra(QQ,-1,-1)
sage: X = random_matrix(Q, n)
sage: Y = random_matrix(Q, n)
- sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
- sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
- sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
+ sage: Xe = QuaternionMatrixEJA.real_embed(X)
+ sage: Ye = QuaternionMatrixEJA.real_embed(Y)
+ sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
sage: Xe*Ye == XYe
True
"""
+ super(QuaternionMatrixEJA,cls).real_embed(M)
quaternions = M.base_ring()
n = M.nrows()
- if M.ncols() != n:
- raise ValueError("the matrix 'M' must be square")
F = QuadraticField(-1, 'I')
i = F.gen()
d = t[3]
cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
[-c + d*i, a - b*i]])
- realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM)
+ realM = ComplexMatrixEJA.real_embed(cplxM)
blocks.append(realM)
# We should have real entries by now, so use the realest field
- @staticmethod
- def real_unembed(M):
+ @classmethod
+ def real_unembed(cls,M):
"""
The inverse of _embed_quaternion_matrix().
SETUP::
- sage: from mjo.eja.eja_algebra import \
- ....: QuaternionMatrixEuclideanJordanAlgebra
+ sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
EXAMPLES::
....: [-2, 1, -4, 3],
....: [-3, 4, 1, -2],
....: [-4, -3, 2, 1]])
- sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
+ sage: QuaternionMatrixEJA.real_unembed(M)
[1 + 2*i + 3*j + 4*k]
TESTS:
sage: set_random_seed()
sage: Q = QuaternionAlgebra(QQ, -1, -1)
sage: M = random_matrix(Q, 3)
- sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
- sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
+ sage: Me = QuaternionMatrixEJA.real_embed(M)
+ sage: QuaternionMatrixEJA.real_unembed(Me) == M
True
"""
+ super(QuaternionMatrixEJA,cls).real_unembed(M)
n = ZZ(M.nrows())
- if M.ncols() != n:
- raise ValueError("the matrix 'M' must be square")
- if not n.mod(4).is_zero():
- raise ValueError("the matrix 'M' must be a quaternion embedding")
+ d = cls.dimension_over_reals()
# Use the base ring of the matrix to ensure that its entries can be
# multiplied by elements of the quaternion algebra.
# 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
# quaternion block.
elements = []
- for l in range(n/4):
- for m in range(n/4):
- submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
- M[4*l:4*l+4,4*m:4*m+4] )
+ for l in range(n/d):
+ for m in range(n/d):
+ submat = ComplexMatrixEJA.real_unembed(
+ M[d*l:d*l+d,d*m:d*m+d] )
if submat[0,0] != submat[1,1].conjugate():
raise ValueError('bad on-diagonal submatrix')
if submat[0,1] != -submat[1,0].conjugate():
z += submat[0,1].imag()*k
elements.append(z)
- return matrix(Q, n/4, elements)
-
-
- @classmethod
- def matrix_inner_product(cls,X,Y):
- """
- Compute a matrix inner product in this algebra directly from
- its real embedding.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
-
- TESTS:
-
- This gives the same answer as the slow, default method implemented
- in :class:`MatrixEuclideanJordanAlgebra`::
-
- sage: set_random_seed()
- sage: J = QuaternionHermitianEJA.random_instance()
- sage: x,y = J.random_elements(2)
- sage: Xe = x.to_matrix()
- sage: Ye = y.to_matrix()
- sage: X = QuaternionHermitianEJA.real_unembed(Xe)
- sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
- sage: expected = (X*Y).trace().coefficient_tuple()[0]
- sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
- sage: actual == expected
- True
-
- """
- return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4
+ return matrix(Q, n/d, elements)
-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
r"""
The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
matrices, the usual symmetric Jordan product, and the
In theory, our "field" can be any subfield of the reals::
- sage: QuaternionHermitianEJA(2, RDF)
+ sage: QuaternionHermitianEJA(2, field=RDF)
Euclidean Jordan algebra of dimension 6 over Real Double Field
- sage: QuaternionHermitianEJA(2, RR)
+ sage: QuaternionHermitianEJA(2, field=RR)
Euclidean Jordan algebra of dimension 6 over Real Field with
53 bits of precision
"""
@classmethod
- def _denormalized_basis(cls, n, field):
+ def _denormalized_basis(cls, n):
"""
Returns a basis for the space of quaternion Hermitian n-by-n matrices.
sage: set_random_seed()
sage: n = ZZ.random_element(1,5)
- sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
+ sage: B = QuaternionHermitianEJA._denormalized_basis(n)
sage: all( M.is_symmetric() for M in B )
True
"""
+ field = ZZ
Q = QuaternionAlgebra(QQ,-1,-1)
I,J,K = Q.gens()
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the quaternion algebra "Q".
- return ( s.change_ring(field) for s in S )
+ return tuple( s.change_ring(field) for s in S )
- def __init__(self, n, field=AA, **kwargs):
- basis = self._denormalized_basis(n,field)
- super(QuaternionHermitianEJA,self).__init__(field,
- basis,
- check_axioms=False,
- **kwargs)
+ def __init__(self, n, **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(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n),
+ self.jordan_product,
+ self.trace_inner_product,
+ **kwargs)
+ # TODO: this could be factored out somehow, but is left here
+ # because the MatrixEJA is not presently a subclass of the
+ # FDEJA class that defines rank() and one().
self.rank.set_cache(n)
+ idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+ self.one.set_cache(self(idV))
+
@staticmethod
def _max_random_instance_size():
return 2 # Dimension 6
@classmethod
- def random_instance(cls, field=AA, **kwargs):
+ def random_instance(cls, **kwargs):
"""
Return a random instance of this type of algebra.
"""
n = ZZ.random_element(cls._max_random_instance_size() + 1)
- return cls(n, field, **kwargs)
+ return cls(n, **kwargs)
-class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class HadamardEJA(ConcreteEJA):
"""
Return the Euclidean Jordan Algebra corresponding to the set
`R^n` under the Hadamard product.
(r0, r1, r2)
"""
- def __init__(self, n, field=AA, **kwargs):
- V = VectorSpace(field, n)
- mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
- for i in range(n) ]
-
- # Inner products are real numbers and not algebra
- # elements, so once we turn the algebra element
- # into a vector in inner_product(), we never go
- # back. As a result -- contrary to what we do with
- # self._multiplication_table -- we store the inner
- # product table as a plain old matrix and not as
- # an algebra operator.
- ip_table = matrix.identity(field,n)
- self._inner_product_matrix = ip_table
-
- super(HadamardEJA, self).__init__(field,
- mult_table,
- check_axioms=False,
+ def __init__(self, n, **kwargs):
+ def jordan_product(x,y):
+ P = x.parent()
+ return P(tuple( xi*yi for (xi,yi) in zip(x,y) ))
+ def inner_product(x,y):
+ return x.inner_product(y)
+
+ # New defaults for keyword arguments. Don't orthonormalize
+ # because our basis is already orthonormal with respect to our
+ # inner-product. Don't check the axioms, because we know this
+ # is a valid EJA... but do double-check if the user passes
+ # check_axioms=True. Note: we DON'T override the "check_field"
+ # default here, because the user can pass in a field!
+ if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+
+ standard_basis = FreeModule(ZZ, n).basis()
+ super(HadamardEJA, self).__init__(standard_basis,
+ jordan_product,
+ inner_product,
**kwargs)
self.rank.set_cache(n)
return 5
@classmethod
- def random_instance(cls, field=AA, **kwargs):
+ def random_instance(cls, **kwargs):
"""
Return a random instance of this type of algebra.
"""
n = ZZ.random_element(cls._max_random_instance_size() + 1)
- return cls(n, field, **kwargs)
+ return cls(n, **kwargs)
-class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class BilinearFormEJA(ConcreteEJA):
r"""
The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
with the half-trace inner product and jordan product ``x*y =
We can check the multiplication condition given in the Jordan, von
Neumann, and Wigner paper (and also discussed on my "On the
symmetry..." paper). Note that this relies heavily on the standard
- choice of basis, as does anything utilizing the bilinear form matrix::
+ choice of basis, as does anything utilizing the bilinear form
+ matrix. We opt not to orthonormalize the basis, because if we
+ did, we would have to normalize the `s_{i}` in a similar manner::
sage: set_random_seed()
sage: n = ZZ.random_element(5)
sage: B22 = M.transpose()*M
sage: B = block_matrix(2,2,[ [B11,0 ],
....: [0, B22 ] ])
- sage: J = BilinearFormEJA(B)
+ sage: J = BilinearFormEJA(B, orthonormalize=False)
sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
sage: V = J.vector_space()
- sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
+ sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
....: for ei in eis ]
sage: actual = [ sis[i]*sis[j]
....: for i in range(n-1)
sage: actual == expected
True
"""
- def __init__(self, B, field=AA, **kwargs):
- n = B.nrows()
-
+ def __init__(self, B, **kwargs):
if not B.is_positive_definite():
raise ValueError("bilinear form is not positive-definite")
- V = VectorSpace(field, n)
- mult_table = [[V.zero() for j in range(n)] for i in range(n)]
- for i in range(n):
- for j in range(n):
- x = V.gen(i)
- y = V.gen(j)
- x0 = x[0]
- xbar = x[1:]
- y0 = y[0]
- ybar = y[1:]
- z0 = (B*x).inner_product(y)
- zbar = y0*xbar + x0*ybar
- z = V([z0] + zbar.list())
- mult_table[i][j] = z
-
- # Inner products are real numbers and not algebra
- # elements, so once we turn the algebra element
- # into a vector in inner_product(), we never go
- # back. As a result -- contrary to what we do with
- # self._multiplication_table -- we store the inner
- # product table as a plain old matrix and not as
- # an algebra operator.
- ip_table = B
- self._inner_product_matrix = ip_table
-
- super(BilinearFormEJA, self).__init__(field,
- mult_table,
- check_axioms=False,
+ def inner_product(x,y):
+ return (B*x).inner_product(y)
+
+ def jordan_product(x,y):
+ P = x.parent()
+ x0 = x[0]
+ xbar = x[1:]
+ y0 = y[0]
+ ybar = y[1:]
+ z0 = inner_product(x,y)
+ zbar = y0*xbar + x0*ybar
+ return P((z0,) + tuple(zbar))
+
+ # 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
+
+ n = B.nrows()
+ standard_basis = FreeModule(ZZ, n).basis()
+ super(BilinearFormEJA, self).__init__(standard_basis,
+ jordan_product,
+ inner_product,
**kwargs)
# The rank of this algebra is two, unless we're in a
return 5
@classmethod
- def random_instance(cls, field=AA, **kwargs):
+ def random_instance(cls, **kwargs):
"""
Return a random instance of this algebra.
"""
n = ZZ.random_element(cls._max_random_instance_size() + 1)
if n.is_zero():
- B = matrix.identity(field, n)
- return cls(B, field, **kwargs)
+ B = matrix.identity(ZZ, n)
+ return cls(B, **kwargs)
- B11 = matrix.identity(field,1)
- M = matrix.random(field, n-1)
- I = matrix.identity(field, n-1)
- alpha = field.zero()
+ B11 = matrix.identity(ZZ, 1)
+ M = matrix.random(ZZ, n-1)
+ I = matrix.identity(ZZ, n-1)
+ alpha = ZZ.zero()
while alpha.is_zero():
- alpha = field.random_element().abs()
+ alpha = ZZ.random_element().abs()
B22 = M.transpose()*M + alpha*I
from sage.matrix.special import block_matrix
B = block_matrix(2,2, [ [B11, ZZ(0) ],
[ZZ(0), B22 ] ])
- return cls(B, field, **kwargs)
+ return cls(B, **kwargs)
class JordanSpinEJA(BilinearFormEJA):
True
"""
- def __init__(self, n, field=AA, **kwargs):
- # This is a special case of the BilinearFormEJA with the identity
- # matrix as its bilinear form.
- B = matrix.identity(field, n)
- super(JordanSpinEJA, self).__init__(B, field, **kwargs)
+ def __init__(self, n, **kwargs):
+ # This is a special case of the BilinearFormEJA with the
+ # identity matrix as its bilinear form.
+ B = matrix.identity(ZZ, n)
+
+ # Don't orthonormalize because our basis is already
+ # orthonormal with respect to our inner-product.
+ if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
+
+ # But also don't pass check_field=False here, because the user
+ # can pass in a field!
+ super(JordanSpinEJA, self).__init__(B, **kwargs)
@staticmethod
def _max_random_instance_size():
return 5
@classmethod
- def random_instance(cls, field=AA, **kwargs):
+ def random_instance(cls, **kwargs):
"""
Return a random instance of this type of algebra.
Needed here to override the implementation for ``BilinearFormEJA``.
"""
n = ZZ.random_element(cls._max_random_instance_size() + 1)
- return cls(n, field, **kwargs)
+ return cls(n, **kwargs)
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class TrivialEJA(ConcreteEJA):
"""
The trivial Euclidean Jordan algebra consisting of only a zero element.
0
"""
- def __init__(self, field=AA, **kwargs):
- mult_table = []
- self._inner_product_matrix = matrix(field,0)
- super(TrivialEJA, self).__init__(field,
- mult_table,
- check_axioms=False,
+ def __init__(self, **kwargs):
+ jordan_product = lambda x,y: x
+ inner_product = lambda x,y: 0
+ basis = ()
+
+ # New defaults for keyword arguments
+ if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+ super(TrivialEJA, self).__init__(basis,
+ jordan_product,
+ inner_product,
**kwargs)
# The rank is zero using my definition, namely the dimension of the
# largest subalgebra generated by any element.
self.one.set_cache( self.zero() )
@classmethod
- def random_instance(cls, field=AA, **kwargs):
+ def random_instance(cls, **kwargs):
# We don't take a "size" argument so the superclass method is
# inappropriate for us.
- return cls(field, **kwargs)
+ return cls(**kwargs)
-class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
+class DirectSumEJA(ConcreteEJA):
r"""
The external (orthogonal) direct sum of two other Euclidean Jordan
algebras. Essentially the Cartesian product of its two factors.
have the same base ring; an error is raised otherwise::
sage: set_random_seed()
- sage: J1 = random_eja(AA)
- sage: J2 = random_eja(QQ)
+ sage: J1 = random_eja(field=AA)
+ sage: J2 = random_eja(field=QQ,orthonormalize=False)
sage: J = DirectSumEJA(J1,J2)
Traceback (most recent call last):
...
n2 = J2.dimension()
n = n1+n2
V = VectorSpace(field, n)
- mult_table = [ [ V.zero() for j in range(n) ]
+ mult_table = [ [ V.zero() for j in range(i+1) ]
for i in range(n) ]
for i in range(n1):
- for j in range(n1):
+ for j in range(i+1):
p = (J1.monomial(i)*J1.monomial(j)).to_vector()
mult_table[i][j] = V(p.list() + [field.zero()]*n2)
for i in range(n2):
- for j in range(n2):
+ for j in range(i+1):
p = (J2.monomial(i)*J2.monomial(j)).to_vector()
mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
+ # TODO: build the IP table here from the two constituent IP
+ # matrices (it'll be block diagonal, I think).
+ ip_table = [ [ field.zero() for j in range(i+1) ]
+ for i in range(n) ]
super(DirectSumEJA, self).__init__(field,
mult_table,
+ ip_table,
check_axioms=False,
**kwargs)
self.rank.set_cache(J1.rank() + J2.rank())
EXAMPLES::
- sage: J1 = HadamardEJA(2,QQ)
- sage: J2 = JordanSpinEJA(3,QQ)
+ sage: J1 = HadamardEJA(2, field=QQ)
+ sage: J2 = JordanSpinEJA(3, field=QQ)
sage: J = DirectSumEJA(J1,J2)
sage: J.factors()
(Euclidean Jordan algebra of dimension 2 over Rational Field,
# zero-by-two matrix (important for composing things).
P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
- pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1)
- pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2)
+ pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
+ pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
return (pi_left, pi_right)
def inclusions(self):
# two-by-zero matrix (important for composing things).
I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
- iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1)
- iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2)
+ iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
+ iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
return (iota_left, iota_right)
def inner_product(self, x, y):
EXAMPLE::
- sage: J1 = HadamardEJA(3,QQ)
- sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
+ sage: J1 = HadamardEJA(3,field=QQ)
+ sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
sage: J = DirectSumEJA(J1,J2)
sage: x1 = J1.one()
sage: x2 = x1
-random_eja = ConcreteEuclideanJordanAlgebra.random_instance
+random_eja = ConcreteEJA.random_instance
projects / sage.d.git / blobdiff