]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: speed up det() for Cartesian product elements.
[sage.d.git] / mjo / eja / eja_algebra.py
index 3361bfaaeb6490f425ac383d0b14904b2e02b76c..2bad32c2f500193e4126b7c5e209c0acb3116ede 100644 (file)
@@ -166,9 +166,10 @@ from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
-from mjo.eja.eja_element import FiniteDimensionalEJAElement
+from mjo.eja.eja_element import (CartesianProductEJAElement,
+                                 FiniteDimensionalEJAElement)
 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
-from mjo.eja.eja_utils import _all2list, _mat2vec
+from mjo.eja.eja_utils import _all2list
 
 def EuclideanJordanAlgebras(field):
     r"""
@@ -1281,7 +1282,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         #
         # Of course, matrices aren't vectors in sage, so we have to
         # appeal to the "long vectors" isometry.
-        oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
+
+        V = VectorSpace(self.base_ring(), self.dimension()**2)
+        oper_vecs = [ V(g.operator().matrix().list()) for g in self.gens() ]
 
         # Now we use basic linear algebra to find the coefficients,
         # of the matrices-as-vectors-linear-combination, which should
@@ -1291,7 +1294,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # We used the isometry on the left-hand side already, but we
         # still need to do it for the right-hand side. Recall that we
         # wanted something that summed to the identity matrix.
-        b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
+        b = V( matrix.identity(self.base_ring(), self.dimension()).list() )
 
         # Now if there's an identity element in the algebra, this
         # should work. We solve on the left to avoid having to
@@ -1743,6 +1746,15 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        check_field=False,
                                        check_axioms=False)
 
+    def rational_algebra(self):
+        # Using None as a flag here (rather than just assigning "self"
+        # to self._rational_algebra by default) feels a little bit
+        # more sane to me in a garbage-collected environment.
+        if self._rational_algebra is None:
+            return self
+        else:
+            return self._rational_algebra
+
     @cached_method
     def _charpoly_coefficients(self):
         r"""
@@ -1767,18 +1779,15 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             Algebraic Real Field
 
         """
-        if self._rational_algebra is None:
-            # There's no need to construct *another* algebra over the
-            # rationals if this one is already over the
-            # rationals. Likewise, if we never orthonormalized our
-            # basis, we might as well just use the given one.
+        if self.rational_algebra() is self:
+            # Bypass the hijinks if they won't benefit us.
             return super()._charpoly_coefficients()
 
         # Do the computation over the rationals. The answer will be
         # the same, because all we've done is a change of basis.
         # Then, change back from QQ to our real base ring
         a = ( a_i.change_ring(self.base_ring())
-              for a_i in self._rational_algebra._charpoly_coefficients() )
+              for a_i in self.rational_algebra()._charpoly_coefficients() )
 
         # Otherwise, convert the coordinate variables back to the
         # deorthonormalized ones.
@@ -2039,7 +2048,6 @@ class MatrixEJA(FiniteDimensionalEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-
         super().__init__(self._denormalized_basis(matrix_space),
                          self.jordan_product,
                          self.trace_inner_product,
@@ -2134,20 +2142,13 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         return cls(n, **kwargs)
 
     def __init__(self, n, field=AA, **kwargs):
-        # We know this is a valid EJA, but will double-check
-        # if the user passes check_axioms=True.
-        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
-
         A = MatrixSpace(field, n)
         super().__init__(A, **kwargs)
 
         from mjo.eja.eja_cache import real_symmetric_eja_coeffs
         a = real_symmetric_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 
@@ -2224,10 +2225,6 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     """
     def __init__(self, n, field=AA, **kwargs):
-        # We know this is a valid EJA, but will double-check
-        # if the user passes check_axioms=True.
-        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
-
         from mjo.hurwitz import ComplexMatrixAlgebra
         A = ComplexMatrixAlgebra(n, scalars=field)
         super().__init__(A, **kwargs)
@@ -2235,10 +2232,7 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
         a = complex_hermitian_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
     @staticmethod
     def _max_random_instance_size(max_dimension):
@@ -2317,10 +2311,6 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     """
     def __init__(self, n, field=AA, **kwargs):
-        # We know this is a valid EJA, but will double-check
-        # if the user passes check_axioms=True.
-        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
-
         from mjo.hurwitz import QuaternionMatrixAlgebra
         A = QuaternionMatrixAlgebra(n, scalars=field)
         super().__init__(A, **kwargs)
@@ -2328,10 +2318,7 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs
         a = quaternion_hermitian_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 
@@ -2487,10 +2474,7 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
         from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs
         a = octonion_hermitian_eja_coeffs(self)
         if a is not None:
-            if self._rational_algebra is None:
-                self._charpoly_coefficients.set_cache(a)
-            else:
-                self._rational_algebra._charpoly_coefficients.set_cache(a)
+            self.rational_algebra()._charpoly_coefficients.set_cache(a)
 
 
 class AlbertEJA(OctonionHermitianEJA):
@@ -3106,6 +3090,7 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         sage: actual == expected             # long time
         True
     """
+    Element = CartesianProductEJAElement
     def __init__(self, factors, **kwargs):
         m = len(factors)
         if m == 0:
@@ -3180,6 +3165,8 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         self._inner_product_matrix = matrix.block_diagonal(
             [J._inner_product_matrix for J in factors]
         )
+        self._inner_product_matrix._cache = {'hermitian': True}
+        self._inner_product_matrix.set_immutable()
 
         # Building the multiplication table is a bit more tricky
         # because we have to embed the entries of the factors'
@@ -3207,6 +3194,34 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         ones = tuple(J.one().to_matrix() for J in factors)
         self.one.set_cache(self(ones))
 
+    def _sets_keys(self):
+        r"""
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  RealSymmetricEJA)
+
+        TESTS:
+
+        The superclass uses ``_sets_keys()`` to implement its
+        ``cartesian_factors()`` method::
+
+            sage: J1 = RealSymmetricEJA(2,
+            ....:                       field=QQ,
+            ....:                       orthonormalize=False,
+            ....:                       prefix="a")
+            sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J = cartesian_product([J1,J2])
+            sage: x = sum(i*J.gens()[i] for i in range(len(J.gens())))
+            sage: x.cartesian_factors()
+            (a1 + 2*a2, 3*b0 + 4*b1 + 5*b2 + 6*b3)
+
+        """
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct,
+        # but returning a tuple instead of a list.
+        return tuple(range(len(self.cartesian_factors())))
+
     def cartesian_factors(self):
         # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
         return self._sets
@@ -3424,9 +3439,9 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+        ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
-        ....:                                  OctonionHermitianEJA,
         ....:                                  RealSymmetricEJA)
 
     EXAMPLES:
@@ -3447,28 +3462,38 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
     The ``cartesian_product()`` function only uses the first factor to
     decide where the result will live; thus we have to be careful to
-    check that all factors do indeed have a `_rational_algebra` member
-    before we try to access it::
-
-        sage: J1 = OctonionHermitianEJA(1) # no rational basis
-        sage: J2 = HadamardEJA(2)
-        sage: cartesian_product([J1,J2])
-        Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
-        (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
-        sage: cartesian_product([J2,J1])
-        Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
-        (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+    check that all factors do indeed have a ``rational_algebra()`` method
+    before we construct an algebra that claims to have a rational basis::
+
+        sage: J1 = HadamardEJA(2)
+        sage: jp = lambda X,Y: X*Y
+        sage: ip = lambda X,Y: X[0,0]*Y[0,0]
+        sage: b1 = matrix(QQ, [[1]])
+        sage: J2 = FiniteDimensionalEJA((b1,), jp, ip)
+        sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA
+        Euclidean Jordan algebra of dimension 1 over Algebraic Real
+        Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic
+        Real Field
+        sage: cartesian_product([J1,J2]) # factor one is RationalBasisEJA
+        Traceback (most recent call last):
+        ...
+        ValueError: factor not a RationalBasisEJA
 
     """
     def __init__(self, algebras, **kwargs):
+        if not all( hasattr(r, "rational_algebra") for r in algebras ):
+            raise ValueError("factor not a RationalBasisEJA")
+
         CartesianProductEJA.__init__(self, algebras, **kwargs)
 
-        self._rational_algebra = None
-        if self.vector_space().base_field() is not QQ:
-            if all( hasattr(r, "_rational_algebra") for r in algebras ):
-                self._rational_algebra = cartesian_product([
-                    r._rational_algebra for r in algebras
-                ])
+    @cached_method
+    def rational_algebra(self):
+        if self.base_ring() is QQ:
+            return self
+
+        return cartesian_product([
+            r.rational_algebra() for r in self.cartesian_factors()
+        ])
 
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA