X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feuclidean_jordan_algebra.py;h=da1145d1851c91cf1f7faef014ba8a519d5d3fcf;hb=90738291f084c62d009a51bba56346a83f88f950;hp=e109b52b6c698108419f241125032c6b759d4652;hpb=cd5ba80d2e11bdce629933e09e7261aa0add2321;p=sage.d.git diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py index e109b52..da1145d 100644 --- a/mjo/eja/euclidean_jordan_algebra.py +++ b/mjo/eja/euclidean_jordan_algebra.py @@ -71,6 +71,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): """ self._rank = rank self._natural_basis = natural_basis + self._multiplication_table = mult_table fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) fda.__init__(field, mult_table, @@ -86,6 +87,54 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): return fmt.format(self.degree(), self.base_ring()) + def characteristic_polynomial(self): + r = self.rank() + n = self.dimension() + + names = ['X' + str(i) for i in range(1,n+1)] + R = PolynomialRing(self.base_ring(), names) + J = FiniteDimensionalEuclideanJordanAlgebra(R, + self._multiplication_table, + rank=r) + x0 = J.zero() + c = 1 + for g in J.gens(): + x0 += c*g + c +=1 + if not x0.is_regular(): + raise ValueError("don't know a regular element") + + # Get the vector space (as opposed to module) so that + # span_of_basis() works. + V = x0.vector().parent().ambient_vector_space() + V1 = V.span_of_basis( (x0**k).vector() for k in range(r) ) + B = V1.basis() + V1.complement().basis() + W = V.span_of_basis(B) + + def e(k): + # The coordinates of e_k with respect to the basis B. + # But, the e_k are elements of B... + return identity_matrix(J.base_ring(), n).column(k-1).column() + + # A matrix implementation 1 + x = J(vector(R, R.gens())) + l1 = [column_matrix(W.coordinates((x**k).vector())) for k in range(r)] + l2 = [e(k) for k in range(r+1, n+1)] + A_of_x = block_matrix(1, n, (l1 + l2)) + xr = W.coordinates((x**r).vector()) + a = [] + for i in range(n): + A_cols = A.columns() + A_cols[i] = xr + numerator = column_matrix(A.base_ring(), A_cols).det() + denominator = A.det() + ai = numerator/denominator + a.append(ai) + + # Note: all entries past the rth should be zero. + return a + + def inner_product(self, x, y): """ The inner product associated with this Euclidean Jordan algebra. @@ -266,6 +315,49 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): return A( (self.operator_matrix()**(n-1))*self.vector() ) + def apply_univariate_polynomial(self, p): + """ + Apply the univariate polynomial ``p`` to this element. + + A priori, SageMath won't allow us to apply a univariate + polynomial to an element of an EJA, because we don't know + that EJAs are rings (they are usually not associative). Of + course, we know that EJAs are power-associative, so the + operation is ultimately kosher. This function sidesteps + the CAS to get the answer we want and expect. + + EXAMPLES:: + + sage: R = PolynomialRing(QQ, 't') + sage: t = R.gen(0) + sage: p = t^4 - t^3 + 5*t - 2 + sage: J = RealCartesianProductEJA(5) + sage: J.one().apply_univariate_polynomial(p) == 3*J.one() + True + + TESTS: + + We should always get back an element of the algebra:: + + sage: set_random_seed() + sage: p = PolynomialRing(QQ, 't').random_element() + sage: J = random_eja() + sage: x = J.random_element() + sage: x.apply_univariate_polynomial(p) in J + True + + """ + if len(p.variables()) > 1: + raise ValueError("not a univariate polynomial") + P = self.parent() + R = P.base_ring() + # Convert the coeficcients to the parent's base ring, + # because a priori they might live in an (unnecessarily) + # larger ring for which P.sum() would fail below. + cs = [ R(c) for c in p.coefficients(sparse=False) ] + return P.sum( cs[k]*(self**k) for k in range(len(cs)) ) + + def characteristic_polynomial(self): """ Return my characteristic polynomial (if I'm a regular @@ -492,8 +584,36 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): We can't use the superclass method because it relies on the algebra being associative. + + ALGORITHM: + + The usual way to do this is to check if the determinant is + zero, but we need the characteristic polynomial for the + determinant. The minimal polynomial is a lot easier to get, + so we use Corollary 2 in Chapter V of Koecher to check + whether or not the paren't algebra's zero element is a root + of this element's minimal polynomial. + + TESTS: + + The identity element is always invertible:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J.one().is_invertible() + True + + The zero element is never invertible:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J.zero().is_invertible() + False + """ - return not self.det().is_zero() + zero = self.parent().zero() + p = self.minimal_polynomial() + return not (p(zero) == zero) def is_nilpotent(self): @@ -598,14 +718,30 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): def minimal_polynomial(self): """ - EXAMPLES:: + Return the minimal polynomial of this element, + as a function of the variable `t`. + + ALGORITHM: + + We restrict ourselves to the associative subalgebra + generated by this element, and then return the minimal + polynomial of this element's operator matrix (in that + subalgebra). This works by Baes Proposition 2.3.16. + + TESTS: + + The minimal polynomial of the identity and zero elements are + always the same:: sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.degree() == x.minimal_polynomial().degree() - True + sage: J = random_eja() + sage: J.one().minimal_polynomial() + t - 1 + sage: J.zero().minimal_polynomial() + t - :: + The degree of an element is (by one definition) the degree + of its minimal polynomial:: sage: set_random_seed() sage: x = random_eja().random_element() @@ -626,31 +762,31 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): sage: y0 = y.vector()[0] sage: y_bar = y.vector()[1:] sage: actual = y.minimal_polynomial() - sage: x = SR.symbol('x', domain='real') - sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2) + sage: t = PolynomialRing(J.base_ring(),'t').gen(0) + sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2) sage: bool(actual == expected) True - """ - # The element we're going to call "minimal_polynomial()" on. - # Either myself, interpreted as an element of a finite- - # dimensional algebra, or an element of an associative - # subalgebra. - elt = None + The minimal polynomial should always kill its element:: - if self.parent().is_associative(): - elt = FiniteDimensionalAlgebraElement(self.parent(), self) - else: - V = self.span_of_powers() - assoc_subalg = self.subalgebra_generated_by() - # Mis-design warning: the basis used for span_of_powers() - # and subalgebra_generated_by() must be the same, and in - # the same order! - elt = assoc_subalg(V.coordinates(self.vector())) + sage: set_random_seed() + sage: x = random_eja().random_element() + sage: p = x.minimal_polynomial() + sage: x.apply_univariate_polynomial(p) + 0 - # Recursive call, but should work since elt lives in an - # associative algebra. - return elt.minimal_polynomial() + """ + V = self.span_of_powers() + assoc_subalg = self.subalgebra_generated_by() + # Mis-design warning: the basis used for span_of_powers() + # and subalgebra_generated_by() must be the same, and in + # the same order! + elt = assoc_subalg(V.coordinates(self.vector())) + + # We get back a symbolic polynomial in 'x' but want a real + # polynomial in 't'. + p_of_x = elt.operator_matrix().minimal_polynomial() + return p_of_x.change_variable_name('t') def natural_representation(self): @@ -853,7 +989,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): # The dimension of the subalgebra can't be greater than # the big algebra, so just put everything into a list # and let span() get rid of the excess. - V = self.vector().parent() + # + # We do the extra ambient_vector_space() in case we're messing + # with polynomials and the direct parent is a module. + V = self.vector().parent().ambient_vector_space() return V.span( (self**d).vector() for d in xrange(V.dimension()) ) @@ -922,7 +1061,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): TESTS:: sage: set_random_seed() - sage: J = eja_rn(5) + sage: J = RealCartesianProductEJA(5) sage: c = J.random_element().subalgebra_idempotent() sage: c^2 == c True @@ -1006,16 +1145,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): return (self*other).trace() -def eja_rn(dimension, field=QQ): +class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra): """ Return the Euclidean Jordan Algebra corresponding to the set `R^n` under the Hadamard 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. + EXAMPLES: This multiplication table can be verified by hand:: - sage: J = eja_rn(3) + sage: J = RealCartesianProductEJA(3) sage: e0,e1,e2 = J.gens() sage: e0*e0 e0 @@ -1031,18 +1174,21 @@ def eja_rn(dimension, field=QQ): e2 """ - # The FiniteDimensionalAlgebra constructor takes a list of - # matrices, the ith representing right multiplication by the ith - # basis element in the vector space. So if e_1 = (1,0,0), then - # right (Hadamard) multiplication of x by e_1 picks out the first - # component of x; and likewise for the ith basis element e_i. - Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i)) - for i in xrange(dimension) ] + @staticmethod + def __classcall_private__(cls, n, field=QQ): + # The FiniteDimensionalAlgebra constructor takes a list of + # matrices, the ith representing right multiplication by the ith + # basis element in the vector space. So if e_1 = (1,0,0), then + # right (Hadamard) multiplication of x by e_1 picks out the first + # component of x; and likewise for the ith basis element e_i. + Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i)) + for i in xrange(n) ] - return FiniteDimensionalEuclideanJordanAlgebra(field, - Qs, - rank=dimension) + fdeja = super(RealCartesianProductEJA, cls) + return fdeja.__classcall_private__(cls, field, Qs, rank=n) + def inner_product(self, x, y): + return _usual_ip(x,y) def random_eja(): @@ -1077,12 +1223,17 @@ def random_eja(): Euclidean Jordan algebra of degree... """ - n = ZZ.random_element(1,5) - constructor = choice([eja_rn, - JordanSpinEJA, - RealSymmetricEJA, - ComplexHermitianEJA, - QuaternionHermitianEJA]) + + # The max_n component lets us choose different upper bounds on the + # value "n" that gets passed to the constructor. This is needed + # because e.g. R^{10} is reasonable to test, while the Hermitian + # 10-by-10 quaternion matrices are not. + (constructor, max_n) = choice([(RealCartesianProductEJA, 6), + (JordanSpinEJA, 6), + (RealSymmetricEJA, 5), + (ComplexHermitianEJA, 4), + (QuaternionHermitianEJA, 3)]) + n = ZZ.random_element(1, max_n) return constructor(n, field=QQ) @@ -1679,18 +1830,11 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra): Qi[0,0] = Qi[0,0] * ~field(2) Qs.append(Qi) + # The rank of the spin algebra is two, unless we're in a + # one-dimensional ambient space (because the rank is bounded by + # the ambient dimension). fdeja = super(JordanSpinEJA, cls) - return fdeja.__classcall_private__(cls, field, Qs) - - def rank(self): - """ - Return the rank of this Jordan Spin Algebra. - - The rank of the spin algebra is two, unless we're in a - one-dimensional ambient space (because the rank is bounded by - the ambient dimension). - """ - return min(self.dimension(),2) + return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2)) def inner_product(self, x, y): return _usual_ip(x,y)