The ``field`` we're given must be real with ``check_field=True``::
- sage: JordanSpinEJA(2,QQbar)
+ sage: JordanSpinEJA(2, field=QQbar)
Traceback (most recent call last):
...
ValueError: scalar field is not real
+ sage: JordanSpinEJA(2, field=QQbar, check_field=False)
+ Euclidean Jordan algebra of dimension 2 over Algebraic Field
The multiplication table must be square with ``check_axioms=True``::
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")
+
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")
+
self._matrix_basis = matrix_basis
if category is None:
# 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.
- self._inner_product_matrix = matrix(field, inner_product_table)
+ 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()
sage: x = J.random_element()
sage: J(x.to_vector().column()) == x
True
+
"""
msg = "not an element of this algebra"
if elt == 0:
# 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))
# 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.
+ # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
epsilon = 1e-16
for i in range(self.dimension()):
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...
"""
def __init__(self,
- field,
basis,
jordan_product,
inner_product,
+ field=AA,
orthonormalize=True,
prefix='e',
category=None,
check_field=True,
check_axioms=True):
+ 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")
+
+ # Temporary(?) hack to ensure that the matrix and vector bases
+ # are over the same ring.
+ basis = tuple( b.change_ring(field) for b in basis )
+
n = len(basis)
vector_basis = basis
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:
# 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
+ self._rational_algebra = None
if orthonormalize:
- # Compute the deorthonormalized tables before we orthonormalize
- # the given basis.
- W = V.span_of_basis( vector_basis )
+ if self.base_ring() is not QQ:
+ # There's no point in constructing the extra algebra if this
+ # one is already rational. If the original basis is rational
+ # but normalization would make it irrational, then this whole
+ # constructor will just fail anyway as it tries to stick an
+ # irrational number into a rational algebra.
+ #
+ # 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 = RationalBasisEuclideanJordanAlgebra(
+ basis,
+ jordan_product,
+ inner_product,
+ field=QQ,
+ orthonormalize=False,
+ prefix=prefix,
+ category=category,
+ check_field=False,
+ check_axioms=False)
- # TODO: use symmetry
- 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) ]
+ # Compute the deorthonormalized tables before we orthonormalize
+ # the given basis. The "check" parameter here guarantees that
+ # the basis is linearly-independent.
+ W = V.span_of_basis( vector_basis, check=check_axioms)
# Note: the Jordan and inner-products are defined in terms
# of the ambient basis. It's important that their arguments
# table is in terms of vectors
elt = _mat2vec(elt)
- # TODO: use symmetry
- elt = W.coordinate_vector(elt)
- self._deortho_multiplication_table[i][j] = elt
- self._deortho_multiplication_table[j][i] = elt
- self._deortho_inner_product_table[i][j] = ip
- 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
- vector_basis = gram_schmidt(vector_basis, inner_product)
- W = V.span_of_basis( vector_basis )
+ 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:
- from mjo.eja.eja_utils import _vec2mat
basis = tuple( map(_vec2mat,basis) )
- W = V.span_of_basis( vector_basis )
+ W = V.span_of_basis( vector_basis, check=check_axioms)
- # TODO: use symmetry
- 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) ]
+ # 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 )
- # Note: the Jordan and inner- products are defined in terms
+ # 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):
# table is in terms of vectors
elt = _mat2vec(elt)
- # TODO: use symmetry
elt = W.coordinate_vector(elt)
mult_table[i][j] = elt
- mult_table[j][i] = elt
ip_table[i][j] = ip
- ip_table[j][i] = ip
+ if 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:
EXAMPLES:
- The returned coefficients should be the same as if we'd never
- orthonormalized the basis to begin with::
-
- sage: B = matrix(QQ, [[1, 0, 0],
- ....: [0, 25, -32],
- ....: [0, -32, 41] ])
- sage: J1 = BilinearFormEJA(B)
- sage: J2 = BilinearFormEJA(B,QQ,orthonormalize=False)
- sage: J1._charpoly_coefficients()
- (X1^2 - 25*X2^2 + 64*X2*X3 - 41*X3^2, -2*X1)
- sage: J2._charpoly_coefficients()
- (X1^2 - 25*X2^2 + 64*X2*X3 - 41*X3^2, -2*X1)
-
The base ring of the resulting polynomial coefficients is what
it should be, and not the rationals (unless the algebra was
already over the rationals)::
Algebraic Real Field
sage: a0.base_ring()
Algebraic Real Field
+
"""
- if self.base_ring() is QQ:
+ if self.base_ring() is QQ or self._rational_algebra is None:
# There's no need to construct *another* algebra over the
- # rationals if this one is already over the rationals.
+ # 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.
superclass = super(RationalBasisEuclideanJordanAlgebra, self)
return superclass._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)
- a = J._charpoly_coefficients()
- return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
+ # 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
+
+ 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:
+class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
r"""
A class for the Euclidean Jordan algebras that we know by name.
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 = [[field.zero() for j in range(algebra_dim)]
- for i in range(algebra_dim)]
- for i in range(algebra_dim):
- for j in range(algebra_dim):
- mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
- mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
-
- try:
- # HACK: ignore the error here if we don't need the
- # inner product (as is the case when we construct
- # a dummy QQ-algebra for fast charpoly coefficients.
- ip_table[i][j] = self.matrix_inner_product(basis[i],
- basis[j])
- except:
- pass
-
- super(MatrixEuclideanJordanAlgebra, self).__init__(field,
- mult_table,
- ip_table,
- matrix_basis=basis,
- **kwargs)
-
- if algebra_dim == 0:
- self.one.set_cache(self.zero())
- 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)))
+ # These all bubble up to the RationalBasisEuclideanJordanAlgebra
+ # superclass constructor, so any (kw)args valid there are also
+ # valid here.
+ return eja_class.random_instance(*args, **kwargs)
- @cached_method
- def _charpoly_coefficients(self):
+
+class MatrixEuclideanJordanAlgebra:
+ @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()
class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@staticmethod
- def real_embed(M):
- """
- The identity function, for embedding real matrices into real
- matrices.
- """
- return M
-
- @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(RealMatrixEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
+ RealMatrixEuclideanJordanAlgebra):
"""
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,
- check_axioms=False,
+ 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 MatrixEuclideanJordanAlgebra 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):
@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 +
True
"""
+ super(ComplexMatrixEuclideanJordanAlgebra,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().
True
"""
+ super(ComplexMatrixEuclideanJordanAlgebra,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
-
- 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
+ return matrix(F, n/d, elements)
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
+ ComplexMatrixEuclideanJordanAlgebra):
"""
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, **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(ComplexHermitianEJA,self).__init__(field,
- basis,
- check_axioms=False,
- **kwargs)
+ 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 MatrixEuclideanJordanAlgebra 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):
@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
True
"""
+ super(QuaternionMatrixEuclideanJordanAlgebra,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()
- @staticmethod
- def real_unembed(M):
+ @classmethod
+ def real_unembed(cls,M):
"""
The inverse of _embed_quaternion_matrix().
True
"""
+ super(QuaternionMatrixEuclideanJordanAlgebra,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):
+ for l in range(n/d):
+ for m in range(n/d):
submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
- M[4*l:4*l+4,4*m:4*m+4] )
+ 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:
+ return matrix(Q, n/d, elements)
- 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
-
-
-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
+ QuaternionMatrixEuclideanJordanAlgebra):
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, **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(QuaternionHermitianEJA,self).__init__(field,
- basis,
- check_axioms=False,
- **kwargs)
+ 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 MatrixEuclideanJordanAlgebra 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(ConcreteEuclideanJordanAlgebra):
"""
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 __init__(self, n, **kwargs):
def jordan_product(x,y):
- return V([ xi*yi for (xi,yi) in zip(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)
- super(HadamardEJA, self).__init__(field,
- basis,
+ # 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)
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(ConcreteEuclideanJordanAlgebra):
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 =
sage: actual == expected
True
"""
- def __init__(self, B, field=AA, **kwargs):
+ def __init__(self, B, **kwargs):
if not B.is_positive_definite():
raise ValueError("bilinear form is not positive-definite")
- n = B.nrows()
- V = VectorSpace(field, n)
-
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 V([z0] + zbar.list())
+ return P((z0,) + tuple(zbar))
- super(BilinearFormEJA, self).__init__(field,
- V.basis(),
+ # 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)
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(ConcreteEuclideanJordanAlgebra):
"""
The trivial Euclidean Jordan algebra consisting of only a zero element.
0
"""
- def __init__(self, field=AA, **kwargs):
- mult_table = []
- ip_table = []
- super(TrivialEJA, self).__init__(field,
- mult_table,
- ip_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):
r"""
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 = None
+ ip_table = [ [ field.zero() for j in range(i+1) ]
+ for i in range(n) ]
super(DirectSumEJA, self).__init__(field,
mult_table,
ip_table,
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,
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
projects / sage.d.git / blobdiff