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.
+ # Save a copy of an algebra with the original, rational basis
+ # and over QQ where computations are fast.
self._rational_algebra = None
- if orthonormalize:
- 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)
+ 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 = RationalBasisEuclideanJordanAlgebra(
+ basis,
+ jordan_product,
+ inner_product,
+ field=QQ,
+ orthonormalize=False,
+ prefix=prefix,
+ category=category,
+ check_field=False,
+ check_axioms=False)
+ if orthonormalize:
# Compute the deorthonormalized tables before we orthonormalize
# the given basis. The "check" parameter here guarantees that
# the basis is linearly-independent.
Algebraic Real Field
"""
- if self.base_ring() is QQ or self._rational_algebra is None:
+ 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. Likewise, if we never orthonormalized our
class MatrixEuclideanJordanAlgebra:
@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()
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(ConcreteEuclideanJordanAlgebra,
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)
- self.one.set_cache(self(matrix.identity(ZZ,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 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,
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)
- # 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():
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 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,
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)
- # TODO: cache one()!
+ idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+ self.one.set_cache(self(idV))
+
@staticmethod
def _max_random_instance_size():
projects / sage.d.git / blobdiff