eja: use izip() instead of zip() in a few places.

[sage.d.git] / mjo / eja / eja_algebra.py

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index ae5baa0eff4a64dac97944bbdaf5a8ec132495bd..79ccc7904830dc739e4af0a532c1bace2a801936 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -5,7 +5,7 @@ are used in optimization, and have some additional nice methods beyond
 what can be supported in a general Jordan Algebra.
 """
 
-from itertools import repeat
+from itertools import izip, repeat
 
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
@@ -61,9 +61,6 @@ 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()
@@ -259,19 +256,6 @@ 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():
@@ -425,7 +409,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             # assign a[r] goes out-of-bounds.
             a.append(1) # corresponds to x^r
 
-        return sum( a[k]*(t**k) for k in range(len(a)) )
+        return sum( a[k]*(t**k) for k in xrange(len(a)) )
 
 
     def inner_product(self, x, y):
@@ -507,7 +491,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         """
         M = list(self._multiplication_table) # copy
-        for i in range(len(M)):
+        for i in xrange(len(M)):
             # M had better be "square"
             M[i] = [self.monomial(i)] + M[i]
         M = [["*"] + list(self.gens())] + M
@@ -694,7 +678,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             # not worry about it.
             raise NotImplementedError
 
-        n = ZZ.random_element(1, cls._max_test_case_size())
+        n = ZZ.random_element(cls._max_test_case_size()) + 1
         return cls(n, field, **kwargs)
 
 
@@ -815,8 +799,8 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     def __init__(self, n, field=QQ, **kwargs):
         V = VectorSpace(field, n)
-        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
-                       for i in range(n) ]
+        mult_table = [ [ V.gen(i)*(i == j) for j in xrange(n) ]
+                       for i in xrange(n) ]
 
         fdeja = super(RealCartesianProductEJA, self)
         return fdeja.__init__(field, mult_table, rank=n, **kwargs)
@@ -899,7 +883,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     def _max_test_case_size():
         # Play it safe, since this will be squared and the underlying
         # field can have dimension 4 (quaternions) too.
-        return 3
+        return 2
 
     @classmethod
     def _denormalized_basis(cls, n, field):
@@ -908,6 +892,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
         S = self._denormalized_basis(n, field)
 
+        # Used in this class's fast _charpoly_coeff() override.
+        self._basis_normalizers = None
+
         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
@@ -932,6 +919,30 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
                               **kwargs)
 
 
+    @cached_method
+    def _charpoly_coeff(self, i):
+        """
+        Override the parent method with something that tries to compute
+        over a faster (non-extension) field.
+        """
+        if self._basis_normalizers is None:
+            # We didn't normalize, so assume that the basis we started
+            # with had entries in a nice field.
+            return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i)
+        else:
+            # If we didn't unembed first, this number would be wrong
+            # by a power-of-two factor for complex/quaternion matrices.
+            n = self.real_unembed(self.natural_basis_space().zero()).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 = izip(x.base_ring().gens(), self._basis_normalizers)
+            substitutions = { v: v*c for (v,c) in pairs }
+            return p.subs(substitutions)
+
+
     @staticmethod
     def multiplication_table_from_matrix_basis(basis):
         """
@@ -953,9 +964,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         V = VectorSpace(field, dimension**2)
         W = V.span_of_basis( _mat2vec(s) for s in basis )
         n = len(basis)
-        mult_table = [[W.zero() for j in range(n)] for i in range(n)]
-        for i in range(n):
-            for j in range(n):
+        mult_table = [[W.zero() for j in xrange(n)] for i in xrange(n)]
+        for i in xrange(n):
+            for j in xrange(n):
                 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
                 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
 
@@ -1138,7 +1149,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
 
     @staticmethod
     def _max_test_case_size():
-        return 5 # Dimension 10
+        return 4 # Dimension 10
 
 
 
@@ -1726,9 +1737,9 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     def __init__(self, n, field=QQ, **kwargs):
         V = VectorSpace(field, n)
-        mult_table = [[V.zero() for j in range(n)] for i in range(n)]
-        for i in range(n):
-            for j in range(n):
+        mult_table = [[V.zero() for j in xrange(n)] for i in xrange(n)]
+        for i in xrange(n):
+            for j in xrange(n):
                 x = V.gen(i)
                 y = V.gen(j)
                 x0 = x[0]