from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
from sage.categories.magmatic_algebras import MagmaticAlgebras
-from sage.combinat.free_module import CombinatorialFreeModule
+from sage.categories.sets_cat import cartesian_product
+from sage.combinat.free_module import (CombinatorialFreeModule,
+ CombinatorialFreeModule_CartesianProduct)
from sage.matrix.constructor import matrix
from sage.matrix.matrix_space import MatrixSpace
from sage.misc.cachefunc import cached_method
QuadraticField)
from mjo.eja.eja_element import FiniteDimensionalEJAElement
from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
-from mjo.eja.eja_utils import _mat2vec
+from mjo.eja.eja_utils import _all2list, _mat2vec
class FiniteDimensionalEJA(CombinatorialFreeModule):
r"""
INPUT:
- - basis -- a tuple of basis elements in their matrix form.
+ - basis -- a tuple of basis elements in "matrix form," which
+ must be the same form as the arguments to ``jordan_product``
+ and ``inner_product``. In reality, "matrix form" can be either
+ vectors, matrices, or a Cartesian product (ordered tuple)
+ of vectors or matrices. All of these would ideally be vector
+ spaces in sage with no special-casing needed; but in reality
+ we turn vectors into column-matrices and Cartesian products
+ `(a,b)` into column matrices `(a,b)^{T}` after converting
+ `a` and `b` themselves.
- jordan_product -- function of two elements (in matrix form)
that returns their jordan product in this algebra; this will
field=AA,
orthonormalize=True,
associative=False,
+ 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):
# Note: this does return true for the real algebraic
# If the basis given to us wasn't over the field that it's
# supposed to be over, fix that. Or, you know, crash.
- basis = tuple( b.change_ring(field) for b in basis )
+ 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 )
+
if check_axioms:
# Check commutativity of the Jordan and inner-products.
if associative:
# Element subalgebras can take advantage of this.
category = category.Associative()
+ if cartesian_product:
+ category = category.CartesianProducts()
# Call the superclass constructor so that we can use its from_vector()
# method to build our multiplication table.
- n = len(basis)
- super().__init__(field,
- range(n),
- prefix=prefix,
- category=category,
- bracket=False)
+ CombinatorialFreeModule.__init__(self,
+ field,
+ range(n),
+ prefix=prefix,
+ category=category,
+ bracket=False)
# Now comes all of the hard work. We'll be constructing an
# ambient vector space V that our (vectorized) basis lives in,
# we see in things like x = 1*e1 + 2*e2.
vector_basis = basis
- from sage.structure.element import is_Matrix
- basis_is_matrices = False
-
degree = 0
if n > 0:
- if is_Matrix(basis[0]):
- if basis[0].is_square():
- # TODO: this ugly is_square() hack works around the problem
- # of passing to_matrix()ed vectors in as the basis from a
- # subalgebra. They aren't REALLY matrices, at least not of
- # the type that we assume here... Ugh.
- basis_is_matrices = True
- from mjo.eja.eja_utils import _vec2mat
- vector_basis = tuple( map(_mat2vec,basis) )
- degree = basis[0].nrows()**2
- else:
- # convert from column matrices to vectors, yuck
- basis = tuple( map(_mat2vec,basis) )
- vector_basis = basis
- degree = basis[0].degree()
- else:
- degree = basis[0].degree()
+ degree = len(_all2list(basis[0]))
- # Build an ambient space that fits...
+ # Build an ambient space that fits our matrix basis when
+ # written out as "long vectors."
V = VectorSpace(field, degree)
- # We overwrite the name "vector_basis" in a second, but never modify it
- # in place, to this effectively makes a copy of it.
- deortho_vector_basis = vector_basis
+ # The matrix that will hole the orthonormal -> unorthonormal
+ # coordinate transformation.
self._deortho_matrix = None
if orthonormalize:
- from mjo.eja.eja_utils import gram_schmidt
- if basis_is_matrices:
- vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y))
- vector_basis = gram_schmidt(vector_basis, vector_ip)
- else:
- vector_basis = gram_schmidt(vector_basis, inner_product)
+ # Save a copy of the un-orthonormalized basis for later.
+ # Convert it to ambient V (vector) coordinates while we're
+ # at it, because we'd have to do it later anyway.
+ deortho_vector_basis = tuple( V(_all2list(b)) for b in basis )
- # Normalize the "matrix" basis, too!
- basis = vector_basis
-
- if basis_is_matrices:
- basis = tuple( map(_vec2mat,basis) )
+ from mjo.eja.eja_utils import gram_schmidt
+ basis = tuple(gram_schmidt(basis, inner_product))
- # Save the matrix "basis" for later... this is the last time we'll
- # reference it in this constructor.
- if basis_is_matrices:
- self._matrix_basis = basis
- else:
- MS = MatrixSpace(self.base_ring(), degree, 1)
- self._matrix_basis = tuple( MS(b) for b in basis )
+ # Save the (possibly orthonormalized) matrix basis for
+ # later...
+ self._matrix_basis = basis
- # Now create the vector space for the algebra...
+ # Now create the vector space for the algebra, which will have
+ # its own set of non-ambient coordinates (in terms of the
+ # supplied basis).
+ vector_basis = tuple( V(_all2list(b)) for b in basis )
W = V.span_of_basis( vector_basis, check=check_axioms)
if orthonormalize:
# Now we actually compute the multiplication and inner-product
# tables/matrices using the possibly-orthonormalized basis.
- self._inner_product_matrix = matrix.zero(field, n)
- self._multiplication_table = [ [0 for j in range(i+1)] for i in range(n) ]
+ self._inner_product_matrix = matrix.identity(field, n)
+ self._multiplication_table = [ [0 for j in range(i+1)]
+ for i in range(n) ]
- print("vector_basis:")
- print(vector_basis)
# Note: the Jordan and inner-products are defined in terms
# of the ambient basis. It's important that their arguments
# are in ambient coordinates as well.
for i in range(n):
for j in range(i+1):
# ortho basis w.r.t. ambient coords
- q_i = vector_basis[i]
- q_j = vector_basis[j]
-
- if basis_is_matrices:
- q_i = _vec2mat(q_i)
- q_j = _vec2mat(q_j)
+ q_i = basis[i]
+ q_j = basis[j]
+ # The jordan product returns a matrixy answer, so we
+ # have to convert it to the algebra coordinates.
elt = jordan_product(q_i, q_j)
- ip = inner_product(q_i, q_j)
-
- if basis_is_matrices:
- # do another mat2vec because the multiplication
- # table is in terms of vectors
- elt = _mat2vec(elt)
-
- # TODO: the jordan product turns things back into
- # matrices here even if they're supposed to be
- # vectors. ugh. Can we get rid of vectors all together
- # please?
- elt = W.coordinate_vector(elt)
+ elt = W.coordinate_vector(V(_all2list(elt)))
self._multiplication_table[i][j] = self.from_vector(elt)
- self._inner_product_matrix[i,j] = ip
- self._inner_product_matrix[j,i] = ip
+
+ if not orthonormalize:
+ # If we're orthonormalizing the basis with respect
+ # to an inner-product, then the inner-product
+ # matrix with respect to the resulting basis is
+ # just going to be the identity.
+ ip = inner_product(q_i, q_j)
+ self._inner_product_matrix[i,j] = ip
+ self._inner_product_matrix[j,i] = ip
self._inner_product_matrix._cache = {'hermitian': True}
self._inner_product_matrix.set_immutable()
sage: y = J.random_element()
sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
True
+
"""
B = self._inner_product_matrix
return (B*x.to_vector()).inner_product(y.to_vector())
- def _is_commutative(self):
+ def is_associative(self):
r"""
- Whether or not this algebra's multiplication table is commutative.
+ Return whether or not this algebra's Jordan product is associative.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+
+ EXAMPLES::
+
+ sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
+ sage: J.is_associative()
+ False
+ sage: x = sum(J.gens())
+ sage: A = x.subalgebra_generated_by(orthonormalize=False)
+ sage: A.is_associative()
+ True
- 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.product_on_basis(i,j) == self.product_on_basis(i,j)
- for i in range(self.dimension())
- for j in range(self.dimension()) )
+ return "Associative" in self.category().axioms()
def _is_jordanian(self):
r"""
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.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
+ return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j])
==
- (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
+ (self.gens()[i])*((self.gens()[i]**2)*self.gens()[j])
for i in range(self.dimension())
for j in range(self.dimension()) )
This method should of course always return ``True``, unless
this algebra was constructed with ``check_axioms=False`` and
- passed an invalid multiplication table.
+ passed an invalid Jordan or inner-product.
"""
# Used to check whether or not something is zero in an inexact
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)
+ x = self.gens()[i]
+ y = self.gens()[j]
+ z = self.gens()[k]
diff = (x*y).inner_product(z) - x.inner_product(y*z)
if self.base_ring().is_exact():
...
ValueError: not an element of this algebra
+ Tuples work as well, provided that the matrix basis for the
+ algebra consists of them::
+
+ sage: J1 = HadamardEJA(3)
+ 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)
+
TESTS:
Ensure that we can convert any element of the two non-matrix
# closure whereas the base ring of the 3-by-3 identity matrix
# could be QQ instead of QQbar.
#
+ # And, we also have to handle Cartesian product bases (when
+ # the matric basis consists of tuples) here. The "good news"
+ # is that we're already converting everything to long vectors,
+ # and that strategy works for tuples as well.
+ #
# We pass check=False because the matrix basis is "guaranteed"
# to be linearly independent... right? Ha ha.
- V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
- W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()),
+ elt = _all2list(elt)
+ V = VectorSpace(self.base_ring(), len(elt))
+ W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()),
check=False)
try:
- coords = W.coordinate_vector(_mat2vec(elt))
+ coords = W.coordinate_vector(V(elt))
except ArithmeticError: # vector is not in free module
raise ValueError(msg)
# And to each subsequent row, prepend an entry that belongs to
# the left-side "header column."
- M += [ [self.monomial(i)] + [ self.product_on_basis(i,j)
- for j in range(n) ]
+ M += [ [self.gens()[i]] + [ self.product_on_basis(i,j)
+ for j in range(n) ]
for i in range(n) ]
return table(M, header_row=True, header_column=True, frame=True)
Why implement this for non-matrix algebras? Avoiding special
cases for the :class:`BilinearFormEJA` pays with simplicity in
its own right. But mainly, we would like to be able to assume
- that elements of a :class:`DirectSumEJA` can be displayed
+ that elements of a :class:`CartesianProductEJA` can be displayed
nicely, without having to have special classes for direct sums
one of whose components was a matrix algebra.
if self.is_trivial():
return MatrixSpace(self.base_ring(), 0)
else:
- return self._matrix_basis[0].matrix_space()
+ return self.matrix_basis()[0].parent()
@cached_method
sage: from mjo.eja.eja_algebra import (HadamardEJA,
....: random_eja)
- EXAMPLES::
+ EXAMPLES:
+
+ We can compute unit element in the Hadamard EJA::
sage: J = HadamardEJA(5)
sage: J.one()
e0 + e1 + e2 + e3 + e4
+ The unit element in the Hadamard EJA is inherited in the
+ subalgebras generated by its elements::
+
+ sage: J = HadamardEJA(5)
+ sage: J.one()
+ e0 + e1 + e2 + e3 + e4
+ sage: x = sum(J.gens())
+ sage: A = x.subalgebra_generated_by(orthonormalize=False)
+ sage: A.one()
+ f0
+ sage: A.one().superalgebra_element()
+ e0 + e1 + e2 + e3 + e4
+
TESTS:
- The identity element acts like the identity::
+ The identity element acts like the identity, regardless of
+ whether or not we orthonormalize::
sage: set_random_seed()
sage: J = random_eja()
sage: x = J.random_element()
sage: J.one()*x == x and x*J.one() == x
True
+ sage: A = x.subalgebra_generated_by()
+ sage: y = A.random_element()
+ sage: A.one()*y == y and y*A.one() == y
+ True
+
+ ::
+
+ sage: set_random_seed()
+ sage: J = random_eja(field=QQ, orthonormalize=False)
+ sage: x = J.random_element()
+ sage: J.one()*x == x and x*J.one() == x
+ True
+ sage: A = x.subalgebra_generated_by(orthonormalize=False)
+ sage: y = A.random_element()
+ sage: A.one()*y == y and y*A.one() == y
+ True
- The matrix of the unit element's operator is the identity::
+ The matrix of the unit element's operator is the identity,
+ regardless of the base field and whether or not we
+ orthonormalize::
sage: set_random_seed()
sage: J = random_eja()
sage: expected = matrix.identity(J.base_ring(), J.dimension())
sage: actual == expected
True
+ sage: x = J.random_element()
+ sage: A = x.subalgebra_generated_by()
+ sage: actual = A.one().operator().matrix()
+ sage: expected = matrix.identity(A.base_ring(), A.dimension())
+ sage: actual == expected
+ True
+
+ ::
+
+ sage: set_random_seed()
+ sage: J = random_eja(field=QQ, orthonormalize=False)
+ sage: actual = J.one().operator().matrix()
+ sage: expected = matrix.identity(J.base_ring(), J.dimension())
+ sage: actual == expected
+ True
+ sage: x = J.random_element()
+ sage: A = x.subalgebra_generated_by(orthonormalize=False)
+ sage: actual = A.one().operator().matrix()
+ sage: expected = matrix.identity(A.base_ring(), A.dimension())
+ sage: actual == expected
+ True
Ensure that the cached unit element (often precomputed by
hand) agrees with the computed one::
sage: J.one() == cached
True
+ ::
+
+ sage: set_random_seed()
+ sage: J = random_eja(field=QQ, orthonormalize=False)
+ sage: cached = J.one()
+ sage: J.one.clear_cache()
+ sage: J.one() == cached
+ True
+
"""
# We can brute-force compute the matrices of the operators
# that correspond to the basis elements of this algebra.
if not c.is_idempotent():
raise ValueError("element is not idempotent: %s" % c)
- from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
-
# Default these to what they should be if they turn out to be
# trivial, because eigenspaces_left() won't return eigenvalues
# 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 = FiniteDimensionalEJASubalgebra(self, ())
+ trivial = self.subalgebra(())
J0 = trivial # eigenvalue zero
J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
J1 = trivial # eigenvalue one
J5 = eigspace
else:
gens = tuple( self.from_vector(b) for b in eigspace.basis() )
- subalg = FiniteDimensionalEJASubalgebra(self,
- gens,
- check_axioms=False)
+ subalg = self.subalgebra(gens, check_axioms=False)
if eigval == 0:
J0 = subalg
elif eigval == 1:
r"""
The `r` polynomial coefficients of the "characteristic polynomial
of" function.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import random_eja
+
+ TESTS:
+
+ The theory shows that these are all homogeneous polynomials of
+ a known degree::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
+ True
+
"""
n = self.dimension()
R = self.coordinate_polynomial_ring()
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.monomial(k).operator().matrix()[i,j]
+ return sum( vars[k]*self.gens()[k].operator().matrix()[i,j]
for k in range(n) )
L_x = matrix(F, n, n, L_x_i_j)
# The theory says that only the first "r" coefficients are
# nonzero, and they actually live in the original polynomial
- # ring and not the fraction field. We negate them because
- # in the actual characteristic polynomial, they get moved
- # to the other side where x^r lives.
- return -A_rref.solve_right(E*b).change_ring(R)[:r]
+ # ring and not the fraction field. We negate them because in
+ # the actual characteristic polynomial, they get moved to the
+ # other side where x^r lives. We don't bother to trim A_rref
+ # down to a square matrix and solve the resulting system,
+ # because the upper-left r-by-r portion of A_rref is
+ # guaranteed to be the identity matrix, so e.g.
+ #
+ # A_rref.solve_right(Y)
+ #
+ # would just be returning Y.
+ return (-E*b)[:r].change_ring(R)
@cached_method
def rank(self):
sage: set_random_seed() # long time
sage: J = random_eja() # long time
- sage: caches = J.rank() # long time
+ sage: cached = J.rank() # long time
sage: J.rank.clear_cache() # long time
sage: J.rank() == cached # long time
True
return len(self._charpoly_coefficients())
+ def subalgebra(self, basis, **kwargs):
+ r"""
+ Create a subalgebra of this algebra from the given basis.
+ """
+ from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+ return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+
+
def vector_space(self):
"""
Return the vector space that underlies this algebra.
return self.zero().to_vector().parent().ambient_vector_space()
- Element = FiniteDimensionalEJAElement
class RationalBasisEJA(FiniteDimensionalEJA):
r"""
jordan_product,
inner_product,
field=AA,
- orthonormalize=True,
check_field=True,
- check_axioms=True,
**kwargs):
if check_field:
if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
raise TypeError("basis not rational")
+ self._rational_algebra = None
if field is not QQ:
# There's no point in constructing the extra algebra if this
# one is already rational.
field=QQ,
orthonormalize=False,
check_field=False,
- check_axioms=False,
- **kwargs)
+ check_axioms=False)
super().__init__(basis,
jordan_product,
inner_product,
field=field,
check_field=check_field,
- check_axioms=check_axioms,
**kwargs)
@cached_method
a = ( a_i.change_ring(self.base_ring())
for a_i in self._rational_algebra._charpoly_coefficients() )
- # Now convert the coordinate variables back to the
+ if self._deortho_matrix is None:
+ # This can happen if our base ring was, say, AA and we
+ # chose not to (or didn't need to) orthonormalize. It's
+ # still faster to do the computations over QQ even if
+ # the numbers in the boxes stay the same.
+ return tuple(a)
+
+ # Otherwise, convert the coordinate variables back to the
# deorthonormalized ones.
R = self.coordinate_polynomial_ring()
from sage.modules.free_module_element import vector
class ComplexMatrixEJA(MatrixEJA):
+ # A manual dictionary-cache for the complex_extension() method,
+ # since apparently @classmethods can't also be @cached_methods.
+ _complex_extension = {}
+
+ @classmethod
+ def complex_extension(cls,field):
+ r"""
+ The complex field that we embed/unembed, as an extension
+ of the given ``field``.
+ """
+ if field in cls._complex_extension:
+ return cls._complex_extension[field]
+
+ # Sage doesn't know how to adjoin the complex "i" (the root of
+ # x^2 + 1) to a field in a general way. Here, we just enumerate
+ # all of the cases that I have cared to support so far.
+ if field is AA:
+ # Sage doesn't know how to embed AA into QQbar, i.e. how
+ # to adjoin sqrt(-1) to AA.
+ F = QQbar
+ elif not field.is_exact():
+ # RDF or RR
+ F = field.complex_field()
+ else:
+ # Works for QQ and... maybe some other fields.
+ R = PolynomialRing(field, 'z')
+ z = R.gen()
+ F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+
+ cls._complex_extension[field] = F
+ return F
+
@staticmethod
def dimension_over_reals():
return 2
blocks = []
for z in M.list():
- a = z.list()[0] # real part, I guess
- b = z.list()[1] # imag part, I guess
- blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
+ a = z.real()
+ b = z.imag()
+ blocks.append(matrix(field, 2, [ [ a, b],
+ [-b, a] ]))
return matrix.block(field, n, blocks)
super(ComplexMatrixEJA,cls).real_unembed(M)
n = ZZ(M.nrows())
d = cls.dimension_over_reals()
-
- # If "M" was normalized, its base ring might have roots
- # adjoined and they can stick around after unembedding.
- field = M.base_ring()
- R = PolynomialRing(field, 'z')
- z = R.gen()
-
- # Sage doesn't know how to adjoin the complex "i" (the root of
- # x^2 + 1) to a field in a general way. Here, we just enumerate
- # all of the cases that I have cared to support so far.
- if field is AA:
- # Sage doesn't know how to embed AA into QQbar, i.e. how
- # to adjoin sqrt(-1) to AA.
- F = QQbar
- elif not field.is_exact():
- # RDF or RR
- F = field.complex_field()
- else:
- # Works for QQ and... maybe some other fields.
- F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+ F = cls.complex_extension(M.base_ring())
i = F.gen()
# Go top-left to bottom-right (reading order), converting every
sage: set_random_seed()
sage: n = ZZ.random_element(1,5)
- sage: field = QuadraticField(2, 'sqrt2')
sage: B = ComplexHermitianEJA._denormalized_basis(n)
sage: all( M.is_symmetric() for M in B)
True
# * The diagonal will (as a result) be real.
#
S = []
+ Eij = matrix.zero(F,n)
for i in range(n):
for j in range(i+1):
- Eij = matrix(F, n, lambda k,l: k==i and l==j)
+ # "build" E_ij
+ Eij[i,j] = 1
if i == j:
Sij = cls.real_embed(Eij)
S.append(Sij)
else:
# The second one has a minus because it's conjugated.
- Sij_real = cls.real_embed(Eij + Eij.transpose())
+ Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
+ Sij_real = cls.real_embed(Eij)
S.append(Sij_real)
- Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
+ # Eij = I*Eij - I*Eij.transpose()
+ Eij[i,j] = I
+ Eij[j,i] = -I
+ Sij_imag = cls.real_embed(Eij)
S.append(Sij_imag)
+ Eij[j,i] = 0
+ # "erase" E_ij
+ Eij[i,j] = 0
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the complex extension "F".
return cls(n, **kwargs)
class QuaternionMatrixEJA(MatrixEJA):
+
+ # A manual dictionary-cache for the quaternion_extension() method,
+ # since apparently @classmethods can't also be @cached_methods.
+ _quaternion_extension = {}
+
+ @classmethod
+ def quaternion_extension(cls,field):
+ r"""
+ The quaternion field that we embed/unembed, as an extension
+ of the given ``field``.
+ """
+ if field in cls._quaternion_extension:
+ return cls._quaternion_extension[field]
+
+ Q = QuaternionAlgebra(field,-1,-1)
+
+ cls._quaternion_extension[field] = Q
+ return Q
+
@staticmethod
def dimension_over_reals():
return 4
# Use the base ring of the matrix to ensure that its entries can be
# multiplied by elements of the quaternion algebra.
- field = M.base_ring()
- Q = QuaternionAlgebra(field,-1,-1)
+ Q = cls.quaternion_extension(M.base_ring())
i,j,k = Q.gens()
# Go top-left to bottom-right (reading order), converting every
# * The diagonal will (as a result) be real.
#
S = []
+ Eij = matrix.zero(Q,n)
for i in range(n):
for j in range(i+1):
- Eij = matrix(Q, n, lambda k,l: k==i and l==j)
+ # "build" E_ij
+ Eij[i,j] = 1
if i == j:
Sij = cls.real_embed(Eij)
S.append(Sij)
else:
# The second, third, and fourth ones have a minus
# because they're conjugated.
- Sij_real = cls.real_embed(Eij + Eij.transpose())
+ # Eij = Eij + Eij.transpose()
+ Eij[j,i] = 1
+ Sij_real = cls.real_embed(Eij)
S.append(Sij_real)
- Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
+ # Eij = I*(Eij - Eij.transpose())
+ Eij[i,j] = I
+ Eij[j,i] = -I
+ Sij_I = cls.real_embed(Eij)
S.append(Sij_I)
- Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
+ # Eij = J*(Eij - Eij.transpose())
+ Eij[i,j] = J
+ Eij[j,i] = -J
+ Sij_J = cls.real_embed(Eij)
S.append(Sij_J)
- Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
+ # Eij = K*(Eij - Eij.transpose())
+ Eij[i,j] = K
+ Eij[j,i] = -K
+ Sij_K = cls.real_embed(Eij)
S.append(Sij_K)
+ Eij[j,i] = 0
+ # "erase" E_ij
+ Eij[i,j] = 0
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the quaternion algebra "Q".
"""
def __init__(self, n, **kwargs):
- def jordan_product(x,y):
- P = x.parent()
- return P(tuple( xi*yi for (xi,yi) in zip(x,y) ))
- def inner_product(x,y):
- return x.inner_product(y)
+ if n == 0:
+ jordan_product = lambda x,y: x
+ inner_product = lambda x,y: x
+ else:
+ def jordan_product(x,y):
+ P = x.parent()
+ return P( xi*yi for (xi,yi) in zip(x,y) )
+
+ def inner_product(x,y):
+ return (x.T*y)[0,0]
# New defaults for keyword arguments. Don't orthonormalize
# because our basis is already orthonormal with respect to our
if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
-
- standard_basis = FreeModule(ZZ, n).basis()
- super(HadamardEJA, self).__init__(standard_basis,
- jordan_product,
- inner_product,
- **kwargs)
+ column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+ super().__init__(column_basis,
+ jordan_product,
+ inner_product,
+ associative=True,
+ **kwargs)
self.rank.set_cache(n)
if n == 0:
....: for j in range(n-1) ]
sage: actual == expected
True
+
"""
def __init__(self, B, **kwargs):
- if not B.is_positive_definite():
- raise ValueError("bilinear form is not positive-definite")
+ # The matrix "B" is supplied by the user in most cases,
+ # so it makes sense to check whether or not its positive-
+ # definite unless we are specifically asked not to...
+ if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
+ if not B.is_positive_definite():
+ raise ValueError("bilinear form is not positive-definite")
+
+ # However, all of the other data for this EJA is computed
+ # by us in manner that guarantees the axioms are
+ # satisfied. So, again, unless we are specifically asked to
+ # verify things, we'll skip the rest of the checks.
+ if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
def inner_product(x,y):
- return (B*x).inner_product(y)
+ return (y.T*B*x)[0,0]
def jordan_product(x,y):
P = x.parent()
- x0 = x[0]
- xbar = x[1:]
- y0 = y[0]
- ybar = y[1:]
- z0 = inner_product(x,y)
+ x0 = x[0,0]
+ xbar = x[1:,0]
+ y0 = y[0,0]
+ ybar = y[1:,0]
+ z0 = inner_product(y,x)
zbar = y0*xbar + x0*ybar
- return P((z0,) + tuple(zbar))
-
- # 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
+ return P([z0] + zbar.list())
n = B.nrows()
- standard_basis = FreeModule(ZZ, n).basis()
- super(BilinearFormEJA, self).__init__(standard_basis,
+ column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
+ super(BilinearFormEJA, self).__init__(column_basis,
jordan_product,
inner_product,
**kwargs)
# inappropriate for us.
return cls(**kwargs)
-class DirectSumEJA(ConcreteEJA):
+
+class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
+ FiniteDimensionalEJA):
r"""
- The external (orthogonal) direct sum of two other Euclidean Jordan
- algebras. Essentially the Cartesian product of its two factors.
- Every Euclidean Jordan algebra decomposes into an orthogonal
- direct sum of simple Euclidean Jordan algebras, so no generality
- is lost by providing only this construction.
+ The external (orthogonal) direct sum of two or more Euclidean
+ Jordan algebras. Every Euclidean Jordan algebra decomposes into an
+ orthogonal direct sum of simple Euclidean Jordan algebras which is
+ then isometric to a Cartesian product, so no generality is lost by
+ providing only this construction.
SETUP::
sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: CartesianProductEJA,
....: HadamardEJA,
- ....: RealSymmetricEJA,
- ....: DirectSumEJA)
+ ....: JordanSpinEJA,
+ ....: RealSymmetricEJA)
- EXAMPLES::
+ EXAMPLES:
+
+ The Jordan product is inherited from our factors and implemented by
+ our CombinatorialFreeModule Cartesian product superclass::
+ sage: set_random_seed()
sage: J1 = HadamardEJA(2)
- sage: J2 = RealSymmetricEJA(3)
- sage: J = DirectSumEJA(J1,J2)
- sage: J.dimension()
- 8
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: x,y = J.random_elements(2)
+ sage: x*y in J
+ True
+
+ The ability to retrieve the original factors is implemented by our
+ CombinatorialFreeModule Cartesian product superclass::
+
+ sage: J1 = HadamardEJA(2, field=QQ)
+ sage: J2 = JordanSpinEJA(3, field=QQ)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.cartesian_factors()
+ (Euclidean Jordan algebra of dimension 2 over Rational Field,
+ Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+ You can provide more than two factors::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = JordanSpinEJA(3)
+ sage: J3 = RealSymmetricEJA(3)
+ sage: cartesian_product([J1,J2,J3])
+ Euclidean Jordan algebra of dimension 2 over Algebraic Real
+ Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
+ Real Field (+) Euclidean Jordan algebra of dimension 6 over
+ Algebraic Real Field
+
+ Rank is additive on a Cartesian product::
+
+ sage: J1 = HadamardEJA(1)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: J1.rank.clear_cache()
+ sage: J2.rank.clear_cache()
+ sage: J.rank.clear_cache()
sage: J.rank()
- 5
+ 3
+ sage: J.rank() == J1.rank() + J2.rank()
+ True
+
+ The same rank computation works over the rationals, with whatever
+ basis you like::
+
+ sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
+ sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
+ sage: J = cartesian_product([J1,J2])
+ sage: J1.rank.clear_cache()
+ sage: J2.rank.clear_cache()
+ sage: J.rank.clear_cache()
+ sage: J.rank()
+ 3
+ sage: J.rank() == J1.rank() + J2.rank()
+ True
+
+ The product algebra will be associative if and only if all of its
+ components are associative::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J1.is_associative()
+ True
+ sage: J2 = HadamardEJA(3)
+ sage: J2.is_associative()
+ True
+ sage: J3 = RealSymmetricEJA(3)
+ sage: J3.is_associative()
+ False
+ sage: CP1 = cartesian_product([J1,J2])
+ sage: CP1.is_associative()
+ True
+ sage: CP2 = cartesian_product([J1,J3])
+ sage: CP2.is_associative()
+ False
TESTS:
- The external direct sum construction is only valid when the two factors
- have the same base ring; an error is raised otherwise::
+ All factors must share the same base field::
- sage: set_random_seed()
- sage: J1 = random_eja(field=AA)
- sage: J2 = random_eja(field=QQ,orthonormalize=False)
- sage: J = DirectSumEJA(J1,J2)
+ sage: J1 = HadamardEJA(2, field=QQ)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: CartesianProductEJA((J1,J2))
Traceback (most recent call last):
...
- ValueError: algebras must share the same base field
+ ValueError: all factors must share the same base field
+
+ The cached unit element is the same one that would be computed::
+
+ sage: set_random_seed() # long time
+ sage: J1 = random_eja() # long time
+ sage: J2 = random_eja() # long time
+ sage: J = cartesian_product([J1,J2]) # long time
+ sage: actual = J.one() # long time
+ sage: J.one.clear_cache() # long time
+ sage: expected = J.one() # long time
+ sage: actual == expected # long time
+ True
"""
- def __init__(self, J1, J2, **kwargs):
- if J1.base_ring() != J2.base_ring():
- raise ValueError("algebras must share the same base field")
- field = J1.base_ring()
-
- self._factors = (J1, J2)
- n1 = J1.dimension()
- n2 = J2.dimension()
- n = n1+n2
- V = VectorSpace(field, n)
- mult_table = [ [ V.zero() for j in range(i+1) ]
- for i in range(n) ]
- for i in range(n1):
- for j in range(i+1):
- p = (J1.monomial(i)*J1.monomial(j)).to_vector()
- mult_table[i][j] = V(p.list() + [field.zero()]*n2)
+ Element = FiniteDimensionalEJAElement
- for i in range(n2):
- for j in range(i+1):
- p = (J2.monomial(i)*J2.monomial(j)).to_vector()
- mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
-
- # TODO: build the IP table here from the two constituent IP
- # matrices (it'll be block diagonal, I think).
- ip_table = [ [ field.zero() for j in range(i+1) ]
- for i in range(n) ]
- super(DirectSumEJA, self).__init__(field,
- mult_table,
- ip_table,
- check_axioms=False,
- **kwargs)
- self.rank.set_cache(J1.rank() + J2.rank())
-
-
- def factors(self):
+
+ 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 ):
+ raise ValueError("all factors must share the same base field")
+
+ associative = all( m.is_associative() for m in algebras )
+
+ # 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()
+ )
+
+ # 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)
+ ))
+
+ def inner_product(x, y):
+ return sum(
+ Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m)
+ )
+
+ # 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,
+ orthonormalize=False,
+ associative=associative,
+ cartesian_product=True,
+ 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))
+
+ def matrix_space(self):
r"""
- Return the pair of this algebra's factors.
+ Return the space that our matrix basis lives in as a Cartesian
+ product.
SETUP::
sage: from mjo.eja.eja_algebra import (HadamardEJA,
- ....: JordanSpinEJA,
- ....: DirectSumEJA)
+ ....: RealSymmetricEJA)
EXAMPLES::
- sage: J1 = HadamardEJA(2, field=QQ)
- sage: J2 = JordanSpinEJA(3, field=QQ)
- sage: J = DirectSumEJA(J1,J2)
- sage: J.factors()
- (Euclidean Jordan algebra of dimension 2 over Rational Field,
- Euclidean Jordan algebra of dimension 3 over Rational Field)
+ 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)
"""
- return self._factors
+ from sage.categories.cartesian_product import cartesian_product
+ return cartesian_product( [J.matrix_space()
+ for J in self.cartesian_factors()] )
- def projections(self):
+ @cached_method
+ def cartesian_projection(self, i):
r"""
- Return a pair of projections onto this algebra's factors.
-
SETUP::
- sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
- ....: ComplexHermitianEJA,
- ....: DirectSumEJA)
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: JordanSpinEJA,
+ ....: HadamardEJA,
+ ....: RealSymmetricEJA,
+ ....: ComplexHermitianEJA)
- EXAMPLES::
+ EXAMPLES:
+
+ The projection morphisms are Euclidean Jordan algebra
+ operators::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.cartesian_projection(0)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [1 0 0 0 0]
+ [0 1 0 0 0]
+ Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
+ Real Field (+) Euclidean Jordan algebra of dimension 3 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
+ Real Field
+ sage: J.cartesian_projection(1)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [0 0 1 0 0]
+ [0 0 0 1 0]
+ [0 0 0 0 1]
+ Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
+ Real Field (+) Euclidean Jordan algebra of dimension 3 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
+ Real Field
+
+ The projections work the way you'd expect on the vector
+ representation of an element::
sage: J1 = JordanSpinEJA(2)
sage: J2 = ComplexHermitianEJA(2)
- sage: J = DirectSumEJA(J1,J2)
- sage: (pi_left, pi_right) = J.projections()
- sage: J.one().to_vector()
- (1, 0, 1, 0, 0, 1)
+ sage: J = cartesian_product([J1,J2])
+ sage: pi_left = J.cartesian_projection(0)
+ sage: pi_right = J.cartesian_projection(1)
sage: pi_left(J.one()).to_vector()
(1, 0)
sage: pi_right(J.one()).to_vector()
(1, 0, 0, 1)
+ sage: J.one().to_vector()
+ (1, 0, 1, 0, 0, 1)
+
+ TESTS:
+
+ The answer never changes::
+
+ sage: set_random_seed()
+ sage: J1 = random_eja()
+ sage: J2 = random_eja()
+ sage: J = cartesian_product([J1,J2])
+ sage: P0 = J.cartesian_projection(0)
+ sage: P1 = J.cartesian_projection(0)
+ sage: P0 == P1
+ True
"""
- (J1,J2) = self.factors()
- m = J1.dimension()
- n = J2.dimension()
- V_basis = self.vector_space().basis()
- # Need to specify the dimensions explicitly so that we don't
- # wind up with a zero-by-zero matrix when we want e.g. a
- # zero-by-two matrix (important for composing things).
- P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
- P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
- pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
- pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
- return (pi_left, pi_right)
-
- def inclusions(self):
- r"""
- Return the pair of inclusion maps from our factors into us.
+ Ji = self.cartesian_factors()[i]
+ # Requires the fix on Trac 31421/31422 to work!
+ Pi = super().cartesian_projection(i)
+ return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
+ @cached_method
+ def cartesian_embedding(self, i):
+ r"""
SETUP::
sage: from mjo.eja.eja_algebra import (random_eja,
....: JordanSpinEJA,
- ....: RealSymmetricEJA,
- ....: DirectSumEJA)
+ ....: HadamardEJA,
+ ....: RealSymmetricEJA)
- EXAMPLES::
+ EXAMPLES:
+
+ The embedding morphisms are Euclidean Jordan algebra
+ operators::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = cartesian_product([J1,J2])
+ sage: J.cartesian_embedding(0)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [1 0]
+ [0 1]
+ [0 0]
+ [0 0]
+ [0 0]
+ Domain: Euclidean Jordan algebra of dimension 2 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 2 over
+ Algebraic Real Field (+) Euclidean Jordan algebra of
+ dimension 3 over Algebraic Real Field
+ sage: J.cartesian_embedding(1)
+ Linear operator between finite-dimensional Euclidean Jordan
+ algebras represented by the matrix:
+ [0 0 0]
+ [0 0 0]
+ [1 0 0]
+ [0 1 0]
+ [0 0 1]
+ Domain: Euclidean Jordan algebra of dimension 3 over
+ Algebraic Real Field
+ Codomain: Euclidean Jordan algebra of dimension 2 over
+ Algebraic Real Field (+) Euclidean Jordan algebra of
+ dimension 3 over Algebraic Real Field
+
+ The embeddings work the way you'd expect on the vector
+ representation of an element::
sage: J1 = JordanSpinEJA(3)
sage: J2 = RealSymmetricEJA(2)
- sage: J = DirectSumEJA(J1,J2)
- sage: (iota_left, iota_right) = J.inclusions()
+ sage: J = cartesian_product([J1,J2])
+ sage: iota_left = J.cartesian_embedding(0)
+ sage: iota_right = J.cartesian_embedding(1)
sage: iota_left(J1.zero()) == J.zero()
True
sage: iota_right(J2.zero()) == J.zero()
TESTS:
+ The answer never changes::
+
+ sage: set_random_seed()
+ sage: J1 = random_eja()
+ sage: J2 = random_eja()
+ sage: J = cartesian_product([J1,J2])
+ sage: E0 = J.cartesian_embedding(0)
+ sage: E1 = J.cartesian_embedding(0)
+ sage: E0 == E1
+ True
+
Composing a projection with the corresponding inclusion should
produce the identity map, and mismatching them should produce
the zero map::
sage: set_random_seed()
sage: J1 = random_eja()
sage: J2 = random_eja()
- sage: J = DirectSumEJA(J1,J2)
- sage: (iota_left, iota_right) = J.inclusions()
- sage: (pi_left, pi_right) = J.projections()
+ sage: J = cartesian_product([J1,J2])
+ sage: iota_left = J.cartesian_embedding(0)
+ sage: iota_right = J.cartesian_embedding(1)
+ sage: pi_left = J.cartesian_projection(0)
+ sage: pi_right = J.cartesian_projection(1)
sage: pi_left*iota_left == J1.one().operator()
True
sage: pi_right*iota_right == J2.one().operator()
True
"""
- (J1,J2) = self.factors()
- m = J1.dimension()
- n = J2.dimension()
- V_basis = self.vector_space().basis()
- # Need to specify the dimensions explicitly so that we don't
- # wind up with a zero-by-zero matrix when we want e.g. a
- # two-by-zero matrix (important for composing things).
- I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
- I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
- iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
- iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
- return (iota_left, iota_right)
-
- def inner_product(self, x, y):
- r"""
- The standard Cartesian inner-product.
-
- We project ``x`` and ``y`` onto our factors, and add up the
- inner-products from the subalgebras.
-
- SETUP::
-
-
- sage: from mjo.eja.eja_algebra import (HadamardEJA,
- ....: QuaternionHermitianEJA,
- ....: DirectSumEJA)
-
- EXAMPLE::
-
- sage: J1 = HadamardEJA(3,field=QQ)
- sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
- sage: J = DirectSumEJA(J1,J2)
- sage: x1 = J1.one()
- sage: x2 = x1
- sage: y1 = J2.one()
- sage: y2 = y1
- sage: x1.inner_product(x2)
- 3
- sage: y1.inner_product(y2)
- 2
- sage: J.one().inner_product(J.one())
- 5
-
- """
- (pi_left, pi_right) = self.projections()
- x1 = pi_left(x)
- x2 = pi_right(x)
- y1 = pi_left(y)
- y2 = pi_right(y)
+ Ji = self.cartesian_factors()[i]
+ # Requires the fix on Trac 31421/31422 to work!
+ Ei = super().cartesian_embedding(i)
+ return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
- return (x1.inner_product(y1) + x2.inner_product(y2))
+FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
random_eja = ConcreteEJA.random_instance