-# -*- coding: utf-8 -*-
-
-from itertools import izip
-
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
-# 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_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.
"""
SETUP::
- sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
+ sage: from mjo.eja.eja_algebra import (HadamardEJA,
....: random_eja)
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 = HadamardEJA(5)
sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
True
We should always get back an element of the algebra::
sage: set_random_seed()
- sage: p = PolynomialRing(QQ, 't').random_element()
+ sage: p = PolynomialRing(AA, 't').random_element()
sage: J = random_eja()
sage: x = J.random_element()
sage: x.apply_univariate_polynomial(p) in J
SETUP::
- sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+ sage: from mjo.eja.eja_algebra import HadamardEJA
EXAMPLES:
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 = HadamardEJA(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 = HadamardEJA(3)
sage: J.zero().characteristic_polynomial()
t^3
to zero on that element::
sage: set_random_seed()
- sage: x = RealCartesianProductEJA(3).random_element()
+ sage: x = HadamardEJA(3).random_element()
sage: p = x.characteristic_polynomial()
sage: x.apply_univariate_polynomial(p)
0
The characteristic polynomials of the zero and unit elements
should be what we think they are in a subalgebra, too::
- sage: J = RealCartesianProductEJA(3)
+ 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
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:
SETUP::
sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+ ....: TrivialEJA,
+ ....: RealSymmetricEJA,
+ ....: ComplexHermitianEJA,
....: random_eja)
EXAMPLES::
sage: J = JordanSpinEJA(2)
- sage: e0,e1 = J.gens()
sage: x = sum( J.gens() )
sage: x.det()
0
::
sage: J = JordanSpinEJA(3)
- sage: e0,e1,e2 = J.gens()
sage: x = sum( J.gens() )
sage: x.det()
-1
+ The determinant of the sole element in the rank-zero trivial
+ algebra is ``1``, by three paths of reasoning. First, its
+ characteristic polynomial is a constant ``1``, so the constant
+ term in that polynomial is ``1``. Second, the characteristic
+ polynomial evaluated at zero is again ``1``. And finally, the
+ (empty) product of its eigenvalues is likewise just unity::
+
+ sage: J = TrivialEJA()
+ sage: J.zero().det()
+ 1
+
TESTS:
An element is invertible if and only if its determinant is
sage: (x*y).det() == x.det()*y.det()
True
+ The determinant in real matrix algebras is 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
+
"""
P = self.parent()
r = P.rank()
- p = P._charpoly_coeff(0)
- # The _charpoly_coeff function already adds the factor of
- # -1 to ensure that _charpoly_coeff(0) is really what
- # appears in front of t^{0} in the charpoly. However,
- # we want (-1)^r times THAT for the determinant.
+
+ if r == 0:
+ # Special case, since we don't get the a0=1
+ # coefficient when the rank of the algebra
+ # is zero.
+ return P.base_ring().one()
+
+ p = P._charpoly_coefficients()[0]
+ # The _charpoly_coeff function already adds the factor of -1
+ # to ensure that _charpoly_coefficients()[0] is really what
+ # appears in front of t^{0} in the charpoly. However, we want
+ # (-1)^r times THAT for the determinant.
return ((-1)**r)*p(*self.to_vector())
+ @cached_method
def inverse(self):
"""
Return the Jordan-multiplicative inverse of this element.
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::
sage: while not x.is_invertible():
....: x = J.random_element()
sage: x_vec = x.to_vector()
- sage: x0 = x_vec[0]
+ sage: x0 = x_vec[:1]
sage: x_bar = x_vec[1:]
- sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
- sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
- sage: x_inverse = coeff*inv_vec
+ sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
+ sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
+ sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
sage: x.inverse() == J.from_vector(x_inverse)
True
sage: JordanSpinEJA(3).zero().inverse()
Traceback (most recent call last):
...
- ValueError: element is not invertible
+ ZeroDivisionError: element is not invertible
TESTS:
....: 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")
-
- return (~self.quadratic_representation())(self)
-
-
+ 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()
+
+ 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 paren't algebra's zero element is a root
- of this element's minimal polynomial.
-
- 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: (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
- # 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)
+ if self.parent()._charpoly_coefficients.is_in_cache():
+ # The determinant will be quicker than inverting the
+ # quadratic representation, most likely.
+ return (not self.det().is_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_minimal_idempotent(self):
+ def is_primitive_idempotent(self):
"""
- Return whether or not this element is a minimal idempotent.
+ Return whether or not this element is a primitive (or minimal)
+ idempotent.
+ A primitive idempotent is a non-zero idempotent that is not
+ the sum of two other non-zero idempotents. Remark 2.7.15 in
+ Baes shows that this is what he refers to as a "minimal
+ idempotent."
An element of a Euclidean Jordan algebra is a minimal idempotent
if it :meth:`is_idempotent` and if its Peirce subalgebra
sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
....: RealSymmetricEJA,
+ ....: TrivialEJA,
....: random_eja)
WARNING::
The spectral decomposition of a non-regular element should always
contain at least one non-minimal idempotent::
- sage: J = RealSymmetricEJA(3, AA)
+ sage: J = RealSymmetricEJA(3)
sage: x = sum(J.gens())
sage: x.is_regular()
False
- sage: [ c.is_minimal_idempotent()
+ sage: [ c.is_primitive_idempotent()
....: for (l,c) in x.spectral_decomposition() ]
[False, True]
On the other hand, the spectral decomposition of a regular
element should always be in terms of minimal idempotents::
- sage: J = JordanSpinEJA(4, AA)
- sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
+ sage: J = JordanSpinEJA(4)
+ sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
sage: x.is_regular()
True
- sage: [ c.is_minimal_idempotent()
+ sage: [ c.is_primitive_idempotent()
....: for (l,c) in x.spectral_decomposition() ]
[True, True]
sage: set_random_seed()
sage: J = random_eja()
- sage: J.rank() == 1 or not J.one().is_minimal_idempotent()
+ sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
True
A non-idempotent cannot be a minimal idempotent::
sage: set_random_seed()
sage: J = JordanSpinEJA(4)
sage: x = J.random_element()
- sage: (not x.is_idempotent()) and x.is_minimal_idempotent()
+ sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
False
Proposition 2.7.19 in Baes says that an element is a minimal
sage: J = JordanSpinEJA(4)
sage: x = J.random_element()
sage: expected = (x.is_idempotent() and x.trace() == 1)
- sage: actual = x.is_minimal_idempotent()
+ sage: actual = x.is_primitive_idempotent()
sage: actual == expected
True
+ Primitive idempotents must be non-zero::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: J.zero().is_idempotent()
+ True
+ sage: J.zero().is_primitive_idempotent()
+ False
+
+ As a consequence of the fact that primitive idempotents must
+ be non-zero, there are no primitive idempotents in a trivial
+ Euclidean Jordan algebra::
+
+ sage: J = TrivialEJA()
+ sage: J.one().is_idempotent()
+ True
+ sage: J.one().is_primitive_idempotent()
+ False
+
"""
- # TODO: when the Peirce decomposition is implemented for real,
- # we can use that instead of finding this eigenspace manually.
- #
- # Trivial eigenspaces don't appear in the list, so we default to the
- # trivial one and override it if there's a nontrivial space in the
- # list.
if not self.is_idempotent():
return False
- J1 = VectorSpace(self.parent().base_ring(), 0)
- for (eigval, eigspace) in self.operator().matrix().left_eigenspaces():
- if eigval == 1:
- J1 = eigspace
+ if self.is_zero():
+ return False
+
+ (_,_,J1) = self.parent().peirce_decomposition(self)
return (J1.dimension() == 1)
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
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::
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
sage: set_random_seed()
sage: J = JordanSpinEJA.random_instance()
+ sage: n = J.dimension()
sage: x = J.random_element()
- sage: x == x.coefficient(0)*J.one() or x.degree() == 2
+ sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
True
TESTS:
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
- return self.subalgebra_generated_by().dimension()
+ 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
+
+ # 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):
TESTS:
The minimal polynomial of the identity and zero elements are
- always the same::
+ always the same, except in trivial algebras where the minimal
+ polynomial of the unit/zero element is ``1``::
sage: set_random_seed()
- sage: J = random_eja(nontrivial=True)
- sage: J.one().minimal_polynomial()
+ sage: J = random_eja()
+ sage: mu = J.one().minimal_polynomial()
+ sage: t = mu.parent().gen()
+ sage: mu + int(J.is_trivial())*(t-2)
t - 1
- sage: J.zero().minimal_polynomial()
+ sage: mu = J.zero().minimal_polynomial()
+ sage: t = mu.parent().gen()
+ sage: mu + int(J.is_trivial())*(t-1)
t
The degree of an element is (by one definition) the degree
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,QQ)
- sage: J2 = RealSymmetricEJA(n,QQ,normalize_basis=False)
- sage: X = random_matrix(QQ,n)
+ sage: J1 = RealSymmetricEJA(n)
+ sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
+ sage: X = random_matrix(AA,n)
sage: X = X*X.transpose()
sage: x1 = J1(X)
sage: x2 = J2(X)
"""
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()
-
- A = self.subalgebra_generated_by()
- return A(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)
+
+ 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()
- 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: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
- ....: QuaternionHermitianEJA)
+ ....: HadamardEJA,
+ ....: QuaternionHermitianEJA,
+ ....: RealSymmetricEJA)
EXAMPLES::
sage: J = ComplexHermitianEJA(3)
sage: J.one()
- e0 + e3 + e8
- sage: J.one().natural_representation()
- [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]
+ b0 + b3 + b8
+ sage: J.one().to_matrix()
+ +---+---+---+
+ | 1 | 0 | 0 |
+ +---+---+---+
+ | 0 | 1 | 0 |
+ +---+---+---+
+ | 0 | 0 | 1 |
+ +---+---+---+
::
- 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]
+ b0 + b5
+ sage: J.one().to_matrix()
+ +---+---+
+ | 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().natural_basis()
- W = self.parent().natural_basis_space()
- return W.linear_combination(izip(B,self.to_vector()))
+ B = self.parent().matrix_basis()
+ W = self.parent().matrix_space()
+
+ 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):
SETUP::
sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
- ....: RealCartesianProductEJA)
+ ....: HadamardEJA)
EXAMPLES::
- sage: J = RealCartesianProductEJA(2)
+ sage: J = HadamardEJA(2)
sage: x = sum(J.gens())
sage: x.norm()
- sqrt(2)
+ 1.414213562373095?
::
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: set_random_seed()
sage: x = JordanSpinEJA.random_instance().random_element()
sage: x_vec = x.to_vector()
+ sage: Q = matrix.identity(x.base_ring(), 0)
sage: n = x_vec.degree()
- sage: x0 = x_vec[0]
- sage: x_bar = x_vec[1:]
- sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
- sage: B = 2*x0*x_bar.row()
- sage: C = 2*x0*x_bar.column()
- sage: D = matrix.identity(QQ, n-1)
- sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
- sage: D = D + 2*x_bar.tensor_product(x_bar)
- sage: Q = matrix.block(2,2,[A,B,C,D])
+ sage: if n > 0:
+ ....: x0 = x_vec[0]
+ ....: x_bar = x_vec[1:]
+ ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
+ ....: B = 2*x0*x_bar.row()
+ ....: C = 2*x0*x_bar.column()
+ ....: D = matrix.identity(x.base_ring(), n-1)
+ ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
+ ....: D = D + 2*x_bar.tensor_product(x_bar)
+ ....: Q = matrix.block(2,2,[A,B,C,D])
sage: Q == x.quadratic_representation().matrix()
True
The spectral decomposition of the identity is ``1`` times itself,
and the spectral decomposition of zero is ``0`` times the identity::
- sage: J = RealSymmetricEJA(3,AA)
+ 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::
- sage: J = RealSymmetricEJA(4,AA)
+ sage: J = RealSymmetricEJA(4)
sage: x = sum(J.gens())
sage: sd = x.spectral_decomposition()
sage: l0 = sd[0][0]
sage: l0*c0 + l1*c1 == x
True
+ The spectral decomposition should work in subalgebras, too::
+
+ sage: J = RealSymmetricEJA(4)
+ 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, c2), (1, c0)]
+
"""
- P = self.parent()
- 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.
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:
True
"""
- 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):
sage: from mjo.eja.eja_algebra import random_eja
- TESTS::
+ TESTS:
+
+ Ensure that we can find an idempotent in a non-trivial algebra
+ where there are non-nilpotent elements, or that we get the dumb
+ solution in the trivial algebra::
sage: set_random_seed()
sage: J = random_eja()
sage: x = J.random_element()
- sage: while x.is_nilpotent():
+ sage: while x.is_nilpotent() and not J.is_trivial():
....: x = J.random_element()
sage: c = x.subalgebra_idempotent()
sage: c^2 == c
True
"""
+ if self.parent().is_trivial():
+ return self
+
if self.is_nilpotent():
raise ValueError("this only works with non-nilpotent elements!")
# will be minimal for some natural number s...
s = 0
minimal_dim = J.dimension()
- for i in xrange(1, minimal_dim):
+ for i in range(1, minimal_dim):
this_dim = (u**i).operator().matrix().image().dimension()
if this_dim < minimal_dim:
minimal_dim = this_dim
SETUP::
sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
- ....: RealCartesianProductEJA,
+ ....: HadamardEJA,
....: TrivialEJA,
....: random_eja)
::
- sage: J = RealCartesianProductEJA(5)
+ sage: J = HadamardEJA(5)
sage: J.one().trace()
5
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()
# the trace is an empty sum.
return P.base_ring().zero()
- p = P._charpoly_coeff(r-1)
+ p = P._charpoly_coefficients()[r-1]
# The _charpoly_coeff function already adds the factor of
# -1 to ensure that _charpoly_coeff(r-1) is really what
# appears in front of t^{r-1} in the charpoly. However,
SETUP::
sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
- ....: RealCartesianProductEJA)
+ ....: HadamardEJA)
EXAMPLES::
- sage: J = RealCartesianProductEJA(2)
+ sage: J = HadamardEJA(2)
sage: x = sum(J.gens())
sage: x.trace_norm()
- sqrt(2)
+ 1.414213562373095?
::
sage: J = JordanSpinEJA(4)
sage: x = sum(J.gens())
sage: x.trace_norm()
- 2*sqrt(2)
+ 2.828427124746190?
"""
return self.trace_inner_product(self).sqrt()
projects / sage.d.git / blobdiff