+class RationalBasisCartesianProductEJA(CartesianProductEJA,
+ RationalBasisEJA):
+ r"""
+ A separate class for products of algebras for which we know a
+ rational basis.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
+ ....: HadamardEJA,
+ ....: JordanSpinEJA,
+ ....: RealSymmetricEJA)
+
+ EXAMPLES:
+
+ This gives us fast characteristic polynomial computations in
+ product algebras, too::
+
+
+ sage: J1 = JordanSpinEJA(2)
+ sage: J2 = RealSymmetricEJA(3)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.characteristic_polynomial_of().degree()
+ 5
+ sage: J.rank()
+ 5
+
+ TESTS:
+
+ 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()`` 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)
+
+ @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
+
+def random_eja(max_dimension=None, *args, **kwargs):
+ r"""
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import random_eja
+
+ TESTS::
+
+ sage: n = ZZ.random_element(1,5)
+ sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
+ sage: J.dimension() <= n
+ True
+
+ """
+ # Use the ConcreteEJA default as the total upper bound (regardless
+ # of any whether or not any individual factors set a lower limit).
+ if max_dimension is None:
+ max_dimension = ConcreteEJA._max_random_instance_dimension()
+ J1 = ConcreteEJA.random_instance(max_dimension, *args, **kwargs)
+
+
+ # Roll the dice to see if we attempt a Cartesian product.
+ dice_roll = ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1)
+ new_max_dimension = max_dimension - J1.dimension()
+ if new_max_dimension == 0 or dice_roll != 0:
+ # If it's already as big as we're willing to tolerate, just
+ # return it and don't worry about Cartesian products.
+ return J1
+ else:
+ # Use random_eja() again so we can get more than two factors
+ # if the sub-call also Decides on a cartesian product.
+ J2 = random_eja(new_max_dimension, *args, **kwargs)
+ return cartesian_product([J1,J2])
+
+
+class ComplexSkewHermitianEJA(RationalBasisEJA):
+ r"""
+ The EJA described in Faraut and Koranyi's Exercise III.1.b.
+ """
+ @staticmethod
+ def _denormalized_basis(A):
+ """
+ SETUP::
+
+ sage: from mjo.hurwitz import ComplexMatrixAlgebra
+ sage: from mjo.eja.eja_algebra import ComplexSkewHermitianEJA
+
+ TESTS::
+
+ sage: n = ZZ.random_element(1,2)
+ sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
+ sage: B = ComplexSkewHermitianEJA._denormalized_basis(A)
+ sage: all( M.is_skew_hermitian() for M in B)
+ True
+
+ """
+ es = A.entry_algebra_gens()
+ gen = lambda A,m: A.monomial(m)
+
+ basis = []
+
+ # The size of the blocks. We're going to treat these thing as
+ # 2x2 block matrices,
+ #
+ # [ x1 x2 ]
+ # [ -x2^* x1-conj ]
+ #
+ # where x1 is skew-Hermitian and x2 is symmetric.
+ #
+ m = A.nrows()/2
+
+ # We only loop through the top half of the matrix, because the
+ # bottom can be constructed from the top.
+ for i in range(m):
+
+ # First do the top-left block, which is skew-Hermitian.
+ # We can compute the bottom-right block in the process.
+ for j in range(i+1):
+ if i == j:
+ # Top-left block's entry.
+ E_ii = gen(A, (i,j,es[1]))
+
+ # Bottom-right block's entry.
+ E_ii += gen(A, (i+m,j+m,es[1])).conjugate()
+ basis.append(E_ii)
+ else:
+ for e in es:
+ # Top-left block's entry.
+ E_ij = gen(A, (i,j,e))
+ E_ij -= E_ij.conjugate_transpose()
+
+ # Bottom-right block's entry.
+ F_ij = gen(A, (i+m,j+m,e)).conjugate()
+ F_ij -= F_ij.conjugate_transpose()
+
+ basis.append(E_ij + F_ij)
+
+ # Now do the top-right block, which is symmetric, and compute
+ # the bottom-left block along the way.
+ for j in range(m,i+m+1):
+ if (i+m) == j:
+ # A symmetric (not Hermitian!) complex matrix can
+ # have both real and complex entries on its
+ # diagonal.
+ for e in es:
+ # Top-right block's entry.
+ E_ii = gen(A, (i,j,e))
+
+ # Bottom-left block's entry.
+ E_ii -= gen(A, (i-m,j-m,e)).conjugate()
+ basis.append(E_ii)
+ else:
+ for e in es:
+ # Top-right block's entry. BEWARE! We're not
+ # reflecting across the main diagonal as in
+ # (i,j)~(j,i). We're only reflecting across
+ # the diagonal for the top-right block.
+ E_ij = gen(A, (i,j,e))
+
+ # Shift it back to non-offset coords, transpose,
+ # and put it back:
+ #
+ # (i,j) -> (i,j-m) -> (j-m, i) -> (j-m, i+m)
+ E_ij += gen(A, (j-m,i+m,e))
+
+ # Bottom-left's block's below-diagonal entry.
+ # Just shift the top-right coords down m and
+ # left m.
+ F_ij = -gen(A, (i+m,j-m,e)).conjugate()
+ F_ij += -gen(A, (j,i,e)).conjugate()
+
+ basis.append(E_ij + F_ij)
+
+ return tuple( basis )
+
+
+ def __init__(self, n, field=AA, **kwargs):
+ # New code; always check the axioms.
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+ from mjo.hurwitz import ComplexMatrixAlgebra
+ A = ComplexMatrixAlgebra(2*n, scalars=field)
+
+ I_n = matrix.identity(ZZ, n)
+ J = matrix.block(ZZ, 2, 2, (0, I_n, -I_n, 0), subdivide=False)
+ J = A.from_list(J.rows())
+
+ def jordan_product(X,Y):
+ return (X*J*Y + Y*J*X)/2
+
+ def inner_product(X,Y):
+ return (X*Y.conjugate_transpose()).trace().real()
+
+ super().__init__(self._denormalized_basis(A),
+ jordan_product,
+ inner_product,
+ field=field,
+ matrix_space=A,
+ **kwargs)