X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=02ed966540132c07af408a1a9fb4e6eebc220c31;hb=b05705ab04692f738a57a6ef387662ba5ea46ceb;hp=ec58b57e8d78b004bcdc0ffaa83fc773c3cfffa8;hpb=ee3262f5130f2e7b38b520d4238cede0cea9b981;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index ec58b57..02ed966 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -5,71 +5,42 @@ 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.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra -from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra -from sage.categories.finite_dimensional_algebras_with_basis import FiniteDimensionalAlgebrasWithBasis -from sage.functions.other import sqrt +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.modules.free_module import VectorSpace -from sage.modules.free_module_element import vector +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.number_field.number_field import NumberField, 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 sage.structure.element import is_Matrix -from sage.structure.category_object import normalize_names - -from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator +from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement +from mjo.eja.eja_utils import _mat2vec -class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): - @staticmethod - def __classcall_private__(cls, - field, - mult_table, - names='e', - assume_associative=False, - category=None, - rank=None, - natural_basis=None): - n = len(mult_table) - mult_table = [b.base_extend(field) for b in mult_table] - for b in mult_table: - b.set_immutable() - if not (is_Matrix(b) and b.dimensions() == (n, n)): - raise ValueError("input is not a multiplication table") - mult_table = tuple(mult_table) - - cat = FiniteDimensionalAlgebrasWithBasis(field) - cat.or_subcategory(category) - if assume_associative: - cat = cat.Associative() - - names = normalize_names(n, names) - - fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls) - return fda.__classcall__(cls, - field, - mult_table, - assume_associative=assume_associative, - names=names, - category=cat, - rank=rank, - natural_basis=natural_basis) - +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, - names='e', - assume_associative=False, + rank, + prefix='e', category=None, - rank=None, natural_basis=None): """ SETUP:: @@ -82,20 +53,104 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): sage: set_random_seed() sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: x*y == y*x True """ self._rank = rank self._natural_basis = natural_basis - self._multiplication_table = mult_table + + if category is None: + category = MagmaticAlgebras(field).FiniteDimensional() + category = category.WithBasis().Unital() + fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) fda.__init__(field, - mult_table, - names=names, + 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): @@ -111,14 +166,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): Ensure that it says what we think it says:: sage: JordanSpinEJA(2, field=QQ) - Euclidean Jordan algebra of degree 2 over Rational Field + Euclidean Jordan algebra of dimension 2 over Rational Field sage: JordanSpinEJA(3, field=RDF) - Euclidean Jordan algebra of degree 3 over Real Double Field + Euclidean Jordan algebra of dimension 3 over Real Double Field """ - fmt = "Euclidean Jordan algebra of degree {} over {}" - return fmt.format(self.degree(), self.base_ring()) + 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): """ @@ -157,12 +214,18 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): determinant). """ z = self._a_regular_element() - V = self.vector_space() - V1 = V.span_of_basis( (z**k).vector() for k in range(self.rank()) ) + # 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): """ @@ -209,14 +272,40 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): r = self.rank() n = self.dimension() - # Construct a new algebra over a multivariate polynomial ring... - names = ['X' + str(i) for i in range(1,n+1)] + # 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) - J = FiniteDimensionalEuclideanJordanAlgebra(R, - self._multiplication_table, - rank=r) - idmat = matrix.identity(J.base_ring(), n) + # 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()) @@ -236,12 +325,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): # 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 = J(W(R.gens())) - l1 = [matrix.column(W.coordinates((x**k).vector())) for k in range(r)] - l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)] - A_of_x = matrix.block(R, 1, n, (l1 + l2)) - xr = W.coordinates((x**r).vector()) - return (A_of_x, x, xr, A_of_x.det()) + 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 @@ -270,7 +357,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): 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().vector() + sage: xvec = J.one().to_vector() sage: p(*xvec) t^2 - 2*t + 1 @@ -302,7 +389,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): # assign a[r] goes out-of-bounds. a.append(1) # corresponds to x^r - return sum( a[k]*(t**k) for k in range(len(a)) ) + return sum( a[k]*(t**k) for k in xrange(len(a)) ) def inner_product(self, x, y): @@ -319,21 +406,76 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): EXAMPLES: - The inner product must satisfy its axiom for this algebra to truly - be a Euclidean Jordan Algebra:: + Our inner product is "associative," which means the following for + a symmetric bilinear form:: sage: set_random_seed() sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() + sage: x,y,z = J.random_elements(3) sage: (x*y).inner_product(z) == y.inner_product(x*z) True """ - if (not x in self) or (not y in self): - raise TypeError("arguments must live in this algebra") - return x.trace_inner_product(y) + 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): @@ -358,18 +500,18 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): sage: J = RealSymmetricEJA(2) sage: J.basis() - Family (e0, e1, e2) + Finite family {0: e0, 1: e1, 2: e2} sage: J.natural_basis() ( - [1 0] [0 1] [0 0] - [0 0], [1 0], [0 1] + [1 0] [ 0 1/2*sqrt2] [0 0] + [0 0], [1/2*sqrt2 0], [0 1] ) :: sage: J = JordanSpinEJA(2) sage: J.basis() - Family (e0, e1) + Finite family {0: e0, 1: e1} sage: J.natural_basis() ( [1] [0] @@ -378,1209 +520,253 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): """ if self._natural_basis is None: - return tuple( b.vector().column() for b in self.basis() ) + M = self.natural_basis_space() + return tuple( M(b.to_vector()) for b in self.basis() ) else: return self._natural_basis - def rank(self): + def natural_basis_space(self): """ - Return the rank of this EJA. + Return the matrix space in which this algebra's natural basis + elements live. """ - if self._rank is None: - raise ValueError("no rank specified at genesis") + if self._natural_basis is None or len(self._natural_basis) == 0: + return MatrixSpace(self.base_ring(), self.dimension(), 1) else: - return self._rank + return self._natural_basis[0].matrix_space() - def vector_space(self): + @staticmethod + def natural_inner_product(X,Y): """ - 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 Rational Field + 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 self.zero().vector().parent().ambient_vector_space() + return (X.conjugate_transpose()*Y).trace() - class Element(FiniteDimensionalAlgebraElement): - """ - An element of a Euclidean Jordan algebra. + @cached_method + def one(self): """ + Return the unit element of this algebra. - def __dir__(self): - """ - Oh man, I should not be doing this. This hides the "disabled" - methods ``left_matrix`` and ``matrix`` from introspection; - in particular it removes them from tab-completion. - """ - return filter(lambda s: s not in ['left_matrix', 'matrix'], - dir(self.__class__) ) - - - def __init__(self, A, elt=None): - """ - - SETUP:: - - sage: from mjo.eja.eja_algebra import 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 - - """ - # Goal: if we're given a matrix, and if it lives in our - # parent algebra's "natural ambient space," convert it - # into an algebra element. - # - # The catch is, we make a recursive call after converting - # the given matrix into a vector that lives in the algebra. - # This we need to try the parent class initializer first, - # to avoid recursing forever if we're given something that - # already fits into the algebra, but also happens to live - # in the parent's "natural ambient space" (this happens with - # vectors in R^n). - try: - FiniteDimensionalAlgebraElement.__init__(self, A, elt) - except ValueError: - natural_basis = A.natural_basis() - if elt in natural_basis[0].matrix_space(): - # Thanks for nothing! Matrix spaces aren't vector - # spaces in Sage, so we have to figure out its - # natural-basis coordinates ourselves. - V = VectorSpace(elt.base_ring(), elt.nrows()**2) - W = V.span( _mat2vec(s) for s in natural_basis ) - coords = W.coordinates(_mat2vec(elt)) - FiniteDimensionalAlgebraElement.__init__(self, A, coords) - - def __pow__(self, n): - """ - Return ``self`` raised to the power ``n``. - - Jordan algebras are always power-associative; see for - example Faraut and Koranyi, Proposition II.1.2 (ii). - - .. WARNING: - - We have to override this because our superclass uses row vectors - instead of column vectors! We, on the other hand, assume column - vectors everywhere. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - EXAMPLES:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.operator()(x) == (x^2) - True - - A few examples of power-associativity:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x*(x*x)*(x*x) == x^5 - True - sage: (x*x)*(x*x*x) == x^5 - True - - We also know that powers operator-commute (Koecher, Chapter - III, Corollary 1):: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: m = ZZ.random_element(0,10) - sage: n = ZZ.random_element(0,10) - sage: Lxm = (x^m).operator() - sage: Lxn = (x^n).operator() - sage: Lxm*Lxn == Lxn*Lxm - True - - """ - if n == 0: - return self.parent().one() - elif n == 1: - return self - else: - return (self.operator()**(n-1))(self) - - - def apply_univariate_polynomial(self, p): - """ - Apply the univariate polynomial ``p`` to this element. - - A priori, SageMath won't allow us to apply a univariate - polynomial to an element of an EJA, because we don't know - that EJAs are rings (they are usually not associative). Of - course, we know that EJAs are power-associative, so the - operation is ultimately kosher. This function sidesteps - the CAS to get the answer we want and expect. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA, - ....: random_eja) - - EXAMPLES:: - - sage: R = PolynomialRing(QQ, 't') - sage: t = R.gen(0) - sage: p = t^4 - t^3 + 5*t - 2 - sage: J = RealCartesianProductEJA(5) - sage: J.one().apply_univariate_polynomial(p) == 3*J.one() - True - - TESTS: - - We should always get back an element of the algebra:: - - sage: set_random_seed() - sage: p = PolynomialRing(QQ, 't').random_element() - sage: J = random_eja() - sage: x = J.random_element() - sage: x.apply_univariate_polynomial(p) in J - True - - """ - if len(p.variables()) > 1: - raise ValueError("not a univariate polynomial") - P = self.parent() - R = P.base_ring() - # Convert the coeficcients to the parent's base ring, - # because a priori they might live in an (unnecessarily) - # larger ring for which P.sum() would fail below. - cs = [ R(c) for c in p.coefficients(sparse=False) ] - return P.sum( cs[k]*(self**k) for k in range(len(cs)) ) - - - def characteristic_polynomial(self): - """ - Return the characteristic polynomial of this element. - - SETUP:: - - sage: from mjo.eja.eja_algebra import RealCartesianProductEJA - - EXAMPLES: - - The rank of `R^3` is three, and the minimal polynomial of - the identity element is `(t-1)` from which it follows that - the characteristic polynomial should be `(t-1)^3`:: - - sage: J = RealCartesianProductEJA(3) - sage: J.one().characteristic_polynomial() - t^3 - 3*t^2 + 3*t - 1 - - Likewise, the characteristic of the zero element in the - rank-three algebra `R^{n}` should be `t^{3}`:: - - sage: J = RealCartesianProductEJA(3) - sage: J.zero().characteristic_polynomial() - t^3 - - The characteristic polynomial of an element should evaluate - to zero on that element:: - - sage: set_random_seed() - sage: x = RealCartesianProductEJA(3).random_element() - sage: p = x.characteristic_polynomial() - sage: x.apply_univariate_polynomial(p) - 0 - - """ - p = self.parent().characteristic_polynomial() - return p(*self.vector()) - - - def inner_product(self, other): - """ - Return the parent algebra's inner product of myself and ``other``. - - SETUP:: - - sage: from mjo.eja.eja_algebra import ( - ....: ComplexHermitianEJA, - ....: JordanSpinEJA, - ....: QuaternionHermitianEJA, - ....: RealSymmetricEJA, - ....: random_eja) - - EXAMPLES: - - The inner product in the Jordan spin algebra is the usual - inner product on `R^n` (this example only works because the - basis for the Jordan algebra is the standard basis in `R^n`):: - - sage: J = JordanSpinEJA(3) - sage: x = vector(QQ,[1,2,3]) - sage: y = vector(QQ,[4,5,6]) - sage: x.inner_product(y) - 32 - sage: J(x).inner_product(J(y)) - 32 - - The inner product on `S^n` is ` = trace(X*Y)`, where - multiplication is the usual matrix multiplication in `S^n`, - so the inner product of the identity matrix with itself - should be the `n`:: - - sage: J = RealSymmetricEJA(3) - sage: J.one().inner_product(J.one()) - 3 - - Likewise, the inner product on `C^n` is ` = - Re(trace(X*Y))`, where we must necessarily take the real - part because the product of Hermitian matrices may not be - Hermitian:: - - sage: J = ComplexHermitianEJA(3) - sage: J.one().inner_product(J.one()) - 3 - - Ditto for the quaternions:: - - sage: J = QuaternionHermitianEJA(3) - sage: J.one().inner_product(J.one()) - 3 - - TESTS: - - Ensure that we can always compute an inner product, and that - it gives us back a real number:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: x.inner_product(y) in RR - True - - """ - P = self.parent() - if not other in P: - raise TypeError("'other' must live in the same algebra") - - return P.inner_product(self, other) - - - def operator_commutes_with(self, other): - """ - Return whether or not this element operator-commutes - with ``other``. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - EXAMPLES: - - The definition of a Jordan algebra says that any element - operator-commutes with its square:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.operator_commutes_with(x^2) - True - - TESTS: - - Test Lemma 1 from Chapter III of Koecher:: - - sage: set_random_seed() - sage: J = random_eja() - sage: u = J.random_element() - sage: v = J.random_element() - sage: lhs = u.operator_commutes_with(u*v) - sage: rhs = v.operator_commutes_with(u^2) - sage: lhs == rhs - True - - Test the first polarization identity from my notes, Koecher Chapter - III, or from Baes (2.3):: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: Lx = x.operator() - sage: Ly = y.operator() - sage: Lxx = (x*x).operator() - sage: Lxy = (x*y).operator() - sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly) - True - - Test the second polarization identity from my notes or from - Baes (2.4):: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: Lx = x.operator() - sage: Ly = y.operator() - sage: Lz = z.operator() - sage: Lzy = (z*y).operator() - sage: Lxy = (x*y).operator() - sage: Lxz = (x*z).operator() - sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly) - True - - Test the third polarization identity from my notes or from - Baes (2.5):: - - sage: set_random_seed() - sage: J = random_eja() - sage: u = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: Lu = u.operator() - sage: Ly = y.operator() - sage: Lz = z.operator() - sage: Lzy = (z*y).operator() - sage: Luy = (u*y).operator() - sage: Luz = (u*z).operator() - sage: Luyz = (u*(y*z)).operator() - sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz - sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly - sage: bool(lhs == rhs) - True - - """ - if not other in self.parent(): - raise TypeError("'other' must live in the same algebra") - - A = self.operator() - B = other.operator() - return (A*B == B*A) - - - def det(self): - """ - Return my determinant, the product of my eigenvalues. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: random_eja) - - EXAMPLES:: - - sage: J = JordanSpinEJA(2) - sage: e0,e1 = J.gens() - sage: x = sum( J.gens() ) - sage: x.det() - 0 - - :: - - sage: J = JordanSpinEJA(3) - sage: e0,e1,e2 = J.gens() - sage: x = sum( J.gens() ) - sage: x.det() - -1 - - TESTS: - - An element is invertible if and only if its determinant is - non-zero:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.is_invertible() == (x.det() != 0) - True - - """ - P = self.parent() - r = P.rank() - p = P._charpoly_coeff(0) - # The _charpoly_coeff function already adds the factor of - # -1 to ensure that _charpoly_coeff(0) is really what - # appears in front of t^{0} in the charpoly. However, - # we want (-1)^r times THAT for the determinant. - return ((-1)**r)*p(*self.vector()) - - - def inverse(self): - """ - Return the Jordan-multiplicative inverse of this element. - - ALGORITHM: - - We appeal to the quadratic representation as in Koecher's - Theorem 12 in Chapter III, Section 5. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: random_eja) - - EXAMPLES: - - The inverse in the spin factor algebra is given in Alizadeh's - Example 11.11:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,10) - sage: J = JordanSpinEJA(n) - sage: x = J.random_element() - sage: while not x.is_invertible(): - ....: x = J.random_element() - sage: x_vec = x.vector() - sage: x0 = x_vec[0] - sage: x_bar = x_vec[1:] - sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar)) - sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list()) - sage: x_inverse = coeff*inv_vec - sage: x.inverse() == J(x_inverse) - True - - TESTS: - - The identity element is its own inverse:: - - sage: set_random_seed() - sage: J = random_eja() - sage: J.one().inverse() == J.one() - True - - If an element has an inverse, it acts like one:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: (not x.is_invertible()) or (x.inverse()*x == J.one()) - True - - The inverse of the inverse is what we started with:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: (not x.is_invertible()) or (x.inverse().inverse() == x) - True - - The zero element is never invertible:: - - sage: set_random_seed() - sage: J = random_eja().zero().inverse() - Traceback (most recent call last): - ... - ValueError: element is not invertible + SETUP:: - """ - if not self.is_invertible(): - raise ValueError("element is not invertible") + sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA, + ....: random_eja) - return (~self.quadratic_representation())(self) + EXAMPLES:: + sage: J = RealCartesianProductEJA(5) + sage: J.one() + e0 + e1 + e2 + e3 + e4 - def is_invertible(self): - """ - Return whether or not this element is invertible. + TESTS: - We can't use the superclass method because it relies on - the algebra being associative. + The identity element acts like the identity:: - ALGORITHM: + 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 usual way to do this is to check if the determinant is - zero, but we need the characteristic polynomial for the - determinant. The minimal polynomial is a lot easier to get, - so we use Corollary 2 in Chapter V of Koecher to check - whether or not the paren't algebra's zero element is a root - of this element's minimal polynomial. + The matrix of the unit element's operator is the identity:: - SETUP:: + 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 - sage: from mjo.eja.eja_algebra import random_eja + """ + # 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() ] - TESTS: + # 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) - The identity element is always invertible:: + # 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()) ) - sage: set_random_seed() - sage: J = random_eja() - sage: J.one().is_invertible() - True + # 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)) - The zero element is never invertible:: - sage: set_random_seed() - sage: J = random_eja() - sage: J.zero().is_invertible() - False + 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() - """ - zero = self.parent().zero() - p = self.minimal_polynomial() - return not (p(zero) == zero) + def random_elements(self, count): + """ + Return ``count`` random elements as a tuple. + SETUP:: - def is_nilpotent(self): - """ - Return whether or not some power of this element is zero. + sage: from mjo.eja.eja_algebra import JordanSpinEJA - The superclass method won't work unless we're in an - associative algebra, and we aren't. However, we generate - an assocoative subalgebra and we're nilpotent there if and - only if we're nilpotent here (probably). + EXAMPLES:: - SETUP:: + 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 - sage: from mjo.eja.eja_algebra import random_eja + """ + return tuple( self.random_element() for idx in xrange(count) ) - TESTS: - The identity element is never nilpotent:: + def rank(self): + """ + Return the rank of this EJA. - sage: set_random_seed() - sage: random_eja().one().is_nilpotent() - False + ALGORITHM: - The additive identity is always nilpotent:: + 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. - sage: set_random_seed() - sage: random_eja().zero().is_nilpotent() - True + SETUP:: - """ - # The element we're going to call "is_nilpotent()" on. - # Either myself, interpreted as an element of a finite- - # dimensional algebra, or an element of an associative - # subalgebra. - elt = None + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA, + ....: random_eja) - if self.parent().is_associative(): - elt = FiniteDimensionalAlgebraElement(self.parent(), self) - else: - V = self.span_of_powers() - assoc_subalg = self.subalgebra_generated_by() - # Mis-design warning: the basis used for span_of_powers() - # and subalgebra_generated_by() must be the same, and in - # the same order! - elt = assoc_subalg(V.coordinates(self.vector())) + EXAMPLES: - # Recursive call, but should work since elt lives in an - # associative algebra. - return elt.is_nilpotent() + 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 - def is_regular(self): - """ - Return whether or not this is a regular element. + The rank of the `n`-by-`n` Hermitian real, complex, or + quaternion matrices is `n`:: - SETUP:: + sage: RealSymmetricEJA(4).rank() + 4 + sage: ComplexHermitianEJA(3).rank() + 3 + sage: QuaternionHermitianEJA(2).rank() + 2 - sage: from mjo.eja.eja_algebra import JordanSpinEJA + TESTS: - EXAMPLES: + Ensure that every EJA that we know how to construct has a + positive integer rank:: - The identity element always has degree one, but any element - linearly-independent from it is regular:: + sage: set_random_seed() + sage: r = random_eja().rank() + sage: r in ZZ and r > 0 + True - sage: J = JordanSpinEJA(5) - sage: J.one().is_regular() - False - sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity - sage: for x in J.gens(): - ....: (J.one() + x).is_regular() - False - True - True - True - True + """ + return self._rank - """ - return self.degree() == self.parent().rank() + def vector_space(self): + """ + Return the vector space that underlies this algebra. - def degree(self): - """ - Compute the degree of this element the straightforward way - according to the definition; by appending powers to a list - and figuring out its dimension (that is, whether or not - they're linearly dependent). + SETUP:: - SETUP:: + sage: from mjo.eja.eja_algebra import RealSymmetricEJA - sage: from mjo.eja.eja_algebra import JordanSpinEJA + EXAMPLES:: - EXAMPLES:: + sage: J = RealSymmetricEJA(2) + sage: J.vector_space() + Vector space of dimension 3 over... - sage: J = JordanSpinEJA(4) - sage: J.one().degree() - 1 - sage: e0,e1,e2,e3 = J.gens() - sage: (e0 - e1).degree() - 2 + """ + return self.zero().to_vector().parent().ambient_vector_space() - In the spin factor algebra (of rank two), all elements that - aren't multiples of the identity are regular:: - sage: set_random_seed() - sage: n = ZZ.random_element(1,10) - sage: J = JordanSpinEJA(n) - sage: x = J.random_element() - sage: x == x.coefficient(0)*J.one() or x.degree() == 2 - True + Element = FiniteDimensionalEuclideanJordanAlgebraElement - """ - return self.span_of_powers().dimension() +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. - def left_matrix(self): - """ - Our parent class defines ``left_matrix`` and ``matrix`` - methods whose names are misleading. We don't want them. - """ - raise NotImplementedError("use operator().matrix() instead") - - matrix = left_matrix - - - def minimal_polynomial(self): - """ - Return the minimal polynomial of this element, - as a function of the variable `t`. - - ALGORITHM: - - We restrict ourselves to the associative subalgebra - generated by this element, and then return the minimal - polynomial of this element's operator matrix (in that - subalgebra). This works by Baes Proposition 2.3.16. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: random_eja) - - TESTS: - - The minimal polynomial of the identity and zero elements are - always the same:: - - sage: set_random_seed() - sage: J = random_eja() - sage: J.one().minimal_polynomial() - t - 1 - sage: J.zero().minimal_polynomial() - t - - The degree of an element is (by one definition) the degree - of its minimal polynomial:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.degree() == x.minimal_polynomial().degree() - True - - The minimal polynomial and the characteristic polynomial coincide - and are known (see Alizadeh, Example 11.11) for all elements of - the spin factor algebra that aren't scalar multiples of the - identity:: - - sage: set_random_seed() - sage: n = ZZ.random_element(2,10) - sage: J = JordanSpinEJA(n) - sage: y = J.random_element() - sage: while y == y.coefficient(0)*J.one(): - ....: y = J.random_element() - sage: y0 = y.vector()[0] - sage: y_bar = y.vector()[1:] - sage: actual = y.minimal_polynomial() - sage: t = PolynomialRing(J.base_ring(),'t').gen(0) - sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2) - sage: bool(actual == expected) - True - - The minimal polynomial should always kill its element:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: p = x.minimal_polynomial() - sage: x.apply_univariate_polynomial(p) - 0 - - """ - V = self.span_of_powers() - assoc_subalg = self.subalgebra_generated_by() - # Mis-design warning: the basis used for span_of_powers() - # and subalgebra_generated_by() must be the same, and in - # the same order! - elt = assoc_subalg(V.coordinates(self.vector())) - return elt.operator().minimal_polynomial() - - - - def natural_representation(self): - """ - Return a more-natural representation of this element. - - 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" representation of this element as a Hermitian - matrix, if it has one. If not, you get the usual representation. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, - ....: QuaternionHermitianEJA) - - EXAMPLES:: - - sage: J = ComplexHermitianEJA(3) - sage: J.one() - e0 + e5 + e8 - sage: J.one().natural_representation() - [1 0 0 0 0 0] - [0 1 0 0 0 0] - [0 0 1 0 0 0] - [0 0 0 1 0 0] - [0 0 0 0 1 0] - [0 0 0 0 0 1] - - :: - - sage: J = QuaternionHermitianEJA(3) - sage: J.one() - e0 + e9 + e14 - sage: J.one().natural_representation() - [1 0 0 0 0 0 0 0 0 0 0 0] - [0 1 0 0 0 0 0 0 0 0 0 0] - [0 0 1 0 0 0 0 0 0 0 0 0] - [0 0 0 1 0 0 0 0 0 0 0 0] - [0 0 0 0 1 0 0 0 0 0 0 0] - [0 0 0 0 0 1 0 0 0 0 0 0] - [0 0 0 0 0 0 1 0 0 0 0 0] - [0 0 0 0 0 0 0 1 0 0 0 0] - [0 0 0 0 0 0 0 0 1 0 0 0] - [0 0 0 0 0 0 0 0 0 1 0 0] - [0 0 0 0 0 0 0 0 0 0 1 0] - [0 0 0 0 0 0 0 0 0 0 0 1] - - """ - B = self.parent().natural_basis() - W = B[0].matrix_space() - return W.linear_combination(zip(self.vector(), B)) - - - def operator(self): - """ - Return the left-multiplication-by-this-element - operator on the ambient algebra. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: x.operator()(y) == x*y - True - sage: y.operator()(x) == x*y - True - - """ - P = self.parent() - fda_elt = FiniteDimensionalAlgebraElement(P, self) - return FiniteDimensionalEuclideanJordanAlgebraOperator( - P, - P, - fda_elt.matrix().transpose() ) - - - def quadratic_representation(self, other=None): - """ - Return the quadratic representation of this element. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: random_eja) - - EXAMPLES: - - The explicit form in the spin factor algebra is given by - Alizadeh's Example 11.12:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,10) - sage: J = JordanSpinEJA(n) - sage: x = J.random_element() - sage: x_vec = x.vector() - sage: x0 = x_vec[0] - sage: x_bar = x_vec[1:] - sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)]) - sage: B = 2*x0*x_bar.row() - sage: C = 2*x0*x_bar.column() - sage: D = matrix.identity(QQ, n-1) - sage: D = (x0^2 - x_bar.inner_product(x_bar))*D - sage: D = D + 2*x_bar.tensor_product(x_bar) - sage: Q = matrix.block(2,2,[A,B,C,D]) - sage: Q == x.quadratic_representation().matrix() - True - - Test all of the properties from Theorem 11.2 in Alizadeh:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: Lx = x.operator() - sage: Lxx = (x*x).operator() - sage: Qx = x.quadratic_representation() - sage: Qy = y.quadratic_representation() - sage: Qxy = x.quadratic_representation(y) - sage: Qex = J.one().quadratic_representation(x) - sage: n = ZZ.random_element(10) - sage: Qxn = (x^n).quadratic_representation() - - Property 1: - - sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy - True - - Property 2 (multiply on the right for :trac:`28272`): - - sage: alpha = QQ.random_element() - sage: (alpha*x).quadratic_representation() == Qx*(alpha^2) - True - - Property 3: - - sage: not x.is_invertible() or ( Qx(x.inverse()) == x ) - True + Not really a subclass of FDEJA because doing that causes method + resolution errors, e.g. - sage: not x.is_invertible() or ( - ....: ~Qx - ....: == - ....: x.inverse().quadratic_representation() ) - True + TypeError: Error when calling the metaclass bases + Cannot create a consistent method resolution + order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra, + KnownRankEJA - sage: Qxy(J.one()) == x*y - True - - Property 4: - - sage: not x.is_invertible() or ( - ....: x.quadratic_representation(x.inverse())*Qx - ....: == Qx*x.quadratic_representation(x.inverse()) ) - True - - sage: not x.is_invertible() or ( - ....: x.quadratic_representation(x.inverse())*Qx - ....: == - ....: 2*x.operator()*Qex - Qx ) - True - - sage: 2*x.operator()*Qex - Qx == Lxx - True - - Property 5: - - sage: Qy(x).quadratic_representation() == Qy*Qx*Qy - True - - Property 6: - - sage: Qxn == (Qx)^n - True - - Property 7: - - sage: not x.is_invertible() or ( - ....: Qx*x.inverse().operator() == Lx ) - True - - Property 8: - - sage: not x.operator_commutes_with(y) or ( - ....: Qx(y)^n == Qxn(y^n) ) - True - - """ - if other is None: - other=self - elif not other in self.parent(): - raise TypeError("'other' must live in the same algebra") - - L = self.operator() - M = other.operator() - return ( L*M + M*L - (self*other).operator() ) - - - def span_of_powers(self): - """ - Return the vector space spanned by successive powers of - this element. - """ - # The dimension of the subalgebra can't be greater than - # the big algebra, so just put everything into a list - # and let span() get rid of the excess. - # - # We do the extra ambient_vector_space() in case we're messing - # with polynomials and the direct parent is a module. - V = self.parent().vector_space() - return V.span( (self**d).vector() for d in xrange(V.dimension()) ) - - - def subalgebra_generated_by(self): - """ - Return the associative subalgebra of the parent EJA generated - by this element. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: x.subalgebra_generated_by().is_associative() - True - - Squaring in the subalgebra should work the same as in - the superalgebra:: - - sage: set_random_seed() - sage: x = random_eja().random_element() - sage: u = x.subalgebra_generated_by().random_element() - sage: u.operator()(u) == u^2 - True - - """ - # First get the subspace spanned by the powers of myself... - V = self.span_of_powers() - F = self.base_ring() - - # Now figure out the entries of the right-multiplication - # matrix for the successive basis elements b0, b1,... of - # that subspace. - mats = [] - for b_right in V.basis(): - eja_b_right = self.parent()(b_right) - b_right_rows = [] - # The first row of the right-multiplication matrix by - # b1 is what we get if we apply that matrix to b1. The - # second row of the right multiplication matrix by b1 - # is what we get when we apply that matrix to b2... - # - # IMPORTANT: this assumes that all vectors are COLUMN - # vectors, unlike our superclass (which uses row vectors). - for b_left in V.basis(): - eja_b_left = self.parent()(b_left) - # Multiply in the original EJA, but then get the - # coordinates from the subalgebra in terms of its - # basis. - this_row = V.coordinates((eja_b_left*eja_b_right).vector()) - b_right_rows.append(this_row) - b_right_matrix = matrix(F, b_right_rows) - mats.append(b_right_matrix) - - # It's an algebra of polynomials in one element, and EJAs - # are power-associative. - # - # TODO: choose generator names intelligently. - return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f') - - - def subalgebra_idempotent(self): - """ - Find an idempotent in the associative subalgebra I generate - using Proposition 2.3.5 in Baes. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: while x.is_nilpotent(): - ....: x = J.random_element() - sage: c = x.subalgebra_idempotent() - sage: c^2 == c - True - - """ - if self.is_nilpotent(): - raise ValueError("this only works with non-nilpotent elements!") - - V = self.span_of_powers() - J = self.subalgebra_generated_by() - # Mis-design warning: the basis used for span_of_powers() - # and subalgebra_generated_by() must be the same, and in - # the same order! - u = J(V.coordinates(self.vector())) - - # The image of the matrix of left-u^m-multiplication - # will be minimal for some natural number s... - s = 0 - minimal_dim = V.dimension() - for i in xrange(1, V.dimension()): - this_dim = (u**i).operator().matrix().image().dimension() - if this_dim < minimal_dim: - minimal_dim = this_dim - s = i - - # Now minimal_matrix should correspond to the smallest - # non-zero subspace in Baes's (or really, Koecher's) - # proposition. - # - # However, we need to restrict the matrix to work on the - # subspace... or do we? Can't we just solve, knowing that - # A(c) = u^(s+1) should have a solution in the big space, - # too? - # - # Beware, solve_right() means that we're using COLUMN vectors. - # Our FiniteDimensionalAlgebraElement superclass uses rows. - u_next = u**(s+1) - A = u_next.operator().matrix() - c_coordinates = A.solve_right(u_next.vector()) - - # Now c_coordinates is the idempotent we want, but it's in - # the coordinate system of the subalgebra. - # - # We need the basis for J, but as elements of the parent algebra. - # - basis = [self.parent(v) for v in V.basis()] - return self.parent().linear_combination(zip(c_coordinates, basis)) - - - def trace(self): - """ - Return my trace, the sum of my eigenvalues. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: RealCartesianProductEJA, - ....: random_eja) - - EXAMPLES:: - - sage: J = JordanSpinEJA(3) - sage: x = sum(J.gens()) - sage: x.trace() - 2 - - :: - - sage: J = RealCartesianProductEJA(5) - sage: J.one().trace() - 5 - - TESTS: - - The trace of an element is a real number:: - - sage: set_random_seed() - sage: J = random_eja() - sage: J.random_element().trace() in J.base_ring() - True - - """ - P = self.parent() - r = P.rank() - p = P._charpoly_coeff(r-1) - # The _charpoly_coeff function already adds the factor of - # -1 to ensure that _charpoly_coeff(r-1) is really what - # appears in front of t^{r-1} in the charpoly. However, - # we want the negative of THAT for the trace. - return -p(*self.vector()) - - - def trace_inner_product(self, other): - """ - Return the trace inner product of myself and ``other``. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS: - - The trace inner product is commutative:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element(); y = J.random_element() - sage: x.trace_inner_product(y) == y.trace_inner_product(x) - True - - The trace inner product is bilinear:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: a = QQ.random_element(); - sage: actual = (a*(x+z)).trace_inner_product(y) - sage: expected = ( a*x.trace_inner_product(y) + - ....: a*z.trace_inner_product(y) ) - sage: actual == expected - True - sage: actual = x.trace_inner_product(a*(y+z)) - sage: expected = ( a*x.trace_inner_product(y) + - ....: a*x.trace_inner_product(z) ) - sage: actual == expected - True - - The trace inner product satisfies the compatibility - condition in the definition of a Euclidean Jordan algebra:: - - sage: set_random_seed() - sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() - sage: z = J.random_element() - sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z) - True + """ + @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 - """ - if not other in self.parent(): - raise TypeError("'other' must live in the same algebra") + @classmethod + def random_instance(cls, field=QQ, **kwargs): + """ + Return a random instance of this type of algebra. - return (self*other).trace() + 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): +class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra, + KnownRankEJA): """ Return the Euclidean Jordan Algebra corresponding to the set `R^n` under the Hadamard product. @@ -1612,22 +798,45 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: e2*e2 e2 + TESTS: + + We can change the generator prefix:: + + sage: RealCartesianProductEJA(3, prefix='r').gens() + (r0, r1, r2) + """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - # The FiniteDimensionalAlgebra constructor takes a list of - # matrices, the ith representing right multiplication by the ith - # basis element in the vector space. So if e_1 = (1,0,0), then - # right (Hadamard) multiplication of x by e_1 picks out the first - # component of x; and likewise for the ith basis element e_i. - Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i)) - for i in xrange(n) ] - - fdeja = super(RealCartesianProductEJA, cls) - return fdeja.__classcall_private__(cls, field, Qs, rank=n) + 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): - return _usual_ip(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(): @@ -1663,487 +872,444 @@ def random_eja(): TESTS:: sage: random_eja() - Euclidean Jordan algebra of degree... + Euclidean Jordan algebra of dimension... """ + classname = choice(KnownRankEJA.__subclasses__()) + return classname.random_instance() - # The max_n component lets us choose different upper bounds on the - # value "n" that gets passed to the constructor. This is needed - # because e.g. R^{10} is reasonable to test, while the Hermitian - # 10-by-10 quaternion matrices are not. - (constructor, max_n) = choice([(RealCartesianProductEJA, 6), - (JordanSpinEJA, 6), - (RealSymmetricEJA, 5), - (ComplexHermitianEJA, 4), - (QuaternionHermitianEJA, 3)]) - n = ZZ.random_element(1, max_n) - return constructor(n, field=QQ) -def _real_symmetric_basis(n, field=QQ): - """ - Return a basis for the space of real symmetric n-by-n matrices. - """ - # 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: - # Beware, orthogonal but not normalized! - Sij = Eij + Eij.transpose() - S.append(Sij) - return tuple(S) - - -def _complex_hermitian_basis(n, field=QQ): - """ - Returns a basis for the space of complex Hermitian n-by-n matrices. - SETUP:: - sage: from mjo.eja.eja_algebra import _complex_hermitian_basis +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 - TESTS:: + 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 = NumberField(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) - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) ) - True - """ - F = QuadraticField(-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(field, n, lambda k,l: k==i and l==j) - if i == j: - Sij = _embed_complex_matrix(Eij) - S.append(Sij) - else: - # Beware, orthogonal but not normalized! The second one - # has a minus because it's conjugated. - Sij_real = _embed_complex_matrix(Eij + Eij.transpose()) - S.append(Sij_real) - Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose()) - S.append(Sij_imag) - return tuple(S) + @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) ) + field = self.base_ring().base_ring() # yeeeaahhhhhhh + J = MatrixEuclideanJordanAlgebra(field, + 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 } + return p.subs(substitutions) -def _quaternion_hermitian_basis(n, field=QQ): - """ - Returns a basis for the space of quaternion Hermitian n-by-n matrices. + @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)) - SETUP:: + return mult_table - sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis - TESTS:: + @staticmethod + def real_embed(M): + """ + Embed the matrix ``M`` into a space of real matrices. - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) ) - True + 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): """ - 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 = _embed_quaternion_matrix(Eij) - S.append(Sij) - else: - # Beware, orthogonal but not normalized! The second, - # third, and fourth ones have a minus because they're - # conjugated. - Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose()) - S.append(Sij_real) - Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose()) - S.append(Sij_I) - Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose()) - S.append(Sij_J) - Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose()) - S.append(Sij_K) - return tuple(S) - - -def _mat2vec(m): - return vector(m.base_ring(), m.list()) - -def _vec2mat(v): - return matrix(v.base_ring(), sqrt(v.degree()), v.list()) - -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. A reordered copy - of the basis is also returned to work around the fact that - the ``span()`` in this function will change the order of the basis - from what we think it is, to... something else. - """ - # 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( _mat2vec(s) for s in basis ) - - # Taking the span above reorders our basis (thanks, jerk!) so we - # need to put our "matrix basis" in the same order as the - # (reordered) vector basis. - S = tuple( _vec2mat(b) for b in W.basis() ) - - Qs = [] - for s in S: - # Brute force the multiplication-by-s matrix by looping - # through all elements of the basis and doing the computation - # to find out what the corresponding row should be. BEWARE: - # these multiplication tables won't be symmetric! It therefore - # becomes REALLY IMPORTANT that the underlying algebra - # constructor uses ROW vectors and not COLUMN vectors. That's - # why we're computing rows here and not columns. - Q_rows = [] - for t in S: - this_row = _mat2vec((s*t + t*s)/2) - Q_rows.append(W.coordinates(this_row)) - Q = matrix(field, W.dimension(), Q_rows) - Qs.append(Q) - - return (Qs, S) - - -def _embed_complex_matrix(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]]``. + 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 _embed_complex_matrix + sage: from mjo.eja.eja_algebra import RealSymmetricEJA 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: _embed_complex_matrix(M) - [ 4 -2| 1 2] - [ 2 4|-2 1] - [-----+-----] - [ 0 -1| 6 0] - [ 1 0| 0 6] + 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 TESTS: - Embedding is a homomorphism (isomorphism, in fact):: + The dimension of this algebra is `(n^2 + n) / 2`:: sage: set_random_seed() - sage: n = ZZ.random_element(5) - sage: F = QuadraticField(-1, 'i') - sage: X = random_matrix(F, n) - sage: Y = random_matrix(F, n) - sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y) - sage: expected = _embed_complex_matrix(X*Y) + 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 - """ - 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.real() - b = z.imag() - blocks.append(matrix(field, 2, [[a,b],[-b,a]])) + 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:: - # We can drop the imaginaries here. - return matrix.block(field.base_ring(), n, blocks) + 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 -def _unembed_complex_matrix(M): """ - The inverse of _embed_complex_matrix(). + @classmethod + def _denormalized_basis(cls, n, field): + """ + Return a basis for the space of real symmetric n-by-n matrices. - SETUP:: + SETUP:: - sage: from mjo.eja.eja_algebra import (_embed_complex_matrix, - ....: _unembed_complex_matrix) + sage: from mjo.eja.eja_algebra import RealSymmetricEJA - EXAMPLES:: + TESTS:: - sage: A = matrix(QQ,[ [ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [ 9, 10, 11, 12], - ....: [-10, 9, -12, 11] ]) - sage: _unembed_complex_matrix(A) - [ 2*i + 1 4*i + 3] - [ 10*i + 9 12*i + 11] + 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 - TESTS: + """ + # 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 - Unembedding is the inverse of embedding:: - sage: set_random_seed() - sage: F = QuadraticField(-1, 'i') - sage: M = random_matrix(F, 3) - sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M - True + @staticmethod + def _max_test_case_size(): + return 4 # Dimension 10 - """ - 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") - - F = QuadraticField(-1, 'i') - 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) - - -def _embed_quaternion_matrix(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:: + def __init__(self, n, field=QQ, **kwargs): + basis = self._denormalized_basis(n, field) + super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs) - sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix - EXAMPLES:: +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]]``. - 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: _embed_quaternion_matrix(M) - [ 1 2 3 4] - [-2 1 -4 3] - [-3 4 1 -2] - [-4 -3 2 1] + SETUP:: - Embedding is a homomorphism (isomorphism, in fact):: + sage: from mjo.eja.eja_algebra import \ + ....: ComplexMatrixEuclideanJordanAlgebra - sage: set_random_seed() - sage: n = ZZ.random_element(5) - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: X = random_matrix(Q, n) - sage: Y = random_matrix(Q, n) - sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y) - sage: expected = _embed_quaternion_matrix(X*Y) - sage: actual == expected - True + EXAMPLES:: - """ - 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] - cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i], - [-c + d*i, a - b*i]]) - blocks.append(_embed_complex_matrix(cplx_matrix)) - - # 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) - - -def _unembed_quaternion_matrix(M): - """ - The inverse of _embed_quaternion_matrix(). + 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] - SETUP:: + TESTS: - sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix, - ....: _unembed_quaternion_matrix) + Embedding is a homomorphism (isomorphism, in fact):: - EXAMPLES:: + 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 - sage: M = matrix(QQ, [[ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [-3, 4, 1, -2], - ....: [-4, -3, 2, 1]]) - sage: _unembed_quaternion_matrix(M) - [1 + 2*i + 3*j + 4*k] + """ + 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]])) - TESTS: + # We can drop the imaginaries here. + return matrix.block(field.base_ring(), n, blocks) - Unembedding is the inverse of embedding:: - sage: set_random_seed() - sage: Q = QuaternionAlgebra(QQ, -1, -1) - sage: M = random_matrix(Q, 3) - sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M - True + @staticmethod + def real_unembed(M): + """ + The inverse of _embed_complex_matrix(). - """ - 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 complex embedding") - - Q = QuaternionAlgebra(QQ,-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 = _unembed_complex_matrix(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].real() + submat[0,0].imag()*i - z += submat[0,1].real()*j + submat[0,1].imag()*k - elements.append(z) - - return matrix(Q, n/4, elements) - - -# The usual inner product on R^n. -def _usual_ip(x,y): - return x.vector().inner_product(y.vector()) - -# The inner product used for the real symmetric simple EJA. -# We keep it as a separate function because e.g. the complex -# algebra uses the same inner product, except divided by 2. -def _matrix_ip(X,Y): - X_mat = X.natural_representation() - Y_mat = Y.natural_representation() - return (X_mat*Y_mat).trace() - - -class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): - """ - 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:: - SETUP:: + sage: from mjo.eja.eja_algebra import \ + ....: ComplexMatrixEuclideanJordanAlgebra - sage: from mjo.eja.eja_algebra import RealSymmetricEJA + EXAMPLES:: - 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] - sage: J = RealSymmetricEJA(2) - sage: e0, e1, e2 = J.gens() - sage: e0*e0 - e0 - sage: e1*e1 - e0 + e2 - sage: e2*e2 - e2 + TESTS: - TESTS: + Unembedding is the inverse of embedding:: - The degree of this algebra is `(n^2 + n) / 2`:: + 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 - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: J = RealSymmetricEJA(n) - sage: J.degree() == (n^2 + n)/2 - 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") + + field = M.base_ring() # This should already have sqrt2 + R = PolynomialRing(field, 'z') + z = R.gen() + F = NumberField(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. - The Jordan multiplication is what we think it is:: + SETUP:: - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: J = RealSymmetricEJA(n) - sage: x = J.random_element() - sage: y = J.random_element() - 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 + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA - """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - S = _real_symmetric_basis(n, field=field) - (Qs, T) = _multiplication_table_from_matrix_basis(S) + TESTS: - fdeja = super(RealSymmetricEJA, cls) - return fdeja.__classcall_private__(cls, - field, - Qs, - rank=n, - natural_basis=T) + This gives the same answer as the slow, default method implemented + in :class:`MatrixEuclideanJordanAlgebra`:: - def inner_product(self, x, y): - return _matrix_ip(x,y) + 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(FiniteDimensionalEuclideanJordanAlgebra): +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, @@ -2156,21 +1322,20 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): TESTS: - The degree of this algebra is `n^2`:: + The dimension of this algebra is `n^2`:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) + sage: n_max = ComplexHermitianEJA._max_test_case_size() + sage: n = ZZ.random_element(1, n_max) sage: J = ComplexHermitianEJA(n) - sage: J.degree() == n^2 + sage: J.dimension() == n^2 True The Jordan multiplication is what we think it is:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: J = ComplexHermitianEJA(n) - sage: x = J.random_element() - sage: y = J.random_element() + 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() @@ -2180,31 +1345,256 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): 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 + """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - S = _complex_hermitian_basis(n) - (Qs, T) = _multiplication_table_from_matrix_basis(S) - fdeja = super(ComplexHermitianEJA, cls) - return fdeja.__classcall_private__(cls, - field, - Qs, - rank=n, - natural_basis=T) + @classmethod + def _denormalized_basis(cls, n, field): + """ + Returns a basis for the space of complex Hermitian n-by-n matrices. - def inner_product(self, x, y): - # Since a+bi on the diagonal is represented as + 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 = NumberField(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) + I = F.gen() + + # This is like the symmetric case, but we need to be careful: # - # a + bi = [ a b ] - # [ -b a ], + # * We want conjugate-symmetry, not just symmetry. + # * The diagonal will (as a result) be real. # - # we'll double-count the "a" entries if we take the trace of - # the embedding. - return _matrix_ip(x,y)/2 + 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) -class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): + + + @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 complex 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 @@ -2217,21 +1607,20 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): TESTS: - The degree of this algebra is `n^2`:: + The dimension of this algebra is `2*n^2 - n`:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) + sage: n_max = QuaternionHermitianEJA._max_test_case_size() + sage: n = ZZ.random_element(1, n_max) sage: J = QuaternionHermitianEJA(n) - sage: J.degree() == 2*(n^2) - n + sage: J.dimension() == 2*(n^2) - n True The Jordan multiplication is what we think it is:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: J = QuaternionHermitianEJA(n) - sage: x = J.random_element() - sage: y = J.random_element() + 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() @@ -2241,33 +1630,90 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): 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 + """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - S = _quaternion_hermitian_basis(n) - (Qs, T) = _multiplication_table_from_matrix_basis(S) + @classmethod + def _denormalized_basis(cls, n, field): + """ + Returns a basis for the space of quaternion Hermitian n-by-n matrices. - fdeja = super(QuaternionHermitianEJA, cls) - return fdeja.__classcall_private__(cls, - field, - Qs, - rank=n, - natural_basis=T) + 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. - def inner_product(self, x, y): - # Since a+bi+cj+dk on the diagonal is represented as - # - # a + bi +cj + dk = [ a b c d] - # [ -b a -d c] - # [ -c d a -b] - # [ -d -c b a], - # - # we'll quadruple-count the "a" entries if we take the trace of - # the embedding. - return _matrix_ip(x,y)/4 + 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 -class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra): + """ + 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) + return 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 = @@ -2299,27 +1745,55 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: e2*e3 0 + We can change the generator prefix:: + + sage: JordanSpinEJA(2, prefix='B').gens() + (B0, B1) + """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - Qs = [] - id_matrix = matrix.identity(field, n) + 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): - ei = id_matrix.column(i) - Qi = matrix.zero(field, n) - Qi.set_row(0, ei) - Qi.set_column(0, ei) - Qi += matrix.diagonal(n, [ei[0]]*n) - # The addition of the diagonal matrix adds an extra ei[0] in the - # upper-left corner of the matrix. - Qi[0,0] = Qi[0,0] * ~field(2) - Qs.append(Qi) + 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, cls) - return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2)) + fdeja = super(JordanSpinEJA, self) + return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs) def inner_product(self, x, y): - return _usual_ip(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())