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"""
- The lowest-level class for representing a Euclidean Jordan algebra.
+ 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
+
+ degree = 0
+ if n > 0:
+ # Works on both column and square matrices...
+ degree = len(basis[0].list())
+
+ # Build an ambient space that fits our matrix basis when
+ # written out as "long vectors."
+ V = VectorSpace(field, degree)
+
+ # The matrix that will hole the orthonormal -> unorthonormal
+ # coordinate transformation.
+ self._deortho_matrix = None
+
+ if orthonormalize:
+ # Save a copy of the un-orthonormalized basis for later.
+ # Convert it to ambient V (vector) coordinates while we're
+ # at it, because we'd have to do it later anyway.
+ deortho_vector_basis = tuple( V(b.list()) for b in basis )
+
+ from mjo.eja.eja_utils import gram_schmidt
+ basis = tuple(gram_schmidt(basis, inner_product))
+
+ # Save the (possibly orthonormalized) matrix basis for
+ # later...
+ self._matrix_basis = basis
+
+ # Now create the vector space for the algebra, which will have
+ # its own set of non-ambient coordinates (in terms of the
+ # supplied basis).
+ vector_basis = tuple( V(b.list()) for b in basis )
+ 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.identity(field, n)
+ self._multiplication_table = [ [0 for j in range(i+1)]
+ for i 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 = basis[i]
+ q_j = basis[j]
+
+ # The jordan product returns a matrixy answer, so we
+ # have to convert it to the algebra coordinates.
+ elt = jordan_product(q_i, q_j)
+ elt = W.coordinate_vector(V(elt.list()))
+ self._multiplication_table[i][j] = self.from_vector(elt)
+
+ if not orthonormalize:
+ # If we're orthonormalizing the basis with respect
+ # to an inner-product, then the inner-product
+ # matrix with respect to the resulting basis is
+ # just going to be the identity.
+ ip = inner_product(q_i, q_j)
+ 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):
"""
Disable the map from the base ring into the algebra.
"""
return None
- def __init__(self,
- field,
- multiplication_table,
- inner_product_table,
- prefix='e',
- category=None,
- matrix_basis=None,
- check_field=True,
- check_axioms=True):
- """
- INPUT:
- * field -- the scalar field for this algebra (must be real)
+ 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]
- * multiplication_table -- the multiplication table for this
- algebra's implicit basis. Only the lower-triangular portion
- of the table is used, since the multiplication is assumed
- to be commutative.
+ 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
- An error is raised if the Jordan product is not commutative::
-
- sage: JP = ((1,2),(0,0))
- sage: IP = ((1,0),(0,1))
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
- Traceback (most recent call last):
- ...
- ValueError: Jordan product is not commutative
-
- An error is raised if the inner-product is not commutative::
+ TESTS:
- sage: JP = ((1,0),(0,1))
- sage: IP = ((1,2),(0,0))
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
- Traceback (most recent call last):
- ...
- ValueError: inner-product is not commutative
+ Ensure that this is the usual inner product for the algebras
+ over `R^n`::
- TESTS:
+ 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 ``field`` we're given must be real with ``check_field=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: JordanSpinEJA(2,QQbar)
- Traceback (most recent call last):
- ...
- ValueError: scalar field is not real
+ 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())
- The multiplication table must be square with ``check_axioms=True``::
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
- Traceback (most recent call last):
- ...
- ValueError: multiplication table is not square
+ def _is_commutative(self):
+ r"""
+ Whether or not this algebra's multiplication table is commutative.
- The multiplication and inner-product tables must be the same
- size (and in particular, the inner-product table must also be
- square) with ``check_axioms=True``::
+ 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()) )
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
- Traceback (most recent call last):
- ...
- ValueError: multiplication and inner-product tables are
- different sizes
- sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
- Traceback (most recent call last):
- ...
- ValueError: multiplication and inner-product tables are
- different sizes
+ 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.
"""
- 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.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)`.
- # The multiplication and inner-product tables should be square
- # if the user wants us to verify them. And we verify them as
- # soon as possible, because we want to exploit their symmetry.
- n = len(multiplication_table)
- if check_axioms:
- if not all( len(l) == n for l in multiplication_table ):
- raise ValueError("multiplication table is not square")
-
- # If the multiplication table is square, we can check if
- # the inner-product table is square by comparing it to the
- # multiplication table's dimensions.
- msg = "multiplication and inner-product tables are different sizes"
- if not len(inner_product_table) == n:
- raise ValueError(msg)
-
- if not all( len(l) == n for l in inner_product_table ):
- raise ValueError(msg)
-
- # Check commutativity of the Jordan product (symmetry of
- # the multiplication table) and the commutativity of the
- # inner-product (symmetry of the inner-product table)
- # first if we're going to check them at all.. This has to
- # be done before we define product_on_basis(), because
- # that method assumes that self._multiplication_table is
- # symmetric. And it has to be done before we build
- # self._inner_product_matrix, because the process used to
- # construct it assumes symmetry as well.
- if not all( multiplication_table[j][i]
- == multiplication_table[i][j]
- for i in range(n)
- for j in range(i+1) ):
- raise ValueError("Jordan product is not commutative")
+ This method should of course always return ``True``, unless
+ this algebra was constructed with ``check_axioms=False`` and
+ passed an invalid Jordan or inner-product.
+ """
- if not all( inner_product_table[j][i]
- == inner_product_table[i][j]
- for i in range(n)
- for j in range(i+1) ):
- raise ValueError("inner-product is not commutative")
+ # 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
- 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.
- #
- # Note: we take advantage of symmetry here, and only store
- # the lower-triangular portion of the table.
- self._multiplication_table = [ [ self.vector_space().zero()
- for j in range(i+1) ]
- for i in range(n) ]
+ 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)
- for i in range(n):
- for j in range(i+1):
- elt = self.from_vector(multiplication_table[i][j])
- self._multiplication_table[i][j] = elt
-
- self._multiplication_table = tuple(map(tuple, self._multiplication_table))
-
- # Save our inner product as a matrix, since the efficiency of
- # matrix multiplication will usually outweigh the fact that we
- # have to store a redundant upper- or lower-triangular part.
- # Pre-cache the fact that these are Hermitian (real symmetric,
- # in fact) in case some e.g. matrix multiplication routine can
- # take advantage of it.
- ip_matrix_constructor = lambda i,j: inner_product_table[i][j] if j <= i else inner_product_table[j][i]
- self._inner_product_matrix = matrix(field, n, ip_matrix_constructor)
- self._inner_product_matrix._cache = {'hermitian': True}
- self._inner_product_matrix.set_immutable()
+ if self.base_ring().is_exact():
+ if diff != 0:
+ return False
+ else:
+ if diff.abs() > epsilon:
+ return False
- 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")
+ 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):
- # 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 _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,orthonormalize=False)
+ sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
sage: J.coordinate_polynomial_ring()
Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
"""
n = self.dimension()
- M = [ [ self.zero() for j in range(n) ]
- for i in range(n) ]
- for i in range(n):
- for j in range(i+1):
- M[i][j] = self._multiplication_table[i][j]
- M[j][i] = M[i][j]
+ # Prepend the header row.
+ M = [["*"] + list(self.gens())]
- for i in range(n):
- # Prepend the left "header" column entry Can't do this in
- # the loop because it messes up the symmetry.
- M[i] = [self.monomial(i)] + M[i]
+ # 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) ]
- # Prepend the header row.
- M = [["*"] + list(self.gens())] + M
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()
+ return self.matrix_basis()[0].parent()
@cached_method
sage: from mjo.eja.eja_algebra import (HadamardEJA,
....: random_eja)
- EXAMPLES::
+ EXAMPLES:
+
+ We can compute unit element in the Hadamard EJA::
+
+ sage: J = HadamardEJA(5)
+ sage: J.one()
+ e0 + e1 + e2 + e3 + e4
+
+ The unit element in the Hadamard EJA is inherited in the
+ subalgebras generated by its elements::
sage: J = HadamardEJA(5)
sage: J.one()
e0 + e1 + e2 + e3 + e4
+ sage: x = sum(J.gens())
+ sage: A = x.subalgebra_generated_by(orthonormalize=False)
+ sage: A.one()
+ f0
+ sage: A.one().superalgebra_element()
+ e0 + e1 + e2 + e3 + e4
TESTS:
- The identity element acts like the identity::
+ The identity element acts like the identity, regardless of
+ whether or not we orthonormalize::
sage: set_random_seed()
sage: J = random_eja()
sage: x = J.random_element()
sage: J.one()*x == x and x*J.one() == x
True
+ sage: A = x.subalgebra_generated_by()
+ sage: y = A.random_element()
+ sage: A.one()*y == y and y*A.one() == y
+ True
- The matrix of the unit element's operator is the identity::
+ ::
+
+ sage: set_random_seed()
+ sage: J = random_eja(field=QQ, orthonormalize=False)
+ sage: x = J.random_element()
+ sage: J.one()*x == x and x*J.one() == x
+ True
+ sage: A = x.subalgebra_generated_by(orthonormalize=False)
+ sage: y = A.random_element()
+ sage: A.one()*y == y and y*A.one() == y
+ True
+
+ The matrix of the unit element's operator is the identity,
+ regardless of the base field and whether or not we
+ orthonormalize::
sage: set_random_seed()
sage: J = random_eja()
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: 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 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:
return self.zero().to_vector().parent().ambient_vector_space()
- Element = FiniteDimensionalEuclideanJordanAlgebraElement
+ Element = FiniteDimensionalEJAElement
-class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+class RationalBasisEJA(FiniteDimensionalEJA):
r"""
New class for algebras whose supplied basis elements have all rational entries.
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,
- field,
- basis,
- jordan_product,
- inner_product,
- orthonormalize=True,
- prefix='e',
- category=None,
- check_field=True,
- check_axioms=True):
-
- n = len(basis)
- vector_basis = basis
-
- 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
- from mjo.eja.eja_utils import _vec2mat
- 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:
- # Compute the deorthonormalized tables before we orthonormalize
- # the given basis.
- W = V.span_of_basis( vector_basis )
-
- if check_axioms:
- # If the superclass constructor is going to verify the
- # symmetry of this table, it has better at least be
- # square...
- self._deortho_multiplication_table = [ [0 for j in range(n)]
- for i in range(n) ]
- self._deortho_inner_product_table = [ [0 for j in range(n)]
- for i in range(n) ]
- else:
- self._deortho_multiplication_table = [ [0 for j in range(i+1)]
- for i in range(n) ]
- self._deortho_inner_product_table = [ [0 for j in range(i+1)]
- for i 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):
- # given 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)
- self._deortho_multiplication_table[i][j] = elt
- self._deortho_inner_product_table[i][j] = ip
- if check_axioms:
- # The tables are square if we're verifying that they
- # are commutative.
- self._deortho_multiplication_table[j][i] = elt
- self._deortho_inner_product_table[j][i] = ip
-
- if self._deortho_multiplication_table is not None:
- self._deortho_multiplication_table = tuple(map(tuple, self._deortho_multiplication_table))
- if self._deortho_inner_product_table is not None:
- self._deortho_inner_product_table = tuple(map(tuple, self._deortho_inner_product_table))
-
- # 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)
-
- W = V.span_of_basis( vector_basis )
-
- # Normalize the "matrix" basis, too!
- basis = vector_basis
-
- if basis_is_matrices:
- basis = tuple( map(_vec2mat,basis) )
-
- W = V.span_of_basis( vector_basis )
-
- # 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 )
- self._deortho_matrix = matrix( U.coordinate_vector(q)
- for q in vector_basis )
-
- # If the superclass constructor is going to verify the
- # symmetry of this table, it has better at least be
- # square...
- if check_axioms:
- mult_table = [ [0 for j in range(n)] for i in range(n) ]
- ip_table = [ [0 for j in range(n)] for i in range(n) ]
- else:
- mult_table = [ [0 for j in range(i+1)] for i in range(n) ]
- ip_table = [ [0 for j in range(i+1)] for i 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)
+ sage: J = BilinearFormEJA(B)
+ sage: J.matrix_basis()
+ (
+ [1] [ 0] [ 0]
+ [0] [1/5] [32/5]
+ [0], [ 0], [ 5]
+ )
- 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
- ip_table[i][j] = ip
- if check_axioms:
- # The tables are square if we're verifying that they
- # are commutative.
- mult_table[j][i] = elt
- 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 )
+ """
+ def __init__(self,
+ basis,
+ jordan_product,
+ inner_product,
+ field=AA,
+ check_field=True,
+ **kwargs):
- super().__init__(field,
- mult_table,
- ip_table,
- prefix,
- category,
- basis, # matrix basis
- check_field,
- check_axioms)
+ 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")
+
+ 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,
+ orthonormalize=False,
+ check_field=False,
+ check_axioms=False)
+
+ super().__init__(basis,
+ jordan_product,
+ inner_product,
+ field=field,
+ check_field=check_field,
+ **kwargs)
@cached_method
def _charpoly_coefficients(self):
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()
+ # 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 all we've done is a change of basis.
- J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
- self._deortho_multiplication_table,
- self._deortho_inner_product_table)
-
- # Change back from QQ to our real base ring
+ # Then, change back from QQ to our real base ring
a = ( a_i.change_ring(self.base_ring())
- for a_i in J._charpoly_coefficients() )
+ for a_i in self._rational_algebra._charpoly_coefficients() )
+
+ if self._deortho_matrix is None:
+ # This can happen if our base ring was, say, AA and we
+ # chose not to (or didn't need to) orthonormalize. It's
+ # still faster to do the computations over QQ even if
+ # the numbers in the boxes stay the same.
+ return tuple(a)
- # Now convert the coordinate variables back to the
+ # Otherwise, convert the coordinate variables back to the
# deorthonormalized ones.
R = self.coordinate_polynomial_ring()
from sage.modules.free_module_element import vector
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 ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
+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
from sage.misc.prandom import choice
eja_class = choice(cls.__subclasses__())
- # These all bubble up to the RationalBasisEuclideanJordanAlgebra
- # superclass constructor, so any (kw)args valid there are also
- # valid here.
+ # 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)
-class MatrixEuclideanJordanAlgebra:
+class MatrixEJA:
@staticmethod
- def real_embed(M):
+ def dimension_over_reals():
+ r"""
+ 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.
+
+ 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):
@classmethod
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):
- @staticmethod
- def real_embed(M):
- """
- The identity function, for embedding real matrices into real
- matrices.
- """
- return M
-
+class RealMatrixEJA(MatrixEJA):
@staticmethod
- def real_unembed(M):
- """
- The identity function, for unembedding real matrices from real
- matrices.
- """
- return M
+ def dimension_over_reals():
+ return 1
-class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
- RealMatrixEuclideanJordanAlgebra):
+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, field=AA, **kwargs):
- basis = self._denormalized_basis(n, field)
- super(RealSymmetricEJA, self).__init__(field,
- basis,
+ 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(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)
- self.one.set_cache(self(matrix.identity(field,n)))
+ idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+ self.one.set_cache(self(idV))
+
+
+
+class ComplexMatrixEJA(MatrixEJA):
+ # A manual dictionary-cache for the complex_extension() method,
+ # since apparently @classmethods can't also be @cached_methods.
+ _complex_extension = {}
+
+ @classmethod
+ def complex_extension(cls,field):
+ r"""
+ The complex field that we embed/unembed, as an extension
+ of the given ``field``.
+ """
+ if field in cls._complex_extension:
+ return cls._complex_extension[field]
+
+ # 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.
+ R = PolynomialRing(field, 'z')
+ z = R.gen()
+ F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+ cls._complex_extension[field] = F
+ return F
-class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@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()
blocks = []
for z in M.list():
- a = z.list()[0] # real part, I guess
- b = z.list()[1] # imag part, I guess
- blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
+ a = z.real()
+ b = z.imag()
+ blocks.append(matrix(field, 2, [ [ a, b],
+ [-b, a] ]))
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")
-
- # 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()
- 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
- else:
- F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+ d = cls.dimension_over_reals()
+ F = cls.complex_extension(M.base_ring())
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 trace_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
-
- 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.trace_inner_product(Xe,Ye)
- sage: actual == expected
- True
-
- """
- return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2
+ return matrix(F, n/d, elements)
-class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
- ComplexMatrixEuclideanJordanAlgebra):
+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:
#
# * The diagonal will (as a result) be real.
#
S = []
+ Eij = matrix.zero(F,n)
for i in range(n):
for j in range(i+1):
- Eij = matrix(F, n, lambda k,l: k==i and l==j)
+ # "build" E_ij
+ Eij[i,j] = 1
if i == j:
Sij = cls.real_embed(Eij)
S.append(Sij)
else:
# The second one has a minus because it's conjugated.
- Sij_real = cls.real_embed(Eij + Eij.transpose())
+ Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
+ Sij_real = cls.real_embed(Eij)
S.append(Sij_real)
- Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
+ # Eij = I*Eij - I*Eij.transpose()
+ Eij[i,j] = I
+ Eij[j,i] = -I
+ Sij_imag = cls.real_embed(Eij)
S.append(Sij_imag)
+ Eij[j,i] = 0
+ # "erase" E_ij
+ Eij[i,j] = 0
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the complex extension "F".
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,
- self.jordan_product,
- self.trace_inner_product,
- **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)
- # TODO: pre-cache the identity!
+ 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 QuaternionMatrixEJA(MatrixEJA):
+
+ # A manual dictionary-cache for the quaternion_extension() method,
+ # since apparently @classmethods can't also be @cached_methods.
+ _quaternion_extension = {}
+
+ @classmethod
+ def quaternion_extension(cls,field):
+ r"""
+ The quaternion field that we embed/unembed, as an extension
+ of the given ``field``.
+ """
+ if field in cls._quaternion_extension:
+ return cls._quaternion_extension[field]
+
+ Q = QuaternionAlgebra(field,-1,-1)
+
+ cls._quaternion_extension[field] = Q
+ return Q
-class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@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.
- field = M.base_ring()
- Q = QuaternionAlgebra(field,-1,-1)
+ Q = cls.quaternion_extension(M.base_ring())
i,j,k = Q.gens()
# Go top-left to bottom-right (reading order), converting every
# 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 trace_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.trace_inner_product(Xe,Ye)
- sage: actual == expected
- True
-
- """
- return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4
+ return matrix(Q, n/d, elements)
-class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
- QuaternionMatrixEuclideanJordanAlgebra):
+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()
# * The diagonal will (as a result) be real.
#
S = []
+ Eij = matrix.zero(Q,n)
for i in range(n):
for j in range(i+1):
- Eij = matrix(Q, n, lambda k,l: k==i and l==j)
+ # "build" E_ij
+ Eij[i,j] = 1
if i == j:
Sij = cls.real_embed(Eij)
S.append(Sij)
else:
# The second, third, and fourth ones have a minus
# because they're conjugated.
- Sij_real = cls.real_embed(Eij + Eij.transpose())
+ # Eij = Eij + Eij.transpose()
+ Eij[j,i] = 1
+ Sij_real = cls.real_embed(Eij)
S.append(Sij_real)
- Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
+ # Eij = I*(Eij - Eij.transpose())
+ Eij[i,j] = I
+ Eij[j,i] = -I
+ Sij_I = cls.real_embed(Eij)
S.append(Sij_I)
- Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
+ # Eij = J*(Eij - Eij.transpose())
+ Eij[i,j] = J
+ Eij[j,i] = -J
+ Sij_J = cls.real_embed(Eij)
S.append(Sij_J)
- Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
+ # Eij = K*(Eij - Eij.transpose())
+ Eij[i,j] = K
+ Eij[j,i] = -K
+ Sij_K = cls.real_embed(Eij)
S.append(Sij_K)
+ Eij[j,i] = 0
+ # "erase" E_ij
+ Eij[i,j] = 0
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the quaternion algebra "Q".
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,
+ 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)
- # TODO: cache one()!
+ 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(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)
- basis = V.basis()
-
- def jordan_product(x,y):
- return V([ xi*yi for (xi,yi) in zip(x,y) ])
- def inner_product(x,y):
- return x.inner_product(y)
-
- super(HadamardEJA, self).__init__(field,
- basis,
- jordan_product,
- inner_product,
- **kwargs)
+ def __init__(self, n, **kwargs):
+ if n == 0:
+ jordan_product = lambda x,y: x
+ inner_product = lambda x,y: x
+ else:
+ def jordan_product(x,y):
+ P = x.parent()
+ return P( xi*yi for (xi,yi) in zip(x,y) )
+
+ def inner_product(x,y):
+ return (x.T*y)[0,0]
+
+ # 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
+
+ column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+ super().__init__(column_basis, jordan_product, inner_product, **kwargs)
self.rank.set_cache(n)
if n == 0:
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(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 =
....: for j in range(n-1) ]
sage: actual == expected
True
- """
- def __init__(self, B, field=AA, **kwargs):
- if not B.is_positive_definite():
- raise ValueError("bilinear form is not positive-definite")
- n = B.nrows()
- V = VectorSpace(field, n)
+ """
+ def __init__(self, B, **kwargs):
+ # The matrix "B" is supplied by the user in most cases,
+ # so it makes sense to check whether or not its positive-
+ # definite unless we are specifically asked not to...
+ if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
+ if not B.is_positive_definite():
+ raise ValueError("bilinear form is not positive-definite")
+
+ # However, all of the other data for this EJA is computed
+ # by us in manner that guarantees the axioms are
+ # satisfied. So, again, unless we are specifically asked to
+ # verify things, we'll skip the rest of the checks.
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
def inner_product(x,y):
- return (B*x).inner_product(y)
+ return (y.T*B*x)[0,0]
def jordan_product(x,y):
- x0 = x[0]
- xbar = x[1:]
- y0 = y[0]
- ybar = y[1:]
- z0 = inner_product(x,y)
+ P = x.parent()
+ x0 = x[0,0]
+ xbar = x[1:,0]
+ y0 = y[0,0]
+ ybar = y[1:,0]
+ z0 = inner_product(y,x)
zbar = y0*xbar + x0*ybar
- return V([z0] + zbar.list())
+ return P([z0] + zbar.list())
- super(BilinearFormEJA, self).__init__(field,
- V.basis(),
+ n = B.nrows()
+ column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+ super(BilinearFormEJA, self).__init__(column_basis,
jordan_product,
inner_product,
**kwargs)
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(ConcreteEuclideanJordanAlgebra):
+class TrivialEJA(ConcreteEJA):
"""
The trivial Euclidean Jordan algebra consisting of only a zero element.
0
"""
- def __init__(self, field=AA, **kwargs):
+ def __init__(self, **kwargs):
jordan_product = lambda x,y: x
- inner_product = lambda x,y: field(0)
+ inner_product = lambda x,y: 0
basis = ()
- super(TrivialEJA, self).__init__(field,
- 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)
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(FiniteDimensionalEJA):
r"""
The external (orthogonal) direct sum of two other Euclidean Jordan
algebras. Essentially the Cartesian product of its two factors.
8
sage: J.rank()
5
+ sage: J.matrix_space()
+ The Cartesian product of (Full MatrixSpace of 2 by 1 dense matrices
+ over Algebraic Real Field, Full MatrixSpace of 3 by 3 dense matrices
+ over Algebraic Real Field)
TESTS:
have the same base ring; an error is raised otherwise::
sage: set_random_seed()
- sage: J1 = random_eja(AA)
- sage: J2 = random_eja(QQ,orthonormalize=False)
+ sage: J1 = random_eja(field=AA)
+ sage: J2 = random_eja(field=QQ,orthonormalize=False)
sage: J = DirectSumEJA(J1,J2)
Traceback (most recent call last):
...
raise ValueError("algebras must share the same base field")
field = J1.base_ring()
- self._factors = (J1, J2)
- n1 = J1.dimension()
- n2 = J2.dimension()
- n = n1+n2
- V = VectorSpace(field, 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(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(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)
+ M = J1.matrix_space().cartesian_product(J2.matrix_space())
+ self._cartprod_algebra = J1.cartesian_product(J2)
+
+ self._matrix_basis = tuple( [M((a,0)) for a in J1.matrix_basis()] +
+ [M((0,b)) for b in J2.matrix_basis()] )
+
+ n = len(self._matrix_basis)
+ self._sets = None
+ CombinatorialFreeModule.__init__(
+ self,
+ field,
+ range(n),
+ category=self._cartprod_algebra.category(),
+ bracket=False,
+ **kwargs)
self.rank.set_cache(J1.rank() + J2.rank())
- def factors(self):
- r"""
- Return the pair of this algebra's factors.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import (HadamardEJA,
- ....: JordanSpinEJA,
- ....: DirectSumEJA)
-
- EXAMPLES::
-
- sage: J1 = HadamardEJA(2,QQ)
- sage: J2 = JordanSpinEJA(3,QQ)
- sage: J = DirectSumEJA(J1,J2)
- sage: J.factors()
- (Euclidean Jordan algebra of dimension 2 over Rational Field,
- Euclidean Jordan algebra of dimension 3 over Rational Field)
-
- """
- return self._factors
- def projections(self):
+ def product(self,x,y):
r"""
- Return a pair of projections onto this algebra's factors.
-
SETUP::
sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
....: ComplexHermitianEJA,
....: DirectSumEJA)
- EXAMPLES::
-
- sage: J1 = JordanSpinEJA(2)
- sage: J2 = ComplexHermitianEJA(2)
- sage: J = DirectSumEJA(J1,J2)
- sage: (pi_left, pi_right) = J.projections()
- sage: J.one().to_vector()
- (1, 0, 1, 0, 0, 1)
- sage: pi_left(J.one()).to_vector()
- (1, 0)
- sage: pi_right(J.one()).to_vector()
- (1, 0, 0, 1)
-
- """
- (J1,J2) = self.factors()
- m = J1.dimension()
- n = J2.dimension()
- V_basis = self.vector_space().basis()
- # Need to specify the dimensions explicitly so that we don't
- # wind up with a zero-by-zero matrix when we want e.g. a
- # 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)
- return (pi_left, pi_right)
-
- def inclusions(self):
- r"""
- Return the pair of inclusion maps from our factors into us.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import (random_eja,
- ....: JordanSpinEJA,
- ....: RealSymmetricEJA,
- ....: DirectSumEJA)
-
- EXAMPLES::
-
- sage: J1 = JordanSpinEJA(3)
- sage: J2 = RealSymmetricEJA(2)
- sage: J = DirectSumEJA(J1,J2)
- sage: (iota_left, iota_right) = J.inclusions()
- sage: iota_left(J1.zero()) == J.zero()
- True
- sage: iota_right(J2.zero()) == J.zero()
- True
- sage: J1.one().to_vector()
- (1, 0, 0)
- sage: iota_left(J1.one()).to_vector()
- (1, 0, 0, 0, 0, 0)
- sage: J2.one().to_vector()
- (1, 0, 1)
- sage: iota_right(J2.one()).to_vector()
- (0, 0, 0, 1, 0, 1)
- sage: J.one().to_vector()
- (1, 0, 0, 1, 0, 1)
-
- TESTS:
-
- Composing a projection with the corresponding inclusion should
- produce the identity map, and mismatching them should produce
- the zero map::
+ TESTS::
sage: set_random_seed()
- sage: J1 = random_eja()
- sage: J2 = random_eja()
+ sage: J1 = JordanSpinEJA(3, field=QQ)
+ sage: J2 = ComplexHermitianEJA(2, field=QQ, orthonormalize=False)
sage: J = DirectSumEJA(J1,J2)
- sage: (iota_left, iota_right) = J.inclusions()
- sage: (pi_left, pi_right) = J.projections()
- sage: pi_left*iota_left == J1.one().operator()
- True
- sage: pi_right*iota_right == J2.one().operator()
- True
- sage: (pi_left*iota_right).is_zero()
- True
- sage: (pi_right*iota_left).is_zero()
+ sage: J.random_element()*J.random_element() in J
True
"""
- (J1,J2) = self.factors()
- m = J1.dimension()
- n = J2.dimension()
- V_basis = self.vector_space().basis()
- # Need to specify the dimensions explicitly so that we don't
- # wind up with a zero-by-zero matrix when we want e.g. a
- # 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)
- return (iota_left, iota_right)
+ xv = self._cartprod_algebra.from_vector(x.to_vector())
+ yv = self._cartprod_algebra.from_vector(y.to_vector())
+ return self.from_vector((xv*yv).to_vector())
- def inner_product(self, x, y):
- r"""
- The standard Cartesian inner-product.
- We project ``x`` and ``y`` onto our factors, and add up the
- inner-products from the subalgebras.
+ def cartesian_factors(self):
+ r"""
+ Return the pair of this algebra's factors.
SETUP::
-
sage: from mjo.eja.eja_algebra import (HadamardEJA,
- ....: QuaternionHermitianEJA,
+ ....: JordanSpinEJA,
....: DirectSumEJA)
- EXAMPLE::
+ EXAMPLES::
- sage: J1 = HadamardEJA(3,QQ)
- sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
+ sage: J1 = HadamardEJA(2, field=QQ)
+ sage: J2 = JordanSpinEJA(3, field=QQ)
sage: J = DirectSumEJA(J1,J2)
- sage: x1 = J1.one()
- sage: x2 = x1
- sage: y1 = J2.one()
- sage: y2 = y1
- sage: x1.inner_product(x2)
- 3
- sage: y1.inner_product(y2)
- 2
- sage: J.one().inner_product(J.one())
- 5
+ sage: J.cartesian_factors()
+ (Euclidean Jordan algebra of dimension 2 over Rational Field,
+ Euclidean Jordan algebra of dimension 3 over Rational Field)
"""
- (pi_left, pi_right) = self.projections()
- x1 = pi_left(x)
- x2 = pi_right(x)
- y1 = pi_left(y)
- y2 = pi_right(y)
-
- return (x1.inner_product(y1) + x2.inner_product(y2))
-
-
-
-random_eja = ConcreteEuclideanJordanAlgebra.random_instance
+ return self._cartprod_algebra.cartesian_factors()
+
+
+# def projections(self):
+# r"""
+# Return a pair of projections onto this algebra's factors.
+
+# SETUP::
+
+# sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+# ....: ComplexHermitianEJA,
+# ....: DirectSumEJA)
+
+# EXAMPLES::
+
+# sage: J1 = JordanSpinEJA(2)
+# sage: J2 = ComplexHermitianEJA(2)
+# sage: J = DirectSumEJA(J1,J2)
+# sage: (pi_left, pi_right) = J.projections()
+# sage: J.one().to_vector()
+# (1, 0, 1, 0, 0, 1)
+# sage: pi_left(J.one()).to_vector()
+# (1, 0)
+# sage: pi_right(J.one()).to_vector()
+# (1, 0, 0, 1)
+
+# """
+# (J1,J2) = self.factors()
+# m = J1.dimension()
+# n = J2.dimension()
+# V_basis = self.vector_space().basis()
+# # Need to specify the dimensions explicitly so that we don't
+# # wind up with a zero-by-zero matrix when we want e.g. a
+# # 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 = FiniteDimensionalEJAOperator(self,J1,P1)
+# pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
+# return (pi_left, pi_right)
+
+# def inclusions(self):
+# r"""
+# Return the pair of inclusion maps from our factors into us.
+
+# SETUP::
+
+# sage: from mjo.eja.eja_algebra import (random_eja,
+# ....: JordanSpinEJA,
+# ....: RealSymmetricEJA,
+# ....: DirectSumEJA)
+
+# EXAMPLES::
+
+# sage: J1 = JordanSpinEJA(3)
+# sage: J2 = RealSymmetricEJA(2)
+# sage: J = DirectSumEJA(J1,J2)
+# sage: (iota_left, iota_right) = J.inclusions()
+# sage: iota_left(J1.zero()) == J.zero()
+# True
+# sage: iota_right(J2.zero()) == J.zero()
+# True
+# sage: J1.one().to_vector()
+# (1, 0, 0)
+# sage: iota_left(J1.one()).to_vector()
+# (1, 0, 0, 0, 0, 0)
+# sage: J2.one().to_vector()
+# (1, 0, 1)
+# sage: iota_right(J2.one()).to_vector()
+# (0, 0, 0, 1, 0, 1)
+# sage: J.one().to_vector()
+# (1, 0, 0, 1, 0, 1)
+
+# TESTS:
+
+# Composing a projection with the corresponding inclusion should
+# produce the identity map, and mismatching them should produce
+# the zero map::
+
+# sage: set_random_seed()
+# sage: J1 = random_eja()
+# sage: J2 = random_eja()
+# sage: J = DirectSumEJA(J1,J2)
+# sage: (iota_left, iota_right) = J.inclusions()
+# sage: (pi_left, pi_right) = J.projections()
+# sage: pi_left*iota_left == J1.one().operator()
+# True
+# sage: pi_right*iota_right == J2.one().operator()
+# True
+# sage: (pi_left*iota_right).is_zero()
+# True
+# sage: (pi_right*iota_left).is_zero()
+# True
+
+# """
+# (J1,J2) = self.factors()
+# m = J1.dimension()
+# n = J2.dimension()
+# V_basis = self.vector_space().basis()
+# # Need to specify the dimensions explicitly so that we don't
+# # wind up with a zero-by-zero matrix when we want e.g. a
+# # 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 = FiniteDimensionalEJAOperator(J1,self,I1)
+# iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
+# return (iota_left, iota_right)
+
+# def inner_product(self, x, y):
+# r"""
+# The standard Cartesian inner-product.
+
+# We project ``x`` and ``y`` onto our factors, and add up the
+# inner-products from the subalgebras.
+
+# SETUP::
+
+
+# sage: from mjo.eja.eja_algebra import (HadamardEJA,
+# ....: QuaternionHermitianEJA,
+# ....: DirectSumEJA)
+
+# EXAMPLE::
+
+# 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
+# sage: y1 = J2.one()
+# sage: y2 = y1
+# sage: x1.inner_product(x2)
+# 3
+# sage: y1.inner_product(y2)
+# 2
+# sage: J.one().inner_product(J.one())
+# 5
+
+# """
+# (pi_left, pi_right) = self.projections()
+# x1 = pi_left(x)
+# x2 = pi_right(x)
+# y1 = pi_left(y)
+# y2 = pi_right(y)
+
+# return (x1.inner_product(y1) + x2.inner_product(y2))
+
+
+
+random_eja = ConcreteEJA.random_instance