X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=c35bf256453cfa9a3758c034a47a043745eb0c2b;hb=7dc61e259969f86d4f31706fcfe411e1677389f2;hp=83ec50e5a215417811fe09d72c61251081772160;hpb=e4d568e25c62d79a2dbe34b81ee4dc21edf09316;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 83ec50e..c35bf25 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -737,7 +737,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): if self.is_trivial(): return MatrixSpace(self.base_ring(), 0) else: - return self._matrix_basis[0].matrix_space() + return self.matrix_basis()[0].parent() @cached_method @@ -1101,6 +1101,21 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): r""" The `r` polynomial coefficients of the "characteristic polynomial of" function. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS: + + 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 + """ n = self.dimension() R = self.coordinate_polynomial_ring() @@ -1136,10 +1151,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # The theory says that only the first "r" coefficients are # nonzero, and they actually live in the original polynomial - # ring and not the fraction field. We negate them because - # in the actual characteristic polynomial, they get moved - # to the other side where x^r lives. - return -A_rref.solve_right(E*b).change_ring(R)[:r] + # ring and not the fraction field. We negate them because in + # the actual characteristic polynomial, they get moved to the + # other side where x^r lives. We don't bother to trim A_rref + # down to a square matrix and solve the resulting system, + # because the upper-left r-by-r portion of A_rref is + # guaranteed to be the identity matrix, so e.g. + # + # A_rref.solve_right(Y) + # + # would just be returning Y. + return (-E*b)[:r].change_ring(R) @cached_method def rank(self): @@ -1200,7 +1222,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: set_random_seed() # long time sage: J = random_eja() # long time - sage: caches = J.rank() # long time + sage: cached = J.rank() # long time sage: J.rank.clear_cache() # long time sage: J.rank() == cached # long time True @@ -1667,6 +1689,38 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): class ComplexMatrixEJA(MatrixEJA): + # 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 @@ -1721,9 +1775,10 @@ class ComplexMatrixEJA(MatrixEJA): blocks = [] for z in M.list(): - a = z.list()[0] # real part, I guess - b = z.list()[1] # imag part, I guess - blocks.append(matrix(field, 2, [[a,b],[-b,a]])) + a = z.real() + b = z.imag() + blocks.append(matrix(field, 2, [ [ a, b], + [-b, a] ])) return matrix.block(field, n, blocks) @@ -1762,26 +1817,7 @@ class ComplexMatrixEJA(MatrixEJA): super(ComplexMatrixEJA,cls).real_unembed(M) n = ZZ(M.nrows()) d = cls.dimension_over_reals() - - # If "M" was normalized, its base ring might have roots - # adjoined and they can stick around after unembedding. - field = M.base_ring() - R = PolynomialRing(field, 'z') - z = R.gen() - - # Sage doesn't know how to adjoin the complex "i" (the root of - # x^2 + 1) to a field in a general way. Here, we just enumerate - # all of the cases that I have cared to support so far. - if field is AA: - # Sage doesn't know how to embed AA into QQbar, i.e. how - # to adjoin sqrt(-1) to AA. - F = QQbar - elif not field.is_exact(): - # RDF or RR - F = field.complex_field() - else: - # Works for QQ and... maybe some other fields. - F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) + F = cls.complex_extension(M.base_ring()) i = F.gen() # Go top-left to bottom-right (reading order), converting every @@ -1877,7 +1913,6 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: field = QuadraticField(2, 'sqrt2') sage: B = ComplexHermitianEJA._denormalized_basis(n) sage: all( M.is_symmetric() for M in B) True @@ -1895,18 +1930,27 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): # * 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): - Eij = matrix(F, n, lambda k,l: k==i and l==j) + # "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. - Sij_real = cls.real_embed(Eij + Eij.transpose()) + Eij[j,i] = 1 # Eij = Eij + Eij.transpose() + Sij_real = cls.real_embed(Eij) S.append(Sij_real) - Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose()) + # 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 these, we can drop back to the "field" that we # started with instead of the complex extension "F". @@ -1942,6 +1986,25 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): return cls(n, **kwargs) class QuaternionMatrixEJA(MatrixEJA): + + # 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 @@ -2046,8 +2109,7 @@ class QuaternionMatrixEJA(MatrixEJA): # Use the base ring of the matrix to ensure that its entries can be # multiplied by elements of the quaternion algebra. - field = M.base_ring() - Q = QuaternionAlgebra(field,-1,-1) + Q = cls.quaternion_extension(M.base_ring()) i,j,k = Q.gens() # Go top-left to bottom-right (reading order), converting every @@ -2162,23 +2224,39 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): # * 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): - Eij = matrix(Q, n, lambda k,l: k==i and l==j) + # "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. - Sij_real = cls.real_embed(Eij + Eij.transpose()) + # Eij = Eij + Eij.transpose() + Eij[j,i] = 1 + Sij_real = cls.real_embed(Eij) S.append(Sij_real) - Sij_I = cls.real_embed(I*Eij - I*Eij.transpose()) + # Eij = I*(Eij - Eij.transpose()) + Eij[i,j] = I + Eij[j,i] = -I + Sij_I = cls.real_embed(Eij) S.append(Sij_I) - Sij_J = cls.real_embed(J*Eij - J*Eij.transpose()) + # Eij = J*(Eij - Eij.transpose()) + Eij[i,j] = J + Eij[j,i] = -J + Sij_J = cls.real_embed(Eij) S.append(Sij_J) - Sij_K = cls.real_embed(K*Eij - K*Eij.transpose()) + # 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 these, we can drop back to the "field" that we # started with instead of the quaternion algebra "Q". @@ -2596,100 +2674,119 @@ class TrivialEJA(ConcreteEJA): # inappropriate for us. return cls(**kwargs) -# class DirectSumEJA(ConcreteEJA): -# r""" -# The external (orthogonal) direct sum of two other Euclidean Jordan -# algebras. Essentially the Cartesian product of its two factors. -# Every Euclidean Jordan algebra decomposes into an orthogonal -# direct sum of simple Euclidean Jordan algebras, so no generality -# is lost by providing only this construction. - -# SETUP:: - -# sage: from mjo.eja.eja_algebra import (random_eja, -# ....: HadamardEJA, -# ....: RealSymmetricEJA, -# ....: DirectSumEJA) - -# EXAMPLES:: - -# sage: J1 = HadamardEJA(2) -# sage: J2 = RealSymmetricEJA(3) -# sage: J = DirectSumEJA(J1,J2) -# sage: J.dimension() -# 8 -# sage: J.rank() -# 5 - -# TESTS: - -# The external direct sum construction is only valid when the two factors -# have the same base ring; an error is raised otherwise:: - -# sage: set_random_seed() -# sage: J1 = random_eja(field=AA) -# sage: J2 = random_eja(field=QQ,orthonormalize=False) -# sage: J = DirectSumEJA(J1,J2) -# Traceback (most recent call last): -# ... -# ValueError: algebras must share the same base field - -# """ -# def __init__(self, J1, J2, **kwargs): -# if J1.base_ring() != J2.base_ring(): -# raise ValueError("algebras must share the same base field") -# field = J1.base_ring() - -# self._factors = (J1, J2) -# n1 = J1.dimension() -# n2 = J2.dimension() -# n = n1+n2 -# V = VectorSpace(field, n) -# mult_table = [ [ V.zero() for j in range(i+1) ] -# for i in range(n) ] -# for i in range(n1): -# for j in range(i+1): -# p = (J1.monomial(i)*J1.monomial(j)).to_vector() -# mult_table[i][j] = V(p.list() + [field.zero()]*n2) - -# for i in range(n2): -# for j in range(i+1): -# p = (J2.monomial(i)*J2.monomial(j)).to_vector() -# mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list()) - -# # TODO: build the IP table here from the two constituent IP -# # matrices (it'll be block diagonal, I think). -# ip_table = [ [ field.zero() for j in range(i+1) ] -# for i in range(n) ] -# super(DirectSumEJA, self).__init__(field, -# mult_table, -# ip_table, -# check_axioms=False, -# **kwargs) -# self.rank.set_cache(J1.rank() + J2.rank()) - - -# def factors(self): -# r""" -# Return the pair of this algebra's factors. -# SETUP:: +class DirectSumEJA(FiniteDimensionalEJA): + r""" + The external (orthogonal) direct sum of two other Euclidean Jordan + algebras. Essentially the Cartesian product of its two factors. + Every Euclidean Jordan algebra decomposes into an orthogonal + direct sum of simple Euclidean Jordan algebras, so no generality + is lost by providing only this construction. -# sage: from mjo.eja.eja_algebra import (HadamardEJA, -# ....: JordanSpinEJA, -# ....: DirectSumEJA) + SETUP:: -# EXAMPLES:: + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: RealSymmetricEJA, + ....: DirectSumEJA) -# sage: J1 = HadamardEJA(2, field=QQ) -# sage: J2 = JordanSpinEJA(3, field=QQ) -# sage: J = DirectSumEJA(J1,J2) -# sage: J.factors() -# (Euclidean Jordan algebra of dimension 2 over Rational Field, -# Euclidean Jordan algebra of dimension 3 over Rational Field) + EXAMPLES:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(3) + sage: J = DirectSumEJA(J1,J2) + sage: J.dimension() + 8 + sage: J.rank() + 5 + sage: J.matrix_space() + The Cartesian product of (Full MatrixSpace of 2 by 1 dense matrices + over Algebraic Real Field, Full MatrixSpace of 3 by 3 dense matrices + over Algebraic Real Field) + + TESTS: + + The external direct sum construction is only valid when the two factors + have the same base ring; an error is raised otherwise:: + + sage: set_random_seed() + sage: J1 = random_eja(field=AA) + sage: J2 = random_eja(field=QQ,orthonormalize=False) + sage: J = DirectSumEJA(J1,J2) + Traceback (most recent call last): + ... + ValueError: algebras must share the same base field + + """ + def __init__(self, J1, J2, **kwargs): + if J1.base_ring() != J2.base_ring(): + raise ValueError("algebras must share the same base field") + field = J1.base_ring() + + M = J1.matrix_space().cartesian_product(J2.matrix_space()) + self._cartprod_algebra = J1.cartesian_product(J2) + + self._matrix_basis = tuple( [M((a,0)) for a in J1.matrix_basis()] + + [M((0,b)) for b in J2.matrix_basis()] ) + + n = len(self._matrix_basis) + self._sets = None + CombinatorialFreeModule.__init__( + self, + field, + range(n), + category=self._cartprod_algebra.category(), + bracket=False, + **kwargs) + self.rank.set_cache(J1.rank() + J2.rank()) + + + + def product(self,x,y): + r""" + SETUP:: + + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: ComplexHermitianEJA, + ....: DirectSumEJA) + + TESTS:: + + sage: set_random_seed() + sage: J1 = JordanSpinEJA(3, field=QQ) + sage: J2 = ComplexHermitianEJA(2, field=QQ, orthonormalize=False) + sage: J = DirectSumEJA(J1,J2) + sage: J.random_element()*J.random_element() in J + True + + """ + xv = self._cartprod_algebra.from_vector(x.to_vector()) + yv = self._cartprod_algebra.from_vector(y.to_vector()) + return self.from_vector((xv*yv).to_vector()) + + + def cartesian_factors(self): + r""" + Return the pair of this algebra's factors. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = JordanSpinEJA(3, field=QQ) + sage: J = DirectSumEJA(J1,J2) + sage: J.cartesian_factors() + (Euclidean Jordan algebra of dimension 2 over Rational Field, + Euclidean Jordan algebra of dimension 3 over Rational Field) + + """ + return self._cartprod_algebra.cartesian_factors() -# """ -# return self._factors # def projections(self): # r"""