X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=106a0cddec06355a952e71d75699677b25dd7da9;hb=fc29add6cf1d9ff4e8a240b0f8f2ca6672d4ea57;hp=f2d7ba70c82ae38375d3496521eef14463566897;hpb=2e4211deee2c4556f2c2b4dccc857106ce6bc89a;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index f2d7ba7..106a0cd 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -32,22 +32,25 @@ for these simple algebras: * :class:`RealSymmetricEJA` * :class:`ComplexHermitianEJA` * :class:`QuaternionHermitianEJA` + * :class:`OctonionHermitianEJA` -Missing from this list is the algebra of three-by-three octononion -Hermitian matrices, as there is (as of yet) no implementation of the -octonions in SageMath. In addition to these, we provide two other -example constructions, +In addition to these, we provide two other example constructions, + * :class:`JordanSpinEJA` * :class:`HadamardEJA` + * :class:`AlbertEJA` * :class:`TrivialEJA` The Jordan spin algebra is a bilinear form algebra where the bilinear form is the identity. The Hadamard EJA is simply a Cartesian product -of one-dimensional spin algebras. And last but not least, the trivial -EJA is exactly what you think. Cartesian products of these are also -supported using the usual ``cartesian_product()`` function; as a -result, we support (up to isomorphism) all Euclidean Jordan algebras -that don't involve octonions. +of one-dimensional spin algebras. The Albert EJA is simply a special +case of the :class:`OctonionHermitianEJA` where the matrices are +three-by-three and the resulting space has dimension 27. And +last/least, the trivial EJA is exactly what you think it is; it could +also be obtained by constructing a dimension-zero instance of any of +the other algebras. Cartesian products of these are also supported +using the usual ``cartesian_product()`` function; as a result, we +support (up to isomorphism) all Euclidean Jordan algebras. SETUP:: @@ -59,8 +62,6 @@ EXAMPLES:: Euclidean Jordan algebra of dimension... """ -from itertools import repeat - from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.magmatic_algebras import MagmaticAlgebras from sage.categories.sets_cat import cartesian_product @@ -171,6 +172,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() @@ -1018,7 +1024,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Full MatrixSpace of 4 by 4 dense matrices over 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(): @@ -1543,7 +1551,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): class RationalBasisEJA(FiniteDimensionalEJA): r""" - New class for algebras whose supplied basis elements have all rational entries. + Algebras whose supplied basis elements have all rational entries. SETUP:: @@ -1574,7 +1582,11 @@ class RationalBasisEJA(FiniteDimensionalEJA): 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 ): + # Use _all2list to get the vector coordinates of octonion + # entries and not the octonions themselves (which are not + # rational). + if not all( all(b_i in QQ for b_i in _all2list(b)) + for b in basis ): raise TypeError("basis not rational") super().__init__(basis, @@ -1726,112 +1738,147 @@ class ConcreteEJA(FiniteDimensionalEJA): class MatrixEJA: @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:: - 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. + 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 - real_embed(M*N) = real_embed(M)*real_embed(N) + :: - """ - 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 = ComplexMatrixAlgebra(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 = 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] + + # Works for real and complex numbers. + return tr.real() -class RealSymmetricEJA(ConcreteEJA, RationalBasisEJA, MatrixEJA): + + +class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1896,38 +1943,6 @@ class RealSymmetricEJA(ConcreteEJA, RationalBasisEJA, 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 @@ -1945,175 +1960,22 @@ class RealSymmetricEJA(ConcreteEJA, RationalBasisEJA, 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), + A = MatrixSpace(field, n) + super().__init__(self._denormalized_basis(A), 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)) - - - -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() + self.one.set_cache(self(A.one())) - # 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(ConcreteEJA, RationalBasisEJA, ComplexMatrixEJA): +class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -2170,90 +2032,24 @@ class ComplexHermitianEJA(ConcreteEJA, RationalBasisEJA, ComplexMatrixEJA): 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 - - super().__init__(self._denormalized_basis(n,field), + from mjo.hurwitz import ComplexMatrixAlgebra + A = ComplexMatrixAlgebra(n, scalars=field) + super().__init__(self._denormalized_basis(A), 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)) + self.one.set_cache(self(A.one())) @staticmethod def _max_random_instance_size(): @@ -2267,157 +2063,8 @@ class ComplexHermitianEJA(ConcreteEJA, RationalBasisEJA, 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(ConcreteEJA, - RationalBasisEJA, - QuaternionMatrixEJA): +class QuaternionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the @@ -2474,100 +2121,24 @@ class QuaternionHermitianEJA(ConcreteEJA, 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), + from mjo.hurwitz import QuaternionMatrixAlgebra + A = QuaternionMatrixAlgebra(n, scalars=field) + super().__init__(self._denormalized_basis(A), 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)) + self.one.set_cache(self(A.one())) @staticmethod @@ -2585,7 +2156,7 @@ class QuaternionHermitianEJA(ConcreteEJA, n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) -class OctonionHermitianEJA(FiniteDimensionalEJA, MatrixEJA): +class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA): r""" SETUP:: @@ -2668,7 +2239,23 @@ class OctonionHermitianEJA(FiniteDimensionalEJA, MatrixEJA): sage: J.rank.clear_cache() # long time sage: J.rank() # long time 2 + """ + @staticmethod + def _max_random_instance_size(): + r""" + The maximum rank of a random QuaternionHermitianEJA. + """ + return 1 # Dimension 1 + + @classmethod + 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, **kwargs) + def __init__(self, n, field=AA, **kwargs): if n > 3: # Otherwise we don't get an EJA. @@ -2678,7 +2265,9 @@ class OctonionHermitianEJA(FiniteDimensionalEJA, 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), + from mjo.hurwitz import OctonionMatrixAlgebra + A = OctonionMatrixAlgebra(n, scalars=field) + super().__init__(self._denormalized_basis(A), self.jordan_product, self.trace_inner_product, field=field, @@ -2688,81 +2277,39 @@ class OctonionHermitianEJA(FiniteDimensionalEJA, MatrixEJA): # 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)) - + self.one.set_cache(self(A.one())) - @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.octonions 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 ) +class AlbertEJA(OctonionHermitianEJA): + r""" + The Albert algebra is the algebra of three-by-three Hermitian + matrices whose entries are octonions. - @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:: - SETUP:: + sage: from mjo.eja.eja_algebra import AlbertEJA - sage: from mjo.eja.eja_algebra import OctonionHermitianEJA + EXAMPLES:: - TESTS:: + sage: AlbertEJA(field=QQ, orthonormalize=False) + Euclidean Jordan algebra of dimension 27 over Rational Field + sage: AlbertEJA() # long time + Euclidean Jordan algebra of dimension 27 over Algebraic Real Field - sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False) - sage: I = J.one().to_matrix() - sage: J.trace_inner_product(I, -I) - -2 + """ + def __init__(self, *args, **kwargs): + super().__init__(3, *args, **kwargs) - """ - return (X*Y).trace().real().coefficient(0) -class HadamardEJA(ConcreteEJA, RationalBasisEJA): +class HadamardEJA(RationalBasisEJA, ConcreteEJA): """ - Return the Euclidean Jordan Algebra corresponding to the set - `R^n` under the Hadamard product. + Return the Euclidean Jordan algebra on `R^n` with the Hadamard + (pointwise real-number multiplication) Jordan product and the + usual inner-product. - Note: this is nothing more than the Cartesian product of ``n`` - copies of the spin algebra. Once Cartesian product algebras - are implemented, this can go. + This is nothing more than the Cartesian product of ``n`` copies of + the one-dimensional Jordan spin algebra, and is the most common + example of a non-simple Euclidean Jordan algebra. SETUP:: @@ -2793,7 +2340,6 @@ class HadamardEJA(ConcreteEJA, RationalBasisEJA): sage: HadamardEJA(3, prefix='r').gens() (r0, r1, r2) - """ def __init__(self, n, field=AA, **kwargs): if n == 0: @@ -2847,7 +2393,7 @@ class HadamardEJA(ConcreteEJA, RationalBasisEJA): return cls(n, **kwargs) -class BilinearFormEJA(ConcreteEJA, RationalBasisEJA): +class BilinearFormEJA(RationalBasisEJA, 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 = @@ -3095,7 +2641,7 @@ class JordanSpinEJA(BilinearFormEJA): return cls(n, **kwargs) -class TrivialEJA(ConcreteEJA, RationalBasisEJA): +class TrivialEJA(RationalBasisEJA, ConcreteEJA): """ The trivial Euclidean Jordan algebra consisting of only a zero element. @@ -3642,7 +3188,9 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA, SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA, + ....: OctonionHermitianEJA, ....: RealSymmetricEJA) EXAMPLES: @@ -3659,15 +3207,32 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA, sage: J.rank() 5 + TESTS: + + The ``cartesian_product()`` function only uses the first factor to + decide where the result will live; thus we have to be careful to + check that all factors do indeed have a `_rational_algebra` member + before we try to access it:: + + sage: J1 = OctonionHermitianEJA(1) # no rational basis + sage: J2 = HadamardEJA(2) + sage: cartesian_product([J1,J2]) + Euclidean Jordan algebra of dimension 1 over Algebraic Real Field + (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field + sage: cartesian_product([J2,J1]) + Euclidean Jordan algebra of dimension 2 over Algebraic Real Field + (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field + """ def __init__(self, algebras, **kwargs): CartesianProductEJA.__init__(self, algebras, **kwargs) self._rational_algebra = None if self.vector_space().base_field() is not QQ: - self._rational_algebra = cartesian_product([ - r._rational_algebra for r in algebras - ]) + if all( hasattr(r, "_rational_algebra") for r in algebras ): + self._rational_algebra = cartesian_product([ + r._rational_algebra for r in algebras + ]) RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA