]> 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 562797318646cb1dd837691f846fd114b44ec131..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
 
@@ -96,9 +113,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             category = MagmaticAlgebras(field).FiniteDimensional()
             category = category.WithBasis().Unital()
 
-        # The multiplication table had better be square
-        n = len(mult_table)
-
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
                      range(n),
@@ -117,7 +131,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             for ls in mult_table
         ]
 
-        if check:
+        if check_axioms:
             if not self._is_commutative():
                 raise ValueError("algebra is not commutative")
             if not self._is_jordanian():
@@ -250,8 +264,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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=False`` and passed
-        an invalid multiplication table.
+        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())
@@ -266,7 +280,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         ``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=False`` and passed an invalid multiplication table.
+        ``check_axioms=False`` and passed an invalid multiplication table.
         """
         return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
                     ==
@@ -280,17 +294,29 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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=False`` and passed
-        an invalid multiplication table.
+        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)
-                    if (x*y).inner_product(z) != x.inner_product(y*z):
-                        return False
+                    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
 
@@ -720,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:
@@ -979,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):
     """
@@ -1078,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():
@@ -1114,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
@@ -1135,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
@@ -1225,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]
@@ -1381,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)
 
 
@@ -1677,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)
 
 
@@ -1978,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 BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
+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(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 =
@@ -2061,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):
@@ -2154,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):
@@ -2188,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)
 
 
@@ -2237,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())