X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=afd0f55c103c885315e99bc90afe647262aadfbd;hb=3844eed972d91ce88b1504818b37ee9428d95c68;hp=79efd3876b6fadd0d2a7494a45856bf46bf3a621;hpb=4d2f874e306baa8d99fd938e978b8c968ca3a82e;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 79efd38..afd0f55 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -166,9 +166,10 @@ from sage.modules.free_module import FreeModule, VectorSpace from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF, PolynomialRing, QuadraticField) -from mjo.eja.eja_element import FiniteDimensionalEJAElement +from mjo.eja.eja_element import (CartesianProductEJAElement, + FiniteDimensionalEJAElement) from mjo.eja.eja_operator import FiniteDimensionalEJAOperator -from mjo.eja.eja_utils import _all2list, _mat2vec +from mjo.eja.eja_utils import _all2list def EuclideanJordanAlgebras(field): r""" @@ -229,7 +230,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): We should compute that an element subalgebra is associative even if we circumvent the element method:: - sage: set_random_seed() sage: J = random_eja(field=QQ,orthonormalize=False) sage: x = J.random_element() sage: A = x.subalgebra_generated_by(orthonormalize=False) @@ -366,7 +366,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): if orthonormalize: # Now "self._matrix_span" is the vector space of our - # algebra coordinates. The variables "X1", "X2",... refer + # algebra coordinates. The variables "X0", "X1",... refer # to the entries of vectors in self._matrix_span. Thus to # convert back and forth between the orthonormal # coordinates and the given ones, we need to stick the @@ -431,7 +431,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): TESTS:: - sage: set_random_seed() sage: J = random_eja() sage: J(1) Traceback (most recent call last): @@ -456,7 +455,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): TESTS:: - sage: set_random_seed() sage: J = random_eja() sage: n = J.dimension() sage: bi = J.zero() @@ -498,7 +496,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Our inner product is "associative," which means the following for a symmetric bilinear form:: - sage: set_random_seed() sage: J = random_eja() sage: x,y,z = J.random_elements(3) sage: (x*y).inner_product(z) == y.inner_product(x*z) @@ -509,7 +506,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Ensure that this is the usual inner product for the algebras over `R^n`:: - sage: set_random_seed() sage: J = HadamardEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = x.inner_product(y) @@ -522,7 +518,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): one). This is in Faraut and Koranyi, and also my "On the symmetry..." paper:: - sage: set_random_seed() sage: J = BilinearFormEJA.random_instance() sage: n = J.dimension() sage: x = J.random_element() @@ -635,7 +630,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): The values we've presupplied to the constructors agree with the computation:: - sage: set_random_seed() sage: J = random_eja() sage: J.is_associative() == J._jordan_product_is_associative() True @@ -757,7 +751,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Ensure that we can convert any element back and forth faithfully between its matrix and algebra representations:: - sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: J(x.to_matrix()) == x @@ -870,7 +863,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J = JordanSpinEJA(3) sage: p = J.characteristic_polynomial_of(); p - X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 + X0^2 - X1^2 - X2^2 + (-2*t)*X0 + t^2 sage: xvec = J.one().to_vector() sage: p(*xvec) t^2 - 2*t + 1 @@ -919,13 +912,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J = HadamardEJA(2) sage: J.coordinate_polynomial_ring() - Multivariate Polynomial Ring in X1, X2... + Multivariate Polynomial Ring in X0, X1... sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False) sage: J.coordinate_polynomial_ring() - Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6... + Multivariate Polynomial Ring in X0, X1, X2, X3, X4, X5... """ - var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) ) + var_names = tuple( "X%d" % z for z in range(self.dimension()) ) return PolynomialRing(self.base_ring(), var_names) def inner_product(self, x, y): @@ -947,7 +940,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Our inner product is "associative," which means the following for a symmetric bilinear form:: - sage: set_random_seed() sage: J = random_eja() sage: x,y,z = J.random_elements(3) sage: (x*y).inner_product(z) == y.inner_product(x*z) @@ -958,7 +950,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Ensure that this is the usual inner product for the algebras over `R^n`:: - sage: set_random_seed() sage: J = HadamardEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = x.inner_product(y) @@ -971,7 +962,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): one). This is in Faraut and Koranyi, and also my "On the symmetry..." paper:: - sage: set_random_seed() sage: J = BilinearFormEJA.random_instance() sage: n = J.dimension() sage: x = J.random_element() @@ -1199,7 +1189,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): The identity element acts like the identity, regardless of whether or not we orthonormalize:: - sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: J.one()*x == x and x*J.one() == x @@ -1211,7 +1200,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): :: - sage: set_random_seed() sage: J = random_eja(field=QQ, orthonormalize=False) sage: x = J.random_element() sage: J.one()*x == x and x*J.one() == x @@ -1225,7 +1213,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): regardless of the base field and whether or not we orthonormalize:: - sage: set_random_seed() sage: J = random_eja() sage: actual = J.one().operator().matrix() sage: expected = matrix.identity(J.base_ring(), J.dimension()) @@ -1240,7 +1227,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): :: - sage: set_random_seed() sage: J = random_eja(field=QQ, orthonormalize=False) sage: actual = J.one().operator().matrix() sage: expected = matrix.identity(J.base_ring(), J.dimension()) @@ -1256,7 +1242,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Ensure that the cached unit element (often precomputed by hand) agrees with the computed one:: - sage: set_random_seed() sage: J = random_eja() sage: cached = J.one() sage: J.one.clear_cache() @@ -1265,7 +1250,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): :: - sage: set_random_seed() sage: J = random_eja(field=QQ, orthonormalize=False) sage: cached = J.one() sage: J.one.clear_cache() @@ -1281,7 +1265,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # # Of course, matrices aren't vectors in sage, so we have to # appeal to the "long vectors" isometry. - oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ] + + V = VectorSpace(self.base_ring(), self.dimension()**2) + oper_vecs = [ V(g.operator().matrix().list()) for g in self.gens() ] # Now we use basic linear algebra to find the coefficients, # of the matrices-as-vectors-linear-combination, which should @@ -1291,7 +1277,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # We used the isometry on the left-hand side already, but we # still need to do it for the right-hand side. Recall that we # wanted something that summed to the identity matrix. - b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) ) + b = V( matrix.identity(self.base_ring(), self.dimension()).list() ) # Now if there's an identity element in the algebra, this # should work. We solve on the left to avoid having to @@ -1376,7 +1362,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Every algebra decomposes trivially with respect to its identity element:: - sage: set_random_seed() sage: J = random_eja() sage: J0,J5,J1 = J.peirce_decomposition(J.one()) sage: J0.dimension() == 0 and J5.dimension() == 0 @@ -1389,7 +1374,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): elements in the two subalgebras are the projections onto their respective subspaces of the superalgebra's identity element:: - sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: if not J.is_trivial(): @@ -1515,6 +1499,64 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): for idx in range(count) ) + def operator_polynomial_matrix(self): + r""" + Return the matrix of polynomials (over this algebra's + :meth:`coordinate_polynomial_ring`) that, when evaluated at + the basis coordinates of an element `x`, produces the basis + representation of `L_{x}`. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA) + + EXAMPLES:: + + sage: J = HadamardEJA(4) + sage: L_x = J.operator_polynomial_matrix() + sage: L_x + [X0 0 0 0] + [ 0 X1 0 0] + [ 0 0 X2 0] + [ 0 0 0 X3] + sage: x = J.one() + sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector()) + sage: L_x.subs(dict(d)) + [1 0 0 0] + [0 1 0 0] + [0 0 1 0] + [0 0 0 1] + + :: + + sage: J = JordanSpinEJA(4) + sage: L_x = J.operator_polynomial_matrix() + sage: L_x + [X0 X1 X2 X3] + [X1 X0 0 0] + [X2 0 X0 0] + [X3 0 0 X0] + sage: x = J.one() + sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector()) + sage: L_x.subs(dict(d)) + [1 0 0 0] + [0 1 0 0] + [0 0 1 0] + [0 0 0 1] + + """ + R = self.coordinate_polynomial_ring() + + def L_x_i_j(i,j): + # From a result in my book, these are the entries of the + # basis representation of L_x. + return sum( v*self.monomial(k).operator().matrix()[i,j] + for (k,v) in enumerate(R.gens()) ) + + n = self.dimension() + return matrix(R, n, n, L_x_i_j) + @cached_method def _charpoly_coefficients(self): r""" @@ -1530,7 +1572,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): The theory shows that these are all homogeneous polynomials of a known degree:: - sage: set_random_seed() sage: J = random_eja() sage: all(p.is_homogeneous() for p in J._charpoly_coefficients()) True @@ -1538,16 +1579,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): """ n = self.dimension() R = self.coordinate_polynomial_ring() - vars = R.gens() F = R.fraction_field() - def L_x_i_j(i,j): - # From a result in my book, these are the entries of the - # basis representation of L_x. - return sum( vars[k]*self.monomial(k).operator().matrix()[i,j] - for k in range(n) ) - - L_x = matrix(F, n, n, L_x_i_j) + L_x = self.operator_polynomial_matrix() r = None if self.rank.is_in_cache(): @@ -1628,7 +1662,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): positive integer rank, unless the algebra is trivial in which case its rank will be zero:: - sage: set_random_seed() sage: J = random_eja() sage: r = J.rank() sage: r in ZZ @@ -1639,7 +1672,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Ensure that computing the rank actually works, since the ranks of all simple algebras are known and will be cached by default:: - sage: set_random_seed() # long time sage: J = random_eja() # long time sage: cached = J.rank() # long time sage: J.rank.clear_cache() # long time @@ -1743,6 +1775,15 @@ class RationalBasisEJA(FiniteDimensionalEJA): check_field=False, check_axioms=False) + def rational_algebra(self): + # Using None as a flag here (rather than just assigning "self" + # to self._rational_algebra by default) feels a little bit + # more sane to me in a garbage-collected environment. + if self._rational_algebra is None: + return self + else: + return self._rational_algebra + @cached_method def _charpoly_coefficients(self): r""" @@ -1759,7 +1800,7 @@ class RationalBasisEJA(FiniteDimensionalEJA): sage: J = JordanSpinEJA(3) sage: J._charpoly_coefficients() - (X1^2 - X2^2 - X3^2, -2*X1) + (X0^2 - X1^2 - X2^2, -2*X0) sage: a0 = J._charpoly_coefficients()[0] sage: J.base_ring() Algebraic Real Field @@ -1767,18 +1808,15 @@ class RationalBasisEJA(FiniteDimensionalEJA): Algebraic Real Field """ - 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 - # basis, we might as well just use the given one. + if self.rational_algebra() is self: + # Bypass the hijinks if they won't benefit us. return super()._charpoly_coefficients() # Do the computation over the rationals. The answer will be # the same, because all we've done is a change of basis. # 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() ) + for a_i in self.rational_algebra()._charpoly_coefficients() ) # Otherwise, convert the coordinate variables back to the # deorthonormalized ones. @@ -1808,7 +1846,6 @@ class ConcreteEJA(FiniteDimensionalEJA): Our basis is normalized with respect to the algebra's inner product, unless we specify otherwise:: - sage: set_random_seed() sage: J = ConcreteEJA.random_instance() sage: all( b.norm() == 1 for b in J.gens() ) True @@ -1819,7 +1856,6 @@ class ConcreteEJA(FiniteDimensionalEJA): natural->EJA basis representation is an isometry and within the EJA the operator is self-adjoint by the Jordan axiom:: - sage: set_random_seed() sage: J = ConcreteEJA.random_instance() sage: x = J.random_element() sage: x.operator().is_self_adjoint() @@ -1893,11 +1929,11 @@ class ConcreteEJA(FiniteDimensionalEJA): return eja_class.random_instance(max_dimension, *args, **kwargs) -class MatrixEJA(FiniteDimensionalEJA): +class HermitianMatrixEJA(FiniteDimensionalEJA): @staticmethod def _denormalized_basis(A): """ - Returns a basis for the space of complex Hermitian n-by-n matrices. + Returns a basis for the given Hermitian matrix space. Why do we embed these? Basically, because all of numerical linear algebra assumes that you're working with vectors consisting of `n` @@ -1910,41 +1946,37 @@ class MatrixEJA(FiniteDimensionalEJA): sage: from mjo.hurwitz import (ComplexMatrixAlgebra, ....: QuaternionMatrixAlgebra, ....: OctonionMatrixAlgebra) - sage: from mjo.eja.eja_algebra import MatrixEJA + sage: from mjo.eja.eja_algebra import HermitianMatrixEJA TESTS:: - sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: A = MatrixSpace(QQ, n) - sage: B = MatrixEJA._denormalized_basis(A) + sage: B = HermitianMatrixEJA._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 = ComplexMatrixAlgebra(n, scalars=QQ) - sage: B = MatrixEJA._denormalized_basis(A) + sage: B = HermitianMatrixEJA._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: B = HermitianMatrixEJA._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: B = HermitianMatrixEJA._denormalized_basis(A) sage: all( M.is_hermitian() for M in B ) True @@ -2039,7 +2071,6 @@ class MatrixEJA(FiniteDimensionalEJA): # 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, @@ -2050,7 +2081,7 @@ class MatrixEJA(FiniteDimensionalEJA): self.rank.set_cache(matrix_space.nrows()) self.one.set_cache( self(matrix_space.one()) ) -class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): +class RealSymmetricEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -2083,7 +2114,6 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): The dimension of this algebra is `(n^2 + n) / 2`:: - sage: set_random_seed() sage: d = RealSymmetricEJA._max_random_instance_dimension() sage: n = RealSymmetricEJA._max_random_instance_size(d) sage: J = RealSymmetricEJA(n) @@ -2092,7 +2122,6 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): The Jordan multiplication is what we think it is:: - sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).to_matrix() @@ -2134,24 +2163,17 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): return cls(n, **kwargs) 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 - A = MatrixSpace(field, n) super().__init__(A, **kwargs) from mjo.eja.eja_cache import real_symmetric_eja_coeffs a = real_symmetric_eja_coeffs(self) if a is not None: - if self._rational_algebra is None: - self._charpoly_coefficients.set_cache(a) - else: - self._rational_algebra._charpoly_coefficients.set_cache(a) + self.rational_algebra()._charpoly_coefficients.set_cache(a) -class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): +class ComplexHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -2191,7 +2213,6 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): The dimension of this algebra is `n^2`:: - sage: set_random_seed() sage: d = ComplexHermitianEJA._max_random_instance_dimension() sage: n = ComplexHermitianEJA._max_random_instance_size(d) sage: J = ComplexHermitianEJA(n) @@ -2200,7 +2221,6 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): The Jordan multiplication is what we think it is:: - sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).to_matrix() @@ -2224,10 +2244,6 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): """ 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 - from mjo.hurwitz import ComplexMatrixAlgebra A = ComplexMatrixAlgebra(n, scalars=field) super().__init__(A, **kwargs) @@ -2235,10 +2251,7 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): from mjo.eja.eja_cache import complex_hermitian_eja_coeffs a = complex_hermitian_eja_coeffs(self) if a is not None: - if self._rational_algebra is None: - self._charpoly_coefficients.set_cache(a) - else: - self._rational_algebra._charpoly_coefficients.set_cache(a) + self.rational_algebra()._charpoly_coefficients.set_cache(a) @staticmethod def _max_random_instance_size(max_dimension): @@ -2259,7 +2272,7 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): return cls(n, **kwargs) -class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): +class QuaternionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA): r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the @@ -2284,7 +2297,6 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): The dimension of this algebra is `2*n^2 - n`:: - sage: set_random_seed() sage: d = QuaternionHermitianEJA._max_random_instance_dimension() sage: n = QuaternionHermitianEJA._max_random_instance_size(d) sage: J = QuaternionHermitianEJA(n) @@ -2293,7 +2305,6 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): The Jordan multiplication is what we think it is:: - sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).to_matrix() @@ -2317,10 +2328,6 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): """ 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 - from mjo.hurwitz import QuaternionMatrixAlgebra A = QuaternionMatrixAlgebra(n, scalars=field) super().__init__(A, **kwargs) @@ -2328,10 +2335,7 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs a = quaternion_hermitian_eja_coeffs(self) if a is not None: - if self._rational_algebra is None: - self._charpoly_coefficients.set_cache(a) - else: - self._rational_algebra._charpoly_coefficients.set_cache(a) + self.rational_algebra()._charpoly_coefficients.set_cache(a) @@ -2356,7 +2360,7 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): n = ZZ.random_element(max_size + 1) return cls(n, **kwargs) -class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): +class OctonionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA): r""" SETUP:: @@ -2446,7 +2450,7 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): @staticmethod def _max_random_instance_size(max_dimension): r""" - The maximum rank of a random QuaternionHermitianEJA. + The maximum rank of a random OctonionHermitianEJA. """ # There's certainly a formula for this, but with only four # cases to worry about, I'm not that motivated to derive it. @@ -2487,10 +2491,7 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs a = octonion_hermitian_eja_coeffs(self) if a is not None: - if self._rational_algebra is None: - self._charpoly_coefficients.set_cache(a) - else: - self._rational_algebra._charpoly_coefficients.set_cache(a) + self.rational_algebra()._charpoly_coefficients.set_cache(a) class AlbertEJA(OctonionHermitianEJA): @@ -2676,7 +2677,6 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): matrix. We opt not to orthonormalize the basis, because if we did, we would have to normalize the `s_{i}` in a similar manner:: - sage: set_random_seed() sage: n = ZZ.random_element(5) sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular') sage: B11 = matrix.identity(QQ,1) @@ -2838,7 +2838,6 @@ class JordanSpinEJA(BilinearFormEJA): Ensure that we have the usual inner product on `R^n`:: - sage: set_random_seed() sage: J = JordanSpinEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = x.inner_product(y) @@ -2949,6 +2948,7 @@ class CartesianProductEJA(FiniteDimensionalEJA): sage: from mjo.eja.eja_algebra import (random_eja, ....: CartesianProductEJA, + ....: ComplexHermitianEJA, ....: HadamardEJA, ....: JordanSpinEJA, ....: RealSymmetricEJA) @@ -2958,7 +2958,6 @@ class CartesianProductEJA(FiniteDimensionalEJA): The Jordan product is inherited from our factors and implemented by our CombinatorialFreeModule Cartesian product superclass:: - sage: set_random_seed() sage: J1 = HadamardEJA(2) sage: J2 = RealSymmetricEJA(2) sage: J = cartesian_product([J1,J2]) @@ -3060,6 +3059,28 @@ class CartesianProductEJA(FiniteDimensionalEJA): | b2 || 0 | 0 | b2 | +----++----+----+----+ + The "matrix space" of a Cartesian product always consists of + ordered pairs (or triples, or...) whose components are the + matrix spaces of its factors:: + + sage: J1 = HadamardEJA(2) + sage: J2 = ComplexHermitianEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.matrix_space() + The Cartesian product of (Full MatrixSpace of 2 by 1 dense + matrices over Algebraic Real Field, Module of 2 by 2 matrices + with entries in Algebraic Field over the scalar ring Algebraic + Real Field) + sage: J.one().to_matrix()[0] + [1] + [1] + sage: J.one().to_matrix()[1] + +---+---+ + | 1 | 0 | + +---+---+ + | 0 | 1 | + +---+---+ + TESTS: All factors must share the same base field:: @@ -3073,7 +3094,6 @@ class CartesianProductEJA(FiniteDimensionalEJA): The cached unit element is the same one that would be computed:: - sage: set_random_seed() # long time sage: J1 = random_eja() # long time sage: J2 = random_eja() # long time sage: J = cartesian_product([J1,J2]) # long time @@ -3082,8 +3102,8 @@ class CartesianProductEJA(FiniteDimensionalEJA): sage: expected = J.one() # long time sage: actual == expected # long time True - """ + Element = CartesianProductEJAElement def __init__(self, factors, **kwargs): m = len(factors) if m == 0: @@ -3101,8 +3121,13 @@ class CartesianProductEJA(FiniteDimensionalEJA): if associative: category = category.Associative() category = category.join([category, category.CartesianProducts()]) - # Compute my matrix space. This category isn't perfect, but - # is good enough for what we need to do. + # Compute my matrix space. We don't simply use the + # ``cartesian_product()`` functor here because it acts + # differently on SageMath MatrixSpaces and our custom + # MatrixAlgebras, which are CombinatorialFreeModules. We + # always want the result to be represented (and indexed) as an + # ordered tuple. 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 ) @@ -3153,6 +3178,8 @@ class CartesianProductEJA(FiniteDimensionalEJA): self._inner_product_matrix = matrix.block_diagonal( [J._inner_product_matrix for J in factors] ) + self._inner_product_matrix._cache = {'hermitian': True} + self._inner_product_matrix.set_immutable() # Building the multiplication table is a bit more tricky # because we have to embed the entries of the factors' @@ -3180,6 +3207,34 @@ class CartesianProductEJA(FiniteDimensionalEJA): ones = tuple(J.one().to_matrix() for J in factors) self.one.set_cache(self(ones)) + def _sets_keys(self): + r""" + + SETUP:: + + sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, + ....: RealSymmetricEJA) + + TESTS: + + The superclass uses ``_sets_keys()`` to implement its + ``cartesian_factors()`` method:: + + sage: J1 = RealSymmetricEJA(2, + ....: field=QQ, + ....: orthonormalize=False, + ....: prefix="a") + sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: J = cartesian_product([J1,J2]) + sage: x = sum(i*J.gens()[i] for i in range(len(J.gens()))) + sage: x.cartesian_factors() + (a1 + 2*a2, 3*b0 + 4*b1 + 5*b2 + 6*b3) + + """ + # Copy/pasted from CombinatorialFreeModule_CartesianProduct, + # but returning a tuple instead of a list. + return tuple(range(len(self.cartesian_factors()))) + def cartesian_factors(self): # Copy/pasted from CombinatorialFreeModule_CartesianProduct. return self._sets @@ -3196,65 +3251,6 @@ class CartesianProductEJA(FiniteDimensionalEJA): return cartesian_product.symbol.join("%s" % factor for factor in self._sets) - def matrix_space(self): - r""" - Return the space that our matrix basis lives in as a Cartesian - product. - - We don't simply use the ``cartesian_product()`` functor here - because it acts differently on SageMath MatrixSpaces and our - custom MatrixAlgebras, which are CombinatorialFreeModules. We - always want the result to be represented (and indexed) as - an ordered tuple. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, - ....: HadamardEJA, - ....: OctonionHermitianEJA, - ....: RealSymmetricEJA) - - EXAMPLES:: - - sage: J1 = HadamardEJA(1) - sage: J2 = RealSymmetricEJA(2) - sage: J = cartesian_product([J1,J2]) - sage: J.matrix_space() - The Cartesian product of (Full MatrixSpace of 1 by 1 dense - matrices over Algebraic Real Field, Full MatrixSpace of 2 - by 2 dense matrices over Algebraic Real Field) - - :: - - sage: J1 = ComplexHermitianEJA(1) - sage: J2 = ComplexHermitianEJA(1) - sage: J = cartesian_product([J1,J2]) - sage: J.one().to_matrix()[0] - +---+ - | 1 | - +---+ - sage: J.one().to_matrix()[1] - +---+ - | 1 | - +---+ - - :: - - sage: J1 = OctonionHermitianEJA(1) - sage: J2 = OctonionHermitianEJA(1) - sage: J = cartesian_product([J1,J2]) - sage: J.one().to_matrix()[0] - +----+ - | e0 | - +----+ - sage: J.one().to_matrix()[1] - +----+ - | e0 | - +----+ - - """ - return super().matrix_space() - @cached_method def cartesian_projection(self, i): @@ -3316,7 +3312,6 @@ class CartesianProductEJA(FiniteDimensionalEJA): The answer never changes:: - sage: set_random_seed() sage: J1 = random_eja() sage: J2 = random_eja() sage: J = cartesian_product([J1,J2]) @@ -3406,7 +3401,6 @@ class CartesianProductEJA(FiniteDimensionalEJA): The answer never changes:: - sage: set_random_seed() sage: J1 = random_eja() sage: J2 = random_eja() sage: J = cartesian_product([J1,J2]) @@ -3419,7 +3413,6 @@ class CartesianProductEJA(FiniteDimensionalEJA): produce the identity map, and mismatching them should produce the zero map:: - sage: set_random_seed() sage: J1 = random_eja() sage: J2 = random_eja() sage: J = cartesian_product([J1,J2]) @@ -3456,9 +3449,9 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA, SETUP:: - sage: from mjo.eja.eja_algebra import (HadamardEJA, + sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA, + ....: HadamardEJA, ....: JordanSpinEJA, - ....: OctonionHermitianEJA, ....: RealSymmetricEJA) EXAMPLES: @@ -3479,28 +3472,38 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA, 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 + check that all factors do indeed have a ``rational_algebra()`` method + before we construct an algebra that claims to have a rational basis:: + + sage: J1 = HadamardEJA(2) + sage: jp = lambda X,Y: X*Y + sage: ip = lambda X,Y: X[0,0]*Y[0,0] + sage: b1 = matrix(QQ, [[1]]) + sage: J2 = FiniteDimensionalEJA((b1,), jp, ip) + sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA + Euclidean Jordan algebra of dimension 1 over Algebraic Real + Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field + sage: cartesian_product([J1,J2]) # factor one is RationalBasisEJA + Traceback (most recent call last): + ... + ValueError: factor not a RationalBasisEJA """ def __init__(self, algebras, **kwargs): + if not all( hasattr(r, "rational_algebra") for r in algebras ): + raise ValueError("factor not a RationalBasisEJA") + CartesianProductEJA.__init__(self, algebras, **kwargs) - self._rational_algebra = None - if self.vector_space().base_field() is not QQ: - if all( hasattr(r, "_rational_algebra") for r in algebras ): - self._rational_algebra = cartesian_product([ - r._rational_algebra for r in algebras - ]) + @cached_method + def rational_algebra(self): + if self.base_ring() is QQ: + return self + + return cartesian_product([ + r.rational_algebra() for r in self.cartesian_factors() + ]) RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA @@ -3514,7 +3517,6 @@ def random_eja(max_dimension=None, *args, **kwargs): TESTS:: - sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False) sage: J.dimension() <= n @@ -3540,3 +3542,234 @@ def random_eja(max_dimension=None, *args, **kwargs): # if the sub-call also Decides on a cartesian product. J2 = random_eja(new_max_dimension, *args, **kwargs) return cartesian_product([J1,J2]) + + +class ComplexSkewSymmetricEJA(RationalBasisEJA, ConcreteEJA): + r""" + The skew-symmetric EJA of order `2n` described in Faraut and + Koranyi's Exercise III.1.b. It has dimension `2n^2 - n`. + + It is (not obviously) isomorphic to the QuaternionHermitianEJA of + order `n`, as can be inferred by comparing rank/dimension or + explicitly from their "characteristic polynomial of" functions, + which just so happen to align nicely. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (ComplexSkewSymmetricEJA, + ....: QuaternionHermitianEJA) + sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator + + EXAMPLES: + + This EJA is isomorphic to the quaternions:: + + sage: J = ComplexSkewSymmetricEJA(2, field=QQ, orthonormalize=False) + sage: K = QuaternionHermitianEJA(2, field=QQ, orthonormalize=False) + sage: jordan_isom_matrix = matrix.diagonal(QQ,[-1,1,1,1,1,-1]) + sage: phi = FiniteDimensionalEJAOperator(J,K,jordan_isom_matrix) + sage: all( phi(x*y) == phi(x)*phi(y) + ....: for x in J.gens() + ....: for y in J.gens() ) + True + sage: x,y = J.random_elements(2) + sage: phi(x*y) == phi(x)*phi(y) + True + + TESTS: + + Random elements should satisfy the same conditions that the basis + elements do:: + + sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ, + ....: orthonormalize=False) + sage: x,y = K.random_elements(2) + sage: z = x*y + sage: x = x.to_matrix() + sage: y = y.to_matrix() + sage: z = z.to_matrix() + sage: all( e.is_skew_symmetric() for e in (x,y,z) ) + True + sage: J = -K.one().to_matrix() + sage: all( e*J == J*e.conjugate() for e in (x,y,z) ) + True + + The power law in Faraut & Koranyi's II.7.a is satisfied. + We're in a subalgebra of theirs, but powers are still + defined the same:: + + sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ, + ....: orthonormalize=False) + sage: x = K.random_element() + sage: k = ZZ.random_element(5) + sage: actual = x^k + sage: J = -K.one().to_matrix() + sage: expected = K(-J*(J*x.to_matrix())^k) + sage: actual == expected + True + + """ + @staticmethod + def _max_random_instance_size(max_dimension): + # Obtained by solving d = 2n^2 - n, which comes from noticing + # that, in 2x2 block form, any element of this algebra has a + # free skew-symmetric top-left block, a Hermitian top-right + # block, and two bottom blocks that are determined by the top. + # The ZZ-int-ZZ thing is just "floor." + return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4)) + + @classmethod + def random_instance(cls, max_dimension=None, *args, **kwargs): + """ + Return a random instance of this type of algebra. + """ + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + return cls(n, **kwargs) + + @staticmethod + def _denormalized_basis(A): + """ + SETUP:: + + sage: from mjo.hurwitz import ComplexMatrixAlgebra + sage: from mjo.eja.eja_algebra import ComplexSkewSymmetricEJA + + TESTS: + + The basis elements are all skew-Hermitian:: + + sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension() + sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max) + sage: n = ZZ.random_element(n_max + 1) + sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ) + sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A) + sage: all( M.is_skew_symmetric() for M in B) + True + + The basis elements ``b`` all satisfy ``b*J == J*b.conjugate()``, + as in the definition of the algebra:: + + sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension() + sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max) + sage: n = ZZ.random_element(n_max + 1) + sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ) + sage: I_n = matrix.identity(ZZ, n) + sage: J = matrix.block(ZZ, 2, 2, (0, I_n, -I_n, 0), subdivide=False) + sage: J = A.from_list(J.rows()) + sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A) + sage: all( b*J == J*b.conjugate() for b in B ) + True + + """ + es = A.entry_algebra_gens() + gen = lambda A,m: A.monomial(m) + + basis = [] + + # The size of the blocks. We're going to treat these thing as + # 2x2 block matrices, + # + # [ x1 x2 ] + # [ -x2-conj x1-conj ] + # + # where x1 is skew-symmetric and x2 is Hermitian. + # + m = A.nrows()/2 + + # We only loop through the top half of the matrix, because the + # bottom can be constructed from the top. + for i in range(m): + # First do the top-left block, which is skew-symmetric. + # We can compute the bottom-right block in the process. + for j in range(i+1): + if i != j: + # Skew-symmetry implies zeros for (i == j). + for e in es: + # Top-left block's entry. + E_ij = gen(A, (i,j,e)) + E_ij -= gen(A, (j,i,e)) + + # Bottom-right block's entry. + F_ij = gen(A, (i+m,j+m,e)).conjugate() + F_ij -= gen(A, (j+m,i+m,e)).conjugate() + + basis.append(E_ij + F_ij) + + # Now do the top-right block, which is Hermitian, and compute + # the bottom-left block along the way. + for j in range(m,i+m+1): + if (i+m) == j: + # Hermitian matrices have real diagonal entries. + # Top-right block's entry. + E_ii = gen(A, (i,j,es[0])) + + # Bottom-left block's entry. Don't conjugate + # 'cause it's real. + E_ii -= gen(A, (i+m,j-m,es[0])) + basis.append(E_ii) + else: + for e in es: + # Top-right block's entry. BEWARE! We're not + # reflecting across the main diagonal as in + # (i,j)~(j,i). We're only reflecting across + # the diagonal for the top-right block. + E_ij = gen(A, (i,j,e)) + + # Shift it back to non-offset coords, transpose, + # conjugate, and put it back: + # + # (i,j) -> (i,j-m) -> (j-m, i) -> (j-m, i+m) + E_ij += gen(A, (j-m,i+m,e)).conjugate() + + # Bottom-left's block's below-diagonal entry. + # Just shift the top-right coords down m and + # left m. + F_ij = -gen(A, (i+m,j-m,e)).conjugate() + F_ij += -gen(A, (j,i,e)) # double-conjugate cancels + + basis.append(E_ij + F_ij) + + return tuple( basis ) + + @staticmethod + @cached_method + def _J_matrix(matrix_space): + n = matrix_space.nrows() // 2 + F = matrix_space.base_ring() + I_n = matrix.identity(F, n) + J = matrix.block(F, 2, 2, (0, I_n, -I_n, 0), subdivide=False) + return matrix_space.from_list(J.rows()) + + def J_matrix(self): + return ComplexSkewSymmetricEJA._J_matrix(self.matrix_space()) + + def __init__(self, n, field=AA, **kwargs): + # New code; always check the axioms. + #if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + + from mjo.hurwitz import ComplexMatrixAlgebra + A = ComplexMatrixAlgebra(2*n, scalars=field) + J = ComplexSkewSymmetricEJA._J_matrix(A) + + def jordan_product(X,Y): + return (X*J*Y + Y*J*X)/2 + + def inner_product(X,Y): + return (X.conjugate_transpose()*Y).trace().real() + + super().__init__(self._denormalized_basis(A), + jordan_product, + inner_product, + field=field, + matrix_space=A, + **kwargs) + + # This algebra is conjectured (by me) to be isomorphic to + # the quaternion Hermitian EJA of size n, and the rank + # would follow from that. + #self.rank.set_cache(n) + self.one.set_cache( self(-J) )