X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=0da3eef9aa6e7835b6a6291d008bed4e31571872;hb=2b95649f57f150fed77ef2e62076eb97f54fa6da;hp=2b769ac447bce5b149782bbf96ada23631f9b2c9;hpb=823fb2a587e26436f46854fe44be0e8df46a6715;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 2b769ac..0da3eef 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -61,7 +61,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): """ SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja) + sage: from mjo.eja.eja_algebra import ( + ....: FiniteDimensionalEuclideanJordanAlgebra, + ....: JordanSpinEJA, + ....: random_eja) EXAMPLES: @@ -75,13 +78,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): TESTS: - The ``field`` we're given must be real:: + The ``field`` we're given must be real with ``check=True``:: sage: JordanSpinEJA(2,QQbar) Traceback (most recent call last): ... ValueError: field is not real + The multiplication table must be square with ``check=True``:: + + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),())) + Traceback (most recent call last): + ... + ValueError: multiplication table is not square + """ if check: if not field.is_subring(RR): @@ -96,9 +106,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): category = MagmaticAlgebras(field).FiniteDimensional() category = category.WithBasis().Unital() + # The multiplication table had better be square + n = len(mult_table) + if check: + if not all( len(l) == n for l in mult_table ): + raise ValueError("multiplication table is not square") + fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) fda.__init__(field, - range(len(mult_table)), + range(n), prefix=prefix, category=category) self.print_options(bracket='') @@ -114,6 +130,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): for ls in mult_table ] + if check: + if not self._is_commutative(): + raise ValueError("algebra is not commutative") + if not self._is_jordanian(): + raise ValueError("Jordan identity does not hold") + if not self._inner_product_is_associative(): + raise ValueError("inner product is not associative") def _element_constructor_(self, elt): """ @@ -235,71 +258,74 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): def product_on_basis(self, i, j): return self._multiplication_table[i][j] - def _a_regular_element(self): - """ - Guess a regular element. Needed to compute the basis for our - characteristic polynomial coefficients. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS: - - Ensure that this hacky method succeeds for every algebra that we - know how to construct:: - - sage: set_random_seed() - sage: J = random_eja() - sage: J._a_regular_element().is_regular() - True + def _is_commutative(self): + r""" + Whether or not this algebra's multiplication table is commutative. + This method should of course always return ``True``, unless + this algebra was constructed with ``check=False`` and passed + an invalid multiplication table. """ - gs = self.gens() - z = self.sum( (i+1)*gs[i] for i in range(len(gs)) ) - if not z.is_regular(): - raise ValueError("don't know a regular element") - return z + return all( self.product_on_basis(i,j) == self.product_on_basis(i,j) + for i in range(self.dimension()) + for j in range(self.dimension()) ) + def _is_jordanian(self): + r""" + Whether or not this algebra's multiplication table respects the + Jordan identity `(x^{2})(xy) = x(x^{2}y)`. + + We only check one arrangement of `x` and `y`, so for a + ``True`` result to be truly true, you should also check + :meth:`_is_commutative`. This method should of course always + return ``True``, unless this algebra was constructed with + ``check=False`` and passed an invalid multiplication table. + """ + return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j)) + == + (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j)) + for i in range(self.dimension()) + for j in range(self.dimension()) ) + + def _inner_product_is_associative(self): + r""" + Return whether or not this algebra's inner product `B` is + associative; that is, whether or not `B(xy,z) = B(x,yz)`. - @cached_method - def _charpoly_basis_space(self): - """ - Return the vector space spanned by the basis used in our - characteristic polynomial coefficients. This is used not only to - compute those coefficients, but also any time we need to - evaluate the coefficients (like when we compute the trace or - determinant). - """ - z = self._a_regular_element() - # Don't use the parent vector space directly here in case this - # happens to be a subalgebra. In that case, we would be e.g. - # two-dimensional but span_of_basis() would expect three - # coordinates. - V = VectorSpace(self.base_ring(), self.vector_space().dimension()) - basis = [ (z**k).to_vector() for k in range(self.rank()) ] - V1 = V.span_of_basis( basis ) - b = (V1.basis() + V1.complement().basis()) - return V.span_of_basis(b) + This method should of course always return ``True``, unless + this algebra was constructed with ``check=False`` and passed + an invalid multiplication table. + """ + # Used to check whether or not something is zero in an inexact + # ring. This number is sufficient to allow the construction of + # QuaternionHermitianEJA(2, RDF) with check=True. + epsilon = 1e-16 - def _charpoly_coeff(self, i): - """ - Return the coefficient polynomial "a_{i}" of this algebra's - general characteristic polynomial. + for i in range(self.dimension()): + for j in range(self.dimension()): + for k in range(self.dimension()): + x = self.monomial(i) + y = self.monomial(j) + z = self.monomial(k) + diff = (x*y).inner_product(z) - x.inner_product(y*z) - Having this be a separate cached method lets us compute and - store the trace/determinant (a_{r-1} and a_{0} respectively) - separate from the entire characteristic polynomial. - """ - return self._charpoly_coefficients()[i] + if self.base_ring().is_exact(): + if diff != 0: + return False + else: + if diff.abs() > epsilon: + return False + return True @cached_method - def characteristic_polynomial(self): + def characteristic_polynomial_of(self): """ - Return a characteristic polynomial that works for all elements - of this algebra. + Return the algebra's "characteristic polynomial of" function, + which is itself a multivariate polynomial that, when evaluated + at the coordinates of some algebra element, returns that + element's characteristic polynomial. The resulting polynomial has `n+1` variables, where `n` is the dimension of this algebra. The first `n` variables correspond to @@ -319,7 +345,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Alizadeh, Example 11.11:: sage: J = JordanSpinEJA(3) - sage: p = J.characteristic_polynomial(); p + sage: p = J.characteristic_polynomial_of(); p X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 sage: xvec = J.one().to_vector() sage: p(*xvec) @@ -332,7 +358,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): any argument:: sage: J = TrivialEJA() - sage: J.characteristic_polynomial() + sage: J.characteristic_polynomial_of() 1 """ @@ -498,8 +524,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): """ Return the matrix space in which this algebra's natural basis elements live. + + Generally this will be an `n`-by-`1` column-vector space, + except when the algebra is trivial. There it's `n`-by-`n` + (where `n` is zero), to ensure that two elements of the + natural basis space (empty matrices) can be multiplied. """ - if self._natural_basis is None or len(self._natural_basis) == 0: + if self.is_trivial(): + return MatrixSpace(self.base_ring(), 0) + elif self._natural_basis is None or len(self._natural_basis) == 0: return MatrixSpace(self.base_ring(), self.dimension(), 1) else: return self._natural_basis[0].matrix_space() @@ -632,6 +665,25 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Vector space of degree 6 and dimension 2... sage: J1 Euclidean Jordan algebra of dimension 3... + sage: J0.one().natural_representation() + [0 0 0] + [0 0 0] + [0 0 1] + sage: orig_df = AA.options.display_format + sage: AA.options.display_format = 'radical' + sage: J.from_vector(J5.basis()[0]).natural_representation() + [ 0 0 1/2*sqrt(2)] + [ 0 0 0] + [1/2*sqrt(2) 0 0] + sage: J.from_vector(J5.basis()[1]).natural_representation() + [ 0 0 0] + [ 0 0 1/2*sqrt(2)] + [ 0 1/2*sqrt(2) 0] + sage: AA.options.display_format = orig_df + sage: J1.one().natural_representation() + [1 0 0] + [0 1 0] + [0 0 0] TESTS: @@ -646,9 +698,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J1.superalgebra() == J and J1.dimension() == J.dimension() True - The identity elements in the two subalgebras are the - projections onto their respective subspaces of the - superalgebra's identity element:: + The decomposition is into eigenspaces, and its components are + therefore necessarily orthogonal. Moreover, the identity + elements in the two subalgebras are the projections onto their + respective subspaces of the superalgebra's identity element:: sage: set_random_seed() sage: J = random_eja() @@ -658,6 +711,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): ....: x = J.random_element() sage: c = x.subalgebra_idempotent() sage: J0,J5,J1 = J.peirce_decomposition(c) + sage: ipsum = 0 + sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()): + ....: w = w.superalgebra_element() + ....: y = J.from_vector(y) + ....: z = z.superalgebra_element() + ....: ipsum += w.inner_product(y).abs() + ....: ipsum += w.inner_product(z).abs() + ....: ipsum += y.inner_product(z).abs() + sage: ipsum + 0 sage: J1(c) == J1.one() True sage: J0(J.one() - c) == J0.one() @@ -693,10 +756,61 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return (J0, J5, J1) - def random_elements(self, count): + def random_element(self, thorough=False): + r""" + Return a random element of this algebra. + + Our algebra superclass method only returns a linear + combination of at most two basis elements. We instead + want the vector space "random element" method that + returns a more diverse selection. + + INPUT: + + - ``thorough`` -- (boolean; default False) whether or not we + should generate irrational coefficients for the random + element when our base ring is irrational; this slows the + algebra operations to a crawl, but any truly random method + should include them + + """ + # For a general base ring... maybe we can trust this to do the + # right thing? Unlikely, but. + V = self.vector_space() + v = V.random_element() + + if self.base_ring() is AA: + # The "random element" method of the algebraic reals is + # stupid at the moment, and only returns integers between + # -2 and 2, inclusive: + # + # https://trac.sagemath.org/ticket/30875 + # + # Instead, we implement our own "random vector" method, + # and then coerce that into the algebra. We use the vector + # space degree here instead of the dimension because a + # subalgebra could (for example) be spanned by only two + # vectors, each with five coordinates. We need to + # generate all five coordinates. + if thorough: + v *= QQbar.random_element().real() + else: + v *= QQ.random_element() + + return self.from_vector(V.coordinate_vector(v)) + + def random_elements(self, count, thorough=False): """ Return ``count`` random elements as a tuple. + INPUT: + + - ``thorough`` -- (boolean; default False) whether or not we + should generate irrational coefficients for the random + elements when our base ring is irrational; this slows the + algebra operations to a crawl, but any truly random method + should include them + SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA @@ -711,7 +825,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): True """ - return tuple( self.random_element() for idx in range(count) ) + return tuple( self.random_element(thorough) + for idx in range(count) ) @classmethod def random_instance(cls, field=AA, **kwargs): @@ -726,7 +841,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # there's only one. return cls(field) - n = ZZ.random_element(cls._max_test_case_size()) + 1 + n = ZZ.random_element(cls._max_test_case_size() + 1) return cls(n, field, **kwargs) @cached_method @@ -748,6 +863,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): for k in range(n) ) L_x = matrix(F, n, n, L_x_i_j) + + r = None + if self.rank.is_in_cache(): + r = self.rank() + # There's no need to pad the system with redundant + # columns if we *know* they'll be redundant. + n = r + # Compute an extra power in case the rank is equal to # the dimension (otherwise, we would stop at x^(r-1)). x_powers = [ (L_x**k)*self.one().to_vector() @@ -756,7 +879,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): AE = A.extended_echelon_form() E = AE[:,n:] A_rref = AE[:,:n] - r = A_rref.rank() + if r is None: + r = A_rref.rank() b = x_powers[r] # The theory says that only the first "r" coefficients are @@ -954,7 +1078,7 @@ class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra): return x.to_vector().inner_product(y.to_vector()) -def random_eja(field=AA, nontrivial=False): +def random_eja(field=AA): """ Return a "random" finite-dimensional Euclidean Jordan Algebra. @@ -968,21 +1092,17 @@ def random_eja(field=AA, nontrivial=False): Euclidean Jordan algebra of dimension... """ - eja_classes = [HadamardEJA, - JordanSpinEJA, - RealSymmetricEJA, - ComplexHermitianEJA, - QuaternionHermitianEJA] - if not nontrivial: - eja_classes.append(TrivialEJA) - classname = choice(eja_classes) + classname = choice([TrivialEJA, + HadamardEJA, + JordanSpinEJA, + RealSymmetricEJA, + ComplexHermitianEJA, + QuaternionHermitianEJA]) return classname.random_instance(field=field) - - class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): @staticmethod def _max_test_case_size(): @@ -1040,11 +1160,33 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): # Do this over the rationals and convert back at the end. # Only works because we know the entries of the basis are - # integers. + # integers. The argument ``check=False`` is required + # because the trace inner-product method for this + # class is a stub and can't actually be checked. J = MatrixEuclideanJordanAlgebra(QQ, basis, - normalize_basis=False) - return J._charpoly_coefficients() + normalize_basis=False, + check=False) + a = J._charpoly_coefficients() + + # Unfortunately, changing the basis does change the + # coefficients of the characteristic polynomial, but since + # these are really the coefficients of the "characteristic + # polynomial of" function, everything is still nice and + # unevaluated. It's therefore "obvious" how scaling the + # basis affects the coordinate variables X1, X2, et + # cetera. Scaling the first basis vector up by "n" adds a + # factor of 1/n into every "X1" term, for example. So here + # we simply undo the basis_normalizer scaling that we + # performed earlier. + # + # The a[0] access here is safe because trivial algebras + # won't have any basis normalizers and therefore won't + # make it to this "else" branch. + XS = a[0].parent().gens() + subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i] + for i in range(len(XS)) } + return tuple( a_i.subs(subs_dict) for a_i in a ) @staticmethod @@ -1062,6 +1204,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): # is supposed to hold the entire long vector, and the subspace W # of V will be spanned by the vectors that arise from symmetric # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3. + if len(basis) == 0: + return [] + field = basis[0].base_ring() dimension = basis[0].nrows() @@ -1219,6 +1364,11 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): sage: x.operator().matrix().is_symmetric() True + We can construct the (trivial) algebra of rank zero:: + + sage: RealSymmetricEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ @classmethod def _denormalized_basis(cls, n, field): @@ -1492,6 +1642,11 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): sage: x.operator().matrix().is_symmetric() True + We can construct the (trivial) algebra of rank zero:: + + sage: ComplexHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ @classmethod @@ -1787,6 +1942,11 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): sage: x.operator().matrix().is_symmetric() True + We can construct the (trivial) algebra of rank zero:: + + sage: QuaternionHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ @classmethod def _denormalized_basis(cls, n, field): @@ -2061,3 +2221,50 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra): # largest subalgebra generated by any element. fdeja.__init__(field, mult_table, **kwargs) self.rank.set_cache(0) + + +class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): + r""" + The external (orthogonal) direct sum of two other Euclidean Jordan + algebras. Essentially the Cartesian product of its two factors. + Every Euclidean Jordan algebra decomposes into an orthogonal + direct sum of simple Euclidean Jordan algebras, so no generality + is lost by providing only this construction. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: RealSymmetricEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(3) + sage: J = DirectSumEJA(J1,J2) + sage: J.dimension() + 8 + sage: J.rank() + 5 + + """ + def __init__(self, J1, J2, field=AA, **kwargs): + n1 = J1.dimension() + n2 = J2.dimension() + n = n1+n2 + V = VectorSpace(field, n) + mult_table = [ [ V.zero() for j in range(n) ] + for i in range(n) ] + for i in range(n1): + for j in range(n1): + p = (J1.monomial(i)*J1.monomial(j)).to_vector() + mult_table[i][j] = V(p.list() + [field.zero()]*n2) + + for i in range(n2): + for j in range(n2): + p = (J2.monomial(i)*J2.monomial(j)).to_vector() + mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list()) + + fdeja = super(DirectSumEJA, self) + fdeja.__init__(field, mult_table, **kwargs) + self.rank.set_cache(J1.rank() + J2.rank())