X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_element.py;h=c98d9a2b64651562928a6489ef88e03fdd2b0dc9;hb=bc02bf48592e22d034310cfffef8fb2a062c0a43;hp=d3e9a33ceba6e1fc4e58b417f3b953e2e2d1c3d7;hpb=95ae8e7b0ddca840da9631603a2f37cca888468b;p=sage.d.git diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py index d3e9a33..c98d9a2 100644 --- a/mjo/eja/eja_element.py +++ b/mjo/eja/eja_element.py @@ -1,9 +1,10 @@ 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 FiniteDimensionalEJAOperator -from mjo.eja.eja_utils import _mat2vec +from mjo.eja.eja_utils import _mat2vec, _scale class FiniteDimensionalEJAElement(IndexedFreeModuleElement): """ @@ -166,8 +167,8 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): sage: J = HadamardEJA(3) sage: p1 = J.one().characteristic_polynomial() sage: q1 = J.zero().characteristic_polynomial() - sage: e0,e1,e2 = J.gens() - sage: A = (e0 + 2*e1 + 3*e2).subalgebra_generated_by() # dim 3 + sage: b0,b1,b2 = J.gens() + sage: A = (b0 + 2*b1 + 3*b2).subalgebra_generated_by() # dim 3 sage: p2 = A.one().characteristic_polynomial() sage: q2 = A.zero().characteristic_polynomial() sage: p1 == p2 @@ -347,7 +348,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): EXAMPLES:: sage: J = JordanSpinEJA(2) - sage: e0,e1 = J.gens() sage: x = sum( J.gens() ) sage: x.det() 0 @@ -355,7 +355,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): :: sage: J = JordanSpinEJA(3) - sage: e0,e1,e2 = J.gens() sage: x = sum( J.gens() ) sage: x.det() -1 @@ -390,7 +389,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): sage: (x*y).det() == x.det()*y.det() True - The determinant in matrix algebras is just the usual determinant:: + The determinant in real matrix algebras is the usual determinant:: sage: set_random_seed() sage: X = matrix.random(QQ,3) @@ -405,21 +404,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): sage: actual2 == expected True - :: - - sage: set_random_seed() - sage: J1 = ComplexHermitianEJA(2) - sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) - sage: X = matrix.random(GaussianIntegers(), 2) - sage: X = X + X.H - sage: expected = AA(X.det()) - sage: actual1 = J1(J1.real_embed(X)).det() - sage: actual2 = J2(J2.real_embed(X)).det() - sage: expected == actual1 - True - sage: expected == actual2 - True - """ P = self.parent() r = P.rank() @@ -438,6 +422,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): return ((-1)**r)*p(*self.to_vector()) + @cached_method def inverse(self): """ Return the Jordan-multiplicative inverse of this element. @@ -482,7 +467,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): sage: JordanSpinEJA(3).zero().inverse() Traceback (most recent call last): ... - ValueError: element is not invertible + ZeroDivisionError: element is not invertible TESTS: @@ -534,40 +519,38 @@ class FiniteDimensionalEJAElement(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 +583,6 @@ class FiniteDimensionalEJAElement(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 +591,17 @@ class FiniteDimensionalEJAElement(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 +647,7 @@ class FiniteDimensionalEJAElement(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() @@ -792,7 +776,9 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): sage: J = JordanSpinEJA(5) sage: J.one().is_regular() False - sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity + sage: b0, b1, b2, b3, b4 = J.gens() + sage: b0 == J.one() + True sage: for x in J.gens(): ....: (J.one() + x).is_regular() False @@ -842,8 +828,8 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): sage: J = JordanSpinEJA(4) sage: J.one().degree() 1 - sage: e0,e1,e2,e3 = J.gens() - sage: (e0 - e1).degree() + sage: b0,b1,b2,b3 = J.gens() + sage: (b0 - b1).degree() 2 In the spin factor algebra (of rank two), all elements that @@ -909,7 +895,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): M = matrix([(self.parent().one()).to_vector()]) old_rank = 1 - # Specifying the row-reduction algorithm can e.g. help over + # 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. # @@ -1046,19 +1032,30 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): """ if self.is_zero(): - # We would generate a zero-dimensional subalgebra - # where the minimal polynomial would be constant. - # That might be correct, but only if *this* algebra - # is trivial too. - if not self.parent().is_trivial(): - # Pretty sure we know what the minimal polynomial of - # the zero operator is going to be. This ensures - # consistency of e.g. the polynomial variable returned - # in the "normal" case without us having to think about it. - return self.operator().minimal_polynomial() - + # Pretty sure we know what the minimal polynomial of + # the zero operator is going to be. This ensures + # consistency of e.g. the polynomial variable returned + # in the "normal" case without us having to think about it. + return self.operator().minimal_polynomial() + + # If we don't orthonormalize the subalgebra's basis, then the + # first two monomials in the subalgebra will be self^0 and + # self^1... assuming that self^1 is not a scalar multiple of + # self^0 (the unit element). We special case these to avoid + # having to solve a system to coerce self into the subalgebra. A = self.subalgebra_generated_by(orthonormalize=False) - return A(self).operator().minimal_polynomial() + + if A.dimension() == 1: + # Does a solve to find the scalar multiple alpha such that + # alpha*unit = self. We have to do this because the basis + # for the subalgebra will be [ self^0 ], and not [ self^1 ]! + unit = self.parent().one() + alpha = self.to_vector() / unit.to_vector() + return (unit.operator()*alpha).minimal_polynomial() + else: + # If the dimension of the subalgebra is >= 2, then we just + # use the second basis element. + return A.monomial(1).operator().minimal_polynomial() @@ -1076,44 +1073,65 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): SETUP:: sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, - ....: QuaternionHermitianEJA) + ....: HadamardEJA, + ....: QuaternionHermitianEJA, + ....: RealSymmetricEJA) EXAMPLES:: sage: J = ComplexHermitianEJA(3) sage: J.one() - e0 + e3 + e8 + b0 + b3 + b8 sage: J.one().to_matrix() - [1 0 0 0 0 0] - [0 1 0 0 0 0] - [0 0 1 0 0 0] - [0 0 0 1 0 0] - [0 0 0 0 1 0] - [0 0 0 0 0 1] + +---+---+---+ + | 1 | 0 | 0 | + +---+---+---+ + | 0 | 1 | 0 | + +---+---+---+ + | 0 | 0 | 1 | + +---+---+---+ :: sage: J = QuaternionHermitianEJA(2) sage: J.one() - e0 + e5 + b0 + b5 sage: J.one().to_matrix() - [1 0 0 0 0 0 0 0] - [0 1 0 0 0 0 0 0] - [0 0 1 0 0 0 0 0] - [0 0 0 1 0 0 0 0] - [0 0 0 0 1 0 0 0] - [0 0 0 0 0 1 0 0] - [0 0 0 0 0 0 1 0] - [0 0 0 0 0 0 0 1] + +---+---+ + | 1 | 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 hasattr(W, 'cartesian_factors'): + # Aaaaand linear combinations don't work in Cartesian + # product spaces, even though they provide a method with + # that name. This is hidden behind an "if" because the + # _scale() function is slow. + pairs = zip(B, self.to_vector()) + return W.sum( _scale(b, alpha) for (b,alpha) in pairs ) + 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): @@ -1320,11 +1338,11 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): sage: J = RealSymmetricEJA(3) sage: J.one() - e0 + e2 + e5 + b0 + b2 + b5 sage: J.one().spectral_decomposition() - [(1, e0 + e2 + e5)] + [(1, b0 + b2 + b5)] sage: J.zero().spectral_decomposition() - [(0, e0 + e2 + e5)] + [(0, b0 + b2 + b5)] TESTS:: @@ -1349,13 +1367,13 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): The spectral decomposition should work in subalgebras, too:: sage: J = RealSymmetricEJA(4) - sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens() - sage: A = 2*e5 - 2*e8 + sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens() + sage: A = 2*b5 - 2*b8 sage: (lambda1, c1) = A.spectral_decomposition()[1] sage: (J0, J5, J1) = J.peirce_decomposition(c1) sage: (f0, f1, f2) = J1.gens() sage: f0.spectral_decomposition() - [(0, f2), (1, f0)] + [(0, c2), (1, c0)] """ A = self.subalgebra_generated_by(orthonormalize=True) @@ -1378,7 +1396,20 @@ class FiniteDimensionalEJAElement(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: @@ -1411,8 +1442,14 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement): True """ - from mjo.eja.eja_element_subalgebra import FiniteDimensionalEJAElementSubalgebra - return FiniteDimensionalEJAElementSubalgebra(self, **kwargs) + 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): @@ -1519,6 +1556,15 @@ class FiniteDimensionalEJAElement(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()