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_operator import FiniteDimensionalEJAOperator
from mjo.eja.eja_utils import _mat2vec
-class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
+class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
"""
An element of a Euclidean Jordan algebra.
"""
return ((-1)**r)*p(*self.to_vector())
+ @cached_method
def inverse(self):
"""
Return the Jordan-multiplicative inverse of this element.
sage: JordanSpinEJA(3).zero().inverse()
Traceback (most recent call last):
...
- ValueError: element is not invertible
+ ZeroDivisionError: element is not invertible
TESTS:
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::
sage: fast = x.is_invertible() # long time
sage: slow == fast # long time
True
-
"""
if self.is_zero():
if self.parent().is_trivial():
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):
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::
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):
# 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()
return W.linear_combination( zip(B, self.to_vector()) )
+
def norm(self):
"""
The norm of this element with respect to :meth:`inner_product`.
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):
[(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.
True
"""
- from mjo.eja.eja_element_subalgebra import FiniteDimensionalEuclideanJordanElementSubalgebra
- return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
+ from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+ powers = tuple( self**k for k in range(self.degree()) )
+ A = FiniteDimensionalEJASubalgebra(self.parent(),
+ powers,
+ associative=True,
+ **kwargs)
+ A.one.set_cache(A(self.parent().one()))
+ return A
def subalgebra_idempotent(self):
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()
"""
return self.trace_inner_product(self).sqrt()
+
+
+
+class CartesianProductEJAElement(FiniteDimensionalEJAElement):
+
+ def to_matrix(self):
+ r"""
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (HadamardEJA,
+ ....: RealSymmetricEJA)
+
+ EXAMPLES::
+
+ 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()
+
+ # Aaaaand linear combinations don't work in Cartesian
+ # product spaces, even though they provide a method
+ # with that name.
+ pairs = zip(B, self.to_vector())
+ sigma = W.zero()
+ for (b,alpha) in pairs:
+ # sum(...) ALSO doesn't work on Cartesian products.
+ sigma += W(tuple(alpha*b_i for b_i in b))
+ return sigma
projects / sage.d.git / blobdiff