]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: use zip() instead of izip(), which doesn't exist in python-3.x.
[sage.d.git] / mjo / eja / eja_algebra.py
index 515e25f9c97de0ae102e8dc6bd7b29db89e39dd0..0f2b655f21795cb1feb0e8b41cfc878261a67962 100644 (file)
@@ -5,7 +5,7 @@ are used in optimization, and have some additional nice methods beyond
 what can be supported in a general Jordan Algebra.
 """
 
-from itertools import izip, repeat
+from itertools import repeat
 
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
@@ -13,16 +13,16 @@ from sage.combinat.free_module import CombinatorialFreeModule
 from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
+from sage.misc.lazy_import import lazy_import
 from sage.misc.prandom import choice
 from sage.misc.table import table
 from sage.modules.free_module import FreeModule, VectorSpace
-from sage.rings.integer_ring import ZZ
-from sage.rings.number_field.number_field import QuadraticField
-from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-from sage.rings.rational_field import QQ
-from sage.rings.real_lazy import CLF, RLF
-
+from sage.rings.all import (ZZ, QQ, RR, RLF, CLF,
+                            PolynomialRing,
+                            QuadraticField)
 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
+lazy_import('mjo.eja.eja_subalgebra',
+            'FiniteDimensionalEuclideanJordanSubalgebra')
 from mjo.eja.eja_utils import _mat2vec
 
 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
@@ -40,11 +40,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                  rank,
                  prefix='e',
                  category=None,
-                 natural_basis=None):
+                 natural_basis=None,
+                 check=True):
         """
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import random_eja
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
 
         EXAMPLES:
 
@@ -56,7 +57,23 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: x*y == y*x
             True
 
+        TESTS:
+
+        The ``field`` we're given must be real::
+
+            sage: JordanSpinEJA(2,QQbar)
+            Traceback (most recent call last):
+            ...
+            ValueError: field is not real
+
         """
+        if check:
+            if not field.is_subring(RR):
+                # Note: this does return true for the real algebraic
+                # field, and any quadratic field where we've specified
+                # a real embedding.
+                raise ValueError('field is not real')
+
         self._rank = rank
         self._natural_basis = natural_basis
 
@@ -237,7 +254,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         """
         (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
         R = A_of_x.base_ring()
-        if i >= self.rank():
+
+        if i == self.rank():
+            return R.one()
+        if i > self.rank():
             # Guaranteed by theory
             return R.zero()
 
@@ -346,7 +366,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
 
         EXAMPLES:
 
@@ -360,12 +380,22 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: p(*xvec)
             t^2 - 2*t + 1
 
+        By definition, the characteristic polynomial is a monic
+        degree-zero polynomial in a rank-zero algebra. Note that
+        Cayley-Hamilton is indeed satisfied since the polynomial
+        ``1`` evaluates to the identity element of the algebra on
+        any argument::
+
+            sage: J = TrivialEJA()
+            sage: J.characteristic_polynomial()
+            1
+
         """
         r = self.rank()
         n = self.dimension()
 
-        # The list of coefficient polynomials a_1, a_2, ..., a_n.
-        a = [ self._charpoly_coeff(i) for i in range(n) ]
+        # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
+        a = [ self._charpoly_coeff(i) for i in range(r+1) ]
 
         # We go to a bit of trouble here to reorder the
         # indeterminates, so that it's easier to evaluate the
@@ -377,17 +407,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         S = PolynomialRing(S, R.variable_names())
         t = S(t)
 
-        # Note: all entries past the rth should be zero. The
-        # coefficient of the highest power (x^r) is 1, but it doesn't
-        # appear in the solution vector which contains coefficients
-        # for the other powers (to make them sum to x^r).
-        if (r < n):
-            a[r] = 1 # corresponds to x^r
-        else:
-            # When the rank is equal to the dimension, trying to
-            # assign a[r] goes out-of-bounds.
-            a.append(1) # corresponds to x^r
-
         return sum( a[k]*(t**k) for k in xrange(len(a)) )
 
 
@@ -428,15 +447,19 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  TrivialEJA)
 
         EXAMPLES::
 
             sage: J = ComplexHermitianEJA(3)
             sage: J.is_trivial()
             False
-            sage: A = J.zero().subalgebra_generated_by()
-            sage: A.is_trivial()
+
+        ::
+
+            sage: J = TrivialEJA()
+            sage: J.is_trivial()
             True
 
         """
@@ -611,13 +634,118 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return self.linear_combination(zip(self.gens(), coeffs))
 
 
-    def random_element(self):
-        # Temporary workaround for https://trac.sagemath.org/ticket/28327
-        if self.is_trivial():
-            return self.zero()
-        else:
-            s = super(FiniteDimensionalEuclideanJordanAlgebra, self)
-            return s.random_element()
+    def peirce_decomposition(self, c):
+        """
+        The Peirce decomposition of this algebra relative to the
+        idempotent ``c``.
+
+        In the future, this can be extended to a complete system of
+        orthogonal idempotents.
+
+        INPUT:
+
+          - ``c`` -- an idempotent of this algebra.
+
+        OUTPUT:
+
+        A triple (J0, J5, J1) containing two subalgebras and one subspace
+        of this algebra,
+
+          - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
+            corresponding to the eigenvalue zero.
+
+          - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
+            corresponding to the eigenvalue one-half.
+
+          - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
+            corresponding to the eigenvalue one.
+
+        These are the only possible eigenspaces for that operator, and this
+        algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
+        orthogonal, and are subalgebras of this algebra with the appropriate
+        restrictions.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
+
+        EXAMPLES:
+
+        The canonical example comes from the symmetric matrices, which
+        decompose into diagonal and off-diagonal parts::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: C = matrix(QQ, [ [1,0,0],
+            ....:                  [0,1,0],
+            ....:                  [0,0,0] ])
+            sage: c = J(C)
+            sage: J0,J5,J1 = J.peirce_decomposition(c)
+            sage: J0
+            Euclidean Jordan algebra of dimension 1...
+            sage: J5
+            Vector space of degree 6 and dimension 2...
+            sage: J1
+            Euclidean Jordan algebra of dimension 3...
+
+        TESTS:
+
+        Every algebra decomposes trivially with respect to its identity
+        element::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J0,J5,J1 = J.peirce_decomposition(J.one())
+            sage: J0.dimension() == 0 and J5.dimension() == 0
+            True
+            sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
+            True
+
+        The identity elements in the two subalgebras are the
+        projections onto their respective subspaces of the
+        superalgebra's identity element::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: if not J.is_trivial():
+            ....:     while x.is_nilpotent():
+            ....:         x = J.random_element()
+            sage: c = x.subalgebra_idempotent()
+            sage: J0,J5,J1 = J.peirce_decomposition(c)
+            sage: J1(c) == J1.one()
+            True
+            sage: J0(J.one() - c) == J0.one()
+            True
+
+        """
+        if not c.is_idempotent():
+            raise ValueError("element is not idempotent: %s" % c)
+
+        # Default these to what they should be if they turn out to be
+        # trivial, because eigenspaces_left() won't return eigenvalues
+        # corresponding to trivial spaces (e.g. it returns only the
+        # eigenspace corresponding to lambda=1 if you take the
+        # decomposition relative to the identity element).
+        trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
+        J0 = trivial                          # eigenvalue zero
+        J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
+        J1 = trivial                          # eigenvalue one
+
+        for (eigval, eigspace) in c.operator().matrix().left_eigenspaces():
+            if eigval == ~(self.base_ring()(2)):
+                J5 = eigspace
+            else:
+                gens = tuple( self.from_vector(b) for b in eigspace.basis() )
+                subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, gens)
+                if eigval == 0:
+                    J0 = subalg
+                elif eigval == 1:
+                    J1 = subalg
+                else:
+                    raise ValueError("unexpected eigenvalue: %s" % eigval)
+
+        return (J0, J5, J1)
+
 
     def random_elements(self, count):
         """
@@ -683,11 +811,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         TESTS:
 
         Ensure that every EJA that we know how to construct has a
-        positive integer rank::
+        positive integer rank, unless the algebra is trivial in
+        which case its rank will be zero::
 
             sage: set_random_seed()
-            sage: r = random_eja().rank()
-            sage: r in ZZ and r > 0
+            sage: J = random_eja()
+            sage: r = J.rank()
+            sage: r in ZZ
+            True
+            sage: r > 0 or (r == 0 and J.is_trivial())
             True
 
         """
@@ -760,6 +892,11 @@ class KnownRankEJA(object):
         Beware, this will crash for "most instances" because the
         constructor below looks wrong.
         """
+        if cls is TrivialEJA:
+            # The TrivialEJA class doesn't take an "n" argument because
+            # there's only one.
+            return cls(field)
+
         n = ZZ.random_element(cls._max_test_case_size()) + 1
         return cls(n, field, **kwargs)
 
@@ -838,32 +975,10 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
         return x.to_vector().inner_product(y.to_vector())
 
 
-def random_eja():
+def random_eja(field=QQ, nontrivial=False):
     """
     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 
-    ALGORITHM:
-
-    For now, we choose a random natural number ``n`` (greater than zero)
-    and then give you back one of the following:
-
-      * The cartesian product of the rational numbers ``n`` times; this is
-        ``QQ^n`` with the Hadamard product.
-
-      * The Jordan spin algebra on ``QQ^n``.
-
-      * The ``n``-by-``n`` rational symmetric matrices with the symmetric
-        product.
-
-      * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
-        in the space of ``2n``-by-``2n`` real symmetric matrices.
-
-      * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
-        in the space of ``4n``-by-``4n`` real symmetric matrices.
-
-    Later this might be extended to return Cartesian products of the
-    EJAs above.
-
     SETUP::
 
         sage: from mjo.eja.eja_algebra import random_eja
@@ -874,8 +989,11 @@ def random_eja():
         Euclidean Jordan algebra of dimension...
 
     """
-    classname = choice(KnownRankEJA.__subclasses__())
-    return classname.random_instance()
+    eja_classes = KnownRankEJA.__subclasses__()
+    if nontrivial:
+        eja_classes.remove(TrivialEJA)
+    classname = choice(eja_classes)
+    return classname.random_instance(field=field)
 
 
 
@@ -914,7 +1032,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
                 basis = tuple( s.change_ring(field) for s in basis )
             self._basis_normalizers = tuple(
                 ~(self.natural_inner_product(s,s).sqrt()) for s in basis )
-            basis = tuple(s*c for (s,c) in izip(basis,self._basis_normalizers))
+            basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
 
         Qs = self.multiplication_table_from_matrix_basis(basis)
 
@@ -937,8 +1055,8 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
             # with had entries in a nice field.
             return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i)
         else:
-            basis = ( (b/n) for (b,n) in izip(self.natural_basis(),
-                                              self._basis_normalizers) )
+            basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
+                                             self._basis_normalizers) )
 
             # Do this over the rationals and convert back at the end.
             J = MatrixEuclideanJordanAlgebra(QQ,
@@ -948,7 +1066,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
             (_,x,_,_) = J._charpoly_matrix_system()
             p = J._charpoly_coeff(i)
             # p might be missing some vars, have to substitute "optionally"
-            pairs = izip(x.base_ring().gens(), self._basis_normalizers)
+            pairs = zip(x.base_ring().gens(), self._basis_normalizers)
             substitutions = { v: v*c for (v,c) in pairs }
             result = p.subs(substitutions)
 
@@ -1223,15 +1341,17 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         n = M.nrows()
         if M.ncols() != n:
             raise ValueError("the matrix 'M' must be square")
-        field = M.base_ring()
+
+        # We don't need any adjoined elements...
+        field = M.base_ring().base_ring()
+
         blocks = []
         for z in M.list():
-            a = z.vector()[0] # real part, I guess
-            b = z.vector()[1] # imag part, I guess
+            a = z.list()[0] # real part, I guess
+            b = z.list()[1] # imag part, I guess
             blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
 
-        # We can drop the imaginaries here.
-        return matrix.block(field.base_ring(), n, blocks)
+        return matrix.block(field, n, blocks)
 
 
     @staticmethod
@@ -1338,6 +1458,16 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
 
         sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
 
+    EXAMPLES:
+
+    In theory, our "field" can be any subfield of the reals::
+
+        sage: ComplexHermitianEJA(2, AA)
+        Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
+        sage: ComplexHermitianEJA(2, RR)
+        Euclidean Jordan algebra of dimension 4 over Real Field with
+        53 bits of precision
+
     TESTS:
 
     The dimension of this algebra is `n^2`::
@@ -1623,6 +1753,16 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
 
         sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
 
+    EXAMPLES:
+
+    In theory, our "field" can be any subfield of the reals::
+
+        sage: QuaternionHermitianEJA(2, AA)
+        Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
+        sage: QuaternionHermitianEJA(2, RR)
+        Euclidean Jordan algebra of dimension 6 over Real Field with
+        53 bits of precision
+
     TESTS:
 
     The dimension of this algebra is `2*n^2 - n`::
@@ -1818,3 +1958,40 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
 
         """
         return x.to_vector().inner_product(y.to_vector())
+
+
+class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+    """
+    The trivial Euclidean Jordan algebra consisting of only a zero element.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import TrivialEJA
+
+    EXAMPLES::
+
+        sage: J = TrivialEJA()
+        sage: J.dimension()
+        0
+        sage: J.zero()
+        0
+        sage: J.one()
+        0
+        sage: 7*J.one()*12*J.one()
+        0
+        sage: J.one().inner_product(J.one())
+        0
+        sage: J.one().norm()
+        0
+        sage: J.one().subalgebra_generated_by()
+        Euclidean Jordan algebra of dimension 0 over Rational Field
+        sage: J.rank()
+        0
+
+    """
+    def __init__(self, field=QQ, **kwargs):
+        mult_table = []
+        fdeja = super(TrivialEJA, self)
+        # The rank is zero using my definition, namely the dimension of the
+        # largest subalgebra generated by any element.
+        return fdeja.__init__(field, mult_table, rank=0, **kwargs)