]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: allow Cartesian products to be returned from random_eja().
[sage.d.git] / mjo / eja / eja_algebra.py
index 48421e357544ef2b5c5613bf9e222827e9a20fec..5bb4cee119e31cde1fd17accc7cc6e0bdee7f041 100644 (file)
@@ -112,6 +112,23 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         generally rules out using the rationals as your ``field``, but
         is required for spectral decompositions.
 
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import random_eja
+
+    TESTS:
+
+    We should compute that an element subalgebra is associative even
+    if we circumvent the element method::
+
+        sage: set_random_seed()
+        sage: J = random_eja(field=QQ,orthonormalize=False)
+        sage: x = J.random_element()
+        sage: A = x.subalgebra_generated_by(orthonormalize=False)
+        sage: basis = tuple(b.superalgebra_element() for b in A.basis())
+        sage: J.subalgebra(basis, orthonormalize=False).is_associative()
+        True
+
     """
     Element = FiniteDimensionalEJAElement
 
@@ -121,23 +138,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  inner_product,
                  field=AA,
                  orthonormalize=True,
-                 associative=False,
+                 associative=None,
                  cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
                  prefix='e'):
 
-        # Keep track of whether or not the matrix basis consists of
-        # tuples, since we need special cases for them damned near
-        # everywhere.  This is INDEPENDENT of whether or not the
-        # algebra is a cartesian product, since a subalgebra of a
-        # cartesian product will have a basis of tuples, but will not
-        # in general itself be a cartesian product algebra.
-        self._matrix_basis_is_cartesian = False
         n = len(basis)
-        if n > 0:
-            if hasattr(basis[0], 'cartesian_factors'):
-                self._matrix_basis_is_cartesian = True
 
         if check_field:
             if not field.is_subring(RR):
@@ -146,20 +153,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 # we've specified a real embedding.
                 raise ValueError("scalar field is not real")
 
+        from mjo.eja.eja_utils import _change_ring
         # If the basis given to us wasn't over the field that it's
         # supposed to be over, fix that. Or, you know, crash.
-        if not cartesian_product:
-            # The field for a cartesian product algebra comes from one
-            # of its factors and is the same for all factors, so
-            # there's no need to "reapply" it on product algebras.
-            if self._matrix_basis_is_cartesian:
-                # OK since if n == 0, the basis does not consist of tuples.
-                P = basis[0].parent()
-                basis = tuple( P(tuple(b_i.change_ring(field) for b_i in b))
-                               for b in basis )
-            else:
-                basis = tuple( b.change_ring(field) for b in basis )
-
+        basis = tuple( _change_ring(b, field) for b in basis )
 
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
@@ -178,12 +175,28 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 
         category = MagmaticAlgebras(field).FiniteDimensional()
-        category = category.WithBasis().Unital()
+        category = category.WithBasis().Unital().Commutative()
+
+        if associative is None:
+            # We should figure it out. As with check_axioms, we have to do
+            # this without the help of the _jordan_product_is_associative()
+            # method because we need to know the category before we
+            # initialize the algebra.
+            associative = all( jordan_product(jordan_product(bi,bj),bk)
+                               ==
+                               jordan_product(bi,jordan_product(bj,bk))
+                               for bi in basis
+                               for bj in basis
+                               for bk in basis)
+
         if associative:
             # Element subalgebras can take advantage of this.
             category = category.Associative()
         if cartesian_product:
-            category = category.CartesianProducts()
+            # 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.
@@ -331,8 +344,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: if n > 0:
             ....:     i = ZZ.random_element(n)
             ....:     j = ZZ.random_element(n)
-            ....:     ei = J.gens()[i]
-            ....:     ej = J.gens()[j]
+            ....:     ei = J.monomial(i)
+            ....:     ej = J.monomial(j)
             ....:     ei_ej = J.product_on_basis(i,j)
             sage: ei*ej == ei_ej
             True
@@ -422,6 +435,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         """
         return "Associative" in self.category().axioms()
 
+    def _is_commutative(self):
+        r"""
+        Whether or not this algebra's multiplication table is commutative.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid multiplication table.
+        """
+        return all( x*y == y*x for x in self.gens() for y in self.gens() )
+
     def _is_jordanian(self):
         r"""
         Whether or not this algebra's multiplication table respects the
@@ -429,16 +452,102 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         We only check one arrangement of `x` and `y`, so for a
         ``True`` result to be truly true, you should also check
-        :meth:`is_commutative`. This method should of course always
+        :meth:`_is_commutative`. This method should of course always
         return ``True``, unless this algebra was constructed with
         ``check_axioms=False`` and passed an invalid multiplication table.
         """
-        return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j])
+        return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
                     ==
-                    (self.gens()[i])*((self.gens()[i]**2)*self.gens()[j])
+                    (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
                     for i in range(self.dimension())
                     for j in range(self.dimension()) )
 
+    def _jordan_product_is_associative(self):
+        r"""
+        Return whether or not this algebra's Jordan product is
+        associative; that is, whether or not `x*(y*z) = (x*y)*z`
+        for all `x,y,x`.
+
+        This method should agree with :meth:`is_associative` unless
+        you lied about the value of the ``associative`` parameter
+        when you constructed the algebra.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA)
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(4, orthonormalize=False)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: A._jordan_product_is_associative()
+            True
+
+        ::
+
+            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: A._jordan_product_is_associative()
+            True
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: A._jordan_product_is_associative()
+            True
+
+        TESTS:
+
+        The values we've presupplied to the constructors agree with
+        the computation::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J.is_associative() == J._jordan_product_is_associative()
+            True
+
+        """
+        R = self.base_ring()
+
+        # Used to check whether or not something is zero.
+        epsilon = R.zero()
+        if not R.is_exact():
+            # I don't know of any examples that make this magnitude
+            # necessary because I don't know how to make an
+            # associative algebra when the element subalgebra
+            # construction is unreliable (as it is over RDF; we can't
+            # find the degree of an element because we can't compute
+            # the rank of a matrix). But even multiplication of floats
+            # is non-associative, so *some* epsilon is needed... let's
+            # just take the one from _inner_product_is_associative?
+            epsilon = 1e-15
+
+        for i in range(self.dimension()):
+            for j in range(self.dimension()):
+                for k in range(self.dimension()):
+                    x = self.monomial(i)
+                    y = self.monomial(j)
+                    z = self.monomial(k)
+                    diff = (x*y)*z - x*(y*z)
+
+                    if diff.norm() > epsilon:
+                        return False
+
+        return True
+
     def _inner_product_is_associative(self):
         r"""
         Return whether or not this algebra's inner product `B` is
@@ -460,9 +569,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         for i in range(self.dimension()):
             for j in range(self.dimension()):
                 for k in range(self.dimension()):
-                    x = self.gens()[i]
-                    y = self.gens()[j]
-                    z = self.gens()[k]
+                    x = self.monomial(i)
+                    y = self.monomial(j)
+                    z = self.monomial(k)
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 
                     if diff.abs() > epsilon:
@@ -480,7 +589,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  JordanSpinEJA,
             ....:                                  HadamardEJA,
             ....:                                  RealSymmetricEJA)
 
@@ -509,22 +619,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J2 = RealSymmetricEJA(2)
             sage: J = cartesian_product([J1,J2])
             sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
-            e(0, 1) + e(1, 2)
+            e1 + e5
 
         TESTS:
 
-        Ensure that we can convert any element of the two non-matrix
-        simple algebras (whose matrix representations are columns)
-        back and forth faithfully::
+        Ensure that we can convert any element back and forth
+        faithfully between its matrix and algebra representations::
 
             sage: set_random_seed()
-            sage: J = HadamardEJA.random_instance()
-            sage: x = J.random_element()
-            sage: J(x.to_vector().column()) == x
-            True
-            sage: J = JordanSpinEJA.random_instance()
+            sage: J = random_eja()
             sage: x = J.random_element()
-            sage: J(x.to_vector().column()) == x
+            sage: J(x.to_matrix()) == x
             True
 
         We cannot coerce elements between algebras just because their
@@ -540,7 +645,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             Traceback (most recent call last):
             ...
             ValueError: not an element of this algebra
-
         """
         msg = "not an element of this algebra"
         if elt in self.base_ring():
@@ -808,7 +912,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # And to each subsequent row, prepend an entry that belongs to
         # the left-side "header column."
-        M += [ [self.gens()[i]] + [ self.product_on_basis(i,j)
+        M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j)
                                     for j in range(n) ]
                for i in range(n) ]
 
@@ -1310,7 +1414,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         def L_x_i_j(i,j):
             # From a result in my book, these are the entries of the
             # basis representation of L_x.
-            return sum( vars[k]*self.gens()[k].operator().matrix()[i,j]
+            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
                         for k in range(n) )
 
         L_x = matrix(F, n, n, L_x_i_j)
@@ -1479,6 +1583,13 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
                 raise TypeError("basis not rational")
 
+        super().__init__(basis,
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         check_field=check_field,
+                         **kwargs)
+
         self._rational_algebra = None
         if field is not QQ:
             # There's no point in constructing the extra algebra if this
@@ -1492,17 +1603,11 @@ class RationalBasisEJA(FiniteDimensionalEJA):
                                        jordan_product,
                                        inner_product,
                                        field=QQ,
+                                       associative=self.is_associative(),
                                        orthonormalize=False,
                                        check_field=False,
                                        check_axioms=False)
 
-        super().__init__(basis,
-                         jordan_product,
-                         inner_product,
-                         field=field,
-                         check_field=check_field,
-                         **kwargs)
-
     @cached_method
     def _charpoly_coefficients(self):
         r"""
@@ -1866,10 +1971,15 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n),
-                                               self.jordan_product,
-                                               self.trace_inner_product,
-                                               **kwargs)
+        associative = False
+        if n <= 1:
+            associative = True
+
+        super().__init__(self._denormalized_basis(n),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         associative=associative,
+                         **kwargs)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
@@ -1959,7 +2069,7 @@ class ComplexMatrixEJA(MatrixEJA):
             True
 
         """
-        super(ComplexMatrixEJA,cls).real_embed(M)
+        super().real_embed(M)
         n = M.nrows()
 
         # We don't need any adjoined elements...
@@ -2006,7 +2116,7 @@ class ComplexMatrixEJA(MatrixEJA):
             True
 
         """
-        super(ComplexMatrixEJA,cls).real_unembed(M)
+        super().real_unembed(M)
         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
         F = cls.complex_extension(M.base_ring())
@@ -2154,10 +2264,15 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n),
-                                                  self.jordan_product,
-                                                  self.trace_inner_product,
-                                                  **kwargs)
+        associative = False
+        if n <= 1:
+            associative = True
+
+        super().__init__(self._denormalized_basis(n),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         associative=associative,
+                         **kwargs)
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
@@ -2240,7 +2355,7 @@ class QuaternionMatrixEJA(MatrixEJA):
             True
 
         """
-        super(QuaternionMatrixEJA,cls).real_embed(M)
+        super().real_embed(M)
         quaternions = M.base_ring()
         n = M.nrows()
 
@@ -2295,7 +2410,7 @@ class QuaternionMatrixEJA(MatrixEJA):
             True
 
         """
-        super(QuaternionMatrixEJA,cls).real_unembed(M)
+        super().real_unembed(M)
         n = ZZ(M.nrows())
         d = cls.dimension_over_reals()
 
@@ -2460,10 +2575,16 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n),
-                                                     self.jordan_product,
-                                                     self.trace_inner_product,
-                                                     **kwargs)
+        associative = False
+        if n <= 1:
+            associative = True
+
+        super().__init__(self._denormalized_basis(n),
+                         self.jordan_product,
+                         self.trace_inner_product,
+                         associative=associative,
+                         **kwargs)
+
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
@@ -2689,10 +2810,17 @@ class BilinearFormEJA(ConcreteEJA):
 
         n = B.nrows()
         column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
-        super(BilinearFormEJA, self).__init__(column_basis,
-                                              jordan_product,
-                                              inner_product,
-                                              **kwargs)
+
+        # TODO: I haven't actually checked this, but it seems legit.
+        associative = False
+        if n <= 2:
+            associative = True
+
+        super().__init__(column_basis,
+                         jordan_product,
+                         inner_product,
+                         associative=associative,
+                         **kwargs)
 
         # The rank of this algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded
@@ -2797,7 +2925,7 @@ class JordanSpinEJA(BilinearFormEJA):
 
         # But also don't pass check_field=False here, because the user
         # can pass in a field!
-        super(JordanSpinEJA, self).__init__(B, **kwargs)
+        super().__init__(B, **kwargs)
 
     @staticmethod
     def _max_random_instance_size():
@@ -2855,10 +2983,12 @@ class TrivialEJA(ConcreteEJA):
         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
-        super(TrivialEJA, self).__init__(basis,
-                                         jordan_product,
-                                         inner_product,
-                                         **kwargs)
+        super().__init__(basis,
+                         jordan_product,
+                         inner_product,
+                         associative=True,
+                         **kwargs)
+
         # The rank is zero using my definition, namely the dimension of the
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)
@@ -2871,8 +3001,7 @@ class TrivialEJA(ConcreteEJA):
         return cls(**kwargs)
 
 
-class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
-                          FiniteDimensionalEJA):
+class CartesianProductEJA(FiniteDimensionalEJA):
     r"""
     The external (orthogonal) direct sum of two or more Euclidean
     Jordan algebras. Every Euclidean Jordan algebra decomposes into an
@@ -2968,6 +3097,33 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         sage: CP2.is_associative()
         False
 
+    Cartesian products of Cartesian products work::
+
+        sage: J1 = JordanSpinEJA(1)
+        sage: J2 = JordanSpinEJA(1)
+        sage: J3 = JordanSpinEJA(1)
+        sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
+        sage: J.multiplication_table()
+        +----++----+----+----+
+        | *  || e0 | e1 | e2 |
+        +====++====+====+====+
+        | e0 || e0 | 0  | 0  |
+        +----++----+----+----+
+        | e1 || 0  | e1 | 0  |
+        +----++----+----+----+
+        | e2 || 0  | 0  | e2 |
+        +----++----+----+----+
+        sage: HadamardEJA(3).multiplication_table()
+        +----++----+----+----+
+        | *  || e0 | e1 | e2 |
+        +====++====+====+====+
+        | e0 || e0 | 0  | 0  |
+        +----++----+----+----+
+        | e1 || 0  | e1 | 0  |
+        +----++----+----+----+
+        | e2 || 0  | 0  | e2 |
+        +----++----+----+----+
+
     TESTS:
 
     All factors must share the same base field::
@@ -2995,37 +3151,41 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
     Element = FiniteDimensionalEJAElement
 
 
-    def __init__(self, algebras, **kwargs):
-        CombinatorialFreeModule_CartesianProduct.__init__(self,
-                                                          algebras,
-                                                          **kwargs)
-        field = algebras[0].base_ring()
-        if not all( J.base_ring() == field for J in algebras ):
+    def __init__(self, factors, **kwargs):
+        m = len(factors)
+        if m == 0:
+            return TrivialEJA()
+
+        self._sets = factors
+
+        field = factors[0].base_ring()
+        if not all( J.base_ring() == field for J in factors ):
             raise ValueError("all factors must share the same base field")
 
-        associative = all( m.is_associative() for m in algebras )
+        associative = all( f.is_associative() for f in factors )
 
-        # The definition of matrix_space() and self.basis() relies
-        # only on the stuff in the CFM_CartesianProduct class, which
-        # we've already initialized.
-        Js = self.cartesian_factors()
-        m = len(Js)
         MS = self.matrix_space()
-        basis = tuple(
-            MS(tuple( self.cartesian_projection(i)(b).to_matrix()
-                      for i in range(m) ))
-            for b in self.basis()
-        )
+        basis = []
+        zero = MS.zero()
+        for i in range(m):
+            for b in factors[i].matrix_basis():
+                z = list(zero)
+                z[i] = b
+                basis.append(z)
+
+        basis = tuple( MS(b) for b in basis )
 
         # Define jordan/inner products that operate on that matrix_basis.
         def jordan_product(x,y):
             return MS(tuple(
-                (Js[i](x[i])*Js[i](y[i])).to_matrix() for i in range(m)
+                (factors[i](x[i])*factors[i](y[i])).to_matrix()
+                for i in range(m)
             ))
 
         def inner_product(x, y):
             return sum(
-                Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m)
+                factors[i](x[i]).inner_product(factors[i](y[i]))
+                for i in range(m)
             )
 
         # There's no need to check the field since it already came
@@ -3045,9 +3205,25 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
                                       check_field=False,
                                       check_axioms=False)
 
-        ones = tuple(J.one() for J in algebras)
-        self.one.set_cache(self._cartesian_product_of_elements(ones))
-        self.rank.set_cache(sum(J.rank() for J in algebras))
+        ones = tuple(J.one().to_matrix() for J in factors)
+        self.one.set_cache(self(ones))
+        self.rank.set_cache(sum(J.rank() for J in factors))
+
+    def cartesian_factors(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        return self._sets
+
+    def cartesian_factor(self, i):
+        r"""
+        Return the ``i``th factor of this algebra.
+        """
+        return self._sets[i]
+
+    def _repr_(self):
+        # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
+        from sage.categories.cartesian_product import cartesian_product
+        return cartesian_product.symbol.join("%s" % factor
+                                             for factor in self._sets)
 
     def matrix_space(self):
         r"""
@@ -3144,9 +3320,12 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             True
 
         """
-        Ji = self.cartesian_factors()[i]
-        # Requires the fix on Trac 31421/31422 to work!
-        Pi = super().cartesian_projection(i)
+        offset = sum( self.cartesian_factor(k).dimension()
+                      for k in range(i) )
+        Ji = self.cartesian_factor(i)
+        Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
+                                   codomain=Ji)
+
         return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
 
     @cached_method
@@ -3252,9 +3431,11 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
             True
 
         """
-        Ji = self.cartesian_factors()[i]
-        # Requires the fix on Trac 31421/31422 to work!
-        Ei = super().cartesian_embedding(i)
+        offset = sum( self.cartesian_factor(k).dimension()
+                      for k in range(i) )
+        Ji = self.cartesian_factor(i)
+        Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
+                                 codomain=self)
         return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
 
 
@@ -3299,4 +3480,15 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
 
-random_eja = ConcreteEJA.random_instance
+def random_eja(*args, **kwargs):
+    J1 = ConcreteEJA.random_instance(*args, **kwargs)
+
+    # This might make Cartesian products appear roughly as often as
+    # any other ConcreteEJA.
+    if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
+        # Use random_eja() again so we can get more than two factors.
+        J2 = random_eja(*args, **kwargs)
+        J = cartesian_product([J1,J2])
+        return J
+    else:
+        return J1