]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: move the TrivialEJA random_instance() into an override method.
[sage.d.git] / mjo / eja / eja_algebra.py
index 26fe1929be872b393fe41926eeda801dbd0a9436..7089c50ab81c222b5be8ddca3c6667da9699cf98 100644 (file)
@@ -217,7 +217,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return self.from_vector(coords)
 
     @staticmethod
-    def _max_test_case_size():
+    def _max_random_instance_size():
         """
         Return an integer "size" that is an upper bound on the size of
         this algebra when it is used in a random test
@@ -233,7 +233,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         interpreted to be far less than the dimension) should override
         with a smaller number.
         """
-        return 5
+        raise NotImplementedError
 
     def _repr_(self):
         """
@@ -839,12 +839,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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)
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
         return cls(n, field, **kwargs)
 
     @cached_method
@@ -1007,81 +1002,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) ]
-
-        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())
-
 
 def random_eja(field=AA):
     """
@@ -1108,9 +1028,68 @@ 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.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+        EXAMPLES:
+
+        The base ring of the resulting polynomial coefficients is what
+        it should be, and not the rationals (unless the algebra was
+        already over the rationals)::
+
+            sage: J = JordanSpinEJA(3)
+            sage: J._charpoly_coefficients()
+            (X1^2 - X2^2 - X3^2, -2*X1)
+            sage: a0 = J._charpoly_coefficients()[0]
+            sage: J.base_ring()
+            Algebraic Real Field
+            sage: a0.base_ring()
+            Algebraic Real 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)
+        a = J._charpoly_coefficients()
+        return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
+
+
 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     @staticmethod
-    def _max_test_case_size():
+    def _max_random_instance_size():
         # Play it safe, since this will be squared and the underlying
         # field can have dimension 4 (quaternions) too.
         return 2
@@ -1156,44 +1135,44 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         Override the parent method with something that tries to compute
         over a faster (non-extension) field.
         """
-        if self._basis_normalizers is None:
-            # We didn't normalize, so assume that the basis we started
-            # with had entries in a nice field.
+        if self._basis_normalizers is None or self.base_ring() is QQ:
+            # We didn't normalize, or the basis we started with had
+            # entries in a nice field already. Just compute the thing.
             return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
-        else:
-            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.
-            # Only works because we know the entries of the basis are
-            # 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,
-                                             check_field=False,
-                                             check_axioms=False)
-            a = J._charpoly_coefficients()
-
-            # Unfortunately, changing the basis does change the
-            # coefficients of the characteristic polynomial, but since
-            # these are really the coefficients of the "characteristic
-            # polynomial of" function, everything is still nice and
-            # unevaluated. It's therefore "obvious" how scaling the
-            # basis affects the coordinate variables X1, X2, et
-            # cetera. Scaling the first basis vector up by "n" adds a
-            # factor of 1/n into every "X1" term, for example. So here
-            # we simply undo the basis_normalizer scaling that we
-            # performed earlier.
-            #
-            # The a[0] access here is safe because trivial algebras
-            # won't have any basis normalizers and therefore won't
-            # make it to this "else" branch.
-            XS = a[0].parent().gens()
-            subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
-                          for i in range(len(XS)) }
-            return tuple( a_i.subs(subs_dict) for a_i in a )
+
+        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.
+        # Only works because we know the entries of the basis are
+        # 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,
+                                         check_field=False,
+                                         check_axioms=False)
+        a = J._charpoly_coefficients()
+
+        # Unfortunately, changing the basis does change the
+        # coefficients of the characteristic polynomial, but since
+        # these are really the coefficients of the "characteristic
+        # polynomial of" function, everything is still nice and
+        # unevaluated. It's therefore "obvious" how scaling the
+        # basis affects the coordinate variables X1, X2, et
+        # cetera. Scaling the first basis vector up by "n" adds a
+        # factor of 1/n into every "X1" term, for example. So here
+        # we simply undo the basis_normalizer scaling that we
+        # performed earlier.
+        #
+        # The a[0] access here is safe because trivial algebras
+        # won't have any basis normalizers and therefore won't
+        # make it to this "else" branch.
+        XS = a[0].parent().gens()
+        subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
+                      for i in range(len(XS)) }
+        return tuple( a_i.subs(subs_dict) for a_i in a )
 
 
     @staticmethod
@@ -1322,7 +1301,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
     The dimension of this algebra is `(n^2 + n) / 2`::
 
         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: J = RealSymmetricEJA(n)
         sage: J.dimension() == (n^2 + n)/2
@@ -1405,7 +1384,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
 
 
     @staticmethod
-    def _max_test_case_size():
+    def _max_random_instance_size():
         return 4 # Dimension 10
 
 
@@ -1451,7 +1430,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         Embedding is a homomorphism (isomorphism, in fact)::
 
             sage: set_random_seed()
-            sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
+            sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_random_instance_size()
             sage: n = ZZ.random_element(n_max)
             sage: F = QuadraticField(-1, 'I')
             sage: X = random_matrix(F, n)
@@ -1603,7 +1582,7 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
     The dimension of this algebra is `n^2`::
 
         sage: set_random_seed()
-        sage: n_max = ComplexHermitianEJA._max_test_case_size()
+        sage: n_max = ComplexHermitianEJA._max_random_instance_size()
         sage: n = ZZ.random_element(1, n_max)
         sage: J = ComplexHermitianEJA(n)
         sage: J.dimension() == n^2
@@ -1747,7 +1726,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         Embedding is a homomorphism (isomorphism, in fact)::
 
             sage: set_random_seed()
-            sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
+            sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_random_instance_size()
             sage: n = ZZ.random_element(n_max)
             sage: Q = QuaternionAlgebra(QQ,-1,-1)
             sage: X = random_matrix(Q, n)
@@ -1906,7 +1885,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
     The dimension of this algebra is `2*n^2 - n`::
 
         sage: set_random_seed()
-        sage: n_max = QuaternionHermitianEJA._max_test_case_size()
+        sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
         sage: n = ZZ.random_element(1, n_max)
         sage: J = QuaternionHermitianEJA(n)
         sage: J.dimension() == 2*(n^2) - n
@@ -2021,7 +2000,87 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
         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)
+
+    @staticmethod
+    def _max_random_instance_size():
+        return 5
+
+    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 =
@@ -2106,6 +2165,10 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
                                               **kwargs)
         self.rank.set_cache(min(n,2))
 
+    @staticmethod
+    def _max_random_instance_size():
+        return 5
+
     def inner_product(self, x, y):
         r"""
         Half of the trace inner product.
@@ -2237,6 +2300,11 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)
 
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        # We don't take a "size" argument so the superclass method is
+        # inappropriate for us.
+        return cls(field, **kwargs)
 
 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
@@ -2264,6 +2332,7 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     """
     def __init__(self, J1, J2, field=AA, **kwargs):
+        self._factors = (J1, J2)
         n1 = J1.dimension()
         n2 = J2.dimension()
         n = n1+n2
@@ -2285,3 +2354,136 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
                                            check_axioms=False,
                                            **kwargs)
         self.rank.set_cache(J1.rank() + J2.rank())
+
+
+    def factors(self):
+        r"""
+        Return the pair of this algebra's factors.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  JordanSpinEJA,
+            ....:                                  DirectSumEJA)
+
+        EXAMPLES::
+
+            sage: J1 = HadamardEJA(2,QQ)
+            sage: J2 = JordanSpinEJA(3,QQ)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: J.factors()
+            (Euclidean Jordan algebra of dimension 2 over Rational Field,
+             Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+        """
+        return self._factors
+
+    def projections(self):
+        r"""
+        Return a pair of projections onto this algebra's factors.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  DirectSumEJA)
+
+        EXAMPLES::
+
+            sage: J1 = JordanSpinEJA(2)
+            sage: J2 = ComplexHermitianEJA(2)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: (pi_left, pi_right) = J.projections()
+            sage: J.one().to_vector()
+            (1, 0, 1, 0, 0, 1)
+            sage: pi_left(J.one()).to_vector()
+            (1, 0)
+            sage: pi_right(J.one()).to_vector()
+            (1, 0, 0, 1)
+
+        """
+        (J1,J2) = self.factors()
+        n = J1.dimension()
+        pi_left  = lambda x: J1.from_vector(x.to_vector()[:n])
+        pi_right = lambda x: J2.from_vector(x.to_vector()[n:])
+        return (pi_left, pi_right)
+
+    def inclusions(self):
+        r"""
+        Return the pair of inclusion maps from our factors into us.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  DirectSumEJA)
+
+        EXAMPLES::
+
+            sage: J1 = JordanSpinEJA(3)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: (iota_left, iota_right) = J.inclusions()
+            sage: iota_left(J1.zero()) == J.zero()
+            True
+            sage: iota_right(J2.zero()) == J.zero()
+            True
+            sage: J1.one().to_vector()
+            (1, 0, 0)
+            sage: iota_left(J1.one()).to_vector()
+            (1, 0, 0, 0, 0, 0)
+            sage: J2.one().to_vector()
+            (1, 0, 1)
+            sage: iota_right(J2.one()).to_vector()
+            (0, 0, 0, 1, 0, 1)
+            sage: J.one().to_vector()
+            (1, 0, 0, 1, 0, 1)
+
+        """
+        (J1,J2) = self.factors()
+        n = J1.dimension()
+        V_basis = self.vector_space().basis()
+        I1 = matrix.column(self.base_ring(), V_basis[:n])
+        I2 = matrix.column(self.base_ring(), V_basis[n:])
+        iota_left = lambda x: self.from_vector(I1*x.to_vector())
+        iota_right = lambda x: self.from_vector(I2*+x.to_vector())
+        return (iota_left, iota_right)
+
+    def inner_product(self, x, y):
+        r"""
+        The standard Cartesian inner-product.
+
+        We project ``x`` and ``y`` onto our factors, and add up the
+        inner-products from the subalgebras.
+
+        SETUP::
+
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  DirectSumEJA)
+
+        EXAMPLE::
+
+            sage: J1 = HadamardEJA(3)
+            sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
+            sage: J = DirectSumEJA(J1,J2)
+            sage: x1 = J1.one()
+            sage: x2 = x1
+            sage: y1 = J2.one()
+            sage: y2 = y1
+            sage: x1.inner_product(x2)
+            3
+            sage: y1.inner_product(y2)
+            2
+            sage: J.one().inner_product(J.one())
+            5
+
+        """
+        (pi_left, pi_right) = self.projections()
+        x1 = pi_left(x)
+        x2 = pi_right(x)
+        y1 = pi_left(y)
+        y2 = pi_right(y)
+
+        return (x1.inner_product(y1) + x2.inner_product(y2))