]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add random_instance() method for algebras.
[sage.d.git] / mjo / eja / eja_algebra.py
index d0e7b074ce41a00f66d6f6dcd487653f3a8b1674..413128c843647e038e8869471daf29d84f1470f2 100644 (file)
@@ -60,6 +60,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         self._rank = rank
         self._natural_basis = natural_basis
 
+        # TODO: HACK for the charpoly.. needs redesign badly.
+        self._basis_normalizers = None
+
         if category is None:
             category = MagmaticAlgebras(field).FiniteDimensional()
             category = category.WithBasis().Unital()
@@ -120,11 +123,11 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         vector representations) back and forth faithfully::
 
             sage: set_random_seed()
-            sage: J = RealCartesianProductEJA(5)
+            sage: J = RealCartesianProductEJA.random_instance()
             sage: x = J.random_element()
             sage: J(x.to_vector().column()) == x
             True
-            sage: J = JordanSpinEJA(5)
+            sage: J = JordanSpinEJA.random_instance()
             sage: x = J.random_element()
             sage: J(x.to_vector().column()) == x
             True
@@ -152,6 +155,26 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return self.from_vector(coords)
 
 
+    @staticmethod
+    def _max_test_case_size():
+        """
+        Return an integer "size" that is an upper bound on the size of
+        this algebra when it is used in a random test
+        case. Unfortunately, the term "size" is quite vague -- when
+        dealing with `R^n` under either the Hadamard or Jordan spin
+        product, the "size" refers to the dimension `n`. When dealing
+        with a matrix algebra (real symmetric or complex/quaternion
+        Hermitian), it refers to the size of the matrix, which is
+        far less than the dimension of the underlying vector space.
+
+        We default to five in this class, which is safe in `R^n`. The
+        matrix algebra subclasses (or any class where the "size" is
+        interpreted to be far less than the dimension) should override
+        with a smaller number.
+        """
+        return 5
+
+
     def _repr_(self):
         """
         Return a string representation of ``self``.
@@ -224,6 +247,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return V.span_of_basis(b)
 
 
+
     @cached_method
     def _charpoly_coeff(self, i):
         """
@@ -234,6 +258,19 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         store the trace/determinant (a_{r-1} and a_{0} respectively)
         separate from the entire characteristic polynomial.
         """
+        if self._basis_normalizers is not None:
+             # Must be a matrix class?
+             # WARNING/TODO: this whole mess is mis-designed.
+             n = self.natural_basis_space().nrows()
+             field = self.base_ring().base_ring() # yeeeeaaaahhh
+             J = self.__class__(n, field, False)
+             (_,x,_,_) = J._charpoly_matrix_system()
+             p = J._charpoly_coeff(i)
+             # p might be missing some vars, have to substitute "optionally"
+             pairs = zip(x.base_ring().gens(), self._basis_normalizers)
+             substitutions = { v: v*c for (v,c) in pairs }
+             return p.subs(substitutions)
+
         (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
         R = A_of_x.base_ring()
         if i >= self.rank():
@@ -418,7 +455,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         """
         X = x.natural_representation()
         Y = y.natural_representation()
-        return self.__class__.natural_inner_product(X,Y)
+        return self.natural_inner_product(X,Y)
 
 
     def is_trivial(self):
@@ -621,6 +658,28 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             return s.random_element()
 
 
+    @classmethod
+    def random_instance(cls, field=QQ, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+
+        In subclasses for algebras that we know how to construct, this
+        is a shortcut for constructing test cases and examples.
+        """
+        if cls is FiniteDimensionalEuclideanJordanAlgebra:
+            # Red flag! But in theory we could do this I guess. The
+            # only finite-dimensional exceptional EJA is the
+            # octononions. So, we could just create an EJA from an
+            # associative matrix algebra (generated by a subset of
+            # elements) with the symmetric product. Or, we could punt
+            # to random_eja() here, override it in our subclasses, and
+            # not worry about it.
+            raise NotImplementedError
+
+        n = ZZ.random_element(1, cls._max_test_case_size())
+        return cls(n, field, **kwargs)
+
+
     def rank(self):
         """
         Return the rank of this EJA.
@@ -744,8 +803,7 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
     Our inner product satisfies the Jordan axiom::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: J = RealCartesianProductEJA(n)
+        sage: J = RealCartesianProductEJA.random_instance()
         sage: x = J.random_element()
         sage: y = J.random_element()
         sage: z = J.random_element()
@@ -775,13 +833,12 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
         over `R^n`::
 
             sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: J = RealCartesianProductEJA(n)
+            sage: J = RealCartesianProductEJA.random_instance()
             sage: x = J.random_element()
             sage: y = J.random_element()
             sage: X = x.natural_representation()
             sage: Y = y.natural_representation()
-            sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y)
+            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
             True
 
         """
@@ -824,18 +881,12 @@ def random_eja():
         Euclidean Jordan algebra of dimension...
 
     """
-
-    # The max_n component lets us choose different upper bounds on the
-    # value "n" that gets passed to the constructor. This is needed
-    # because e.g. R^{10} is reasonable to test, while the Hermitian
-    # 10-by-10 quaternion matrices are not.
-    (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
-                                   (JordanSpinEJA, 6),
-                                   (RealSymmetricEJA, 5),
-                                   (ComplexHermitianEJA, 4),
-                                   (QuaternionHermitianEJA, 3)])
-    n = ZZ.random_element(1, max_n)
-    return constructor(n, field=QQ)
+    classname = choice([RealCartesianProductEJA,
+                        JordanSpinEJA,
+                        RealSymmetricEJA,
+                        ComplexHermitianEJA,
+                        QuaternionHermitianEJA])
+    return classname.random_instance()
 
 
 
@@ -924,7 +975,7 @@ def _complex_hermitian_basis(n, field):
 
 
 
-def _quaternion_hermitian_basis(n, field, normalize):
+def _quaternion_hermitian_basis(n, field):
     """
     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
 
@@ -942,7 +993,7 @@ def _quaternion_hermitian_basis(n, field, normalize):
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
-        sage: B = _quaternion_hermitian_basis(n, QQ, False)
+        sage: B = _quaternion_hermitian_basis(n, QQ)
         sage: all( M.is_symmetric() for M in B )
         True
 
@@ -1015,7 +1066,8 @@ def _embed_complex_matrix(M):
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import _embed_complex_matrix
+        sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
+        ....:                                  ComplexHermitianEJA)
 
     EXAMPLES::
 
@@ -1037,7 +1089,8 @@ def _embed_complex_matrix(M):
     Embedding is a homomorphism (isomorphism, in fact)::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(5)
+        sage: n_max = ComplexHermitianEJA._max_test_case_size()
+        sage: n = ZZ.random_element(n_max)
         sage: F = QuadraticField(-1, 'i')
         sage: X = random_matrix(F, n)
         sage: Y = random_matrix(F, n)
@@ -1129,7 +1182,8 @@ def _embed_quaternion_matrix(M):
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
+        sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
+        ....:                                  QuaternionHermitianEJA)
 
     EXAMPLES::
 
@@ -1146,7 +1200,8 @@ def _embed_quaternion_matrix(M):
     Embedding is a homomorphism (isomorphism, in fact)::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(5)
+        sage: n_max = QuaternionHermitianEJA._max_test_case_size()
+        sage: n = ZZ.random_element(n_max)
         sage: Q = QuaternionAlgebra(QQ,-1,-1)
         sage: X = random_matrix(Q, n)
         sage: Y = random_matrix(Q, n)
@@ -1215,7 +1270,10 @@ def _unembed_quaternion_matrix(M):
     if not n.mod(4).is_zero():
         raise ValueError("the matrix 'M' must be a complex embedding")
 
-    Q = QuaternionAlgebra(QQ,-1,-1)
+    # Use the base ring of the matrix to ensure that its entries can be
+    # multiplied by elements of the quaternion algebra.
+    field = M.base_ring()
+    Q = QuaternionAlgebra(field,-1,-1)
     i,j,k = Q.gens()
 
     # Go top-left to bottom-right (reading order), converting every
@@ -1229,8 +1287,10 @@ def _unembed_quaternion_matrix(M):
                 raise ValueError('bad on-diagonal submatrix')
             if submat[0,1] != -submat[1,0].conjugate():
                 raise ValueError('bad off-diagonal submatrix')
-            z  = submat[0,0].real() + submat[0,0].imag()*i
-            z += submat[0,1].real()*j + submat[0,1].imag()*k
+            z  = submat[0,0].vector()[0]   # real part
+            z += submat[0,0].vector()[1]*i # imag part
+            z += submat[0,1].vector()[0]*j # real part
+            z += submat[0,1].vector()[1]*k # imag part
             elements.append(z)
 
     return matrix(Q, n/4, elements)
@@ -1271,7 +1331,8 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
     The dimension of this algebra is `(n^2 + n) / 2`::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
+        sage: n_max = RealSymmetricEJA._max_test_case_size()
+        sage: n = ZZ.random_element(1, n_max)
         sage: J = RealSymmetricEJA(n)
         sage: J.dimension() == (n^2 + n)/2
         True
@@ -1279,8 +1340,7 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
     The Jordan multiplication is what we think it is::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: J = RealSymmetricEJA(n)
+        sage: J = RealSymmetricEJA.random_instance()
         sage: x = J.random_element()
         sage: y = J.random_element()
         sage: actual = (x*y).natural_representation()
@@ -1300,28 +1360,29 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
     Our inner product satisfies the Jordan axiom::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: J = RealSymmetricEJA(n)
+        sage: J = RealSymmetricEJA.random_instance()
         sage: x = J.random_element()
         sage: y = J.random_element()
         sage: z = J.random_element()
         sage: (x*y).inner_product(z) == y.inner_product(x*z)
         True
 
-    Our basis is normalized with respect to the natural inner product::
+    Our natural basis is normalized with respect to the natural inner
+    product unless we specify otherwise::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: J = RealSymmetricEJA(n)
+        sage: J = RealSymmetricEJA.random_instance()
         sage: all( b.norm() == 1 for b in J.gens() )
         True
 
-    Left-multiplication operators are symmetric because they satisfy
-    the Jordan axiom::
+    Since our natural basis is normalized with respect to the natural
+    inner product, and since we know that this algebra is an EJA, any
+    left-multiplication operator's matrix will be symmetric because
+    natural->EJA basis representation is an isometry and within the EJA
+    the operator is self-adjoint by the Jordan axiom::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: x = RealSymmetricEJA(n).random_element()
+        sage: x = RealSymmetricEJA.random_instance().random_element()
         sage: x.operator().matrix().is_symmetric()
         True
 
@@ -1338,11 +1399,10 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
             p = z**2 - 2
             if p.is_irreducible():
                 field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
-            S = [ s.change_ring(field) for s in S ]
-            self._basis_denormalizers = tuple(
-                self.__class__.natural_inner_product(s,s).sqrt()
-                for s in S )
-            S = tuple( s/c for (s,c) in zip(S,self._basis_denormalizers) )
+                S = [ s.change_ring(field) for s in S ]
+            self._basis_normalizers = tuple(
+                ~(self.natural_inner_product(s,s).sqrt()) for s in S )
+            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )
 
         Qs = _multiplication_table_from_matrix_basis(S)
 
@@ -1353,6 +1413,9 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
                               natural_basis=S,
                               **kwargs)
 
+    @staticmethod
+    def _max_test_case_size():
+        return 5
 
 
 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
@@ -1371,7 +1434,8 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     The dimension of this algebra is `n^2`::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
+        sage: n_max = ComplexHermitianEJA._max_test_case_size()
+        sage: n = ZZ.random_element(1, n_max)
         sage: J = ComplexHermitianEJA(n)
         sage: J.dimension() == n^2
         True
@@ -1379,8 +1443,7 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     The Jordan multiplication is what we think it is::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: J = ComplexHermitianEJA(n)
+        sage: J = ComplexHermitianEJA.random_instance()
         sage: x = J.random_element()
         sage: y = J.random_element()
         sage: actual = (x*y).natural_representation()
@@ -1400,28 +1463,29 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     Our inner product satisfies the Jordan axiom::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: J = ComplexHermitianEJA(n)
+        sage: J = ComplexHermitianEJA.random_instance()
         sage: x = J.random_element()
         sage: y = J.random_element()
         sage: z = J.random_element()
         sage: (x*y).inner_product(z) == y.inner_product(x*z)
         True
 
-    Our basis is normalized with respect to the natural inner product::
+    Our natural basis is normalized with respect to the natural inner
+    product unless we specify otherwise::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,4)
-        sage: J = ComplexHermitianEJA(n)
+        sage: J = ComplexHermitianEJA.random_instance()
         sage: all( b.norm() == 1 for b in J.gens() )
         True
 
-    Left-multiplication operators are symmetric because they satisfy
-    the Jordan axiom::
+    Since our natural basis is normalized with respect to the natural
+    inner product, and since we know that this algebra is an EJA, any
+    left-multiplication operator's matrix will be symmetric because
+    natural->EJA basis representation is an isometry and within the EJA
+    the operator is self-adjoint by the Jordan axiom::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: x = ComplexHermitianEJA(n).random_element()
+        sage: x = ComplexHermitianEJA.random_instance().random_element()
         sage: x.operator().matrix().is_symmetric()
         True
 
@@ -1438,11 +1502,10 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
             p = z**2 - 2
             if p.is_irreducible():
                 field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
-            S = [ s.change_ring(field) for s in S ]
-            self._basis_denormalizers = tuple(
-                self.__class__.natural_inner_product(s,s).sqrt()
-                for s in S )
-            S = tuple( s/c for (s,c) in zip(S,self._basis_denormalizers) )
+                S = [ s.change_ring(field) for s in S ]
+            self._basis_normalizers = tuple(
+                ~(self.natural_inner_product(s,s).sqrt()) for s in S )
+            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )
 
         Qs = _multiplication_table_from_matrix_basis(S)
 
@@ -1454,6 +1517,10 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
                               **kwargs)
 
 
+    @staticmethod
+    def _max_test_case_size():
+        return 4
+
     @staticmethod
     def natural_inner_product(X,Y):
         Xu = _unembed_complex_matrix(X)
@@ -1461,6 +1528,8 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         # The trace need not be real; consider Xu = (i*I) and Yu = I.
         return ((Xu*Yu).trace()).vector()[0] # real part, I guess
 
+
+
 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
@@ -1474,10 +1543,11 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     TESTS:
 
-    The dimension of this algebra is `n^2`::
+    The dimension of this algebra is `2*n^2 - n`::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
+        sage: n_max = QuaternionHermitianEJA._max_test_case_size()
+        sage: n = ZZ.random_element(1, n_max)
         sage: J = QuaternionHermitianEJA(n)
         sage: J.dimension() == 2*(n^2) - n
         True
@@ -1485,8 +1555,7 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     The Jordan multiplication is what we think it is::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: J = QuaternionHermitianEJA(n)
+        sage: J = QuaternionHermitianEJA.random_instance()
         sage: x = J.random_element()
         sage: y = J.random_element()
         sage: actual = (x*y).natural_representation()
@@ -1506,17 +1575,50 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     Our inner product satisfies the Jordan axiom::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: J = QuaternionHermitianEJA(n)
+        sage: J = QuaternionHermitianEJA.random_instance()
         sage: x = J.random_element()
         sage: y = J.random_element()
         sage: z = J.random_element()
         sage: (x*y).inner_product(z) == y.inner_product(x*z)
         True
 
+    Our natural basis is normalized with respect to the natural inner
+    product unless we specify otherwise::
+
+        sage: set_random_seed()
+        sage: J = QuaternionHermitianEJA.random_instance()
+        sage: all( b.norm() == 1 for b in J.gens() )
+        True
+
+    Since our natural basis is normalized with respect to the natural
+    inner product, and since we know that this algebra is an EJA, any
+    left-multiplication operator's matrix will be symmetric because
+    natural->EJA basis representation is an isometry and within the EJA
+    the operator is self-adjoint by the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: x = QuaternionHermitianEJA.random_instance().random_element()
+        sage: x.operator().matrix().is_symmetric()
+        True
+
     """
     def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
-        S = _quaternion_hermitian_basis(n, field, normalize_basis)
+        S = _quaternion_hermitian_basis(n, field)
+
+        if n > 1 and normalize_basis:
+            # We'll need sqrt(2) to normalize the basis, and this
+            # winds up in the multiplication table, so the whole
+            # algebra needs to be over the field extension.
+            R = PolynomialRing(field, 'z')
+            z = R.gen()
+            p = z**2 - 2
+            if p.is_irreducible():
+                field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
+                S = [ s.change_ring(field) for s in S ]
+            self._basis_normalizers = tuple(
+                ~(self.natural_inner_product(s,s).sqrt()) for s in S )
+            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )
+
         Qs = _multiplication_table_from_matrix_basis(S)
 
         fdeja = super(QuaternionHermitianEJA, self)
@@ -1526,17 +1628,20 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
                               natural_basis=S,
                               **kwargs)
 
-    def inner_product(self, x, y):
-        # Since a+bi+cj+dk on the diagonal is represented as
-        #
-        #   a + bi +cj + dk = [  a  b  c  d]
-        #                     [ -b  a -d  c]
-        #                     [ -c  d  a -b]
-        #                     [ -d -c  b  a],
-        #
-        # we'll quadruple-count the "a" entries if we take the trace of
-        # the embedding.
-        return _matrix_ip(x,y)/4
+    @staticmethod
+    def _max_test_case_size():
+        return 3
+
+    @staticmethod
+    def natural_inner_product(X,Y):
+        Xu = _unembed_quaternion_matrix(X)
+        Yu = _unembed_quaternion_matrix(Y)
+        # The trace need not be real; consider Xu = (i*I) and Yu = I.
+        # The result will be a quaternion algebra element, which doesn't
+        # have a "vector" method, but does have coefficient_tuple() method
+        # that returns the coefficients of 1, i, j, and k -- in that order.
+        return ((Xu*Yu).trace()).coefficient_tuple()[0]
+
 
 
 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
@@ -1579,8 +1684,7 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
     Our inner product satisfies the Jordan axiom::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
-        sage: J = JordanSpinEJA(n)
+        sage: J = JordanSpinEJA.random_instance()
         sage: x = J.random_element()
         sage: y = J.random_element()
         sage: z = J.random_element()
@@ -1625,13 +1729,12 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
         over `R^n`::
 
             sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: J = JordanSpinEJA(n)
+            sage: J = JordanSpinEJA.random_instance()
             sage: x = J.random_element()
             sage: y = J.random_element()
             sage: X = x.natural_representation()
             sage: Y = y.natural_representation()
-            sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y)
+            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
             True
 
         """