from sage.modules.free_module import VectorSpace
from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
-# TODO: make this unnecessary somehow.
-from sage.misc.lazy_import import lazy_import
-lazy_import('mjo.eja.eja_algebra', 'FiniteDimensionalEuclideanJordanAlgebra')
-lazy_import('mjo.eja.eja_element_subalgebra',
- 'FiniteDimensionalEuclideanJordanElementSubalgebra')
-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.
"""
True
"""
- p = self.parent().characteristic_polynomial()
+ p = self.parent().characteristic_polynomial_of()
return p(*self.to_vector())
Ditto for the quaternions::
- sage: J = QuaternionHermitianEJA(3)
+ sage: J = QuaternionHermitianEJA(2)
sage: J.one().inner_product(J.one())
- 3
+ 2
TESTS:
sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
....: TrivialEJA,
+ ....: RealSymmetricEJA,
+ ....: ComplexHermitianEJA,
....: random_eja)
EXAMPLES::
sage: x,y = J.random_elements(2)
sage: (x*y).det() == x.det()*y.det()
True
+
+ The determinant in matrix algebras is just the usual determinant::
+
+ sage: set_random_seed()
+ sage: X = matrix.random(QQ,3)
+ sage: X = X + X.T
+ sage: J1 = RealSymmetricEJA(3)
+ sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
+ sage: expected = X.det()
+ sage: actual1 = J1(X).det()
+ sage: actual2 = J2(X).det()
+ sage: actual1 == expected
+ True
+ 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()
ALGORITHM:
- We appeal to the quadratic representation as in Koecher's
- Theorem 12 in Chapter III, Section 5.
+ In general we appeal to the quadratic representation as in
+ Koecher's Theorem 12 in Chapter III, Section 5. But if the
+ parent algebra's "characteristic polynomial of" coefficients
+ happen to be cached, then we use Proposition II.2.4 in Faraut
+ and Korányi which gives a formula for the inverse based on the
+ characteristic polynomial and the Cayley-Hamilton theorem for
+ Euclidean Jordan algebras::
SETUP::
....: x.operator().inverse()(J.one()) == x.inverse() )
True
- Proposition II.2.4 in Faraut and Korányi gives a formula for
- the inverse based on the characteristic polynomial and the
- Cayley-Hamilton theorem for Euclidean Jordan algebras::
+ Check that the fast (cached) and slow algorithms give the same
+ answer::
- sage: set_random_seed()
- sage: J = ComplexHermitianEJA(3)
- sage: x = J.random_element()
- sage: while not x.is_invertible():
- ....: x = J.random_element()
- sage: r = J.rank()
- sage: a = x.characteristic_polynomial().coefficients(sparse=False)
- sage: expected = (-1)^(r+1)/x.det()
- sage: expected *= sum( a[i+1]*x^i for i in range(r) )
- sage: x.inverse() == expected
+ sage: set_random_seed() # long time
+ sage: J = random_eja(field=QQ, orthonormalize=False) # long time
+ sage: x = J.random_element() # long time
+ sage: while not x.is_invertible(): # long time
+ ....: x = J.random_element() # long time
+ sage: slow = x.inverse() # long time
+ sage: _ = J._charpoly_coefficients() # long time
+ sage: fast = x.inverse() # long time
+ sage: slow == fast # long time
True
-
"""
if not self.is_invertible():
raise ValueError("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!
+ 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)
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
+ 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.
sage: (not J.is_trivial()) and J.zero().is_invertible()
False
+ Test that the fast (cached) and slow algorithms give the same
+ answer::
+
+ sage: set_random_seed() # long time
+ sage: J = random_eja(field=QQ, orthonormalize=False) # long time
+ sage: x = J.random_element() # long time
+ sage: slow = x.is_invertible() # long time
+ sage: _ = J._charpoly_coefficients() # long time
+ sage: fast = x.is_invertible() # long time
+ sage: slow == fast # long time
+ True
+
"""
if self.is_zero():
if self.parent().is_trivial():
else:
return False
+ if self.parent()._charpoly_coefficients.is_in_cache():
+ # The determinant will be quicker than computing the minimal
+ # polynomial from scratch, 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()
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):
two here so that said elements actually exist::
sage: set_random_seed()
- sage: n_max = max(2, JordanSpinEJA._max_test_case_size())
+ sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())
sage: n = ZZ.random_element(2, n_max)
sage: J = JordanSpinEJA(n)
sage: y = J.random_element()
and in particular, a re-scaling of the basis::
sage: set_random_seed()
- sage: n_max = RealSymmetricEJA._max_test_case_size()
+ sage: n_max = RealSymmetricEJA._max_random_instance_size()
sage: n = ZZ.random_element(1, n_max)
sage: J1 = RealSymmetricEJA(n)
- sage: J2 = RealSymmetricEJA(n,normalize_basis=False)
+ sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
sage: X = random_matrix(AA,n)
sage: X = X*X.transpose()
sage: x1 = J1(X)
# in the "normal" case without us having to think about it.
return self.operator().minimal_polynomial()
- A = self.subalgebra_generated_by()
+ A = self.subalgebra_generated_by(orthonormalize_basis=False)
return A(self).operator().minimal_polynomial()
- def natural_representation(self):
+ def to_matrix(self):
"""
- Return a more-natural representation of this element.
+ Return an (often more natural) representation of this element as a
+ matrix.
- Every finite-dimensional Euclidean Jordan Algebra is a
- direct sum of five simple algebras, four of which comprise
- Hermitian matrices. This method returns the original
- "natural" representation of this element as a Hermitian
- matrix, if it has one. If not, you get the usual representation.
+ Every finite-dimensional Euclidean Jordan Algebra is a direct
+ sum of five simple algebras, four of which comprise Hermitian
+ matrices. This method returns a "natural" matrix
+ representation of this element as either a Hermitian matrix or
+ column vector.
SETUP::
sage: J = ComplexHermitianEJA(3)
sage: J.one()
e0 + e3 + e8
- sage: J.one().natural_representation()
+ sage: J.one().to_matrix()
[1 0 0 0 0 0]
[0 1 0 0 0 0]
[0 0 1 0 0 0]
::
- sage: J = QuaternionHermitianEJA(3)
+ sage: J = QuaternionHermitianEJA(2)
sage: J.one()
- e0 + e5 + e14
- sage: J.one().natural_representation()
- [1 0 0 0 0 0 0 0 0 0 0 0]
- [0 1 0 0 0 0 0 0 0 0 0 0]
- [0 0 1 0 0 0 0 0 0 0 0 0]
- [0 0 0 1 0 0 0 0 0 0 0 0]
- [0 0 0 0 1 0 0 0 0 0 0 0]
- [0 0 0 0 0 1 0 0 0 0 0 0]
- [0 0 0 0 0 0 1 0 0 0 0 0]
- [0 0 0 0 0 0 0 1 0 0 0 0]
- [0 0 0 0 0 0 0 0 1 0 0 0]
- [0 0 0 0 0 0 0 0 0 1 0 0]
- [0 0 0 0 0 0 0 0 0 0 1 0]
- [0 0 0 0 0 0 0 0 0 0 0 1]
+ e0 + e5
+ 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]
"""
- B = self.parent().natural_basis()
- W = self.parent().natural_basis_space()
- return W.linear_combination(zip(B,self.to_vector()))
+ 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()) )
def norm(self):
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):
sage: (J0, J5, J1) = J.peirce_decomposition(c1)
sage: (f0, f1, f2) = J1.gens()
sage: f0.spectral_decomposition()
- [(0, 1.000000000000000?*f2), (1, 1.000000000000000?*f0)]
+ [(0, f2), (1, f0)]
"""
A = self.subalgebra_generated_by(orthonormalize_basis=True)
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
"""
- return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
+ from mjo.eja.eja_element_subalgebra import FiniteDimensionalEJAElementSubalgebra
+ return FiniteDimensionalEJAElementSubalgebra(self, **kwargs)
def subalgebra_idempotent(self):