]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/euclidean_jordan_algebra.py
eja: add apply_univariate_polynomial() for elements and use it in a test.
[sage.d.git] / mjo / eja / euclidean_jordan_algebra.py
index 4a2979248fc5e7b90938a0c18b8a7545bb780e97..da1145d1851c91cf1f7faef014ba8a519d5d3fcf 100644 (file)
@@ -71,6 +71,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         """
         self._rank = rank
         self._natural_basis = natural_basis
+        self._multiplication_table = mult_table
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
                      mult_table,
@@ -86,6 +87,54 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         return fmt.format(self.degree(), self.base_ring())
 
 
+    def characteristic_polynomial(self):
+        r = self.rank()
+        n = self.dimension()
+
+        names = ['X' + str(i) for i in range(1,n+1)]
+        R = PolynomialRing(self.base_ring(), names)
+        J = FiniteDimensionalEuclideanJordanAlgebra(R,
+                                                    self._multiplication_table,
+                                                    rank=r)
+        x0 = J.zero()
+        c = 1
+        for g in J.gens():
+            x0 += c*g
+            c +=1
+        if not x0.is_regular():
+            raise ValueError("don't know a regular element")
+
+        # Get the vector space (as opposed to module) so that
+        # span_of_basis() works.
+        V = x0.vector().parent().ambient_vector_space()
+        V1 = V.span_of_basis( (x0**k).vector() for k in range(r) )
+        B = V1.basis() + V1.complement().basis()
+        W = V.span_of_basis(B)
+
+        def e(k):
+            # The coordinates of e_k with respect to the basis B.
+            # But, the e_k are elements of B...
+            return identity_matrix(J.base_ring(), n).column(k-1).column()
+
+        # A matrix implementation 1
+        x = J(vector(R, R.gens()))
+        l1 = [column_matrix(W.coordinates((x**k).vector())) for k in range(r)]
+        l2 = [e(k) for k in range(r+1, n+1)]
+        A_of_x = block_matrix(1, n, (l1 + l2))
+        xr = W.coordinates((x**r).vector())
+        a = []
+        for i in range(n):
+            A_cols = A.columns()
+            A_cols[i] = xr
+            numerator = column_matrix(A.base_ring(), A_cols).det()
+            denominator = A.det()
+            ai = numerator/denominator
+            a.append(ai)
+
+        # Note: all entries past the rth should be zero.
+        return a
+
+
     def inner_product(self, x, y):
         """
         The inner product associated with this Euclidean Jordan algebra.
@@ -266,6 +315,49 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 return A( (self.operator_matrix()**(n-1))*self.vector() )
 
 
+        def apply_univariate_polynomial(self, p):
+            """
+            Apply the univariate polynomial ``p`` to this element.
+
+            A priori, SageMath won't allow us to apply a univariate
+            polynomial to an element of an EJA, because we don't know
+            that EJAs are rings (they are usually not associative). Of
+            course, we know that EJAs are power-associative, so the
+            operation is ultimately kosher. This function sidesteps
+            the CAS to get the answer we want and expect.
+
+            EXAMPLES::
+
+                sage: R = PolynomialRing(QQ, 't')
+                sage: t = R.gen(0)
+                sage: p = t^4 - t^3 + 5*t - 2
+                sage: J = RealCartesianProductEJA(5)
+                sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
+                True
+
+            TESTS:
+
+            We should always get back an element of the algebra::
+
+                sage: set_random_seed()
+                sage: p = PolynomialRing(QQ, 't').random_element()
+                sage: J = random_eja()
+                sage: x = J.random_element()
+                sage: x.apply_univariate_polynomial(p) in J
+                True
+
+            """
+            if len(p.variables()) > 1:
+                raise ValueError("not a univariate polynomial")
+            P = self.parent()
+            R = P.base_ring()
+            # Convert the coeficcients to the parent's base ring,
+            # because a priori they might live in an (unnecessarily)
+            # larger ring for which P.sum() would fail below.
+            cs = [ R(c) for c in p.coefficients(sparse=False) ]
+            return P.sum( cs[k]*(self**k) for k in range(len(cs)) )
+
+
         def characteristic_polynomial(self):
             """
             Return my characteristic polynomial (if I'm a regular
@@ -626,6 +718,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
         def minimal_polynomial(self):
             """
+            Return the minimal polynomial of this element,
+            as a function of the variable `t`.
+
             ALGORITHM:
 
             We restrict ourselves to the associative subalgebra
@@ -633,14 +728,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             polynomial of this element's operator matrix (in that
             subalgebra). This works by Baes Proposition 2.3.16.
 
-            EXAMPLES::
+            TESTS:
+
+            The minimal polynomial of the identity and zero elements are
+            always the same::
 
                 sage: set_random_seed()
-                sage: x = random_eja().random_element()
-                sage: x.degree() == x.minimal_polynomial().degree()
-                True
+                sage: J = random_eja()
+                sage: J.one().minimal_polynomial()
+                t - 1
+                sage: J.zero().minimal_polynomial()
+                t
 
-            ::
+            The degree of an element is (by one definition) the degree
+            of its minimal polynomial::
 
                 sage: set_random_seed()
                 sage: x = random_eja().random_element()
@@ -661,11 +762,19 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: y0 = y.vector()[0]
                 sage: y_bar = y.vector()[1:]
                 sage: actual = y.minimal_polynomial()
-                sage: x = SR.symbol('x', domain='real')
-                sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
+                sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
+                sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
                 sage: bool(actual == expected)
                 True
 
+            The minimal polynomial should always kill its element::
+
+                sage: set_random_seed()
+                sage: x = random_eja().random_element()
+                sage: p = x.minimal_polynomial()
+                sage: x.apply_univariate_polynomial(p)
+                0
+
             """
             V = self.span_of_powers()
             assoc_subalg = self.subalgebra_generated_by()
@@ -673,7 +782,11 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             # and subalgebra_generated_by() must be the same, and in
             # the same order!
             elt = assoc_subalg(V.coordinates(self.vector()))
-            return elt.operator_matrix().minimal_polynomial()
+
+            # We get back a symbolic polynomial in 'x' but want a real
+            # polynomial in 't'.
+            p_of_x = elt.operator_matrix().minimal_polynomial()
+            return p_of_x.change_variable_name('t')
 
 
         def natural_representation(self):
@@ -1110,12 +1223,17 @@ def random_eja():
         Euclidean Jordan algebra of degree...
 
     """
-    n = ZZ.random_element(1,5)
-    constructor = choice([RealCartesianProductEJA,
-                          JordanSpinEJA,
-                          RealSymmetricEJA,
-                          ComplexHermitianEJA,
-                          QuaternionHermitianEJA])
+
+    # 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)