X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;ds=sidebyside;f=mjo%2Feja%2Feuclidean_jordan_algebra.py;h=2c496cd3ed893b3fd34f82f6e3ce730c757e0926;hb=ab8ae4abd5b5414d57059b985259337b58528793;hp=29248fcfaabfad33567070eaa65efffef7df87ff;hpb=119af99174a49d71d1b11fdfd1b0de451fee733f;p=sage.d.git diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py index 29248fc..2c496cd 100644 --- a/mjo/eja/euclidean_jordan_algebra.py +++ b/mjo/eja/euclidean_jordan_algebra.py @@ -69,6 +69,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): True """ + self._charpoly = None # for caching self._rank = rank self._natural_basis = natural_basis self._multiplication_table = mult_table @@ -87,29 +88,34 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): return fmt.format(self.degree(), self.base_ring()) + def characteristic_polynomial(self): + """ + EXAMPLES: + + The characteristic polynomial in the spin algebra is given in + Alizadeh, Example 11.11:: + + sage: J = JordanSpinEJA(3) + sage: p = J.characteristic_polynomial(); p + X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 + sage: xvec = J.one().vector() + sage: p(*xvec) + t^2 - 2*t + 1 + + """ + if self._charpoly is not None: + return self._charpoly + r = self.rank() n = self.dimension() + # Now switch to the polynomial rings. 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. @@ -118,21 +124,46 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): # A matrix implementation 1 x = J(vector(R, R.gens())) - l1 = [column_matrix(W.coordinates((x**k).vector())) for k in range(r)] + l1 = [column_matrix((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()) + xr = (x**r).vector() a = [] + denominator = A_of_x.det() # This is constant for i in range(n): - A_cols = A.columns() + A_cols = A_of_x.columns() A_cols[i] = xr - numerator = column_matrix(A.base_ring(), A_cols).det() - denominator = A.det() + numerator = column_matrix(A_of_x.base_ring(), A_cols).det() ai = numerator/denominator a.append(ai) - # Note: all entries past the rth should be zero. - return a + # We go to a bit of trouble here to reorder the + # indeterminates, so that it's easier to evaluate the + # characteristic polynomial at x's coordinates and get back + # something in terms of t, which is what we want. + S = PolynomialRing(self.base_ring(),'t') + t = S.gen(0) + S = PolynomialRing(S, R.variable_names()) + t = S(t) + + # We're relying on the theory here to ensure that each entry + # a[i] is indeed back in R, and the added negative signs are + # to make the whole expression sum to zero. + a = [R(-ai) for ai in a] # corresponds to powerx x^0 through x^(r-1) + + # Note: all entries past the rth should be zero. The + # coefficient of the highest power (x^r) is 1, but it doesn't + # appear in the solution vector which contains coefficients + # for the other powers (to make them sum to x^r). + if (r < n): + a[r] = 1 # corresponds to x^r + else: + # When the rank is equal to the dimension, trying to + # assign a[r] goes out-of-bounds. + a.append(1) # corresponds to x^r + + self._charpoly = sum( a[k]*(t**k) for k in range(len(a)) ) + return self._charpoly def inner_product(self, x, y): @@ -315,19 +346,82 @@ 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 - element). + Return the characteristic polynomial of this element. + + EXAMPLES: + + The rank of `R^3` is three, and the minimal polynomial of + the identity element is `(t-1)` from which it follows that + the characteristic polynomial should be `(t-1)^3`:: + + sage: J = RealCartesianProductEJA(3) + sage: J.one().characteristic_polynomial() + t^3 - 3*t^2 + 3*t - 1 + + Likewise, the characteristic of the zero element in the + rank-three algebra `R^{n}` should be `t^{3}`:: + + sage: J = RealCartesianProductEJA(3) + sage: J.zero().characteristic_polynomial() + t^3 + + The characteristic polynomial of an element should evaluate + to zero on that element:: + + sage: set_random_seed() + sage: x = RealCartesianProductEJA(3).random_element() + sage: p = x.characteristic_polynomial() + sage: x.apply_univariate_polynomial(p) + 0 - Eventually this should be implemented in terms of the parent - algebra's characteristic polynomial that works for ALL - elements. """ - if self.is_regular(): - return self.minimal_polynomial() - else: - raise NotImplementedError('irregular element') + p = self.parent().characteristic_polynomial() + return p(*self.vector()) def inner_product(self, other): @@ -675,6 +769,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): def minimal_polynomial(self): """ + Return the minimal polynomial of this element, + as a function of the variable `t`. + ALGORITHM: We restrict ourselves to the associative subalgebra @@ -682,14 +779,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): polynomial of this element's operator matrix (in that subalgebra). This works by Baes Proposition 2.3.16. - EXAMPLES:: + 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() @@ -710,11 +813,19 @@ 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 minimal polynomial should always kill its element:: + + sage: set_random_seed() + sage: x = random_eja().random_element() + sage: p = x.minimal_polynomial() + sage: x.apply_univariate_polynomial(p) + 0 + """ V = self.span_of_powers() assoc_subalg = self.subalgebra_generated_by() @@ -722,7 +833,11 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): # and subalgebra_generated_by() must be the same, and in # the same order! elt = assoc_subalg(V.coordinates(self.vector())) - return elt.operator_matrix().minimal_polynomial() + + # 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): @@ -1159,12 +1274,17 @@ def random_eja(): Euclidean Jordan algebra of degree... """ - n = ZZ.random_element(1,5) - constructor = choice([RealCartesianProductEJA, - 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)