]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add an ungodly hack to get fast charpolys back.
[sage.d.git] / mjo / eja / eja_algebra.py
index adf9581c541bf81e02c612d281a4d671c3b9df74..698aa37e2d6698c680ec51d3db1b272e16b6df4d 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()
@@ -224,6 +227,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return V.span_of_basis(b)
 
 
+
     @cached_method
     def _charpoly_coeff(self, i):
         """
@@ -234,6 +238,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():
@@ -416,9 +433,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
-        if (not x in self) or (not y in self):
-            raise TypeError("arguments must live in this algebra")
-        return x.trace_inner_product(y)
+        X = x.natural_representation()
+        Y = y.natural_representation()
+        return self.__class__.natural_inner_product(X,Y)
 
 
     def is_trivial(self):
@@ -537,6 +554,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             return self._natural_basis[0].matrix_space()
 
 
+    @staticmethod
+    def natural_inner_product(X,Y):
+        """
+        Compute the inner product of two naturally-represented elements.
+
+        For example in the real symmetric matrix EJA, this will compute
+        the trace inner-product of two n-by-n symmetric matrices. The
+        default should work for the real cartesian product EJA, the
+        Jordan spin EJA, and the real symmetric matrices. The others
+        will have to be overridden.
+        """
+        return (X.conjugate_transpose()*Y).trace()
+
+
     @cached_method
     def one(self):
         """
@@ -748,7 +779,30 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
         return fdeja.__init__(field, mult_table, rank=n, **kwargs)
 
     def inner_product(self, x, y):
-        return _usual_ip(x,y)
+        """
+        Faster to reimplement than to use natural representations.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+
+        TESTS:
+
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: J = RealCartesianProductEJA(n)
+            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)
+            True
+
+        """
+        return x.to_vector().inner_product(y.to_vector())
 
 
 def random_eja():
@@ -802,7 +856,7 @@ def random_eja():
 
 
 
-def _real_symmetric_basis(n, field, normalize):
+def _real_symmetric_basis(n, field):
     """
     Return a basis for the space of real symmetric n-by-n matrices.
 
@@ -814,7 +868,7 @@ def _real_symmetric_basis(n, field, normalize):
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
-        sage: B = _real_symmetric_basis(n, QQbar, False)
+        sage: B = _real_symmetric_basis(n, QQ)
         sage: all( M.is_symmetric() for M in  B)
         True
 
@@ -829,13 +883,11 @@ def _real_symmetric_basis(n, field, normalize):
                 Sij = Eij
             else:
                 Sij = Eij + Eij.transpose()
-            if normalize:
-                Sij = Sij / _real_symmetric_matrix_ip(Sij,Sij).sqrt()
             S.append(Sij)
     return tuple(S)
 
 
-def _complex_hermitian_basis(n, field, normalize):
+def _complex_hermitian_basis(n, field):
     """
     Returns a basis for the space of complex Hermitian n-by-n matrices.
 
@@ -854,7 +906,7 @@ def _complex_hermitian_basis(n, field, normalize):
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
         sage: field = QuadraticField(2, 'sqrt2')
-        sage: B = _complex_hermitian_basis(n, field, False)
+        sage: B = _complex_hermitian_basis(n, field)
         sage: all( M.is_symmetric() for M in  B)
         True
 
@@ -877,8 +929,7 @@ def _complex_hermitian_basis(n, field, normalize):
                 Sij = _embed_complex_matrix(Eij)
                 S.append(Sij)
             else:
-                # Beware, orthogonal but not normalized! The second one
-                # has a minus because it's conjugated.
+                # The second one has a minus because it's conjugated.
                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
                 S.append(Sij_real)
                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
@@ -886,11 +937,7 @@ def _complex_hermitian_basis(n, field, normalize):
 
     # Since we embedded these, we can drop back to the "field" that we
     # started with instead of the complex extension "F".
-    S = [ s.change_ring(field) for s in S ]
-    if normalize:
-        S = [ s / _complex_hermitian_matrix_ip(s,s).sqrt() for s in S ]
-
-    return tuple(S)
+    return tuple( s.change_ring(field) for s in S )
 
 
 
@@ -1206,10 +1253,6 @@ def _unembed_quaternion_matrix(M):
     return matrix(Q, n/4, elements)
 
 
-# The usual inner product on R^n.
-def _usual_ip(x,y):
-    return x.to_vector().inner_product(y.to_vector())
-
 # The inner product used for the real symmetric simple EJA.
 # We keep it as a separate function because e.g. the complex
 # algebra uses the same inner product, except divided by 2.
@@ -1218,15 +1261,6 @@ def _matrix_ip(X,Y):
     Y_mat = Y.natural_representation()
     return (X_mat*Y_mat).trace()
 
-def _real_symmetric_matrix_ip(X,Y):
-    return (X*Y).trace()
-
-def _complex_hermitian_matrix_ip(X,Y):
-    # This takes EMBEDDED matrices.
-    Xu = _unembed_complex_matrix(X)
-    Yu = _unembed_complex_matrix(Y)
-    # The trace need not be real; consider Xu = (i*I) and Yu = I.
-    return ((Xu*Yu).trace()).vector()[0] # real part, I guess
 
 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
@@ -1310,7 +1344,9 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     """
     def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
-        if n > 1:
+        S = _real_symmetric_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.
@@ -1319,8 +1355,12 @@ 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_normalizers = 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_normalizers) )
 
-        S = _real_symmetric_basis(n, field, normalize_basis)
         Qs = _multiplication_table_from_matrix_basis(S)
 
         fdeja = super(RealSymmetricEJA, self)
@@ -1330,10 +1370,6 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
                               natural_basis=S,
                               **kwargs)
 
-    def inner_product(self, x, y):
-        X = x.natural_representation()
-        Y = y.natural_representation()
-        return _real_symmetric_matrix_ip(X,Y)
 
 
 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
@@ -1408,7 +1444,9 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     """
     def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
-        if n > 1:
+        S = _complex_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.
@@ -1417,8 +1455,12 @@ 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_normalizers = 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_normalizers) )
 
-        S = _complex_hermitian_basis(n, field, normalize_basis)
         Qs = _multiplication_table_from_matrix_basis(S)
 
         fdeja = super(ComplexHermitianEJA, self)
@@ -1429,11 +1471,12 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
                               **kwargs)
 
 
-    def inner_product(self, x, y):
-        X = x.natural_representation()
-        Y = y.natural_representation()
-        return _complex_hermitian_matrix_ip(X,Y)
-
+    @staticmethod
+    def natural_inner_product(X,Y):
+        Xu = _unembed_complex_matrix(X)
+        Yu = _unembed_complex_matrix(Y)
+        # 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):
     """
@@ -1586,4 +1629,27 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
         return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
 
     def inner_product(self, x, y):
-        return _usual_ip(x,y)
+        """
+        Faster to reimplement than to use natural representations.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+        TESTS:
+
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: J = JordanSpinEJA(n)
+            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)
+            True
+
+        """
+        return x.to_vector().inner_product(y.to_vector())