X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=b2891e577e9aa0c5e36b64b6a537807e9330d0e6;hb=bc02bf48592e22d034310cfffef8fb2a062c0a43;hp=2ecf8cd274bbee0be8c31e25c8a35a6e9561123a;hpb=e28bd3518185e3a87866c61d973876f84fdeea66;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 2ecf8cd..b2891e5 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -104,6 +104,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): product. This will be applied to ``basis`` to compute an inner-product table (basically a matrix) for this algebra. + - ``matrix_space`` -- the space that your matrix basis lives in, + or ``None`` (the default). So long as your basis does not have + length zero you can omit this. But in trivial algebras, it is + required. + - ``field`` -- a subfield of the reals (default: ``AA``); the scalar field for the algebra. @@ -128,7 +133,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: basis = tuple(b.superalgebra_element() for b in A.basis()) sage: J.subalgebra(basis, orthonormalize=False).is_associative() True - """ Element = FiniteDimensionalEJAElement @@ -137,6 +141,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): jordan_product, inner_product, field=AA, + matrix_space=None, orthonormalize=True, associative=None, cartesian_product=False, @@ -172,6 +177,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): category = MagmaticAlgebras(field).FiniteDimensional() category = category.WithBasis().Unital().Commutative() + if n <= 1: + # All zero- and one-dimensional algebras are just the real + # numbers with (some positive multiples of) the usual + # multiplication as its Jordan and inner-product. + associative = True if associative is None: # We should figure it out. As with check_axioms, we have to do # this without the help of the _jordan_product_is_associative() @@ -231,8 +241,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): basis = tuple(gram_schmidt(basis, inner_product)) # Save the (possibly orthonormalized) matrix basis for - # later... + # later, as well as the space that its elements live in. + # In most cases we can deduce the matrix space, but when + # n == 0 (that is, there are no basis elements) we cannot. self._matrix_basis = basis + if matrix_space is None: + self._matrix_space = self._matrix_basis[0].parent() + else: + self._matrix_space = matrix_space # Now create the vector space for the algebra, which will have # its own set of non-ambient coordinates (in terms of the @@ -1016,16 +1032,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) sage: J.matrix_space() - Full MatrixSpace of 4 by 4 dense matrices over Rational Field + Module of 2 by 2 matrices with entries in Algebraic Field over + the scalar ring Rational Field sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False) sage: J.matrix_space() - Full MatrixSpace of 4 by 4 dense matrices over Rational Field + Module of 1 by 1 matrices with entries in Quaternion + Algebra (-1, -1) with base ring Rational Field over + the scalar ring Rational Field """ - if self.is_trivial(): - return MatrixSpace(self.base_ring(), 0) - else: - return self.matrix_basis()[0].parent() + return self._matrix_space @cached_method @@ -1602,6 +1618,7 @@ class RationalBasisEJA(FiniteDimensionalEJA): jordan_product, inner_product, field=QQ, + matrix_space=self.matrix_space(), associative=self.is_associative(), orthonormalize=False, check_field=False, @@ -1729,114 +1746,164 @@ class ConcreteEJA(FiniteDimensionalEJA): return eja_class.random_instance(*args, **kwargs) -class MatrixEJA: +class MatrixEJA(FiniteDimensionalEJA): @staticmethod - def jordan_product(X,Y): - return (X*Y + Y*X)/2 - - @staticmethod - def trace_inner_product(X,Y): - r""" - A trace inner-product for matrices that aren't embedded in the - reals. It takes MATRICES as arguments, not EJA elements. + def _denormalized_basis(A): """ - return (X*Y).trace().real() + Returns a basis for the space of complex Hermitian n-by-n matrices. -class RealEmbeddedMatrixEJA(MatrixEJA): - @staticmethod - def dimension_over_reals(): - r""" - The dimension of this matrix's base ring over the reals. + Why do we embed these? Basically, because all of numerical linear + algebra assumes that you're working with vectors consisting of `n` + entries from a field and scalars from the same field. There's no way + to tell SageMath that (for example) the vectors contain complex + numbers, while the scalar field is real. - 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. + SETUP:: - This is used to determine the size of the matrix returned from - :meth:`real_embed`, among other things. - """ - raise NotImplementedError + sage: from mjo.hurwitz import (ComplexMatrixAlgebra, + ....: QuaternionMatrixAlgebra, + ....: OctonionMatrixAlgebra) + sage: from mjo.eja.eja_algebra import MatrixEJA - @classmethod - def real_embed(cls,M): - """ - Embed the matrix ``M`` into a space of real matrices. + TESTS:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = MatrixSpace(QQ, n) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B) + True - The matrix ``M`` can have entries in any field at the moment: - the real numbers, complex numbers, or quaternions. And although - they are not a field, we can probably support octonions at some - point, too. This function returns a real matrix that "acts like" - the original with respect to matrix multiplication; i.e. + :: - real_embed(M*N) = real_embed(M)*real_embed(N) + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = ComplexMatrixAlgebra(n, scalars=QQ) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B) + True - """ - if M.ncols() != M.nrows(): - raise ValueError("the matrix 'M' must be square") - return M + :: + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = QuaternionMatrixAlgebra(n, scalars=QQ) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B ) + True + + :: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: A = OctonionMatrixAlgebra(n, scalars=QQ) + sage: B = MatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B ) + True - @classmethod - def real_unembed(cls,M): - """ - The inverse of :meth:`real_embed`. """ - 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 + # These work for real MatrixSpace, whose monomials only have + # two coordinates (because the last one would always be "1"). + es = A.base_ring().gens() + gen = lambda A,m: A.monomial(m[:2]) + if hasattr(A, 'entry_algebra_gens'): + # We've got a MatrixAlgebra, and its monomials will have + # three coordinates. + es = A.entry_algebra_gens() + gen = lambda A,m: A.monomial(m) - @classmethod - def trace_inner_product(cls,X,Y): + basis = [] + for i in range(A.nrows()): + for j in range(i+1): + if i == j: + E_ii = gen(A, (i,j,es[0])) + basis.append(E_ii) + else: + for e in es: + E_ij = gen(A, (i,j,e)) + E_ij += E_ij.conjugate_transpose() + basis.append(E_ij) + + return tuple( basis ) + + @staticmethod + def jordan_product(X,Y): + return (X*Y + Y*X)/2 + + @staticmethod + def trace_inner_product(X,Y): r""" - Compute the trace inner-product of two real-embeddings. + A trace inner-product for matrices that aren't embedded in the + reals. It takes MATRICES as arguments, not EJA elements. SETUP:: - sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, - ....: QuaternionHermitianEJA) + sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA, + ....: OctonionHermitianEJA) EXAMPLES:: - 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: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 :: - 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 + sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 + + :: + + sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 + + :: + + sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 """ - # This does in fact compute the real part of the trace. - # If we compute the trace of e.g. a complex matrix M, - # then we do so by adding up its diagonal entries -- - # call them z_1 through z_n. The real embedding of z_1 - # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace - # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth. - return (X*Y).trace()/cls.dimension_over_reals() + tr = (X*Y).trace() + if hasattr(tr, 'coefficient'): + # Works for octonions, and has to come first because they + # also have a "real()" method that doesn't return an + # element of the scalar ring. + return tr.coefficient(0) + elif hasattr(tr, 'coefficient_tuple'): + # Works for quaternions. + return tr.coefficient_tuple()[0] -class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): + # Works for real and complex numbers. + return tr.real() + + + def __init__(self, matrix_space, **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().__init__(self._denormalized_basis(matrix_space), + self.jordan_product, + self.trace_inner_product, + field=matrix_space.base_ring(), + matrix_space=matrix_space, + **kwargs) + + self.rank.set_cache(matrix_space.nrows()) + self.one.set_cache( self(matrix_space.one()) ) + +class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1901,38 +1968,6 @@ class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ - @classmethod - def _denormalized_basis(cls, n, field): - """ - Return a basis for the space of real symmetric n-by-n matrices. - - SETUP:: - - sage: from mjo.eja.eja_algebra import RealSymmetricEJA - - TESTS:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ) - sage: all( M.is_symmetric() for M in B) - True - - """ - # The basis of symmetric matrices, as matrices, in their R^(n-by-n) - # coordinates. - S = [] - for i in range(n): - for j in range(i+1): - Eij = matrix(field, n, lambda k,l: k==i and l==j) - if i == j: - Sij = Eij - else: - Sij = Eij + Eij.transpose() - S.append(Sij) - return tuple(S) - - @staticmethod def _max_random_instance_size(): return 4 # Dimension 10 @@ -1950,175 +1985,12 @@ class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - associative = False - if n <= 1: - associative = True - - super().__init__(self._denormalized_basis(n,field), - self.jordan_product, - self.trace_inner_product, - field=field, - associative=associative, - **kwargs) - - # TODO: this could be factored out somehow, but is left here - # because the MatrixEJA is not presently a subclass of the - # FDEJA class that defines rank() and one(). - self.rank.set_cache(n) - idV = self.matrix_space().one() - self.one.set_cache(self(idV)) - - + A = MatrixSpace(field, n) + super().__init__(A, **kwargs) -class ComplexMatrixEJA(RealEmbeddedMatrixEJA): - # 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 - - @staticmethod - 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 + - bi` to the block matrix ``[[a,b],[-b,a]]``. - - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexMatrixEJA - - EXAMPLES:: - - sage: F = QuadraticField(-1, 'I') - sage: x1 = F(4 - 2*i) - sage: x2 = F(1 + 2*i) - sage: x3 = F(-i) - sage: x4 = F(6) - sage: M = matrix(F,2,[[x1,x2],[x3,x4]]) - sage: ComplexMatrixEJA.real_embed(M) - [ 4 -2| 1 2] - [ 2 4|-2 1] - [-----+-----] - [ 0 -1| 6 0] - [ 1 0| 0 6] - - TESTS: - - Embedding is a homomorphism (isomorphism, in fact):: - - sage: set_random_seed() - sage: n = ZZ.random_element(3) - sage: F = QuadraticField(-1, 'I') - sage: X = random_matrix(F, n) - sage: Y = random_matrix(F, n) - 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().real_embed(M) - n = M.nrows() - - # We don't need any adjoined elements... - field = M.base_ring().base_ring() - - blocks = [] - for z in M.list(): - a = z.real() - b = z.imag() - blocks.append(matrix(field, 2, [ [ a, b], - [-b, a] ])) - return matrix.block(field, n, blocks) - - - @classmethod - def real_unembed(cls,M): - """ - The inverse of _embed_complex_matrix(). - - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexMatrixEJA - - EXAMPLES:: - - sage: A = matrix(QQ,[ [ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [ 9, 10, 11, 12], - ....: [-10, 9, -12, 11] ]) - sage: ComplexMatrixEJA.real_unembed(A) - [ 2*I + 1 4*I + 3] - [ 10*I + 9 12*I + 11] - - TESTS: - - Unembedding is the inverse of embedding:: - - sage: set_random_seed() - sage: F = QuadraticField(-1, 'I') - sage: M = random_matrix(F, 3) - sage: Me = ComplexMatrixEJA.real_embed(M) - sage: ComplexMatrixEJA.real_unembed(Me) == M - True - - """ - super().real_unembed(M) - n = ZZ(M.nrows()) - 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/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]: - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0] + submat[0,1]*i - elements.append(z) - - return matrix(F, n/d, elements) - - -class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, ComplexMatrixEJA): +class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -2131,13 +2003,28 @@ class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, ComplexMatrixEJA): EXAMPLES: - In theory, our "field" can be any subfield of the reals:: + In theory, our "field" can be any subfield of the reals, but we + can't use inexact real fields at the moment because SageMath + doesn't know how to convert their elements into complex numbers, + or even into algebraic reals:: - sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True) - Euclidean Jordan algebra of dimension 4 over Real Double Field - sage: ComplexHermitianEJA(2, field=RR, check_axioms=True) - Euclidean Jordan algebra of dimension 4 over Real Field with - 53 bits of precision + sage: QQbar(RDF(1)) + Traceback (most recent call last): + ... + TypeError: Illegal initializer for algebraic number + sage: AA(RR(1)) + Traceback (most recent call last): + ... + TypeError: Illegal initializer for algebraic number + + This causes the following error when we try to scale a matrix of + complex numbers by an inexact real number:: + + sage: ComplexHermitianEJA(2,field=RR) + Traceback (most recent call last): + ... + TypeError: Unable to coerce entries (=(1.00000000000000, + -0.000000000000000)) to coefficients in Algebraic Real Field TESTS: @@ -2173,92 +2060,16 @@ class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, ComplexMatrixEJA): sage: ComplexHermitianEJA(0) Euclidean Jordan algebra of dimension 0 over Algebraic Real Field - """ - - @classmethod - def _denormalized_basis(cls, n, field): - """ - Returns a basis for the space of complex Hermitian n-by-n matrices. - - Why do we embed these? Basically, because all of numerical linear - algebra assumes that you're working with vectors consisting of `n` - entries from a field and scalars from the same field. There's no way - to tell SageMath that (for example) the vectors contain complex - numbers, while the scalar field is real. - - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA - - TESTS:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ) - sage: all( M.is_symmetric() for M in B) - True - - """ - R = PolynomialRing(ZZ, 'z') - z = R.gen() - F = ZZ.extension(z**2 + 1, 'I') - I = F.gen(1) - - # This is like the symmetric case, but we need to be careful: - # - # * We want conjugate-symmetry, not just symmetry. - # * 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): - # "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. - Eij[j,i] = 1 # Eij = Eij + Eij.transpose() - Sij_real = cls.real_embed(Eij) - S.append(Sij_real) - # 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 the entries, we can drop back to the - # desired real "field" instead of the extension "F". - return tuple( s.change_ring(field) for s in S ) - - def __init__(self, n, field=AA, **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 - associative = False - if n <= 1: - associative = True + from mjo.hurwitz import ComplexMatrixAlgebra + A = ComplexMatrixAlgebra(n, scalars=field) + super().__init__(A, **kwargs) - super().__init__(self._denormalized_basis(n,field), - self.jordan_product, - self.trace_inner_product, - field=field, - associative=associative, - **kwargs) - # TODO: this could be factored out somehow, but is left here - # because the MatrixEJA is not presently a subclass of the - # FDEJA class that defines rank() and one(). - self.rank.set_cache(n) - idV = matrix.identity(ZZ, self.dimension_over_reals()*n) - self.one.set_cache(self(idV)) @staticmethod def _max_random_instance_size(): @@ -2272,157 +2083,8 @@ class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, ComplexMatrixEJA): n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) -class QuaternionMatrixEJA(RealEmbeddedMatrixEJA): - # 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 - - @staticmethod - 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 - = a + bi + cj + dk` to the block-complex matrix ``[[a + bi, - c+di],[-c + di, a-bi]]`, and then embedding those into a real - matrix. - - SETUP:: - - sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA - - EXAMPLES:: - - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: i,j,k = Q.gens() - sage: x = 1 + 2*i + 3*j + 4*k - sage: M = matrix(Q, 1, [[x]]) - sage: QuaternionMatrixEJA.real_embed(M) - [ 1 2 3 4] - [-2 1 -4 3] - [-3 4 1 -2] - [-4 -3 2 1] - - Embedding is a homomorphism (isomorphism, in fact):: - - sage: set_random_seed() - sage: n = ZZ.random_element(2) - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: X = random_matrix(Q, n) - sage: Y = random_matrix(Q, n) - 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().real_embed(M) - quaternions = M.base_ring() - n = M.nrows() - - F = QuadraticField(-1, 'I') - i = F.gen() - - blocks = [] - for z in M.list(): - t = z.coefficient_tuple() - a = t[0] - b = t[1] - c = t[2] - d = t[3] - cplxM = matrix(F, 2, [[ a + b*i, c + d*i], - [-c + d*i, a - b*i]]) - realM = ComplexMatrixEJA.real_embed(cplxM) - blocks.append(realM) - - # We should have real entries by now, so use the realest field - # we've got for the return value. - return matrix.block(quaternions.base_ring(), n, blocks) - - - - @classmethod - def real_unembed(cls,M): - """ - The inverse of _embed_quaternion_matrix(). - - SETUP:: - - sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA - - EXAMPLES:: - - sage: M = matrix(QQ, [[ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [-3, 4, 1, -2], - ....: [-4, -3, 2, 1]]) - sage: QuaternionMatrixEJA.real_unembed(M) - [1 + 2*i + 3*j + 4*k] - - TESTS: - - Unembedding is the inverse of embedding:: - - sage: set_random_seed() - sage: Q = QuaternionAlgebra(QQ, -1, -1) - sage: M = random_matrix(Q, 3) - sage: Me = QuaternionMatrixEJA.real_embed(M) - sage: QuaternionMatrixEJA.real_unembed(Me) == M - True - - """ - super().real_unembed(M) - n = ZZ(M.nrows()) - 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. - 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/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(): - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0].real() - z += submat[0,0].imag()*i - z += submat[0,1].real()*j - z += submat[0,1].imag()*k - elements.append(z) - - return matrix(Q, n/d, elements) - - -class QuaternionHermitianEJA(RationalBasisEJA, - ConcreteEJA, - QuaternionMatrixEJA): +class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the @@ -2479,100 +2141,14 @@ class QuaternionHermitianEJA(RationalBasisEJA, Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ - @classmethod - def _denormalized_basis(cls, n, field): - """ - Returns a basis for the space of quaternion Hermitian n-by-n matrices. - - Why do we embed these? Basically, because all of numerical - linear algebra assumes that you're working with vectors consisting - of `n` entries from a field and scalars from the same field. There's - no way to tell SageMath that (for example) the vectors contain - complex numbers, while the scalar field is real. - - SETUP:: - - sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA - - TESTS:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ) - sage: all( M.is_symmetric() for M in B ) - True - - """ - Q = QuaternionAlgebra(QQ,-1,-1) - I,J,K = Q.gens() - - # This is like the symmetric case, but we need to be careful: - # - # * We want conjugate-symmetry, not just symmetry. - # * 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): - # "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. - # Eij = Eij + Eij.transpose() - Eij[j,i] = 1 - Sij_real = cls.real_embed(Eij) - S.append(Sij_real) - # Eij = I*(Eij - Eij.transpose()) - Eij[i,j] = I - Eij[j,i] = -I - Sij_I = cls.real_embed(Eij) - S.append(Sij_I) - # Eij = J*(Eij - Eij.transpose()) - Eij[i,j] = J - Eij[j,i] = -J - Sij_J = cls.real_embed(Eij) - S.append(Sij_J) - # 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 the entries, we can drop back to the - # desired real "field" instead of the quaternion algebra "Q". - return tuple( s.change_ring(field) for s in S ) - - def __init__(self, n, field=AA, **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 - associative = False - if n <= 1: - associative = True - - super().__init__(self._denormalized_basis(n,field), - self.jordan_product, - self.trace_inner_product, - field=field, - associative=associative, - **kwargs) - - # TODO: this could be factored out somehow, but is left here - # because the MatrixEJA is not presently a subclass of the - # FDEJA class that defines rank() and one(). - self.rank.set_cache(n) - idV = matrix.identity(ZZ, self.dimension_over_reals()*n) - self.one.set_cache(self(idV)) + from mjo.hurwitz import QuaternionMatrixAlgebra + A = QuaternionMatrixAlgebra(n, scalars=field) + super().__init__(A, **kwargs) @staticmethod @@ -2590,12 +2166,13 @@ class QuaternionHermitianEJA(RationalBasisEJA, n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) -class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): +class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): r""" SETUP:: sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA, ....: OctonionHermitianEJA) + sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra EXAMPLES: @@ -2610,7 +2187,8 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): After a change-of-basis, the 2-by-2 algebra has the same multiplication table as the ten-dimensional Jordan spin algebra:: - sage: b = OctonionHermitianEJA._denormalized_basis(2,QQ) + sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ) + sage: b = OctonionHermitianEJA._denormalized_basis(A) sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],) sage: jp = OctonionHermitianEJA.jordan_product sage: ip = OctonionHermitianEJA.trace_inner_product @@ -2699,82 +2277,9 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - super().__init__(self._denormalized_basis(n,field), - self.jordan_product, - self.trace_inner_product, - field=field, - **kwargs) - - # TODO: this could be factored out somehow, but is left here - # because the MatrixEJA is not presently a subclass of the - # FDEJA class that defines rank() and one(). - self.rank.set_cache(n) - idV = self.matrix_space().one() - self.one.set_cache(self(idV)) - - - @classmethod - def _denormalized_basis(cls, n, field): - """ - Returns a basis for the space of octonion Hermitian n-by-n - matrices. - - SETUP:: - - sage: from mjo.eja.eja_algebra import OctonionHermitianEJA - - EXAMPLES:: - - sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ) - sage: all( M.is_hermitian() for M in B ) - True - sage: len(B) - 27 - - """ from mjo.hurwitz import OctonionMatrixAlgebra - MS = OctonionMatrixAlgebra(n, scalars=field) - es = MS.entry_algebra().gens() - - basis = [] - for i in range(n): - for j in range(i+1): - if i == j: - E_ii = MS.monomial( (i,j,es[0]) ) - basis.append(E_ii) - else: - for e in es: - E_ij = MS.monomial( (i,j,e) ) - ec = e.conjugate() - # If the conjugate has a negative sign in front - # of it, (j,i,ec) won't be a monomial! - if (j,i,ec) in MS.indices(): - E_ij += MS.monomial( (j,i,ec) ) - else: - E_ij -= MS.monomial( (j,i,-ec) ) - basis.append(E_ij) - - return tuple( basis ) - - @staticmethod - def trace_inner_product(X,Y): - r""" - The octonions don't know that the reals are embedded in them, - so we have to take the e0 component ourselves. - - SETUP:: - - sage: from mjo.eja.eja_algebra import OctonionHermitianEJA - - TESTS:: - - sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False) - sage: I = J.one().to_matrix() - sage: J.trace_inner_product(I, -I) - -2 - - """ - return (X*Y).trace().real().coefficient(0) + A = OctonionMatrixAlgebra(n, scalars=field) + super().__init__(A, **kwargs) class AlbertEJA(OctonionHermitianEJA): @@ -2839,13 +2344,14 @@ class HadamardEJA(RationalBasisEJA, ConcreteEJA): (r0, r1, r2) """ def __init__(self, n, field=AA, **kwargs): + MS = MatrixSpace(field, n, 1) + 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) ) + return MS( xi*yi for (xi,yi) in zip(x,y) ) def inner_product(x,y): return (x.T*y)[0,0] @@ -2859,20 +2365,17 @@ class HadamardEJA(RationalBasisEJA, ConcreteEJA): 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(field, n).basis() ) + column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() ) super().__init__(column_basis, jordan_product, inner_product, field=field, + matrix_space=MS, associative=True, **kwargs) self.rank.set_cache(n) - if n == 0: - self.one.set_cache( self.zero() ) - else: - self.one.set_cache( sum(self.gens()) ) + self.one.set_cache( self.sum(self.gens()) ) @staticmethod def _max_random_instance_size(): @@ -2986,22 +2489,22 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): # verify things, we'll skip the rest of the checks. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + n = B.nrows() + MS = MatrixSpace(field, n, 1) + def inner_product(x,y): return (y.T*B*x)[0,0] def jordan_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 P([z0] + zbar.list()) + return MS([z0] + zbar.list()) - n = B.nrows() - column_basis = tuple( b.column() - for b in FreeModule(field, n).basis() ) + column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() ) # TODO: I haven't actually checked this, but it seems legit. associative = False @@ -3012,6 +2515,7 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): jordan_product, inner_product, field=field, + matrix_space=MS, associative=associative, **kwargs) @@ -3019,7 +2523,6 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): # one-dimensional ambient space (because the rank is bounded # by the ambient dimension). self.rank.set_cache(min(n,2)) - if n == 0: self.one.set_cache( self.zero() ) else: @@ -3167,10 +2670,11 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA): 0 """ - def __init__(self, **kwargs): + def __init__(self, field=AA, **kwargs): jordan_product = lambda x,y: x - inner_product = lambda x,y: 0 + inner_product = lambda x,y: field.zero() basis = () + MS = MatrixSpace(field,0) # New defaults for keyword arguments if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False @@ -3180,6 +2684,8 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA): jordan_product, inner_product, associative=True, + field=field, + matrix_space=MS, **kwargs) # The rank is zero using my definition, namely the dimension of the @@ -3357,7 +2863,14 @@ class CartesianProductEJA(FiniteDimensionalEJA): associative = all( f.is_associative() for f in factors ) - MS = self.matrix_space() + # Compute my matrix space. This category isn't perfect, but + # is good enough for what we need to do. + MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis() + MS_cat = MS_cat.Unital().CartesianProducts() + MS_factors = tuple( J.matrix_space() for J in factors ) + from sage.sets.cartesian_product import CartesianProduct + MS = CartesianProduct(MS_factors, MS_cat) + basis = [] zero = MS.zero() for i in range(m): @@ -3392,15 +2905,16 @@ class CartesianProductEJA(FiniteDimensionalEJA): jordan_product, inner_product, field=field, + matrix_space=MS, orthonormalize=False, associative=associative, cartesian_product=True, check_field=False, check_axioms=False) + self.rank.set_cache(sum(J.rank() for J in factors)) ones = tuple(J.one().to_matrix() for J in factors) self.one.set_cache(self(ones)) - self.rank.set_cache(sum(J.rank() for J in factors)) def cartesian_factors(self): # Copy/pasted from CombinatorialFreeModule_CartesianProduct. @@ -3452,11 +2966,13 @@ class CartesianProductEJA(FiniteDimensionalEJA): sage: J2 = ComplexHermitianEJA(1) sage: J = cartesian_product([J1,J2]) sage: J.one().to_matrix()[0] - [1 0] - [0 1] + +---+ + | 1 | + +---+ sage: J.one().to_matrix()[1] - [1 0] - [0 1] + +---+ + | 1 | + +---+ :: @@ -3473,16 +2989,7 @@ class CartesianProductEJA(FiniteDimensionalEJA): +----+ """ - scalars = self.cartesian_factor(0).base_ring() - - # This category isn't perfect, but is good enough for what we - # need to do. - cat = MagmaticAlgebras(scalars).FiniteDimensional().WithBasis() - cat = cat.Unital().CartesianProducts() - factors = tuple( J.matrix_space() for J in self.cartesian_factors() ) - - from sage.sets.cartesian_product import CartesianProduct - return CartesianProduct(factors, cat) + return super().matrix_space() @cached_method