X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_element.py;h=e8c183b51d522efd14aba57b779b4b5ebfa4b896;hb=45f207f28a8396426469fedb026b4da82e30fbf5;hp=6812e2807a7ff3c3a4f7253c7af9750f5de8fbb2;hpb=0ecada5937cb3c3aa7307029f165996d83a16133;p=sage.d.git diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py index 6812e28..e8c183b 100644 --- a/mjo/eja/eja_element.py +++ b/mjo/eja/eja_element.py @@ -1,11 +1,12 @@ from sage.matrix.constructor import matrix +from sage.misc.cachefunc import cached_method from sage.modules.free_module import VectorSpace from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement -from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator -from mjo.eja.eja_utils import _mat2vec +from mjo.eja.eja_operator import FiniteDimensionalEJAOperator +from mjo.eja.eja_utils import _mat2vec, _scale -class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): +class FiniteDimensionalEJAElement(IndexedFreeModuleElement): """ An element of a Euclidean Jordan algebra. """ @@ -438,6 +439,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): return ((-1)**r)*p(*self.to_vector()) + @cached_method def inverse(self): """ Return the Jordan-multiplicative inverse of this element. @@ -482,7 +484,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: JordanSpinEJA(3).zero().inverse() Traceback (most recent call last): ... - ValueError: element is not invertible + ZeroDivisionError: element is not invertible TESTS: @@ -534,40 +536,38 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: slow == fast # long time True """ - if not self.is_invertible(): - raise ValueError("element is not invertible") - + not_invertible_msg = "element is not invertible" if self.parent()._charpoly_coefficients.is_in_cache(): # We can invert using our charpoly if it will be fast to # compute. If the coefficients are cached, our rank had # better be too! + if self.det().is_zero(): + raise ZeroDivisionError(not_invertible_msg) r = self.parent().rank() a = self.characteristic_polynomial().coefficients(sparse=False) return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det() - return (~self.quadratic_representation())(self) + try: + inv = (~self.quadratic_representation())(self) + self.is_invertible.set_cache(True) + return inv + except ZeroDivisionError: + self.is_invertible.set_cache(False) + raise ZeroDivisionError(not_invertible_msg) + @cached_method def is_invertible(self): """ Return whether or not this element is invertible. 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 parent algebra's zero element is a root - of this element's minimal polynomial. - - That is... unless the coefficients of our algebra's - "characteristic polynomial of" function are already cached! - In that case, we just use the determinant (which will be fast - as a result). - - Beware that we can't use the superclass method, because it - relies on the algebra being associative. + If computing my determinant will be fast, we do so and compare + with zero (Proposition II.2.4 in Faraut and + Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi + reduces the problem to the invertibility of my quadratic + representation. SETUP:: @@ -600,7 +600,6 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: fast = x.is_invertible() # long time sage: slow == fast # long time True - """ if self.is_zero(): if self.parent().is_trivial(): @@ -609,15 +608,17 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): return False if self.parent()._charpoly_coefficients.is_in_cache(): - # The determinant will be quicker than computing the minimal - # polynomial from scratch, most likely. + # The determinant will be quicker than inverting the + # quadratic representation, most likely. return (not self.det().is_zero()) - # In fact, we only need to know if the constant term is non-zero, - # so we can pass in the field's zero element instead. - zero = self.base_ring().zero() - p = self.minimal_polynomial() - return not (p(zero) == zero) + # The easiest way to determine if I'm invertible is to try. + try: + inv = (~self.quadratic_representation())(self) + self.inverse.set_cache(inv) + return True + except ZeroDivisionError: + return False def is_primitive_idempotent(self): @@ -663,7 +664,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): element should always be in terms of minimal idempotents:: sage: J = JordanSpinEJA(4) - sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) ) + sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) ) sage: x.is_regular() True sage: [ c.is_primitive_idempotent() @@ -830,10 +831,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): ALGORITHM: - For now, we skip the messy minimal polynomial computation - and instead return the dimension of the vector space spanned - by the powers of this element. The latter is a bit more - straightforward to compute. + ......... SETUP:: @@ -881,12 +879,59 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): True """ - if self.is_zero() and not self.parent().is_trivial(): + n = self.parent().dimension() + + if n == 0: + # The minimal polynomial is an empty product, i.e. the + # constant polynomial "1" having degree zero. + return 0 + elif self.is_zero(): # The minimal polynomial of zero in a nontrivial algebra - # is "t"; in a trivial algebra it's "1" by convention - # (it's an empty product). + # is "t", and is of degree one. + return 1 + elif n == 1: + # If this is a nonzero element of a nontrivial algebra, it + # has degree at least one. It follows that, in an algebra + # of dimension one, the degree must be actually one. return 1 - return self.subalgebra_generated_by().dimension() + + # BEWARE: The subalgebra_generated_by() method uses the result + # of this method to construct a basis for the subalgebra. That + # means, in particular, that we cannot implement this method + # as ``self.subalgebra_generated_by().dimension()``. + + # Algorithm: keep appending (vector representations of) powers + # self as rows to a matrix and echelonizing it. When its rank + # stops increasing, we've reached a redundancy. + + # Given the special cases above, we can assume that "self" is + # nonzero, the algebra is nontrivial, and that its dimension + # is at least two. + M = matrix([(self.parent().one()).to_vector()]) + old_rank = 1 + + # Specifying the row-reduction algorithm can e.g. help over + # AA because it avoids the RecursionError that gets thrown + # when we have to look too hard for a root. + # + # Beware: QQ supports an entirely different set of "algorithm" + # keywords than do AA and RR. + algo = None + from sage.rings.all import QQ + if self.parent().base_ring() is not QQ: + algo = "scaled_partial_pivoting" + + for d in range(1,n): + M = matrix(M.rows() + [(self**d).to_vector()]) + M.echelonize(algo) + new_rank = M.rank() + if new_rank == old_rank: + return new_rank + else: + old_rank = new_rank + + return n + def left_matrix(self): @@ -1013,7 +1058,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): # in the "normal" case without us having to think about it. return self.operator().minimal_polynomial() - A = self.subalgebra_generated_by(orthonormalize_basis=False) + A = self.subalgebra_generated_by(orthonormalize=False) return A(self).operator().minimal_polynomial() @@ -1032,7 +1077,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): SETUP:: sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, - ....: QuaternionHermitianEJA) + ....: HadamardEJA, + ....: QuaternionHermitianEJA, + ....: RealSymmetricEJA) EXAMPLES:: @@ -1062,14 +1109,36 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): [0 0 0 0 0 0 1 0] [0 0 0 0 0 0 0 1] + This also works in Cartesian product algebras:: + + sage: J1 = HadamardEJA(1) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: x = sum(J.gens()) + sage: x.to_matrix()[0] + [1] + sage: x.to_matrix()[1] + [ 1 0.7071067811865475?] + [0.7071067811865475? 1] + """ B = self.parent().matrix_basis() W = self.parent().matrix_space() - # This is just a manual "from_vector()", but of course - # matrix spaces aren't vector spaces in sage, so they - # don't have a from_vector() method. - return W.linear_combination( zip(B, self.to_vector()) ) + if self.parent()._matrix_basis_is_cartesian: + # Aaaaand linear combinations don't work in Cartesian + # product spaces, even though they provide a method + # with that name. This is special-cased because the + # _scale() function is slow. + pairs = zip(B, self.to_vector()) + return sum( ( _scale(b, alpha) for (b,alpha) in pairs ), + W.zero()) + else: + # This is just a manual "from_vector()", but of course + # matrix spaces aren't vector spaces in sage, so they + # don't have a from_vector() method. + return W.linear_combination( zip(B, self.to_vector()) ) + def norm(self): @@ -1122,10 +1191,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): P = self.parent() left_mult_by_self = lambda y: self*y L = P.module_morphism(function=left_mult_by_self, codomain=P) - return FiniteDimensionalEuclideanJordanAlgebraOperator( - P, - P, - L.matrix() ) + return FiniteDimensionalEJAOperator(P, P, L.matrix() ) def quadratic_representation(self, other=None): @@ -1317,13 +1383,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): [(0, f2), (1, f0)] """ - A = self.subalgebra_generated_by(orthonormalize_basis=True) + A = self.subalgebra_generated_by(orthonormalize=True) result = [] for (evalue, proj) in A(self).operator().spectral_decomposition(): result.append( (evalue, proj(A.one()).superalgebra_element()) ) return result - def subalgebra_generated_by(self, orthonormalize_basis=False): + def subalgebra_generated_by(self, **kwargs): """ Return the associative subalgebra of the parent EJA generated by this element. @@ -1337,7 +1403,20 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): SETUP:: - sage: from mjo.eja.eja_algebra import random_eja + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: RealSymmetricEJA) + + EXAMPLES: + + We can create subalgebras of Cartesian product EJAs that are not + themselves Cartesian product EJAs (they're just "regular" EJAs):: + + sage: J1 = HadamardEJA(3) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.one().subalgebra_generated_by() + Euclidean Jordan algebra of dimension 1 over Algebraic Real Field TESTS: @@ -1370,8 +1449,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): True """ - from mjo.eja.eja_element_subalgebra import FiniteDimensionalEuclideanJordanElementSubalgebra - return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis) + powers = tuple( self**k for k in range(self.degree()) ) + A = self.parent().subalgebra(powers, + associative=True, + check_field=False, + check_axioms=False, + **kwargs) + A.one.set_cache(A(self.parent().one())) + return A def subalgebra_idempotent(self): @@ -1478,6 +1563,15 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: J.random_element().trace() in RLF True + The trace is linear:: + + sage: set_random_seed() + sage: J = random_eja() + sage: x,y = J.random_elements(2) + sage: alpha = J.base_ring().random_element() + sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace() + True + """ P = self.parent() r = P.rank()