]> 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 4a1c6f94f589c80b9021b71ee01693765aabcaeb..2bad32c2f500193e4126b7c5e209c0acb3116ede 100644 (file)
@@ -166,9 +166,21 @@ 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"""
+    The category of Euclidean Jordan algebras over ``field``, which
+    must be a subfield of the real numbers. For now this is just a
+    convenient wrapper around all of the other category axioms that
+    apply to all EJAs.
+    """
+    category = MagmaticAlgebras(field).FiniteDimensional()
+    category = category.WithBasis().Unital().Commutative()
+    return category
 
 class FiniteDimensionalEJA(CombinatorialFreeModule):
     r"""
@@ -228,6 +240,26 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
     """
     Element = FiniteDimensionalEJAElement
 
+    @staticmethod
+    def _check_input_field(field):
+        if not field.is_subring(RR):
+            # Note: this does return true for the real algebraic
+            # field, the rationals, and any quadratic field where
+            # we've specified a real embedding.
+            raise ValueError("scalar field is not real")
+
+    @staticmethod
+    def _check_input_axioms(basis, jordan_product, inner_product):
+        if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
+                    for bi in basis
+                    for bj in basis ):
+            raise ValueError("Jordan product is not commutative")
+
+        if not all( inner_product(bi,bj) == inner_product(bj,bi)
+                    for bi in basis
+                    for bj in basis ):
+            raise ValueError("inner-product is not commutative")
+
     def __init__(self,
                  basis,
                  jordan_product,
@@ -236,7 +268,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  matrix_space=None,
                  orthonormalize=True,
                  associative=None,
-                 cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
                  prefix="b"):
@@ -244,30 +275,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         n = len(basis)
 
         if check_field:
-            if not field.is_subring(RR):
-                # Note: this does return true for the real algebraic
-                # field, the rationals, and any quadratic field where
-                # we've specified a real embedding.
-                raise ValueError("scalar field is not real")
+            self._check_input_field(field)
 
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
             # This has to be done before we build the multiplication
             # and inner-product tables/matrices, because we take
             # advantage of symmetry in the process.
-            if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
-                        for bi in basis
-                        for bj in basis ):
-                raise ValueError("Jordan product is not commutative")
-
-            if not all( inner_product(bi,bj) == inner_product(bj,bi)
-                        for bi in basis
-                        for bj in basis ):
-                raise ValueError("inner-product is not commutative")
-
-
-        category = MagmaticAlgebras(field).FiniteDimensional()
-        category = category.WithBasis().Unital().Commutative()
+            self._check_input_axioms(basis, jordan_product, inner_product)
 
         if n <= 1:
             # All zero- and one-dimensional algebras are just the real
@@ -286,14 +301,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                                for bj in basis
                                for bk in basis)
 
+        category = EuclideanJordanAlgebras(field)
+
         if associative:
             # Element subalgebras can take advantage of this.
             category = category.Associative()
-        if cartesian_product:
-            # Use join() here because otherwise we only get the
-            # "Cartesian product of..." and not the things themselves.
-            category = category.join([category,
-                                      category.CartesianProducts()])
 
         # Call the superclass constructor so that we can use its from_vector()
         # method to build our multiplication table.
@@ -368,7 +380,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # Now we actually compute the multiplication and inner-product
         # tables/matrices using the possibly-orthonormalized basis.
         self._inner_product_matrix = matrix.identity(field, n)
-        self._multiplication_table = [ [0 for j in range(i+1)]
+        zed = self.zero()
+        self._multiplication_table = [ [zed for j in range(i+1)]
                                        for i in range(n) ]
 
         # Note: the Jordan and inner-products are defined in terms
@@ -779,8 +792,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             if elt.parent().superalgebra() == self:
                 return elt.superalgebra_element()
 
-        if hasattr(elt, 'column'):
-            # Convert a vector into a column-matrix...
+        if hasattr(elt, 'sparse_vector'):
+            # Convert a vector into a column-matrix. We check for
+            # "sparse_vector" and not "column" because matrices also
+            # have a "column" method.
             elt = elt.column()
 
         if elt not in self.matrix_space():
@@ -1267,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
@@ -1277,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
@@ -1407,7 +1424,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # corresponding to trivial spaces (e.g. it returns only the
         # eigenspace corresponding to lambda=1 if you take the
         # decomposition relative to the identity element).
-        trivial = self.subalgebra(())
+        trivial = self.subalgebra((), check_axioms=False)
         J0 = trivial                          # eigenvalue zero
         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
         J1 = trivial                          # eigenvalue one
@@ -1729,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"""
@@ -1753,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.
@@ -2025,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,
@@ -2120,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)
 
 
 
@@ -2210,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)
@@ -2221,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):
@@ -2303,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)
@@ -2314,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)
 
 
 
@@ -2473,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):
@@ -2935,6 +2933,7 @@ class CartesianProductEJA(FiniteDimensionalEJA):
 
         sage: from mjo.eja.eja_algebra import (random_eja,
         ....:                                  CartesianProductEJA,
+        ....:                                  ComplexHermitianEJA,
         ....:                                  HadamardEJA,
         ....:                                  JordanSpinEJA,
         ....:                                  RealSymmetricEJA)
@@ -3046,6 +3045,28 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         | b2 || 0  | 0  | b2 |
         +----++----+----+----+
 
+    The "matrix space" of a Cartesian product always consists of
+    ordered pairs (or triples, or...) whose components are the
+    matrix spaces of its factors::
+
+            sage: J1 = HadamardEJA(2)
+            sage: J2 = ComplexHermitianEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.matrix_space()
+            The Cartesian product of (Full MatrixSpace of 2 by 1 dense
+            matrices over Algebraic Real Field, Module of 2 by 2 matrices
+            with entries in Algebraic Field over the scalar ring Algebraic
+            Real Field)
+            sage: J.one().to_matrix()[0]
+            [1]
+            [1]
+            sage: J.one().to_matrix()[1]
+            +---+---+
+            | 1 | 0 |
+            +---+---+
+            | 0 | 1 |
+            +---+---+
+
     TESTS:
 
     All factors must share the same base field::
@@ -3068,11 +3089,8 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         sage: expected = J.one()             # long time
         sage: actual == expected             # long time
         True
-
     """
-    Element = FiniteDimensionalEJAElement
-
-
+    Element = CartesianProductEJAElement
     def __init__(self, factors, **kwargs):
         m = len(factors)
         if m == 0:
@@ -3084,68 +3102,126 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         if not all( J.base_ring() == field for J in factors ):
             raise ValueError("all factors must share the same base field")
 
+        # Figure out the category to use.
         associative = all( f.is_associative() for f in factors )
-
-        # Compute my matrix space. This category isn't perfect, but
-        # is good enough for what we need to do.
+        category = EuclideanJordanAlgebras(field)
+        if associative: category = category.Associative()
+        category = category.join([category, category.CartesianProducts()])
+
+        # Compute my matrix space.  We don't simply use the
+        # ``cartesian_product()`` functor here because it acts
+        # differently on SageMath MatrixSpaces and our custom
+        # MatrixAlgebras, which are CombinatorialFreeModules. We
+        # always want the result to be represented (and indexed) as an
+        # ordered tuple. This category isn't perfect, but is good
+        # enough for what we need to do.
         MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
         MS_cat = MS_cat.Unital().CartesianProducts()
         MS_factors = tuple( J.matrix_space() for J in factors )
         from sage.sets.cartesian_product import CartesianProduct
-        MS = CartesianProduct(MS_factors, MS_cat)
+        self._matrix_space = CartesianProduct(MS_factors, MS_cat)
 
-        basis = []
-        zero = MS.zero()
+        self._matrix_basis = []
+        zero = self._matrix_space.zero()
         for i in range(m):
             for b in factors[i].matrix_basis():
                 z = list(zero)
                 z[i] = b
-                basis.append(z)
+                self._matrix_basis.append(z)
 
-        basis = tuple( MS(b) for b in basis )
+        self._matrix_basis = tuple( self._matrix_space(b)
+                                    for b in self._matrix_basis )
+        n = len(self._matrix_basis)
 
-        # Define jordan/inner products that operate on that matrix_basis.
-        def jordan_product(x,y):
-            return MS(tuple(
-                (factors[i](x[i])*factors[i](y[i])).to_matrix()
-                for i in range(m)
-            ))
-
-        def inner_product(x, y):
-            return sum(
-                factors[i](x[i]).inner_product(factors[i](y[i]))
-                for i in range(m)
-            )
+        # We already have what we need for the super-superclass constructor.
+        CombinatorialFreeModule.__init__(self,
+                                         field,
+                                         range(n),
+                                         prefix="b",
+                                         category=category,
+                                         bracket=False)
 
-        # There's no need to check the field since it already came
-        # from an EJA. Likewise the axioms are guaranteed to be
-        # satisfied, unless the guy writing this class sucks.
-        #
-        # If you want the basis to be orthonormalized, orthonormalize
-        # the factors.
-        FiniteDimensionalEJA.__init__(self,
-                                      basis,
-                                      jordan_product,
-                                      inner_product,
-                                      field=field,
-                                      matrix_space=MS,
-                                      orthonormalize=False,
-                                      associative=associative,
-                                      cartesian_product=True,
-                                      check_field=False,
-                                      check_axioms=False)
+        # Now create the vector space for the algebra, which will have
+        # its own set of non-ambient coordinates (in terms of the
+        # supplied basis).
+        degree = sum( f._matrix_span.ambient_vector_space().degree()
+                      for f in factors )
+        V = VectorSpace(field, degree)
+        vector_basis = tuple( V(_all2list(b)) for b in self._matrix_basis )
+
+        # Save the span of our matrix basis (when written out as long
+        # vectors) because otherwise we'll have to reconstruct it
+        # every time we want to coerce a matrix into the algebra.
+        self._matrix_span = V.span_of_basis( vector_basis, check=False)
 
         # Since we don't (re)orthonormalize the basis, the FDEJA
         # constructor is going to set self._deortho_matrix to the
         # identity matrix. Here we set it to the correct value using
         # the deortho matrices from our factors.
-        self._deortho_matrix = matrix.block_diagonal( [J._deortho_matrix
-                                                       for J in factors] )
+        self._deortho_matrix = matrix.block_diagonal(
+            [J._deortho_matrix for J in factors]
+        )
+
+        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'
+        # multiplication tables into the product EJA.
+        zed = self.zero()
+        self._multiplication_table = [ [zed for j in range(i+1)]
+                                       for i in range(n) ]
+
+        # Keep track of an offset that tallies the dimensions of all
+        # previous factors. If the second factor is dim=2 and if the
+        # first one is dim=3, then we want to skip the first 3x3 block
+        # when copying the multiplication table for the second factor.
+        offset = 0
+        for f in range(m):
+            phi_f = self.cartesian_embedding(f)
+            factor_dim = factors[f].dimension()
+            for i in range(factor_dim):
+                for j in range(i+1):
+                    f_ij = factors[f]._multiplication_table[i][j]
+                    e = phi_f(f_ij)
+                    self._multiplication_table[offset+i][offset+j] = e
+            offset += factor_dim
 
         self.rank.set_cache(sum(J.rank() for J in factors))
         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
@@ -3162,65 +3238,6 @@ class CartesianProductEJA(FiniteDimensionalEJA):
         return cartesian_product.symbol.join("%s" % factor
                                              for factor in self._sets)
 
-    def matrix_space(self):
-        r"""
-        Return the space that our matrix basis lives in as a Cartesian
-        product.
-
-        We don't simply use the ``cartesian_product()`` functor here
-        because it acts differently on SageMath MatrixSpaces and our
-        custom MatrixAlgebras, which are CombinatorialFreeModules. We
-        always want the result to be represented (and indexed) as
-        an ordered tuple.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
-            ....:                                  HadamardEJA,
-            ....:                                  OctonionHermitianEJA,
-            ....:                                  RealSymmetricEJA)
-
-        EXAMPLES::
-
-            sage: J1 = HadamardEJA(1)
-            sage: J2 = RealSymmetricEJA(2)
-            sage: J = cartesian_product([J1,J2])
-            sage: J.matrix_space()
-            The Cartesian product of (Full MatrixSpace of 1 by 1 dense
-            matrices over Algebraic Real Field, Full MatrixSpace of 2
-            by 2 dense matrices over Algebraic Real Field)
-
-        ::
-
-            sage: J1 = ComplexHermitianEJA(1)
-            sage: J2 = ComplexHermitianEJA(1)
-            sage: J = cartesian_product([J1,J2])
-            sage: J.one().to_matrix()[0]
-            +---+
-            | 1 |
-            +---+
-            sage: J.one().to_matrix()[1]
-            +---+
-            | 1 |
-            +---+
-
-        ::
-
-            sage: J1 = OctonionHermitianEJA(1)
-            sage: J2 = OctonionHermitianEJA(1)
-            sage: J = cartesian_product([J1,J2])
-            sage: J.one().to_matrix()[0]
-            +----+
-            | e0 |
-            +----+
-            sage: J.one().to_matrix()[1]
-            +----+
-            | e0 |
-            +----+
-
-        """
-        return super().matrix_space()
-
 
     @cached_method
     def cartesian_projection(self, i):
@@ -3422,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:
@@ -3445,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