"""
SETUP::
- sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
+ sage: from mjo.eja.eja_algebra import (
+ ....: FiniteDimensionalEuclideanJordanAlgebra,
+ ....: JordanSpinEJA,
+ ....: random_eja)
EXAMPLES:
TESTS:
- The ``field`` we're given must be real::
+ The ``field`` we're given must be real with ``check=True``::
sage: JordanSpinEJA(2,QQbar)
Traceback (most recent call last):
...
ValueError: field is not real
+ The multiplication table must be square with ``check=True``::
+
+ sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
+ Traceback (most recent call last):
+ ...
+ ValueError: multiplication table is not square
+
"""
if check:
if not field.is_subring(RR):
category = MagmaticAlgebras(field).FiniteDimensional()
category = category.WithBasis().Unital()
+ # The multiplication table had better be square
+ n = len(mult_table)
+ if check:
+ if not all( len(l) == n for l in mult_table ):
+ raise ValueError("multiplication table is not square")
+
fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
fda.__init__(field,
- range(len(mult_table)),
+ range(n),
prefix=prefix,
category=category)
self.print_options(bracket='')
for ls in mult_table
]
+ if check:
+ if not self._is_commutative():
+ raise ValueError("algebra is not commutative")
+ if not self._is_jordanian():
+ raise ValueError("Jordan identity does not hold")
+ if not self._inner_product_is_associative():
+ raise ValueError("inner product is not associative")
def _element_constructor_(self, elt):
"""
def product_on_basis(self, i, j):
return self._multiplication_table[i][j]
- def _a_regular_element(self):
- """
- Guess a regular element. Needed to compute the basis for our
- characteristic polynomial coefficients.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import random_eja
-
- TESTS:
-
- Ensure that this hacky method succeeds for every algebra that we
- know how to construct::
-
- sage: set_random_seed()
- sage: J = random_eja()
- sage: J._a_regular_element().is_regular()
- True
+ 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=False`` and passed
+ an invalid multiplication table.
"""
- gs = self.gens()
- z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
- if not z.is_regular():
- raise ValueError("don't know a regular element")
- return z
+ 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()) )
+ def _is_jordanian(self):
+ r"""
+ Whether or not this algebra's multiplication table respects the
+ Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
+
+ 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
+ return ``True``, unless this algebra was constructed with
+ ``check=False`` and passed an invalid multiplication table.
+ """
+ return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
+ ==
+ (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
+ for i in range(self.dimension())
+ for j in range(self.dimension()) )
+
+ def _inner_product_is_associative(self):
+ r"""
+ Return whether or not this algebra's inner product `B` is
+ associative; that is, whether or not `B(xy,z) = B(x,yz)`.
- @cached_method
- def _charpoly_basis_space(self):
- """
- Return the vector space spanned by the basis used in our
- characteristic polynomial coefficients. This is used not only to
- compute those coefficients, but also any time we need to
- evaluate the coefficients (like when we compute the trace or
- determinant).
- """
- z = self._a_regular_element()
- # Don't use the parent vector space directly here in case this
- # happens to be a subalgebra. In that case, we would be e.g.
- # two-dimensional but span_of_basis() would expect three
- # coordinates.
- V = VectorSpace(self.base_ring(), self.vector_space().dimension())
- basis = [ (z**k).to_vector() for k in range(self.rank()) ]
- V1 = V.span_of_basis( basis )
- b = (V1.basis() + V1.complement().basis())
- return V.span_of_basis(b)
+ This method should of course always return ``True``, unless
+ this algebra was constructed with ``check=False`` and passed
+ an invalid multiplication table.
+ """
+ # Used to check whether or not something is zero in an inexact
+ # ring. This number is sufficient to allow the construction of
+ # QuaternionHermitianEJA(2, RDF) with check=True.
+ epsilon = 1e-16
+ 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).inner_product(z) - x.inner_product(y*z)
- @cached_method
- def _charpoly_coeff(self, i):
- """
- Return the coefficient polynomial "a_{i}" of this algebra's
- general characteristic polynomial.
-
- Having this be a separate cached method lets us compute and
- store the trace/determinant (a_{r-1} and a_{0} respectively)
- separate from the entire characteristic polynomial.
- """
- (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
- R = A_of_x.base_ring()
-
- if i == self.rank():
- return R.one()
- if i > self.rank():
- # Guaranteed by theory
- return R.zero()
-
- # Danger: the in-place modification is done for performance
- # reasons (reconstructing a matrix with huge polynomial
- # entries is slow), but I don't know how cached_method works,
- # so it's highly possible that we're modifying some global
- # list variable by reference, here. In other words, you
- # probably shouldn't call this method twice on the same
- # algebra, at the same time, in two threads
- Ai_orig = A_of_x.column(i)
- A_of_x.set_column(i,xr)
- numerator = A_of_x.det()
- A_of_x.set_column(i,Ai_orig)
-
- # We're relying on the theory here to ensure that each a_i is
- # indeed back in R, and the added negative signs are to make
- # the whole charpoly expression sum to zero.
- return R(-numerator/detA)
+ if self.base_ring().is_exact():
+ if diff != 0:
+ return False
+ else:
+ if diff.abs() > epsilon:
+ return False
+ return True
@cached_method
- def _charpoly_matrix_system(self):
- """
- Compute the matrix whose entries A_ij are polynomials in
- X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
- corresponding to `x^r` and the determinent of the matrix A =
- [A_ij]. In other words, all of the fixed (cachable) data needed
- to compute the coefficients of the characteristic polynomial.
+ def characteristic_polynomial_of(self):
"""
- r = self.rank()
- n = self.dimension()
-
- # Turn my vector space into a module so that "vectors" can
- # have multivatiate polynomial entries.
- names = tuple('X' + str(i) for i in range(1,n+1))
- R = PolynomialRing(self.base_ring(), names)
-
- # Using change_ring() on the parent's vector space doesn't work
- # here because, in a subalgebra, that vector space has a basis
- # and change_ring() tries to bring the basis along with it. And
- # that doesn't work unless the new ring is a PID, which it usually
- # won't be.
- V = FreeModule(R,n)
-
- # Now let x = (X1,X2,...,Xn) be the vector whose entries are
- # indeterminates...
- x = V(names)
-
- # And figure out the "left multiplication by x" matrix in
- # that setting.
- lmbx_cols = []
- monomial_matrices = [ self.monomial(i).operator().matrix()
- for i in range(n) ] # don't recompute these!
- for k in range(n):
- ek = self.monomial(k).to_vector()
- lmbx_cols.append(
- sum( x[i]*(monomial_matrices[i]*ek)
- for i in range(n) ) )
- Lx = matrix.column(R, lmbx_cols)
-
- # Now we can compute powers of x "symbolically"
- x_powers = [self.one().to_vector(), x]
- for d in range(2, r+1):
- x_powers.append( Lx*(x_powers[-1]) )
-
- idmat = matrix.identity(R, n)
-
- W = self._charpoly_basis_space()
- W = W.change_ring(R.fraction_field())
-
- # Starting with the standard coordinates x = (X1,X2,...,Xn)
- # and then converting the entries to W-coordinates allows us
- # to pass in the standard coordinates to the charpoly and get
- # back the right answer. Specifically, with x = (X1,X2,...,Xn),
- # we have
- #
- # W.coordinates(x^2) eval'd at (standard z-coords)
- # =
- # W-coords of (z^2)
- # =
- # W-coords of (standard coords of x^2 eval'd at std-coords of z)
- #
- # We want the middle equivalent thing in our matrix, but use
- # the first equivalent thing instead so that we can pass in
- # standard coordinates.
- x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
- l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
- A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
- return (A_of_x, x, x_powers[r], A_of_x.det())
-
-
- @cached_method
- def characteristic_polynomial(self):
- """
- Return a characteristic polynomial that works for all elements
- of this algebra.
+ Return the algebra's "characteristic polynomial of" function,
+ which is itself a multivariate polynomial that, when evaluated
+ at the coordinates of some algebra element, returns that
+ element's characteristic polynomial.
The resulting polynomial has `n+1` variables, where `n` is the
dimension of this algebra. The first `n` variables correspond to
Alizadeh, Example 11.11::
sage: J = JordanSpinEJA(3)
- sage: p = J.characteristic_polynomial(); p
+ sage: p = J.characteristic_polynomial_of(); p
X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
sage: xvec = J.one().to_vector()
sage: p(*xvec)
any argument::
sage: J = TrivialEJA()
- sage: J.characteristic_polynomial()
+ sage: J.characteristic_polynomial_of()
1
"""
r = self.rank()
n = self.dimension()
- # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
- a = [ self._charpoly_coeff(i) for i in range(r+1) ]
+ # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
+ a = self._charpoly_coefficients()
# We go to a bit of trouble here to reorder the
# indeterminates, so that it's easier to evaluate the
# characteristic polynomial at x's coordinates and get back
# something in terms of t, which is what we want.
- R = a[0].parent()
S = PolynomialRing(self.base_ring(),'t')
t = S.gen(0)
- S = PolynomialRing(S, R.variable_names())
- t = S(t)
+ if r > 0:
+ R = a[0].parent()
+ S = PolynomialRing(S, R.variable_names())
+ t = S(t)
- return sum( a[k]*(t**k) for k in range(len(a)) )
+ return (t**r + sum( a[k]*(t**k) for k in range(r) ))
def inner_product(self, x, y):
"""
Return the matrix space in which this algebra's natural basis
elements live.
+
+ Generally this will be an `n`-by-`1` column-vector space,
+ except when the algebra is trivial. There it's `n`-by-`n`
+ (where `n` is zero), to ensure that two elements of the
+ natural basis space (empty matrices) can be multiplied.
"""
- if self._natural_basis is None or len(self._natural_basis) == 0:
+ if self.is_trivial():
+ return MatrixSpace(self.base_ring(), 0)
+ elif self._natural_basis is None or len(self._natural_basis) == 0:
return MatrixSpace(self.base_ring(), self.dimension(), 1)
else:
return self._natural_basis[0].matrix_space()
Vector space of degree 6 and dimension 2...
sage: J1
Euclidean Jordan algebra of dimension 3...
+ sage: J0.one().natural_representation()
+ [0 0 0]
+ [0 0 0]
+ [0 0 1]
+ sage: orig_df = AA.options.display_format
+ sage: AA.options.display_format = 'radical'
+ sage: J.from_vector(J5.basis()[0]).natural_representation()
+ [ 0 0 1/2*sqrt(2)]
+ [ 0 0 0]
+ [1/2*sqrt(2) 0 0]
+ sage: J.from_vector(J5.basis()[1]).natural_representation()
+ [ 0 0 0]
+ [ 0 0 1/2*sqrt(2)]
+ [ 0 1/2*sqrt(2) 0]
+ sage: AA.options.display_format = orig_df
+ sage: J1.one().natural_representation()
+ [1 0 0]
+ [0 1 0]
+ [0 0 0]
TESTS:
sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
True
- The identity elements in the two subalgebras are the
- projections onto their respective subspaces of the
- superalgebra's identity element::
+ The decomposition is into eigenspaces, and its components are
+ therefore necessarily orthogonal. Moreover, the identity
+ elements in the two subalgebras are the projections onto their
+ respective subspaces of the superalgebra's identity element::
sage: set_random_seed()
sage: J = random_eja()
....: x = J.random_element()
sage: c = x.subalgebra_idempotent()
sage: J0,J5,J1 = J.peirce_decomposition(c)
+ sage: ipsum = 0
+ sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
+ ....: w = w.superalgebra_element()
+ ....: y = J.from_vector(y)
+ ....: z = z.superalgebra_element()
+ ....: ipsum += w.inner_product(y).abs()
+ ....: ipsum += w.inner_product(z).abs()
+ ....: ipsum += y.inner_product(z).abs()
+ sage: ipsum
+ 0
sage: J1(c) == J1.one()
True
sage: J0(J.one() - c) == J0.one()
return (J0, J5, J1)
- def random_elements(self, count):
+ def random_element(self, thorough=False):
+ r"""
+ Return a random element of this algebra.
+
+ Our algebra superclass method only returns a linear
+ combination of at most two basis elements. We instead
+ want the vector space "random element" method that
+ returns a more diverse selection.
+
+ INPUT:
+
+ - ``thorough`` -- (boolean; default False) whether or not we
+ should generate irrational coefficients for the random
+ element when our base ring is irrational; this slows the
+ algebra operations to a crawl, but any truly random method
+ should include them
+
+ """
+ # For a general base ring... maybe we can trust this to do the
+ # right thing? Unlikely, but.
+ V = self.vector_space()
+ v = V.random_element()
+
+ if self.base_ring() is AA:
+ # The "random element" method of the algebraic reals is
+ # stupid at the moment, and only returns integers between
+ # -2 and 2, inclusive:
+ #
+ # https://trac.sagemath.org/ticket/30875
+ #
+ # Instead, we implement our own "random vector" method,
+ # and then coerce that into the algebra. We use the vector
+ # space degree here instead of the dimension because a
+ # subalgebra could (for example) be spanned by only two
+ # vectors, each with five coordinates. We need to
+ # generate all five coordinates.
+ if thorough:
+ v *= QQbar.random_element().real()
+ else:
+ v *= QQ.random_element()
+
+ return self.from_vector(V.coordinate_vector(v))
+
+ def random_elements(self, count, thorough=False):
"""
Return ``count`` random elements as a tuple.
+ INPUT:
+
+ - ``thorough`` -- (boolean; default False) whether or not we
+ should generate irrational coefficients for the random
+ elements when our base ring is irrational; this slows the
+ algebra operations to a crawl, but any truly random method
+ should include them
+
SETUP::
sage: from mjo.eja.eja_algebra import JordanSpinEJA
True
"""
- return tuple( self.random_element() for idx in range(count) )
+ return tuple( self.random_element(thorough)
+ for idx in range(count) )
@classmethod
def random_instance(cls, field=AA, **kwargs):
# there's only one.
return cls(field)
- n = ZZ.random_element(cls._max_test_case_size()) + 1
+ n = ZZ.random_element(cls._max_test_case_size() + 1)
return cls(n, field, **kwargs)
@cached_method
- def rank(self):
+ def _charpoly_coefficients(self):
+ r"""
+ The `r` polynomial coefficients of the "characteristic polynomial
+ of" function.
"""
- Return the rank of this EJA.
+ n = self.dimension()
+ var_names = [ "X" + str(z) for z in range(1,n+1) ]
+ R = PolynomialRing(self.base_ring(), var_names)
+ vars = R.gens()
+ F = R.fraction_field()
+
+ 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]
+ for k in range(n) )
- ALGORITHM:
+ L_x = matrix(F, n, n, L_x_i_j)
+
+ r = None
+ if self.rank.is_in_cache():
+ r = self.rank()
+ # There's no need to pad the system with redundant
+ # columns if we *know* they'll be redundant.
+ n = r
+
+ # Compute an extra power in case the rank is equal to
+ # the dimension (otherwise, we would stop at x^(r-1)).
+ x_powers = [ (L_x**k)*self.one().to_vector()
+ for k in range(n+1) ]
+ A = matrix.column(F, x_powers[:n])
+ AE = A.extended_echelon_form()
+ E = AE[:,n:]
+ A_rref = AE[:,:n]
+ if r is None:
+ r = A_rref.rank()
+ b = x_powers[r]
+
+ # 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]
- We first compute the polynomial "column matrices" `p_{k}` that
- evaluate to `x^k` on the coordinates of `x`. Then, we begin
- adding them to a matrix one at a time, and trying to solve the
- system that makes `p_{0}`,`p_{1}`,..., `p_{s-1}` add up to
- `p_{s}`. This will succeed only when `s` is the rank of the
- algebra, as proven in a recent draft paper of mine.
+ @cached_method
+ def rank(self):
+ r"""
+ Return the rank of this EJA.
+
+ This is a cached method because we know the rank a priori for
+ all of the algebras we can construct. Thus we can avoid the
+ expensive ``_charpoly_coefficients()`` call unless we truly
+ need to compute the whole characteristic polynomial.
SETUP::
sage: J.rank.clear_cache()
sage: J.rank()
2
-
"""
- n = self.dimension()
- if n == 0:
- return 0
- elif n == 1:
- return 1
-
- var_names = [ "X" + str(z) for z in range(1,n+1) ]
- R = PolynomialRing(self.base_ring(), var_names)
- vars = R.gens()
-
- 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]
- for k in range(n) )
-
- L_x = matrix(R, n, n, L_x_i_j)
- x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
- for k in range(n) ]
-
- # Can assume n >= 2
- M = matrix([x_powers[0]])
- old_rank = 1
-
- for d in range(1,n):
- M = matrix(M.rows() + [x_powers[d]])
- M.echelonize()
- # TODO: we've basically solved the system here.
- # We should save the echelonized matrix somehow
- # so that it can be reused in the charpoly method.
- new_rank = M.rank()
- if new_rank == old_rank:
- return new_rank
- else:
- old_rank = new_rank
-
- return n
+ return len(self._charpoly_coefficients())
def vector_space(self):
return x.to_vector().inner_product(y.to_vector())
-def random_eja(field=AA, nontrivial=False):
+def random_eja(field=AA):
"""
Return a "random" finite-dimensional Euclidean Jordan Algebra.
Euclidean Jordan algebra of dimension...
"""
- eja_classes = [HadamardEJA,
- JordanSpinEJA,
- RealSymmetricEJA,
- ComplexHermitianEJA,
- QuaternionHermitianEJA]
- if not nontrivial:
- eja_classes.append(TrivialEJA)
- classname = choice(eja_classes)
+ classname = choice([TrivialEJA,
+ HadamardEJA,
+ JordanSpinEJA,
+ RealSymmetricEJA,
+ ComplexHermitianEJA,
+ QuaternionHermitianEJA])
return classname.random_instance(field=field)
-
-
class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
@staticmethod
def _max_test_case_size():
a multiplication table because the latter can be computed in terms
of the former when the product is known (like it is here).
"""
- # Used in this class's fast _charpoly_coeff() override.
+ # Used in this class's fast _charpoly_coefficients() override.
self._basis_normalizers = None
# We're going to loop through this a few times, so now's a good
@cached_method
- def rank(self):
+ def _charpoly_coefficients(self):
r"""
Override the parent method with something that tries to compute
over a faster (non-extension) field.
if self._basis_normalizers is None:
# We didn't normalize, so assume that the basis we started
# with had entries in a nice field.
- return super(MatrixEuclideanJordanAlgebra, self).rank()
+ return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
else:
basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
self._basis_normalizers) )
# Do this over the rationals and convert back at the end.
# Only works because we know the entries of the basis are
- # integers.
- J = MatrixEuclideanJordanAlgebra(QQ,
- basis,
- normalize_basis=False)
- return J.rank()
-
- @cached_method
- def _charpoly_coeff(self, i):
- """
- Override the parent method with something that tries to compute
- over a faster (non-extension) field.
- """
- if self._basis_normalizers is None:
- # We didn't normalize, so assume that the basis we started
- # with had entries in a nice field.
- return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i)
- else:
- basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
- self._basis_normalizers) )
-
- # Do this over the rationals and convert back at the end.
+ # integers. The argument ``check=False`` is required
+ # because the trace inner-product method for this
+ # class is a stub and can't actually be checked.
J = MatrixEuclideanJordanAlgebra(QQ,
basis,
- normalize_basis=False)
- (_,x,_,_) = J._charpoly_matrix_system()
- p = J._charpoly_coeff(i)
- # p might be missing some vars, have to substitute "optionally"
- pairs = zip(x.base_ring().gens(), self._basis_normalizers)
- substitutions = { v: v*c for (v,c) in pairs }
- result = p.subs(substitutions)
-
- # The result of "subs" can be either a coefficient-ring
- # element or a polynomial. Gotta handle both cases.
- if result in QQ:
- return self.base_ring()(result)
- else:
- return result.change_ring(self.base_ring())
+ normalize_basis=False,
+ check=False)
+ a = J._charpoly_coefficients()
+
+ # Unfortunately, changing the basis does change the
+ # coefficients of the characteristic polynomial, but since
+ # these are really the coefficients of the "characteristic
+ # polynomial of" function, everything is still nice and
+ # unevaluated. It's therefore "obvious" how scaling the
+ # basis affects the coordinate variables X1, X2, et
+ # cetera. Scaling the first basis vector up by "n" adds a
+ # factor of 1/n into every "X1" term, for example. So here
+ # we simply undo the basis_normalizer scaling that we
+ # performed earlier.
+ #
+ # The a[0] access here is safe because trivial algebras
+ # won't have any basis normalizers and therefore won't
+ # make it to this "else" branch.
+ XS = a[0].parent().gens()
+ subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
+ for i in range(len(XS)) }
+ return tuple( a_i.subs(subs_dict) for a_i in a )
@staticmethod
# is supposed to hold the entire long vector, and the subspace W
# of V will be spanned by the vectors that arise from symmetric
# matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
+ if len(basis) == 0:
+ return []
+
field = basis[0].base_ring()
dimension = basis[0].nrows()
sage: x.operator().matrix().is_symmetric()
True
+ We can construct the (trivial) algebra of rank zero::
+
+ sage: RealSymmetricEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
"""
@classmethod
def _denormalized_basis(cls, n, field):
sage: x.operator().matrix().is_symmetric()
True
+ We can construct the (trivial) algebra of rank zero::
+
+ sage: ComplexHermitianEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
"""
@classmethod
sage: x.operator().matrix().is_symmetric()
True
+ We can construct the (trivial) algebra of rank zero::
+
+ sage: QuaternionHermitianEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
"""
@classmethod
def _denormalized_basis(cls, n, field):
# largest subalgebra generated by any element.
fdeja.__init__(field, mult_table, **kwargs)
self.rank.set_cache(0)
+
+
+class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
+ 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.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (HadamardEJA,
+ ....: RealSymmetricEJA,
+ ....: DirectSumEJA)
+
+ EXAMPLES::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = RealSymmetricEJA(3)
+ sage: J = DirectSumEJA(J1,J2)
+ sage: J.dimension()
+ 8
+ sage: J.rank()
+ 5
+
+ """
+ def __init__(self, J1, J2, field=AA, **kwargs):
+ n1 = J1.dimension()
+ n2 = J2.dimension()
+ n = n1+n2
+ V = VectorSpace(field, n)
+ mult_table = [ [ V.zero() for j in range(n) ]
+ for i in range(n) ]
+ for i in range(n1):
+ for j in range(n1):
+ p = (J1.monomial(i)*J1.monomial(j)).to_vector()
+ mult_table[i][j] = V(p.list() + [field.zero()]*n2)
+
+ for i in range(n2):
+ for j in range(n2):
+ p = (J2.monomial(i)*J2.monomial(j)).to_vector()
+ mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
+
+ fdeja = super(DirectSumEJA, self)
+ fdeja.__init__(field, mult_table, **kwargs)
+ self.rank.set_cache(J1.rank() + J2.rank())