]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: require the rank argument for an EJA, because we can't compute it.
[sage.d.git] / mjo / eja / eja_algebra.py
index 8ab6afa33e8b8ccb9241099f5420fd035bd7b65c..cd39704407deda564f857677a3439bb382f2bcb7 100644 (file)
@@ -5,11 +5,22 @@ are used in optimization, and have some additional nice methods beyond
 what can be supported in a general Jordan Algebra.
 """
 
-from sage.all import *
+
 
 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra
 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement
+from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.finite_dimensional_algebras_with_basis import FiniteDimensionalAlgebrasWithBasis
+from sage.functions.other import sqrt
+from sage.matrix.constructor import matrix
+from sage.misc.cachefunc import cached_method
+from sage.misc.prandom import choice
+from sage.modules.free_module import VectorSpace
+from sage.modules.free_module_element import vector
+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.structure.element import is_Matrix
 from sage.structure.category_object import normalize_names
 
@@ -21,10 +32,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
     def __classcall_private__(cls,
                               field,
                               mult_table,
+                              rank,
                               names='e',
                               assume_associative=False,
                               category=None,
-                              rank=None,
                               natural_basis=None):
         n = len(mult_table)
         mult_table = [b.base_extend(field) for b in mult_table]
@@ -45,20 +56,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         return fda.__classcall__(cls,
                                  field,
                                  mult_table,
+                                 rank=rank,
                                  assume_associative=assume_associative,
                                  names=names,
                                  category=cat,
-                                 rank=rank,
                                  natural_basis=natural_basis)
 
 
     def __init__(self,
                  field,
                  mult_table,
+                 rank,
                  names='e',
                  assume_associative=False,
                  category=None,
-                 rank=None,
                  natural_basis=None):
         """
         SETUP::
@@ -90,6 +101,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
     def _repr_(self):
         """
         Return a string representation of ``self``.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+        TESTS:
+
+        Ensure that it says what we think it says::
+
+            sage: JordanSpinEJA(2, field=QQ)
+            Euclidean Jordan algebra of degree 2 over Rational Field
+            sage: JordanSpinEJA(3, field=RDF)
+            Euclidean Jordan algebra of degree 3 over Real Double Field
+
         """
         fmt = "Euclidean Jordan algebra of degree {} over {}"
         return fmt.format(self.degree(), self.base_ring())
@@ -99,6 +124,21 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         """
         Guess a regular element. Needed to compute the basis for our
         characteristic polynomial coefficients.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja
+
+        TESTS:
+
+        Ensure that this hacky method succeeds for every algebra that we
+        know how to construct::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J._a_regular_element().is_regular()
+            True
+
         """
         gs = self.gens()
         z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
@@ -176,7 +216,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                                                     self._multiplication_table,
                                                     rank=r)
 
-        idmat = identity_matrix(J.base_ring(), n)
+        idmat = matrix.identity(J.base_ring(), n)
 
         W = self._charpoly_basis_space()
         W = W.change_ring(R.fraction_field())
@@ -196,10 +236,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         # We want the middle equivalent thing in our matrix, but use
         # the first equivalent thing instead so that we can pass in
         # standard coordinates.
-        x = J(vector(R, R.gens()))
-        l1 = [column_matrix(W.coordinates((x**k).vector())) for k in range(r)]
+        x = J(W(R.gens()))
+        l1 = [matrix.column(W.coordinates((x**k).vector())) for k in range(r)]
         l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
-        A_of_x = block_matrix(R, 1, n, (l1 + l2))
+        A_of_x = matrix.block(R, 1, n, (l1 + l2))
         xr = W.coordinates((x**r).vector())
         return (A_of_x, x, xr, A_of_x.det())
 
@@ -207,19 +247,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
     @cached_method
     def characteristic_polynomial(self):
         """
+        Return a characteristic polynomial that works for all elements
+        of this algebra.
 
-        .. WARNING::
-
-            This implementation doesn't guarantee that the polynomial
-            denominator in the coefficients is not identically zero, so
-            theoretically it could crash. The way that this is handled
-            in e.g. Faraut and Koranyi is to use a basis that guarantees
-            the denominator is non-zero. But, doing so requires knowledge
-            of at least one regular element, and we don't even know how
-            to do that. The trade-off is that, if we use the standard basis,
-            the resulting polynomial will accept the "usual" coordinates. In
-            other words, we don't have to do a change of basis before e.g.
-            computing the trace or determinant.
+        The resulting polynomial has `n+1` variables, where `n` is the
+        dimension of this algebra. The first `n` variables correspond to
+        the coordinates of an algebra element: when evaluated at the
+        coordinates of an algebra element with respect to a certain
+        basis, the result is a univariate polynomial (in the one
+        remaining variable ``t``), namely the characteristic polynomial
+        of that element.
 
         SETUP::
 
@@ -349,11 +386,61 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
     def rank(self):
         """
         Return the rank of this EJA.
+
+        ALGORITHM:
+
+        The author knows of no algorithm to compute the rank of an EJA
+        where only the multiplication table is known. In lieu of one, we
+        require the rank to be specified when the algebra is created,
+        and simply pass along that number here.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  random_eja)
+
+        EXAMPLES:
+
+        The rank of the Jordan spin algebra is always two::
+
+            sage: JordanSpinEJA(2).rank()
+            2
+            sage: JordanSpinEJA(3).rank()
+            2
+            sage: JordanSpinEJA(4).rank()
+            2
+
+        The rank of the `n`-by-`n` Hermitian real, complex, or
+        quaternion matrices is `n`::
+
+            sage: RealSymmetricEJA(2).rank()
+            2
+            sage: ComplexHermitianEJA(2).rank()
+            2
+            sage: QuaternionHermitianEJA(2).rank()
+            2
+            sage: RealSymmetricEJA(5).rank()
+            5
+            sage: ComplexHermitianEJA(5).rank()
+            5
+            sage: QuaternionHermitianEJA(5).rank()
+            5
+
+        TESTS:
+
+        Ensure that every EJA that we know how to construct has a
+        positive integer rank::
+
+            sage: set_random_seed()
+            sage: r = random_eja().rank()
+            sage: r in ZZ and r > 0
+            True
+
         """
-        if self._rank is None:
-            raise ValueError("no rank specified at genesis")
-        else:
-            return self._rank
+        return self._rank
 
 
     def vector_space(self):
@@ -401,7 +488,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             The identity in `S^n` is converted to the identity in the EJA::
 
                 sage: J = RealSymmetricEJA(3)
-                sage: I = identity_matrix(QQ,3)
+                sage: I = matrix.identity(QQ,3)
                 sage: J(I) == J.one()
                 True
 
@@ -1208,10 +1295,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
                 sage: B = 2*x0*x_bar.row()
                 sage: C = 2*x0*x_bar.column()
-                sage: D = identity_matrix(QQ, n-1)
+                sage: D = matrix.identity(QQ, n-1)
                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
                 sage: D = D + 2*x_bar.tensor_product(x_bar)
-                sage: Q = block_matrix(2,2,[A,B,C,D])
+                sage: Q = matrix.block(2,2,[A,B,C,D])
                 sage: Q == x.quadratic_representation().matrix()
                 True
 
@@ -1377,7 +1464,17 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             # are power-associative.
             #
             # TODO: choose generator names intelligently.
-            return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
+            #
+            # The rank is the highest possible degree of a minimal polynomial,
+            # and is bounded above by the dimension. We know in this case that
+            # there's an element whose minimal polynomial has the same degree
+            # as the space's dimension, so that must be its rank too.
+            return FiniteDimensionalEuclideanJordanAlgebra(
+                     F,
+                     mats,
+                     V.dimension(),
+                     assume_associative=True,
+                     names='f')
 
 
         def subalgebra_idempotent(self):
@@ -1854,7 +1951,7 @@ def _embed_complex_matrix(M):
         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
 
     # We can drop the imaginaries here.
-    return block_matrix(field.base_ring(), n, blocks)
+    return matrix.block(field.base_ring(), n, blocks)
 
 
 def _unembed_complex_matrix(M):
@@ -1970,7 +2067,7 @@ def _embed_quaternion_matrix(M):
 
     # We should have real entries by now, so use the realest field
     # we've got for the return value.
-    return block_matrix(quaternions.base_ring(), n, blocks)
+    return matrix.block(quaternions.base_ring(), n, blocks)
 
 
 def _unembed_quaternion_matrix(M):
@@ -2266,13 +2363,13 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
     @staticmethod
     def __classcall_private__(cls, n, field=QQ):
         Qs = []
-        id_matrix = identity_matrix(field, n)
+        id_matrix = matrix.identity(field, n)
         for i in xrange(n):
             ei = id_matrix.column(i)
-            Qi = zero_matrix(field, n)
+            Qi = matrix.zero(field, n)
             Qi.set_row(0, ei)
             Qi.set_column(0, ei)
-            Qi += diagonal_matrix(n, [ei[0]]*n)
+            Qi += matrix.diagonal(n, [ei[0]]*n)
             # The addition of the diagonal matrix adds an extra ei[0] in the
             # upper-left corner of the matrix.
             Qi[0,0] = Qi[0,0] * ~field(2)