""" Euclidean Jordan Algebras. These are formally-real Jordan Algebras; specifically those where u^2 + v^2 = 0 implies that u = v = 0. They are used in optimization, and have some additional nice methods beyond what can be supported in a general Jordan Algebra. """ from itertools import izip, repeat from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.magmatic_algebras import MagmaticAlgebras from sage.combinat.free_module import CombinatorialFreeModule from sage.matrix.constructor import matrix from sage.matrix.matrix_space import MatrixSpace from sage.misc.cachefunc import cached_method from sage.misc.prandom import choice from sage.misc.table import table from sage.modules.free_module import FreeModule, VectorSpace from sage.rings.integer_ring import ZZ from sage.rings.number_field.number_field import QuadraticField from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.rational_field import QQ from sage.rings.real_lazy import CLF, RLF from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement from mjo.eja.eja_utils import _mat2vec class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # This is an ugly hack needed to prevent the category framework # from implementing a coercion from our base ring (e.g. the # rationals) into the algebra. First of all -- such a coercion is # nonsense to begin with. But more importantly, it tries to do so # in the category of rings, and since our algebras aren't # associative they generally won't be rings. _no_generic_basering_coercion = True def __init__(self, field, mult_table, rank, prefix='e', category=None, natural_basis=None): """ SETUP:: sage: from mjo.eja.eja_algebra import random_eja EXAMPLES: By definition, Jordan multiplication commutes:: sage: set_random_seed() sage: J = random_eja() sage: x,y = J.random_elements(2) sage: x*y == y*x True """ self._rank = rank self._natural_basis = natural_basis if category is None: category = MagmaticAlgebras(field).FiniteDimensional() category = category.WithBasis().Unital() fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) fda.__init__(field, range(len(mult_table)), prefix=prefix, category=category) self.print_options(bracket='') # The multiplication table we're given is necessarily in terms # of vectors, because we don't have an algebra yet for # anything to be an element of. However, it's faster in the # long run to have the multiplication table be in terms of # algebra elements. We do this after calling the superclass # constructor so that from_vector() knows what to do. self._multiplication_table = [ map(lambda x: self.from_vector(x), ls) for ls in mult_table ] def _element_constructor_(self, elt): """ Construct an element of this algebra from its natural representation. This gets called only after the parent element _call_ method fails to find a coercion for the argument. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: RealCartesianProductEJA, ....: RealSymmetricEJA) EXAMPLES: The identity in `S^n` is converted to the identity in the EJA:: sage: J = RealSymmetricEJA(3) sage: I = matrix.identity(QQ,3) sage: J(I) == J.one() True This skew-symmetric matrix can't be represented in the EJA:: sage: J = RealSymmetricEJA(3) sage: A = matrix(QQ,3, lambda i,j: i-j) sage: J(A) Traceback (most recent call last): ... ArithmeticError: vector is not in free module TESTS: Ensure that we can convert any element of the two non-matrix simple algebras (whose natural representations are their usual vector representations) back and forth faithfully:: sage: set_random_seed() sage: J = RealCartesianProductEJA.random_instance() sage: x = J.random_element() sage: J(x.to_vector().column()) == x True sage: J = JordanSpinEJA.random_instance() sage: x = J.random_element() sage: J(x.to_vector().column()) == x True """ if elt == 0: # The superclass implementation of random_element() # needs to be able to coerce "0" into the algebra. return self.zero() natural_basis = self.natural_basis() basis_space = natural_basis[0].matrix_space() if elt not in basis_space: raise ValueError("not a naturally-represented algebra element") # Thanks for nothing! Matrix spaces aren't vector spaces in # Sage, so we have to figure out its natural-basis coordinates # ourselves. We use the basis space's ring instead of the # element's ring because the basis space might be an algebraic # closure whereas the base ring of the 3-by-3 identity matrix # could be QQ instead of QQbar. V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols()) W = V.span_of_basis( _mat2vec(s) for s in natural_basis ) coords = W.coordinate_vector(_mat2vec(elt)) return self.from_vector(coords) def _repr_(self): """ Return a string representation of ``self``. SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA TESTS: Ensure that it says what we think it says:: sage: JordanSpinEJA(2, field=QQ) Euclidean Jordan algebra of dimension 2 over Rational Field sage: JordanSpinEJA(3, field=RDF) Euclidean Jordan algebra of dimension 3 over Real Double Field """ fmt = "Euclidean Jordan algebra of dimension {} over {}" return fmt.format(self.dimension(), self.base_ring()) 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 """ 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 @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) @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(): # 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) @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. """ 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. The resulting polynomial has `n+1` variables, where `n` is the dimension of this algebra. The first `n` variables correspond to the coordinates of an algebra element: when evaluated at the coordinates of an algebra element with respect to a certain basis, the result is a univariate polynomial (in the one remaining variable ``t``), namely the characteristic polynomial of that element. SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA EXAMPLES: The characteristic polynomial in the spin algebra is given in Alizadeh, Example 11.11:: sage: J = JordanSpinEJA(3) sage: p = J.characteristic_polynomial(); p X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 sage: xvec = J.one().to_vector() sage: p(*xvec) t^2 - 2*t + 1 """ r = self.rank() n = self.dimension() # The list of coefficient polynomials a_1, a_2, ..., a_n. a = [ self._charpoly_coeff(i) for i in range(n) ] # 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) # Note: all entries past the rth should be zero. The # coefficient of the highest power (x^r) is 1, but it doesn't # appear in the solution vector which contains coefficients # for the other powers (to make them sum to x^r). if (r < n): a[r] = 1 # corresponds to x^r else: # When the rank is equal to the dimension, trying to # assign a[r] goes out-of-bounds. a.append(1) # corresponds to x^r return sum( a[k]*(t**k) for k in xrange(len(a)) ) def inner_product(self, x, y): """ The inner product associated with this Euclidean Jordan algebra. Defaults to the trace inner product, but can be overridden by subclasses if they are sure that the necessary properties are satisfied. SETUP:: sage: from mjo.eja.eja_algebra import random_eja EXAMPLES: Our inner product is "associative," which means the following for a symmetric bilinear form:: sage: set_random_seed() sage: J = random_eja() sage: x,y,z = J.random_elements(3) sage: (x*y).inner_product(z) == y.inner_product(x*z) True """ X = x.natural_representation() Y = y.natural_representation() return self.natural_inner_product(X,Y) def is_trivial(self): """ Return whether or not this algebra is trivial. A trivial algebra contains only the zero element. SETUP:: sage: from mjo.eja.eja_algebra import ComplexHermitianEJA EXAMPLES:: sage: J = ComplexHermitianEJA(3) sage: J.is_trivial() False sage: A = J.zero().subalgebra_generated_by() sage: A.is_trivial() True """ return self.dimension() == 0 def multiplication_table(self): """ Return a visual representation of this algebra's multiplication table (on basis elements). SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA EXAMPLES:: sage: J = JordanSpinEJA(4) sage: J.multiplication_table() +----++----+----+----+----+ | * || e0 | e1 | e2 | e3 | +====++====+====+====+====+ | e0 || e0 | e1 | e2 | e3 | +----++----+----+----+----+ | e1 || e1 | e0 | 0 | 0 | +----++----+----+----+----+ | e2 || e2 | 0 | e0 | 0 | +----++----+----+----+----+ | e3 || e3 | 0 | 0 | e0 | +----++----+----+----+----+ """ M = list(self._multiplication_table) # copy for i in xrange(len(M)): # M had better be "square" M[i] = [self.monomial(i)] + M[i] M = [["*"] + list(self.gens())] + M return table(M, header_row=True, header_column=True, frame=True) def natural_basis(self): """ Return a more-natural representation of this algebra's basis. Every finite-dimensional Euclidean Jordan Algebra is a direct sum of five simple algebras, four of which comprise Hermitian matrices. This method returns the original "natural" basis for our underlying vector space. (Typically, the natural basis is used to construct the multiplication table in the first place.) Note that this will always return a matrix. The standard basis in `R^n` will be returned as `n`-by-`1` column matrices. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: RealSymmetricEJA) EXAMPLES:: sage: J = RealSymmetricEJA(2) sage: J.basis() Finite family {0: e0, 1: e1, 2: e2} sage: J.natural_basis() ( [1 0] [ 0 1/2*sqrt2] [0 0] [0 0], [1/2*sqrt2 0], [0 1] ) :: sage: J = JordanSpinEJA(2) sage: J.basis() Finite family {0: e0, 1: e1} sage: J.natural_basis() ( [1] [0] [0], [1] ) """ if self._natural_basis is None: M = self.natural_basis_space() return tuple( M(b.to_vector()) for b in self.basis() ) else: return self._natural_basis def natural_basis_space(self): """ Return the matrix space in which this algebra's natural basis elements live. """ if 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() @staticmethod def natural_inner_product(X,Y): """ Compute the inner product of two naturally-represented elements. For example in the real symmetric matrix EJA, this will compute the trace inner-product of two n-by-n symmetric matrices. The default should work for the real cartesian product EJA, the Jordan spin EJA, and the real symmetric matrices. The others will have to be overridden. """ return (X.conjugate_transpose()*Y).trace() @cached_method def one(self): """ Return the unit element of this algebra. SETUP:: sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA, ....: random_eja) EXAMPLES:: sage: J = RealCartesianProductEJA(5) sage: J.one() e0 + e1 + e2 + e3 + e4 TESTS: The identity element acts like the identity:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: J.one()*x == x and x*J.one() == x True The matrix of the unit element's operator is the identity:: sage: set_random_seed() sage: J = random_eja() sage: actual = J.one().operator().matrix() sage: expected = matrix.identity(J.base_ring(), J.dimension()) sage: actual == expected True """ # We can brute-force compute the matrices of the operators # that correspond to the basis elements of this algebra. # If some linear combination of those basis elements is the # algebra identity, then the same linear combination of # their matrices has to be the identity matrix. # # Of course, matrices aren't vectors in sage, so we have to # appeal to the "long vectors" isometry. oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ] # Now we use basis linear algebra to find the coefficients, # of the matrices-as-vectors-linear-combination, which should # work for the original algebra basis too. A = matrix.column(self.base_ring(), oper_vecs) # We used the isometry on the left-hand side already, but we # still need to do it for the right-hand side. Recall that we # wanted something that summed to the identity matrix. b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) ) # Now if there's an identity element in the algebra, this should work. coeffs = A.solve_right(b) return self.linear_combination(zip(self.gens(), coeffs)) def random_element(self): # Temporary workaround for https://trac.sagemath.org/ticket/28327 if self.is_trivial(): return self.zero() else: s = super(FiniteDimensionalEuclideanJordanAlgebra, self) return s.random_element() def random_elements(self, count): """ Return ``count`` random elements as a tuple. SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA EXAMPLES:: sage: J = JordanSpinEJA(3) sage: x,y,z = J.random_elements(3) sage: all( [ x in J, y in J, z in J ]) True sage: len( J.random_elements(10) ) == 10 True """ return tuple( self.random_element() for idx in xrange(count) ) def rank(self): """ Return the rank of this EJA. ALGORITHM: The author knows of no algorithm to compute the rank of an EJA where only the multiplication table is known. In lieu of one, we require the rank to be specified when the algebra is created, and simply pass along that number here. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: RealSymmetricEJA, ....: ComplexHermitianEJA, ....: QuaternionHermitianEJA, ....: random_eja) EXAMPLES: The rank of the Jordan spin algebra is always two:: sage: JordanSpinEJA(2).rank() 2 sage: JordanSpinEJA(3).rank() 2 sage: JordanSpinEJA(4).rank() 2 The rank of the `n`-by-`n` Hermitian real, complex, or quaternion matrices is `n`:: sage: RealSymmetricEJA(4).rank() 4 sage: ComplexHermitianEJA(3).rank() 3 sage: QuaternionHermitianEJA(2).rank() 2 TESTS: Ensure that every EJA that we know how to construct has a positive integer rank:: sage: set_random_seed() sage: r = random_eja().rank() sage: r in ZZ and r > 0 True """ return self._rank def vector_space(self): """ Return the vector space that underlies this algebra. SETUP:: sage: from mjo.eja.eja_algebra import RealSymmetricEJA EXAMPLES:: sage: J = RealSymmetricEJA(2) sage: J.vector_space() Vector space of dimension 3 over... """ return self.zero().to_vector().parent().ambient_vector_space() Element = FiniteDimensionalEuclideanJordanAlgebraElement class KnownRankEJA(object): """ A class for algebras that we actually know we can construct. The main issue is that, for most of our methods to make sense, we need to know the rank of our algebra. Thus we can't simply generate a "random" algebra, or even check that a given basis and product satisfy the axioms; because even if everything looks OK, we wouldn't know the rank we need to actuallty build the thing. Not really a subclass of FDEJA because doing that causes method resolution errors, e.g. TypeError: Error when calling the metaclass bases Cannot create a consistent method resolution order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA """ @staticmethod def _max_test_case_size(): """ Return an integer "size" that is an upper bound on the size of this algebra when it is used in a random test case. Unfortunately, the term "size" is quite vague -- when dealing with `R^n` under either the Hadamard or Jordan spin product, the "size" refers to the dimension `n`. When dealing with a matrix algebra (real symmetric or complex/quaternion Hermitian), it refers to the size of the matrix, which is far less than the dimension of the underlying vector space. We default to five in this class, which is safe in `R^n`. The matrix algebra subclasses (or any class where the "size" is interpreted to be far less than the dimension) should override with a smaller number. """ return 5 @classmethod def random_instance(cls, field=QQ, **kwargs): """ Return a random instance of this type of algebra. Beware, this will crash for "most instances" because the constructor below looks wrong. """ n = ZZ.random_element(cls._max_test_case_size()) + 1 return cls(n, field, **kwargs) class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): """ Return the Euclidean Jordan Algebra corresponding to the set `R^n` under the Hadamard product. Note: this is nothing more than the Cartesian product of ``n`` copies of the spin algebra. Once Cartesian product algebras are implemented, this can go. SETUP:: sage: from mjo.eja.eja_algebra import RealCartesianProductEJA EXAMPLES: This multiplication table can be verified by hand:: sage: J = RealCartesianProductEJA(3) sage: e0,e1,e2 = J.gens() sage: e0*e0 e0 sage: e0*e1 0 sage: e0*e2 0 sage: e1*e1 e1 sage: e1*e2 0 sage: e2*e2 e2 TESTS: We can change the generator prefix:: sage: RealCartesianProductEJA(3, prefix='r').gens() (r0, r1, r2) """ def __init__(self, n, field=QQ, **kwargs): V = VectorSpace(field, n) mult_table = [ [ V.gen(i)*(i == j) for j in xrange(n) ] for i in xrange(n) ] fdeja = super(RealCartesianProductEJA, self) return fdeja.__init__(field, mult_table, rank=n, **kwargs) def inner_product(self, x, y): """ Faster to reimplement than to use natural representations. SETUP:: sage: from mjo.eja.eja_algebra import RealCartesianProductEJA TESTS: Ensure that this is the usual inner product for the algebras over `R^n`:: sage: set_random_seed() sage: J = RealCartesianProductEJA.random_instance() sage: x,y = J.random_elements(2) sage: X = x.natural_representation() sage: Y = y.natural_representation() sage: x.inner_product(y) == J.natural_inner_product(X,Y) True """ return x.to_vector().inner_product(y.to_vector()) def random_eja(): """ Return a "random" finite-dimensional Euclidean Jordan Algebra. ALGORITHM: For now, we choose a random natural number ``n`` (greater than zero) and then give you back one of the following: * The cartesian product of the rational numbers ``n`` times; this is ``QQ^n`` with the Hadamard product. * The Jordan spin algebra on ``QQ^n``. * The ``n``-by-``n`` rational symmetric matrices with the symmetric product. * The ``n``-by-``n`` complex-rational Hermitian matrices embedded in the space of ``2n``-by-``2n`` real symmetric matrices. * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded in the space of ``4n``-by-``4n`` real symmetric matrices. Later this might be extended to return Cartesian products of the EJAs above. SETUP:: sage: from mjo.eja.eja_algebra import random_eja TESTS:: sage: random_eja() Euclidean Jordan algebra of dimension... """ classname = choice(KnownRankEJA.__subclasses__()) return classname.random_instance() class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): @staticmethod def _max_test_case_size(): # Play it safe, since this will be squared and the underlying # field can have dimension 4 (quaternions) too. return 2 def __init__(self, field, basis, rank, normalize_basis=True, **kwargs): """ Compared to the superclass constructor, we take a basis instead of 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. self._basis_normalizers = None # We're going to loop through this a few times, so now's a good # time to ensure that it isn't a generator expression. basis = tuple(basis) if rank > 1 and normalize_basis: # We'll need sqrt(2) to normalize the basis, and this # winds up in the multiplication table, so the whole # algebra needs to be over the field extension. R = PolynomialRing(field, 'z') z = R.gen() p = z**2 - 2 if p.is_irreducible(): field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt()) basis = tuple( s.change_ring(field) for s in basis ) self._basis_normalizers = tuple( ~(self.natural_inner_product(s,s).sqrt()) for s in basis ) basis = tuple(s*c for (s,c) in izip(basis,self._basis_normalizers)) Qs = self.multiplication_table_from_matrix_basis(basis) fdeja = super(MatrixEuclideanJordanAlgebra, self) return fdeja.__init__(field, Qs, rank=rank, natural_basis=basis, **kwargs) @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 izip(self.natural_basis(), self._basis_normalizers) ) # Do this over the rationals and convert back at the end. J = MatrixEuclideanJordanAlgebra(QQ, basis, self.rank(), normalize_basis=False) (_,x,_,_) = J._charpoly_matrix_system() p = J._charpoly_coeff(i) # p might be missing some vars, have to substitute "optionally" pairs = izip(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()) @staticmethod def multiplication_table_from_matrix_basis(basis): """ At least three of the five simple Euclidean Jordan algebras have the symmetric multiplication (A,B) |-> (AB + BA)/2, where the multiplication on the right is matrix multiplication. Given a basis for the underlying matrix space, this function returns a multiplication table (obtained by looping through the basis elements) for an algebra of those matrices. """ # In S^2, for example, we nominally have four coordinates even # though the space is of dimension three only. The vector space V # 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. field = basis[0].base_ring() dimension = basis[0].nrows() V = VectorSpace(field, dimension**2) W = V.span_of_basis( _mat2vec(s) for s in basis ) n = len(basis) mult_table = [[W.zero() for j in xrange(n)] for i in xrange(n)] for i in xrange(n): for j in xrange(n): mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) return mult_table @staticmethod def real_embed(M): """ Embed the matrix ``M`` into a space of real matrices. The matrix ``M`` can have entries in any field at the moment: the real numbers, complex numbers, or quaternions. And although they are not a field, we can probably support octonions at some point, too. This function returns a real matrix that "acts like" the original with respect to matrix multiplication; i.e. real_embed(M*N) = real_embed(M)*real_embed(N) """ raise NotImplementedError @staticmethod def real_unembed(M): """ The inverse of :meth:`real_embed`. """ raise NotImplementedError @classmethod def natural_inner_product(cls,X,Y): Xu = cls.real_unembed(X) Yu = cls.real_unembed(Y) tr = (Xu*Yu).trace() if tr in RLF: # It's real already. return tr # Otherwise, try the thing that works for complex numbers; and # if that doesn't work, the thing that works for quaternions. try: return tr.vector()[0] # real part, imag part is index 1 except AttributeError: # A quaternions doesn't have a vector() method, but does # have coefficient_tuple() method that returns the # coefficients of 1, i, j, and k -- in that order. return tr.coefficient_tuple()[0] class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @staticmethod def real_embed(M): """ The identity function, for embedding real matrices into real matrices. """ return M @staticmethod def real_unembed(M): """ The identity function, for unembedding real matrices from real matrices. """ return M class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner product. It has dimension `(n^2 + n)/2` over the reals. SETUP:: sage: from mjo.eja.eja_algebra import RealSymmetricEJA EXAMPLES:: sage: J = RealSymmetricEJA(2) sage: e0, e1, e2 = J.gens() sage: e0*e0 e0 sage: e1*e1 1/2*e0 + 1/2*e2 sage: e2*e2 e2 In theory, our "field" can be any subfield of the reals:: sage: RealSymmetricEJA(2, AA) Euclidean Jordan algebra of dimension 3 over Algebraic Real Field sage: RealSymmetricEJA(2, RR) Euclidean Jordan algebra of dimension 3 over Real Field with 53 bits of precision TESTS: The dimension of this algebra is `(n^2 + n) / 2`:: sage: set_random_seed() sage: n_max = RealSymmetricEJA._max_test_case_size() sage: n = ZZ.random_element(1, n_max) sage: J = RealSymmetricEJA(n) sage: J.dimension() == (n^2 + n)/2 True The Jordan multiplication is what we think it is:: sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).natural_representation() sage: X = x.natural_representation() sage: Y = y.natural_representation() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True sage: J(expected) == x*y True We can change the generator prefix:: sage: RealSymmetricEJA(3, prefix='q').gens() (q0, q1, q2, q3, q4, q5) Our natural basis is normalized with respect to the natural inner product unless we specify otherwise:: sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() sage: all( b.norm() == 1 for b in J.gens() ) True Since our natural basis is normalized with respect to the natural inner product, and since we know that this algebra is an EJA, any left-multiplication operator's matrix will be symmetric because natural->EJA basis representation is an isometry and within the EJA the operator is self-adjoint by the Jordan axiom:: sage: set_random_seed() sage: x = RealSymmetricEJA.random_instance().random_element() sage: x.operator().matrix().is_symmetric() True """ @classmethod def _denormalized_basis(cls, n, field): """ Return a basis for the space of real symmetric n-by-n matrices. SETUP:: sage: from mjo.eja.eja_algebra import RealSymmetricEJA TESTS:: sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: B = RealSymmetricEJA._denormalized_basis(n,QQ) sage: all( M.is_symmetric() for M in B) True """ # The basis of symmetric matrices, as matrices, in their R^(n-by-n) # coordinates. S = [] for i in xrange(n): for j in xrange(i+1): Eij = matrix(field, n, lambda k,l: k==i and l==j) if i == j: Sij = Eij else: Sij = Eij + Eij.transpose() S.append(Sij) return S @staticmethod def _max_test_case_size(): return 4 # Dimension 10 def __init__(self, n, field=QQ, **kwargs): basis = self._denormalized_basis(n, field) super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs) class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @staticmethod def real_embed(M): """ Embed the n-by-n complex matrix ``M`` into the space of real matrices of size 2n-by-2n via the map the sends each entry `z = a + bi` to the block matrix ``[[a,b],[-b,a]]``. SETUP:: sage: from mjo.eja.eja_algebra import \ ....: ComplexMatrixEuclideanJordanAlgebra EXAMPLES:: sage: F = QuadraticField(-1, 'i') sage: x1 = F(4 - 2*i) sage: x2 = F(1 + 2*i) sage: x3 = F(-i) sage: x4 = F(6) sage: M = matrix(F,2,[[x1,x2],[x3,x4]]) sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M) [ 4 -2| 1 2] [ 2 4|-2 1] [-----+-----] [ 0 -1| 6 0] [ 1 0| 0 6] TESTS: Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size() sage: n = ZZ.random_element(n_max) sage: F = QuadraticField(-1, 'i') sage: X = random_matrix(F, n) sage: Y = random_matrix(F, n) sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X) sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y) sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y) sage: Xe*Ye == XYe True """ n = M.nrows() if M.ncols() != n: raise ValueError("the matrix 'M' must be square") field = M.base_ring() blocks = [] for z in M.list(): a = z.vector()[0] # real part, I guess b = z.vector()[1] # imag part, I guess blocks.append(matrix(field, 2, [[a,b],[-b,a]])) # We can drop the imaginaries here. return matrix.block(field.base_ring(), n, blocks) @staticmethod def real_unembed(M): """ The inverse of _embed_complex_matrix(). SETUP:: sage: from mjo.eja.eja_algebra import \ ....: ComplexMatrixEuclideanJordanAlgebra EXAMPLES:: sage: A = matrix(QQ,[ [ 1, 2, 3, 4], ....: [-2, 1, -4, 3], ....: [ 9, 10, 11, 12], ....: [-10, 9, -12, 11] ]) sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A) [ 2*i + 1 4*i + 3] [ 10*i + 9 12*i + 11] TESTS: Unembedding is the inverse of embedding:: sage: set_random_seed() sage: F = QuadraticField(-1, 'i') sage: M = random_matrix(F, 3) sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M) sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M True """ n = ZZ(M.nrows()) if M.ncols() != n: raise ValueError("the matrix 'M' must be square") if not n.mod(2).is_zero(): raise ValueError("the matrix 'M' must be a complex embedding") # 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() F = field.extension(z**2 + 1, 'i', embedding=CLF(-1).sqrt()) i = F.gen() # Go top-left to bottom-right (reading order), converting every # 2-by-2 block we see to a single complex element. elements = [] for k in xrange(n/2): for j in xrange(n/2): submat = M[2*k:2*k+2,2*j:2*j+2] if submat[0,0] != submat[1,1]: raise ValueError('bad on-diagonal submatrix') if submat[0,1] != -submat[1,0]: raise ValueError('bad off-diagonal submatrix') z = submat[0,0] + submat[0,1]*i elements.append(z) return matrix(F, n/2, elements) @classmethod def natural_inner_product(cls,X,Y): """ Compute a natural inner product in this algebra directly from its real embedding. SETUP:: sage: from mjo.eja.eja_algebra import ComplexHermitianEJA TESTS: This gives the same answer as the slow, default method implemented in :class:`MatrixEuclideanJordanAlgebra`:: sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: Xe = x.natural_representation() sage: Ye = y.natural_representation() sage: X = ComplexHermitianEJA.real_unembed(Xe) sage: Y = ComplexHermitianEJA.real_unembed(Ye) sage: expected = (X*Y).trace().vector()[0] sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye) sage: actual == expected True """ return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, and the real-part-of-trace inner product. It has dimension `n^2` over the reals. SETUP:: sage: from mjo.eja.eja_algebra import ComplexHermitianEJA TESTS: The dimension of this algebra is `n^2`:: sage: set_random_seed() sage: n_max = ComplexHermitianEJA._max_test_case_size() sage: n = ZZ.random_element(1, n_max) sage: J = ComplexHermitianEJA(n) sage: J.dimension() == n^2 True The Jordan multiplication is what we think it is:: sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).natural_representation() sage: X = x.natural_representation() sage: Y = y.natural_representation() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True sage: J(expected) == x*y True We can change the generator prefix:: sage: ComplexHermitianEJA(2, prefix='z').gens() (z0, z1, z2, z3) Our natural basis is normalized with respect to the natural inner product unless we specify otherwise:: sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: all( b.norm() == 1 for b in J.gens() ) True Since our natural basis is normalized with respect to the natural inner product, and since we know that this algebra is an EJA, any left-multiplication operator's matrix will be symmetric because natural->EJA basis representation is an isometry and within the EJA the operator is self-adjoint by the Jordan axiom:: sage: set_random_seed() sage: x = ComplexHermitianEJA.random_instance().random_element() sage: x.operator().matrix().is_symmetric() True """ @classmethod def _denormalized_basis(cls, n, field): """ Returns a basis for the space of complex Hermitian n-by-n matrices. Why do we embed these? Basically, because all of numerical linear algebra assumes that you're working with vectors consisting of `n` entries from a field and scalars from the same field. There's no way to tell SageMath that (for example) the vectors contain complex numbers, while the scalar field is real. SETUP:: sage: from mjo.eja.eja_algebra import ComplexHermitianEJA TESTS:: sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: field = QuadraticField(2, 'sqrt2') sage: B = ComplexHermitianEJA._denormalized_basis(n, field) sage: all( M.is_symmetric() for M in B) True """ R = PolynomialRing(field, 'z') z = R.gen() F = field.extension(z**2 + 1, 'I') I = F.gen() # This is like the symmetric case, but we need to be careful: # # * We want conjugate-symmetry, not just symmetry. # * The diagonal will (as a result) be real. # S = [] for i in xrange(n): for j in xrange(i+1): Eij = matrix(F, n, lambda k,l: k==i and l==j) 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()) S.append(Sij_real) Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose()) S.append(Sij_imag) # Since we embedded these, we can drop back to the "field" that we # started with instead of the complex extension "F". return ( s.change_ring(field) for s in S ) def __init__(self, n, field=QQ, **kwargs): basis = self._denormalized_basis(n,field) super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs) class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @staticmethod def real_embed(M): """ Embed the n-by-n quaternion matrix ``M`` into the space of real matrices of size 4n-by-4n by first sending each quaternion entry `z = a + bi + cj + dk` to the block-complex matrix ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into a real matrix. SETUP:: sage: from mjo.eja.eja_algebra import \ ....: QuaternionMatrixEuclideanJordanAlgebra EXAMPLES:: sage: Q = QuaternionAlgebra(QQ,-1,-1) sage: i,j,k = Q.gens() sage: x = 1 + 2*i + 3*j + 4*k sage: M = matrix(Q, 1, [[x]]) sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M) [ 1 2 3 4] [-2 1 -4 3] [-3 4 1 -2] [-4 -3 2 1] Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size() sage: n = ZZ.random_element(n_max) sage: Q = QuaternionAlgebra(QQ,-1,-1) sage: X = random_matrix(Q, n) sage: Y = random_matrix(Q, n) sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X) sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y) sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y) sage: Xe*Ye == XYe True """ quaternions = M.base_ring() n = M.nrows() if M.ncols() != n: raise ValueError("the matrix 'M' must be square") F = QuadraticField(-1, 'i') i = F.gen() blocks = [] for z in M.list(): t = z.coefficient_tuple() a = t[0] b = t[1] c = t[2] d = t[3] cplxM = matrix(F, 2, [[ a + b*i, c + d*i], [-c + d*i, a - b*i]]) realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM) blocks.append(realM) # We should have real entries by now, so use the realest field # we've got for the return value. return matrix.block(quaternions.base_ring(), n, blocks) @staticmethod def real_unembed(M): """ The inverse of _embed_quaternion_matrix(). SETUP:: sage: from mjo.eja.eja_algebra import \ ....: QuaternionMatrixEuclideanJordanAlgebra EXAMPLES:: sage: M = matrix(QQ, [[ 1, 2, 3, 4], ....: [-2, 1, -4, 3], ....: [-3, 4, 1, -2], ....: [-4, -3, 2, 1]]) sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M) [1 + 2*i + 3*j + 4*k] TESTS: Unembedding is the inverse of embedding:: sage: set_random_seed() sage: Q = QuaternionAlgebra(QQ, -1, -1) sage: M = random_matrix(Q, 3) sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M) sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M True """ n = ZZ(M.nrows()) if M.ncols() != n: raise ValueError("the matrix 'M' must be square") if not n.mod(4).is_zero(): raise ValueError("the matrix 'M' must be a quaternion embedding") # 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) i,j,k = Q.gens() # Go top-left to bottom-right (reading order), converting every # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1 # quaternion block. elements = [] for l in xrange(n/4): for m in xrange(n/4): submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed( M[4*l:4*l+4,4*m:4*m+4] ) if submat[0,0] != submat[1,1].conjugate(): raise ValueError('bad on-diagonal submatrix') if submat[0,1] != -submat[1,0].conjugate(): raise ValueError('bad off-diagonal submatrix') z = submat[0,0].vector()[0] # real part z += submat[0,0].vector()[1]*i # imag part z += submat[0,1].vector()[0]*j # real part z += submat[0,1].vector()[1]*k # imag part elements.append(z) return matrix(Q, n/4, elements) @classmethod def natural_inner_product(cls,X,Y): """ Compute a natural inner product in this algebra directly from its real embedding. SETUP:: sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA TESTS: This gives the same answer as the slow, default method implemented in :class:`MatrixEuclideanJordanAlgebra`:: sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: Xe = x.natural_representation() sage: Ye = y.natural_representation() sage: X = QuaternionHermitianEJA.real_unembed(Xe) sage: Y = QuaternionHermitianEJA.real_unembed(Ye) sage: expected = (X*Y).trace().coefficient_tuple()[0] sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye) sage: actual == expected True """ return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, KnownRankEJA): """ The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the real-part-of-trace inner product. It has dimension `2n^2 - n` over the reals. SETUP:: sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA TESTS: The dimension of this algebra is `2*n^2 - n`:: sage: set_random_seed() sage: n_max = QuaternionHermitianEJA._max_test_case_size() sage: n = ZZ.random_element(1, n_max) sage: J = QuaternionHermitianEJA(n) sage: J.dimension() == 2*(n^2) - n True The Jordan multiplication is what we think it is:: sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).natural_representation() sage: X = x.natural_representation() sage: Y = y.natural_representation() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True sage: J(expected) == x*y True We can change the generator prefix:: sage: QuaternionHermitianEJA(2, prefix='a').gens() (a0, a1, a2, a3, a4, a5) Our natural basis is normalized with respect to the natural inner product unless we specify otherwise:: sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: all( b.norm() == 1 for b in J.gens() ) True Since our natural basis is normalized with respect to the natural inner product, and since we know that this algebra is an EJA, any left-multiplication operator's matrix will be symmetric because natural->EJA basis representation is an isometry and within the EJA the operator is self-adjoint by the Jordan axiom:: sage: set_random_seed() sage: x = QuaternionHermitianEJA.random_instance().random_element() sage: x.operator().matrix().is_symmetric() True """ @classmethod def _denormalized_basis(cls, n, field): """ Returns a basis for the space of quaternion Hermitian n-by-n matrices. Why do we embed these? Basically, because all of numerical linear algebra assumes that you're working with vectors consisting of `n` entries from a field and scalars from the same field. There's no way to tell SageMath that (for example) the vectors contain complex numbers, while the scalar field is real. SETUP:: sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA TESTS:: sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ) sage: all( M.is_symmetric() for M in B ) True """ Q = QuaternionAlgebra(QQ,-1,-1) I,J,K = Q.gens() # This is like the symmetric case, but we need to be careful: # # * We want conjugate-symmetry, not just symmetry. # * The diagonal will (as a result) be real. # S = [] for i in xrange(n): for j in xrange(i+1): Eij = matrix(Q, n, lambda k,l: k==i and l==j) 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()) S.append(Sij_real) Sij_I = cls.real_embed(I*Eij - I*Eij.transpose()) S.append(Sij_I) Sij_J = cls.real_embed(J*Eij - J*Eij.transpose()) S.append(Sij_J) Sij_K = cls.real_embed(K*Eij - K*Eij.transpose()) S.append(Sij_K) # Since we embedded these, we can drop back to the "field" that we # started with instead of the quaternion algebra "Q". return ( s.change_ring(field) for s in S ) def __init__(self, n, field=QQ, **kwargs): basis = self._denormalized_basis(n,field) super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs) class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): """ The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` with the usual inner product and jordan product ``x*y = (, x0*y_bar + y0*x_bar)``. It has dimension `n` over the reals. SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA EXAMPLES: This multiplication table can be verified by hand:: sage: J = JordanSpinEJA(4) sage: e0,e1,e2,e3 = J.gens() sage: e0*e0 e0 sage: e0*e1 e1 sage: e0*e2 e2 sage: e0*e3 e3 sage: e1*e2 0 sage: e1*e3 0 sage: e2*e3 0 We can change the generator prefix:: sage: JordanSpinEJA(2, prefix='B').gens() (B0, B1) """ def __init__(self, n, field=QQ, **kwargs): V = VectorSpace(field, n) mult_table = [[V.zero() for j in xrange(n)] for i in xrange(n)] for i in xrange(n): for j in xrange(n): x = V.gen(i) y = V.gen(j) x0 = x[0] xbar = x[1:] y0 = y[0] ybar = y[1:] # z = x*y z0 = x.inner_product(y) zbar = y0*xbar + x0*ybar z = V([z0] + zbar.list()) mult_table[i][j] = z # The rank of the spin algebra is two, unless we're in a # one-dimensional ambient space (because the rank is bounded by # the ambient dimension). fdeja = super(JordanSpinEJA, self) return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs) def inner_product(self, x, y): """ Faster to reimplement than to use natural representations. SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA TESTS: Ensure that this is the usual inner product for the algebras over `R^n`:: sage: set_random_seed() sage: J = JordanSpinEJA.random_instance() sage: x,y = J.random_elements(2) sage: X = x.natural_representation() sage: Y = y.natural_representation() sage: x.inner_product(y) == J.natural_inner_product(X,Y) True """ return x.to_vector().inner_product(y.to_vector())