]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: begin generalizing the charpoly-over-QQ optimizations.
[sage.d.git] / mjo / eja / eja_algebra.py
index f327bf51aada40b33fd05fb0d1c38c124d4b545e..1fc3618005eb0ac8be12e0e3e09849ab6f983819 100644 (file)
@@ -57,11 +57,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                  prefix='e',
                  category=None,
                  natural_basis=None,
-                 check=True):
+                 check_field=True,
+                 check_axioms=True):
         """
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
+            sage: from mjo.eja.eja_algebra import (
+            ....:   FiniteDimensionalEuclideanJordanAlgebra,
+            ....:   JordanSpinEJA,
+            ....:   random_eja)
 
         EXAMPLES:
 
@@ -75,20 +79,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         TESTS:
 
-        The ``field`` we're given must be real::
+        The ``field`` we're given must be real with ``check_field=True``::
 
             sage: JordanSpinEJA(2,QQbar)
             Traceback (most recent call last):
             ...
-            ValueError: field is not real
+            ValueError: scalar field is not real
+
+        The multiplication table must be square with ``check_axioms=True``::
+
+            sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
+            Traceback (most recent call last):
+            ...
+            ValueError: multiplication table is not square
 
         """
-        if check:
+        if check_field:
             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')
+                # field, the rationals, and any quadratic field where
+                # we've specified a real embedding.
+                raise ValueError("scalar field is not real")
+
+        # The multiplication table had better be square
+        n = len(mult_table)
+        if check_axioms:
+            if not all( len(l) == n for l in mult_table ):
+                raise ValueError("multiplication table is not square")
 
         self._natural_basis = natural_basis
 
@@ -98,7 +115,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
-                     range(len(mult_table)),
+                     range(n),
                      prefix=prefix,
                      category=category)
         self.print_options(bracket='')
@@ -114,6 +131,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             for ls in mult_table
         ]
 
+        if check_axioms:
+            if not self._is_commutative():
+                raise ValueError("algebra is not commutative")
+            if not self._is_jordanian():
+                raise ValueError("Jordan identity does not hold")
+            if not self._inner_product_is_associative():
+                raise ValueError("inner product is not associative")
 
     def _element_constructor_(self, elt):
         """
@@ -235,6 +259,67 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     def product_on_basis(self, i, j):
         return self._multiplication_table[i][j]
 
+    def _is_commutative(self):
+        r"""
+        Whether or not this algebra's multiplication table is commutative.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid multiplication table.
+        """
+        return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
+                    for i in range(self.dimension())
+                    for j in range(self.dimension()) )
+
+    def _is_jordanian(self):
+        r"""
+        Whether or not this algebra's multiplication table respects the
+        Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
+
+        We only check one arrangement of `x` and `y`, so for a
+        ``True`` result to be truly true, you should also check
+        :meth:`_is_commutative`. This method should of course always
+        return ``True``, unless this algebra was constructed with
+        ``check_axioms=False`` and passed an invalid multiplication table.
+        """
+        return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
+                    ==
+                    (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
+                    for i in range(self.dimension())
+                    for j in range(self.dimension()) )
+
+    def _inner_product_is_associative(self):
+        r"""
+        Return whether or not this algebra's inner product `B` is
+        associative; that is, whether or not `B(xy,z) = B(x,yz)`.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid multiplication table.
+        """
+
+        # Used to check whether or not something is zero in an inexact
+        # ring. This number is sufficient to allow the construction of
+        # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
+        epsilon = 1e-16
+
+        for i in range(self.dimension()):
+            for j in range(self.dimension()):
+                for k in range(self.dimension()):
+                    x = self.monomial(i)
+                    y = self.monomial(j)
+                    z = self.monomial(k)
+                    diff = (x*y).inner_product(z) - x.inner_product(y*z)
+
+                    if self.base_ring().is_exact():
+                        if diff != 0:
+                            return False
+                    else:
+                        if diff.abs() > epsilon:
+                            return False
+
+        return True
+
     @cached_method
     def characteristic_polynomial_of(self):
         """
@@ -661,7 +746,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                 J5 = eigspace
             else:
                 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
-                subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, gens)
+                subalg = FiniteDimensionalEuclideanJordanSubalgebra(self,
+                                                                    gens,
+                                                                    check_axioms=False)
                 if eigval == 0:
                     J0 = subalg
                 elif eigval == 1:
@@ -920,79 +1007,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
 
-class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra):
-    """
-    Return the Euclidean Jordan Algebra corresponding to the set
-    `R^n` under the Hadamard product.
-
-    Note: this is nothing more than the Cartesian product of ``n``
-    copies of the spin algebra. Once Cartesian product algebras
-    are implemented, this can go.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import HadamardEJA
-
-    EXAMPLES:
-
-    This multiplication table can be verified by hand::
-
-        sage: J = HadamardEJA(3)
-        sage: e0,e1,e2 = J.gens()
-        sage: e0*e0
-        e0
-        sage: e0*e1
-        0
-        sage: e0*e2
-        0
-        sage: e1*e1
-        e1
-        sage: e1*e2
-        0
-        sage: e2*e2
-        e2
-
-    TESTS:
-
-    We can change the generator prefix::
-
-        sage: HadamardEJA(3, prefix='r').gens()
-        (r0, r1, r2)
-
-    """
-    def __init__(self, n, field=AA, **kwargs):
-        V = VectorSpace(field, n)
-        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
-                       for i in range(n) ]
-
-        fdeja = super(HadamardEJA, self)
-        fdeja.__init__(field, mult_table, **kwargs)
-        self.rank.set_cache(n)
-
-    def inner_product(self, x, y):
-        """
-        Faster to reimplement than to use natural representations.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import HadamardEJA
-
-        TESTS:
-
-        Ensure that this is the usual inner product for the algebras
-        over `R^n`::
-
-            sage: set_random_seed()
-            sage: J = HadamardEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: X = x.natural_representation()
-            sage: Y = y.natural_representation()
-            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
-            True
-
-        """
-        return x.to_vector().inner_product(y.to_vector())
-
 
 def random_eja(field=AA):
     """
@@ -1019,6 +1033,44 @@ def random_eja(field=AA):
 
 
 
+class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+    r"""
+    Algebras whose basis consists of vectors with rational
+    entries. Equivalently, algebras whose multiplication tables
+    contain only rational coefficients.
+
+    When an EJA has a basis that can be made rational, we can speed up
+    the computation of its characteristic polynomial by doing it over
+    ``QQ``. All of the named EJA constructors that we provide fall
+    into this category.
+    """
+    @cached_method
+    def _charpoly_coefficients(self):
+        r"""
+        Override the parent method with something that tries to compute
+        over a faster (non-extension) field.
+        """
+        if self.base_ring() is QQ:
+            # There's no need to construct *another* algebra over the
+            # rationals if this one is already over the rationals.
+            superclass = super(RationalBasisEuclideanJordanAlgebra, self)
+            return superclass._charpoly_coefficients()
+
+        mult_table = tuple(
+            map(lambda x: x.to_vector(), ls)
+            for ls in self._multiplication_table
+        )
+
+        # Do the computation over the rationals. The answer will be
+        # the same, because our basis coordinates are (essentially)
+        # rational.
+        J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
+                                                    mult_table,
+                                                    check_field=False,
+                                                    check_axioms=False)
+        return J._charpoly_coefficients()
+
+
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     @staticmethod
     def _max_test_case_size():
@@ -1055,9 +1107,10 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
 
         Qs = self.multiplication_table_from_matrix_basis(basis)
 
-        fdeja = super(MatrixEuclideanJordanAlgebra, self)
-        fdeja.__init__(field, Qs, natural_basis=basis, **kwargs)
-        return
+        super(MatrixEuclideanJordanAlgebra, self).__init__(field,
+                                                           Qs,
+                                                           natural_basis=basis,
+                                                           **kwargs)
 
 
     @cached_method
@@ -1076,10 +1129,14 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
 
             # Do this over the rationals and convert back at the end.
             # Only works because we know the entries of the basis are
-            # integers.
+            # integers. The argument ``check_axioms=False`` is required
+            # because the trace inner-product method for this
+            # class is a stub and can't actually be checked.
             J = MatrixEuclideanJordanAlgebra(QQ,
                                              basis,
-                                             normalize_basis=False)
+                                             normalize_basis=False,
+                                             check_field=False,
+                                             check_axioms=False)
             a = J._charpoly_coefficients()
 
             # Unfortunately, changing the basis does change the
@@ -1166,16 +1223,11 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()
 
-        if tr in RLF:
-            # It's real already.
-            return tr
-
-        # Otherwise, try the thing that works for complex numbers; and
-        # if that doesn't work, the thing that works for quaternions.
         try:
-            return tr.vector()[0] # real part, imag part is index 1
+            # Works in QQ, AA, RDF, et cetera.
+            return tr.real()
         except AttributeError:
-            # A quaternions doesn't have a vector() method, but does
+            # A quaternion doesn't have a real() method, but does
             # have coefficient_tuple() method that returns the
             # coefficients of 1, i, j, and k -- in that order.
             return tr.coefficient_tuple()[0]
@@ -1322,7 +1374,10 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n, field)
-        super(RealSymmetricEJA, self).__init__(field, basis, **kwargs)
+        super(RealSymmetricEJA, self).__init__(field,
+                                               basis,
+                                               check_axioms=False,
+                                               **kwargs)
         self.rank.set_cache(n)
 
 
@@ -1618,7 +1673,10 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(ComplexHermitianEJA,self).__init__(field, basis, **kwargs)
+        super(ComplexHermitianEJA,self).__init__(field,
+                                                 basis,
+                                                 check_axioms=False,
+                                                 **kwargs)
         self.rank.set_cache(n)
 
 
@@ -1919,11 +1977,90 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n,field)
-        super(QuaternionHermitianEJA,self).__init__(field, basis, **kwargs)
+        super(QuaternionHermitianEJA,self).__init__(field,
+                                                    basis,
+                                                    check_axioms=False,
+                                                    **kwargs)
+        self.rank.set_cache(n)
+
+
+class HadamardEJA(RationalBasisEuclideanJordanAlgebra):
+    """
+    Return the Euclidean Jordan Algebra corresponding to the set
+    `R^n` under the Hadamard product.
+
+    Note: this is nothing more than the Cartesian product of ``n``
+    copies of the spin algebra. Once Cartesian product algebras
+    are implemented, this can go.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import HadamardEJA
+
+    EXAMPLES:
+
+    This multiplication table can be verified by hand::
+
+        sage: J = HadamardEJA(3)
+        sage: e0,e1,e2 = J.gens()
+        sage: e0*e0
+        e0
+        sage: e0*e1
+        0
+        sage: e0*e2
+        0
+        sage: e1*e1
+        e1
+        sage: e1*e2
+        0
+        sage: e2*e2
+        e2
+
+    TESTS:
+
+    We can change the generator prefix::
+
+        sage: HadamardEJA(3, prefix='r').gens()
+        (r0, r1, r2)
+
+    """
+    def __init__(self, n, field=AA, **kwargs):
+        V = VectorSpace(field, n)
+        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
+                       for i in range(n) ]
+
+        super(HadamardEJA, self).__init__(field,
+                                          mult_table,
+                                          check_axioms=False,
+                                          **kwargs)
         self.rank.set_cache(n)
 
+    def inner_product(self, x, y):
+        """
+        Faster to reimplement than to use natural representations.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import HadamardEJA
+
+        TESTS:
+
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::
+
+            sage: set_random_seed()
+            sage: J = HadamardEJA.random_instance()
+            sage: x,y = J.random_elements(2)
+            sage: X = x.natural_representation()
+            sage: Y = y.natural_representation()
+            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
+            True
+
+        """
+        return x.to_vector().inner_product(y.to_vector())
+
 
-class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
+class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
     r"""
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the half-trace inner product and jordan product ``x*y =
@@ -2002,8 +2139,10 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded
         # by the ambient dimension).
-        fdeja = super(BilinearFormEJA, self)
-        fdeja.__init__(field, mult_table, **kwargs)
+        super(BilinearFormEJA, self).__init__(field,
+                                              mult_table,
+                                              check_axioms=False,
+                                              **kwargs)
         self.rank.set_cache(min(n,2))
 
     def inner_product(self, x, y):
@@ -2095,7 +2234,7 @@ class JordanSpinEJA(BilinearFormEJA):
     def __init__(self, n, field=AA, **kwargs):
         # This is a special case of the BilinearFormEJA with the identity
         # matrix as its bilinear form.
-        return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
+        super(JordanSpinEJA, self).__init__(n, field, **kwargs)
 
 
 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
@@ -2129,10 +2268,12 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     def __init__(self, field=AA, **kwargs):
         mult_table = []
-        fdeja = super(TrivialEJA, self)
+        super(TrivialEJA, self).__init__(field,
+                                         mult_table,
+                                         check_axioms=False,
+                                         **kwargs)
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
-        fdeja.__init__(field, mult_table, **kwargs)
         self.rank.set_cache(0)
 
 
@@ -2178,6 +2319,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
                 p = (J2.monomial(i)*J2.monomial(j)).to_vector()
                 mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
 
-        fdeja = super(DirectSumEJA, self)
-        fdeja.__init__(field, mult_table, **kwargs)
+        super(DirectSumEJA, self).__init__(field,
+                                           mult_table,
+                                           check_axioms=False,
+                                           **kwargs)
         self.rank.set_cache(J1.rank() + J2.rank())