eja: add an ungodly hack to get fast charpolys back.

diff --git a/mjo/eja/TODO b/mjo/eja/TODO
index 67f390ba24f06e6d0f419609ddb8654afd0a65bf..750cdc87a3f204f422aa53ec5fb9affe72fa9d1b 100644 (file)
--- a/mjo/eja/TODO
+++ b/mjo/eja/TODO
@@ -9,9 +9,5 @@
 5. Factor out the unit-norm basis (and operator symmetry) tests once
    all of the algebras pass.
 
-6. Create Element subclasses for the matrix EJAs, and then override
-   their characteristic_polynomial() method to create a new algebra
-   over the rationals (with a non-normalized basis). We can then
-   compute the charpoly quickly by passing the natural representation
-   of the given element into the new algebra and computing its charpoly
-   there. (Relies on the theory to ensure that the charpolys are equal.)
\ No newline at end of file
+6. Refactor the current ungodly fast charpoly hack (relies on the
+   theory to ensure that the charpolys are equal.)

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index d0e7b074ce41a00f66d6f6dcd487653f3a8b1674..698aa37e2d6698c680ec51d3db1b272e16b6df4d 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -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():
@@ -1339,10 +1356,10 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
             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()
+            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_denormalizers) )
+            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )
 
         Qs = _multiplication_table_from_matrix_basis(S)
 
@@ -1439,10 +1456,10 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
             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()
+            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_denormalizers) )
+            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )
 
         Qs = _multiplication_table_from_matrix_basis(S)