X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=51ff79054eab8cdb83fe48c3d566131598b46278;hb=83d718190836138f62989d68c7b44494ed52c9fd;hp=e075ed21ab287fcefcfbf7c9674a39281b48b69c;hpb=e224054da4441a27efe2301ce18e688f17408ddb;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index e075ed2..51ff790 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -20,7 +20,9 @@ from itertools import repeat from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.magmatic_algebras import MagmaticAlgebras -from sage.combinat.free_module import CombinatorialFreeModule +from sage.categories.sets_cat import cartesian_product +from sage.combinat.free_module import (CombinatorialFreeModule, + CombinatorialFreeModule_CartesianProduct) from sage.matrix.constructor import matrix from sage.matrix.matrix_space import MatrixSpace from sage.misc.cachefunc import cached_method @@ -29,14 +31,197 @@ from sage.modules.free_module import FreeModule, VectorSpace from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF, PolynomialRing, QuadraticField) -from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement -from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator +from mjo.eja.eja_element import (CartesianProductEJAElement, + FiniteDimensionalEJAElement) +from mjo.eja.eja_operator import FiniteDimensionalEJAOperator from mjo.eja.eja_utils import _mat2vec -class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): +class FiniteDimensionalEJA(CombinatorialFreeModule): r""" - The lowest-level class for representing a Euclidean Jordan algebra. + A finite-dimensional Euclidean Jordan algebra. + + INPUT: + + - basis -- a tuple of basis elements in "matrix form," which + must be the same form as the arguments to ``jordan_product`` + and ``inner_product``. In reality, "matrix form" can be either + vectors, matrices, or a Cartesian product (ordered tuple) + of vectors or matrices. All of these would ideally be vector + spaces in sage with no special-casing needed; but in reality + we turn vectors into column-matrices and Cartesian products + `(a,b)` into column matrices `(a,b)^{T}` after converting + `a` and `b` themselves. + + - jordan_product -- function of two elements (in matrix form) + that returns their jordan product in this algebra; this will + be applied to ``basis`` to compute a multiplication table for + the algebra. + + - inner_product -- function of two elements (in matrix form) that + returns their inner product. This will be applied to ``basis`` to + compute an inner-product table (basically a matrix) for this algebra. + """ + Element = FiniteDimensionalEJAElement + + def __init__(self, + basis, + jordan_product, + inner_product, + field=AA, + orthonormalize=True, + associative=False, + cartesian_product=False, + check_field=True, + check_axioms=True, + prefix='e'): + + if check_field: + if not field.is_subring(RR): + # Note: this does return true for the real algebraic + # field, the rationals, and any quadratic field where + # we've specified a real embedding. + raise ValueError("scalar field is not real") + + # If the basis given to us wasn't over the field that it's + # supposed to be over, fix that. Or, you know, crash. + if not cartesian_product: + # The field for a cartesian product algebra comes from one + # of its factors and is the same for all factors, so + # there's no need to "reapply" it on product algebras. + basis = tuple( b.change_ring(field) for b in basis ) + + + if check_axioms: + # Check commutativity of the Jordan and inner-products. + # This has to be done before we build the multiplication + # and inner-product tables/matrices, because we take + # advantage of symmetry in the process. + if not all( jordan_product(bi,bj) == jordan_product(bj,bi) + for bi in basis + for bj in basis ): + raise ValueError("Jordan product is not commutative") + + if not all( inner_product(bi,bj) == inner_product(bj,bi) + for bi in basis + for bj in basis ): + raise ValueError("inner-product is not commutative") + + + category = MagmaticAlgebras(field).FiniteDimensional() + category = category.WithBasis().Unital() + if associative: + # Element subalgebras can take advantage of this. + category = category.Associative() + if cartesian_product: + category = category.CartesianProducts() + + # Call the superclass constructor so that we can use its from_vector() + # method to build our multiplication table. + n = len(basis) + super().__init__(field, + range(n), + prefix=prefix, + category=category, + bracket=False) + + # Now comes all of the hard work. We'll be constructing an + # ambient vector space V that our (vectorized) basis lives in, + # as well as a subspace W of V spanned by those (vectorized) + # basis elements. The W-coordinates are the coefficients that + # we see in things like x = 1*e1 + 2*e2. + vector_basis = basis + + def flatten(b): + # flatten a vector, matrix, or cartesian product of those + # things into a long list. + if cartesian_product: + return sum(( b_i.list() for b_i in b ), []) + else: + return b.list() + + degree = 0 + if n > 0: + degree = len(flatten(basis[0])) + + # Build an ambient space that fits our matrix basis when + # written out as "long vectors." + V = VectorSpace(field, degree) + + # The matrix that will hole the orthonormal -> unorthonormal + # coordinate transformation. + self._deortho_matrix = None + + if orthonormalize: + # Save a copy of the un-orthonormalized basis for later. + # Convert it to ambient V (vector) coordinates while we're + # at it, because we'd have to do it later anyway. + deortho_vector_basis = tuple( V(flatten(b)) for b in basis ) + + from mjo.eja.eja_utils import gram_schmidt + basis = tuple(gram_schmidt(basis, inner_product)) + + # Save the (possibly orthonormalized) matrix basis for + # later... + self._matrix_basis = basis + + # Now create the vector space for the algebra, which will have + # its own set of non-ambient coordinates (in terms of the + # supplied basis). + vector_basis = tuple( V(flatten(b)) for b in basis ) + W = V.span_of_basis( vector_basis, check=check_axioms) + + if orthonormalize: + # Now "W" is the vector space of our algebra coordinates. The + # variables "X1", "X2",... refer to the entries of vectors in + # W. Thus to convert back and forth between the orthonormal + # coordinates and the given ones, we need to stick the original + # basis in W. + U = V.span_of_basis( deortho_vector_basis, check=check_axioms) + self._deortho_matrix = matrix( U.coordinate_vector(q) + for q in vector_basis ) + + + # Now we actually compute the multiplication and inner-product + # tables/matrices using the possibly-orthonormalized basis. + self._inner_product_matrix = matrix.identity(field, n) + self._multiplication_table = [ [0 for j in range(i+1)] + for i in range(n) ] + + # Note: the Jordan and inner-products are defined in terms + # of the ambient basis. It's important that their arguments + # are in ambient coordinates as well. + for i in range(n): + for j in range(i+1): + # ortho basis w.r.t. ambient coords + q_i = basis[i] + q_j = basis[j] + + # The jordan product returns a matrixy answer, so we + # have to convert it to the algebra coordinates. + elt = jordan_product(q_i, q_j) + elt = W.coordinate_vector(V(flatten(elt))) + self._multiplication_table[i][j] = self.from_vector(elt) + + if not orthonormalize: + # If we're orthonormalizing the basis with respect + # to an inner-product, then the inner-product + # matrix with respect to the resulting basis is + # just going to be the identity. + ip = inner_product(q_i, q_j) + self._inner_product_matrix[i,j] = ip + self._inner_product_matrix[j,i] = ip + + self._inner_product_matrix._cache = {'hermitian': True} + self._inner_product_matrix.set_immutable() + + if check_axioms: + 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 _coerce_map_from_base_ring(self): """ Disable the map from the base ring into the algebra. @@ -61,190 +246,140 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): """ return None - def __init__(self, - field, - multiplication_table, - inner_product_table, - prefix='e', - category=None, - matrix_basis=None, - check_field=True, - check_axioms=True): - """ - INPUT: - * field -- the scalar field for this algebra (must be real) + def product_on_basis(self, i, j): + # We only stored the lower-triangular portion of the + # multiplication table. + if j <= i: + return self._multiplication_table[i][j] + else: + return self._multiplication_table[j][i] + + def inner_product(self, x, y): + """ + The inner product associated with this Euclidean Jordan algebra. - * multiplication_table -- the multiplication table for this - algebra's implicit basis. Only the lower-triangular portion - of the table is used, since the multiplication is assumed - to be commutative. + Defaults to the trace inner product, but can be overridden by + subclasses if they are sure that the necessary properties are + satisfied. SETUP:: - sage: from mjo.eja.eja_algebra import ( - ....: FiniteDimensionalEuclideanJordanAlgebra, - ....: JordanSpinEJA, - ....: random_eja) + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: BilinearFormEJA) EXAMPLES: - By definition, Jordan multiplication commutes:: + Our inner product is "associative," which means the following for + a symmetric bilinear form:: sage: set_random_seed() sage: J = random_eja() + sage: x,y,z = J.random_elements(3) + sage: (x*y).inner_product(z) == y.inner_product(x*z) + True + + TESTS: + + Ensure that this is the usual inner product for the algebras + over `R^n`:: + + sage: set_random_seed() + sage: J = HadamardEJA.random_instance() sage: x,y = J.random_elements(2) - sage: x*y == y*x + sage: actual = x.inner_product(y) + sage: expected = x.to_vector().inner_product(y.to_vector()) + sage: actual == expected True - An error is raised if the Jordan product is not commutative:: + Ensure that this is one-half of the trace inner-product in a + BilinearFormEJA that isn't just the reals (when ``n`` isn't + one). This is in Faraut and Koranyi, and also my "On the + symmetry..." paper:: - sage: JP = ((1,2),(0,0)) - sage: IP = ((1,0),(0,1)) - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP) - Traceback (most recent call last): - ... - ValueError: Jordan product is not commutative + sage: set_random_seed() + sage: J = BilinearFormEJA.random_instance() + sage: n = J.dimension() + sage: x = J.random_element() + sage: y = J.random_element() + sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) + True - An error is raised if the inner-product is not commutative:: + """ + B = self._inner_product_matrix + return (B*x.to_vector()).inner_product(y.to_vector()) - sage: JP = ((1,0),(0,1)) - sage: IP = ((1,2),(0,0)) - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP) - Traceback (most recent call last): - ... - ValueError: inner-product is not commutative - TESTS: + def is_associative(self): + r""" + Return whether or not this algebra's Jordan product is associative. - The ``field`` we're given must be real with ``check_field=True``:: + SETUP:: - sage: JordanSpinEJA(2, field=QQbar) - Traceback (most recent call last): - ... - ValueError: scalar field is not real - sage: JordanSpinEJA(2, field=QQbar, check_field=False) - Euclidean Jordan algebra of dimension 2 over Algebraic Field + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA - The multiplication table must be square with ``check_axioms=True``:: + EXAMPLES:: - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),)) - Traceback (most recent call last): - ... - ValueError: multiplication table is not square + sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False) + sage: J.is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: A.is_associative() + True - The multiplication and inner-product tables must be the same - size (and in particular, the inner-product table must also be - square) with ``check_axioms=True``:: + """ + return "Associative" in self.category().axioms() - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(())) - Traceback (most recent call last): - ... - ValueError: multiplication and inner-product tables are - different sizes - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),)) - Traceback (most recent call last): - ... - ValueError: multiplication and inner-product tables are - different sizes + 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_axioms=False`` and passed an invalid multiplication table. """ - if check_field: - if not field.is_subring(RR): - # Note: this does return true for the real algebraic - # field, the rationals, and any quadratic field where - # we've specified a real embedding. - raise ValueError("scalar field is not real") + return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j]) + == + (self.gens()[i])*((self.gens()[i]**2)*self.gens()[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)`. - # The multiplication and inner-product tables should be square - # if the user wants us to verify them. And we verify them as - # soon as possible, because we want to exploit their symmetry. - n = len(multiplication_table) - if check_axioms: - if not all( len(l) == n for l in multiplication_table ): - raise ValueError("multiplication table is not square") - - # If the multiplication table is square, we can check if - # the inner-product table is square by comparing it to the - # multiplication table's dimensions. - msg = "multiplication and inner-product tables are different sizes" - if not len(inner_product_table) == n: - raise ValueError(msg) - - if not all( len(l) == n for l in inner_product_table ): - raise ValueError(msg) - - # Check commutativity of the Jordan product (symmetry of - # the multiplication table) and the commutativity of the - # inner-product (symmetry of the inner-product table) - # first if we're going to check them at all.. This has to - # be done before we define product_on_basis(), because - # that method assumes that self._multiplication_table is - # symmetric. And it has to be done before we build - # self._inner_product_matrix, because the process used to - # construct it assumes symmetry as well. - if not all( multiplication_table[j][i] - == multiplication_table[i][j] - for i in range(n) - for j in range(i+1) ): - raise ValueError("Jordan product is not commutative") + This method should of course always return ``True``, unless + this algebra was constructed with ``check_axioms=False`` and + passed an invalid Jordan or inner-product. + """ - if not all( inner_product_table[j][i] - == inner_product_table[i][j] - for i in range(n) - for j in range(i+1) ): - raise ValueError("inner-product is not commutative") + # Used to check whether or not something is zero in an inexact + # ring. This number is sufficient to allow the construction of + # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True. + epsilon = 1e-16 - self._matrix_basis = matrix_basis - - if category is None: - category = MagmaticAlgebras(field).FiniteDimensional() - category = category.WithBasis().Unital() - - fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) - fda.__init__(field, - range(n), - 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. - # - # Note: we take advantage of symmetry here, and only store - # the lower-triangular portion of the table. - self._multiplication_table = [ [ self.vector_space().zero() - for j in range(i+1) ] - for i in range(n) ] + for i in range(self.dimension()): + for j in range(self.dimension()): + for k in range(self.dimension()): + x = self.gens()[i] + y = self.gens()[j] + z = self.gens()[k] + diff = (x*y).inner_product(z) - x.inner_product(y*z) - for i in range(n): - for j in range(i+1): - elt = self.from_vector(multiplication_table[i][j]) - self._multiplication_table[i][j] = elt - - self._multiplication_table = tuple(map(tuple, self._multiplication_table)) - - # Save our inner product as a matrix, since the efficiency of - # matrix multiplication will usually outweigh the fact that we - # have to store a redundant upper- or lower-triangular part. - # Pre-cache the fact that these are Hermitian (real symmetric, - # in fact) in case some e.g. matrix multiplication routine can - # take advantage of it. - ip_matrix_constructor = lambda i,j: inner_product_table[i][j] if j <= i else inner_product_table[j][i] - self._inner_product_matrix = matrix(field, n, ip_matrix_constructor) - self._inner_product_matrix._cache = {'hermitian': True} - self._inner_product_matrix.set_immutable() + if self.base_ring().is_exact(): + if diff != 0: + return False + else: + if diff.abs() > epsilon: + return False - if check_axioms: - 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") + return True def _element_constructor_(self, elt): """ @@ -307,6 +442,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # that the integer 3 belongs to the space of 2-by-2 matrices. raise ValueError(msg) + try: + elt = elt.column() + except (AttributeError, TypeError): + # Try to convert a vector into a column-matrix + pass + if elt not in self.matrix_space(): raise ValueError(msg) @@ -351,74 +492,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): fmt = "Euclidean Jordan algebra of dimension {} over {}" return fmt.format(self.dimension(), self.base_ring()) - def product_on_basis(self, i, j): - # We only stored the lower-triangular portion of the - # multiplication table. - if j <= i: - return self._multiplication_table[i][j] - else: - return self._multiplication_table[j][i] - - 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_axioms=False`` and - passed an invalid multiplication table. - """ - return all( self.product_on_basis(i,j) == self.product_on_basis(i,j) - for i in range(self.dimension()) - for j in range(self.dimension()) ) - - 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_axioms=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)`. - - This method should of course always return ``True``, unless - this algebra was constructed with ``check_axioms=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, field=RDF) with check_axioms=True. - epsilon = 1e-16 - - for i in range(self.dimension()): - for j in range(self.dimension()): - for k in range(self.dimension()): - x = self.monomial(i) - y = self.monomial(j) - z = self.monomial(k) - diff = (x*y).inner_product(z) - x.inner_product(y*z) - - 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_of(self): @@ -614,20 +687,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): """ n = self.dimension() - M = [ [ self.zero() for j in range(n) ] - for i in range(n) ] - for i in range(n): - for j in range(i+1): - M[i][j] = self._multiplication_table[i][j] - M[j][i] = M[i][j] + # Prepend the header row. + M = [["*"] + list(self.gens())] - for i in range(n): - # Prepend the left "header" column entry Can't do this in - # the loop because it messes up the symmetry. - M[i] = [self.monomial(i)] + M[i] + # And to each subsequent row, prepend an entry that belongs to + # the left-side "header column." + M += [ [self.gens()[i]] + [ self.product_on_basis(i,j) + for j in range(n) ] + for i in range(n) ] - # Prepend the header row. - M = [["*"] + list(self.gens())] + M return table(M, header_row=True, header_column=True, frame=True) @@ -655,7 +723,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Why implement this for non-matrix algebras? Avoiding special cases for the :class:`BilinearFormEJA` pays with simplicity in its own right. But mainly, we would like to be able to assume - that elements of a :class:`DirectSumEJA` can be displayed + that elements of a :class:`CartesianProductEJA` can be displayed nicely, without having to have special classes for direct sums one of whose components was a matrix algebra. @@ -686,11 +754,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): [0], [1] ) """ - if self._matrix_basis is None: - M = self.matrix_space() - return tuple( M(b.to_vector()) for b in self.basis() ) - else: - return self._matrix_basis + return self._matrix_basis def matrix_space(self): @@ -708,10 +772,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): """ if self.is_trivial(): return MatrixSpace(self.base_ring(), 0) - elif self._matrix_basis is None or len(self._matrix_basis) == 0: - return MatrixSpace(self.base_ring(), self.dimension(), 1) else: - return self._matrix_basis[0].matrix_space() + return self.matrix_basis()[0].parent() @cached_method @@ -724,23 +786,57 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: random_eja) - EXAMPLES:: + EXAMPLES: + + We can compute unit element in the Hadamard EJA:: + + sage: J = HadamardEJA(5) + sage: J.one() + e0 + e1 + e2 + e3 + e4 + + The unit element in the Hadamard EJA is inherited in the + subalgebras generated by its elements:: sage: J = HadamardEJA(5) sage: J.one() e0 + e1 + e2 + e3 + e4 + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: A.one() + f0 + sage: A.one().superalgebra_element() + e0 + e1 + e2 + e3 + e4 TESTS: - The identity element acts like the identity:: + The identity element acts like the identity, regardless of + whether or not we orthonormalize:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: J.one()*x == x and x*J.one() == x True + sage: A = x.subalgebra_generated_by() + sage: y = A.random_element() + sage: A.one()*y == y and y*A.one() == y + True - The matrix of the unit element's operator is the identity:: + :: + + sage: set_random_seed() + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: x = J.random_element() + sage: J.one()*x == x and x*J.one() == x + True + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: y = A.random_element() + sage: A.one()*y == y and y*A.one() == y + True + + The matrix of the unit element's operator is the identity, + regardless of the base field and whether or not we + orthonormalize:: sage: set_random_seed() sage: J = random_eja() @@ -748,6 +844,27 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: expected = matrix.identity(J.base_ring(), J.dimension()) sage: actual == expected True + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by() + sage: actual = A.one().operator().matrix() + sage: expected = matrix.identity(A.base_ring(), A.dimension()) + sage: actual == expected + True + + :: + + sage: set_random_seed() + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: actual = J.one().operator().matrix() + sage: expected = matrix.identity(J.base_ring(), J.dimension()) + sage: actual == expected + True + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: actual = A.one().operator().matrix() + sage: expected = matrix.identity(A.base_ring(), A.dimension()) + sage: actual == expected + True Ensure that the cached unit element (often precomputed by hand) agrees with the computed one:: @@ -759,6 +876,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J.one() == cached True + :: + + sage: set_random_seed() + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: cached = J.one() + sage: J.one.clear_cache() + sage: J.one() == cached + True + """ # We can brute-force compute the matrices of the operators # that correspond to the basis elements of this algebra. @@ -903,14 +1029,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): if not c.is_idempotent(): raise ValueError("element is not idempotent: %s" % c) - from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra + from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra # Default these to what they should be if they turn out to be # trivial, because eigenspaces_left() won't return eigenvalues # corresponding to trivial spaces (e.g. it returns only the # eigenspace corresponding to lambda=1 if you take the # decomposition relative to the identity element). - trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ()) + trivial = FiniteDimensionalEJASubalgebra(self, ()) J0 = trivial # eigenvalue zero J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half J1 = trivial # eigenvalue one @@ -920,9 +1046,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): J5 = eigspace else: gens = tuple( self.from_vector(b) for b in eigspace.basis() ) - subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, - gens, - check_axioms=False) + subalg = FiniteDimensionalEJASubalgebra(self, + gens, + check_axioms=False) if eigval == 0: J0 = subalg elif eigval == 1: @@ -1011,6 +1137,21 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): r""" The `r` polynomial coefficients of the "characteristic polynomial of" function. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS: + + The theory shows that these are all homogeneous polynomials of + a known degree:: + + sage: set_random_seed() + sage: J = random_eja() + sage: all(p.is_homogeneous() for p in J._charpoly_coefficients()) + True + """ n = self.dimension() R = self.coordinate_polynomial_ring() @@ -1020,7 +1161,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): def L_x_i_j(i,j): # From a result in my book, these are the entries of the # basis representation of L_x. - return sum( vars[k]*self.monomial(k).operator().matrix()[i,j] + return sum( vars[k]*self.gens()[k].operator().matrix()[i,j] for k in range(n) ) L_x = matrix(F, n, n, L_x_i_j) @@ -1046,10 +1187,17 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # The theory says that only the first "r" coefficients are # nonzero, and they actually live in the original polynomial - # ring and not the fraction field. We negate them because - # in the actual characteristic polynomial, they get moved - # to the other side where x^r lives. - return -A_rref.solve_right(E*b).change_ring(R)[:r] + # ring and not the fraction field. We negate them because in + # the actual characteristic polynomial, they get moved to the + # other side where x^r lives. We don't bother to trim A_rref + # down to a square matrix and solve the resulting system, + # because the upper-left r-by-r portion of A_rref is + # guaranteed to be the identity matrix, so e.g. + # + # A_rref.solve_right(Y) + # + # would just be returning Y. + return (-E*b)[:r].change_ring(R) @cached_method def rank(self): @@ -1110,7 +1258,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: set_random_seed() # long time sage: J = random_eja() # long time - sage: caches = J.rank() # long time + sage: cached = J.rank() # long time sage: J.rank.clear_cache() # long time sage: J.rank() == cached # long time True @@ -1137,9 +1285,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return self.zero().to_vector().parent().ambient_vector_space() - Element = FiniteDimensionalEuclideanJordanAlgebraElement + Element = FiniteDimensionalEJAElement -class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): +class RationalBasisEJA(FiniteDimensionalEJA): r""" New class for algebras whose supplied basis elements have all rational entries. @@ -1166,11 +1314,8 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr jordan_product, inner_product, field=AA, - orthonormalize=True, - prefix='e', - category=None, check_field=True, - check_axioms=True): + **kwargs): if check_field: # Abuse the check_field parameter to check that the entries of @@ -1178,168 +1323,29 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr if not all( all(b_i in QQ for b_i in b.list()) for b in basis ): raise TypeError("basis not rational") - # Temporary(?) hack to ensure that the matrix and vector bases - # are over the same ring. - basis = tuple( b.change_ring(field) for b in basis ) - - n = len(basis) - vector_basis = basis - - from sage.structure.element import is_Matrix - basis_is_matrices = False - - degree = 0 - if n > 0: - if is_Matrix(basis[0]): - basis_is_matrices = True - from mjo.eja.eja_utils import _vec2mat - vector_basis = tuple( map(_mat2vec,basis) ) - degree = basis[0].nrows()**2 - else: - degree = basis[0].degree() - - V = VectorSpace(field, degree) - - # If we were asked to orthonormalize, and if the orthonormal - # basis is different from the given one, then we also want to - # compute multiplication and inner-product tables for the - # deorthonormalized basis. These can be used later to - # construct a deorthonormalized copy of this algebra over QQ - # in which several operations are much faster. self._rational_algebra = None - - if orthonormalize: - if self.base_ring() is not QQ: - # There's no point in constructing the extra algebra if this - # one is already rational. If the original basis is rational - # but normalization would make it irrational, then this whole - # constructor will just fail anyway as it tries to stick an - # irrational number into a rational algebra. - # - # Note: the same Jordan and inner-products work here, - # because they are necessarily defined with respect to - # ambient coordinates and not any particular basis. - self._rational_algebra = RationalBasisEuclideanJordanAlgebra( - basis, - jordan_product, - inner_product, - field=QQ, - orthonormalize=False, - prefix=prefix, - category=category, - check_field=False, - check_axioms=False) - - # Compute the deorthonormalized tables before we orthonormalize - # the given basis. The "check" parameter here guarantees that - # the basis is linearly-independent. - W = V.span_of_basis( vector_basis, check=check_axioms) - - # Note: the Jordan and inner-products are defined in terms - # of the ambient basis. It's important that their arguments - # are in ambient coordinates as well. - for i in range(n): - for j in range(i+1): - # given basis w.r.t. ambient coords - q_i = vector_basis[i] - q_j = vector_basis[j] - - if basis_is_matrices: - q_i = _vec2mat(q_i) - q_j = _vec2mat(q_j) - - elt = jordan_product(q_i, q_j) - ip = inner_product(q_i, q_j) - - if basis_is_matrices: - # do another mat2vec because the multiplication - # table is in terms of vectors - elt = _mat2vec(elt) - - # We overwrite the name "vector_basis" in a second, but never modify it - # in place, to this effectively makes a copy of it. - deortho_vector_basis = vector_basis - self._deortho_matrix = None - - if orthonormalize: - from mjo.eja.eja_utils import gram_schmidt - if basis_is_matrices: - vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y)) - vector_basis = gram_schmidt(vector_basis, vector_ip) - else: - vector_basis = gram_schmidt(vector_basis, inner_product) - - # Normalize the "matrix" basis, too! - basis = vector_basis - - if basis_is_matrices: - basis = tuple( map(_vec2mat,basis) ) - - W = V.span_of_basis( vector_basis, check=check_axioms) - - # Now "W" is the vector space of our algebra coordinates. The - # variables "X1", "X2",... refer to the entries of vectors in - # W. Thus to convert back and forth between the orthonormal - # coordinates and the given ones, we need to stick the original - # basis in W. - U = V.span_of_basis( deortho_vector_basis, check=check_axioms) - self._deortho_matrix = matrix( U.coordinate_vector(q) - for q in vector_basis ) - - # If the superclass constructor is going to verify the - # symmetry of this table, it has better at least be - # square... - if check_axioms: - mult_table = [ [0 for j in range(n)] for i in range(n) ] - ip_table = [ [0 for j in range(n)] for i in range(n) ] - else: - mult_table = [ [0 for j in range(i+1)] for i in range(n) ] - ip_table = [ [0 for j in range(i+1)] for i in range(n) ] - - # Note: the Jordan and inner-products are defined in terms - # of the ambient basis. It's important that their arguments - # are in ambient coordinates as well. - for i in range(n): - for j in range(i+1): - # ortho basis w.r.t. ambient coords - q_i = vector_basis[i] - q_j = vector_basis[j] - - if basis_is_matrices: - q_i = _vec2mat(q_i) - q_j = _vec2mat(q_j) - - elt = jordan_product(q_i, q_j) - ip = inner_product(q_i, q_j) - - if basis_is_matrices: - # do another mat2vec because the multiplication - # table is in terms of vectors - elt = _mat2vec(elt) - - elt = W.coordinate_vector(elt) - mult_table[i][j] = elt - ip_table[i][j] = ip - if check_axioms: - # The tables are square if we're verifying that they - # are commutative. - mult_table[j][i] = elt - ip_table[j][i] = ip - - if basis_is_matrices: - for m in basis: - m.set_immutable() - else: - basis = tuple( x.column() for x in basis ) - - super().__init__(field, - mult_table, - ip_table, - prefix, - category, - basis, # matrix basis - check_field, - check_axioms) + if field is not QQ: + # There's no point in constructing the extra algebra if this + # one is already rational. + # + # Note: the same Jordan and inner-products work here, + # because they are necessarily defined with respect to + # ambient coordinates and not any particular basis. + self._rational_algebra = FiniteDimensionalEJA( + basis, + jordan_product, + inner_product, + field=QQ, + orthonormalize=False, + check_field=False, + check_axioms=False) + + super().__init__(basis, + jordan_product, + inner_product, + field=field, + check_field=check_field, + **kwargs) @cached_method def _charpoly_coefficients(self): @@ -1365,13 +1371,12 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr Algebraic Real Field """ - if self.base_ring() is QQ or self._rational_algebra is None: + if self._rational_algebra is None: # There's no need to construct *another* algebra over the # rationals if this one is already over the # rationals. Likewise, if we never orthonormalized our # basis, we might as well just use the given one. - superclass = super(RationalBasisEuclideanJordanAlgebra, self) - return superclass._charpoly_coefficients() + return super()._charpoly_coefficients() # Do the computation over the rationals. The answer will be # the same, because all we've done is a change of basis. @@ -1379,7 +1384,14 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr a = ( a_i.change_ring(self.base_ring()) for a_i in self._rational_algebra._charpoly_coefficients() ) - # Now convert the coordinate variables back to the + if self._deortho_matrix is None: + # This can happen if our base ring was, say, AA and we + # chose not to (or didn't need to) orthonormalize. It's + # still faster to do the computations over QQ even if + # the numbers in the boxes stay the same. + return tuple(a) + + # Otherwise, convert the coordinate variables back to the # deorthonormalized ones. R = self.coordinate_polynomial_ring() from sage.modules.free_module_element import vector @@ -1389,7 +1401,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr subs_dict = { X[i]: BX[i] for i in range(len(X)) } return tuple( a_i.subs(subs_dict) for a_i in a ) -class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): +class ConcreteEJA(RationalBasisEJA): r""" A class for the Euclidean Jordan algebras that we know by name. @@ -1400,7 +1412,7 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): SETUP:: - sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra + sage: from mjo.eja.eja_algebra import ConcreteEJA TESTS: @@ -1408,7 +1420,7 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): product, unless we specify otherwise:: sage: set_random_seed() - sage: J = ConcreteEuclideanJordanAlgebra.random_instance() + sage: J = ConcreteEJA.random_instance() sage: all( b.norm() == 1 for b in J.gens() ) True @@ -1419,7 +1431,7 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): EJA the operator is self-adjoint by the Jordan axiom:: sage: set_random_seed() - sage: J = ConcreteEuclideanJordanAlgebra.random_instance() + sage: J = ConcreteEJA.random_instance() sage: x = J.random_element() sage: x.operator().is_self_adjoint() True @@ -1451,15 +1463,30 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): from sage.misc.prandom import choice eja_class = choice(cls.__subclasses__()) - # These all bubble up to the RationalBasisEuclideanJordanAlgebra - # superclass constructor, so any (kw)args valid there are also - # valid here. + # These all bubble up to the RationalBasisEJA superclass + # constructor, so any (kw)args valid there are also valid + # here. return eja_class.random_instance(*args, **kwargs) -class MatrixEuclideanJordanAlgebra: +class MatrixEJA: @staticmethod - def real_embed(M): + def dimension_over_reals(): + r""" + The dimension of this matrix's base ring over the reals. + + The reals are dimension one over themselves, obviously; that's + just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi` + have dimension two. Finally, the quaternions have dimension + four over the reals. + + This is used to determine the size of the matrix returned from + :meth:`real_embed`, among other things. + """ + raise NotImplementedError + + @classmethod + def real_embed(cls,M): """ Embed the matrix ``M`` into a space of real matrices. @@ -1469,25 +1496,86 @@ class MatrixEuclideanJordanAlgebra: 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) + real_embed(M*N) = real_embed(M)*real_embed(N) + + """ + if M.ncols() != M.nrows(): + raise ValueError("the matrix 'M' must be square") + return M + + + @classmethod + def real_unembed(cls,M): + """ + The inverse of :meth:`real_embed`. + """ + if M.ncols() != M.nrows(): + raise ValueError("the matrix 'M' must be square") + if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero(): + raise ValueError("the matrix 'M' must be a real embedding") + return M + + @staticmethod + def jordan_product(X,Y): + return (X*Y + Y*X)/2 + + @classmethod + def trace_inner_product(cls,X,Y): + r""" + Compute the trace inner-product of two real-embeddings. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA) + + EXAMPLES:: + + This gives the same answer as it would if we computed the trace + from the unembedded (original) matrices:: + + sage: set_random_seed() + sage: J = RealSymmetricEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: Xe = x.to_matrix() + sage: Ye = y.to_matrix() + sage: X = J.real_unembed(Xe) + sage: Y = J.real_unembed(Ye) + sage: expected = (X*Y).trace() + sage: actual = J.trace_inner_product(Xe,Ye) + sage: actual == expected + True - """ - raise NotImplementedError + :: + sage: set_random_seed() + sage: J = ComplexHermitianEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: Xe = x.to_matrix() + sage: Ye = y.to_matrix() + sage: X = J.real_unembed(Xe) + sage: Y = J.real_unembed(Ye) + sage: expected = (X*Y).trace().real() + sage: actual = J.trace_inner_product(Xe,Ye) + sage: actual == expected + True - @staticmethod - def real_unembed(M): - """ - The inverse of :meth:`real_embed`. - """ - raise NotImplementedError + :: - @staticmethod - def jordan_product(X,Y): - return (X*Y + Y*X)/2 + sage: set_random_seed() + sage: J = QuaternionHermitianEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: Xe = x.to_matrix() + sage: Ye = y.to_matrix() + sage: X = J.real_unembed(Xe) + sage: Y = J.real_unembed(Ye) + sage: expected = (X*Y).trace().coefficient_tuple()[0] + sage: actual = J.trace_inner_product(Xe,Ye) + sage: actual == expected + True - @classmethod - def trace_inner_product(cls,X,Y): + """ Xu = cls.real_unembed(X) Yu = cls.real_unembed(Y) tr = (Xu*Yu).trace() @@ -1502,26 +1590,13 @@ class MatrixEuclideanJordanAlgebra: 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 - +class RealMatrixEJA(MatrixEJA): @staticmethod - def real_unembed(M): - """ - The identity function, for unembedding real matrices from real - matrices. - """ - return M + def dimension_over_reals(): + return 1 -class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, - RealMatrixEuclideanJordanAlgebra): +class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1639,13 +1714,55 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, self.jordan_product, self.trace_inner_product, **kwargs) + + # TODO: this could be factored out somehow, but is left here + # because the MatrixEJA is not presently a subclass of the + # FDEJA class that defines rank() and one(). self.rank.set_cache(n) - self.one.set_cache(self(matrix.identity(ZZ,n))) + idV = matrix.identity(ZZ, self.dimension_over_reals()*n) + self.one.set_cache(self(idV)) + -class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): +class ComplexMatrixEJA(MatrixEJA): + # A manual dictionary-cache for the complex_extension() method, + # since apparently @classmethods can't also be @cached_methods. + _complex_extension = {} + + @classmethod + def complex_extension(cls,field): + r""" + The complex field that we embed/unembed, as an extension + of the given ``field``. + """ + if field in cls._complex_extension: + return cls._complex_extension[field] + + # Sage doesn't know how to adjoin the complex "i" (the root of + # x^2 + 1) to a field in a general way. Here, we just enumerate + # all of the cases that I have cared to support so far. + if field is AA: + # Sage doesn't know how to embed AA into QQbar, i.e. how + # to adjoin sqrt(-1) to AA. + F = QQbar + elif not field.is_exact(): + # RDF or RR + F = field.complex_field() + else: + # Works for QQ and... maybe some other fields. + R = PolynomialRing(field, 'z') + z = R.gen() + F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) + + cls._complex_extension[field] = F + return F + @staticmethod - def real_embed(M): + def dimension_over_reals(): + return 2 + + @classmethod + def real_embed(cls,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 + @@ -1653,8 +1770,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): SETUP:: - sage: from mjo.eja.eja_algebra import \ - ....: ComplexMatrixEuclideanJordanAlgebra + sage: from mjo.eja.eja_algebra import ComplexMatrixEJA EXAMPLES:: @@ -1664,7 +1780,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: x3 = F(-i) sage: x4 = F(6) sage: M = matrix(F,2,[[x1,x2],[x3,x4]]) - sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M) + sage: ComplexMatrixEJA.real_embed(M) [ 4 -2| 1 2] [ 2 4|-2 1] [-----+-----] @@ -1680,38 +1796,37 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): 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 = ComplexMatrixEJA.real_embed(X) + sage: Ye = ComplexMatrixEJA.real_embed(Y) + sage: XYe = ComplexMatrixEJA.real_embed(X*Y) sage: Xe*Ye == XYe True """ + super(ComplexMatrixEJA,cls).real_embed(M) n = M.nrows() - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") # We don't need any adjoined elements... field = M.base_ring().base_ring() blocks = [] for z in M.list(): - a = z.list()[0] # real part, I guess - b = z.list()[1] # imag part, I guess - blocks.append(matrix(field, 2, [[a,b],[-b,a]])) + a = z.real() + b = z.imag() + blocks.append(matrix(field, 2, [ [ a, b], + [-b, a] ])) return matrix.block(field, n, blocks) - @staticmethod - def real_unembed(M): + @classmethod + def real_unembed(cls,M): """ The inverse of _embed_complex_matrix(). SETUP:: - sage: from mjo.eja.eja_algebra import \ - ....: ComplexMatrixEuclideanJordanAlgebra + sage: from mjo.eja.eja_algebra import ComplexMatrixEJA EXAMPLES:: @@ -1719,7 +1834,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): ....: [-2, 1, -4, 3], ....: [ 9, 10, 11, 12], ....: [-10, 9, -12, 11] ]) - sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A) + sage: ComplexMatrixEJA.real_unembed(A) [ 2*I + 1 4*I + 3] [ 10*I + 9 12*I + 11] @@ -1730,36 +1845,23 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): 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 + sage: Me = ComplexMatrixEJA.real_embed(M) + sage: ComplexMatrixEJA.real_unembed(Me) == M True """ + super(ComplexMatrixEJA,cls).real_unembed(M) n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(2).is_zero(): - raise ValueError("the matrix 'M' must be a complex embedding") - - # If "M" was normalized, its base ring might have roots - # adjoined and they can stick around after unembedding. - field = M.base_ring() - R = PolynomialRing(field, 'z') - z = R.gen() - if field is AA: - # Sage doesn't know how to embed AA into QQbar, i.e. how - # to adjoin sqrt(-1) to AA. - F = QQbar - else: - F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) + d = cls.dimension_over_reals() + F = cls.complex_extension(M.base_ring()) 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 range(n/2): - for j in range(n/2): - submat = M[2*k:2*k+2,2*j:2*j+2] + for k in range(n/d): + for j in range(n/d): + submat = M[d*k:d*k+d,d*j:d*j+d] if submat[0,0] != submat[1,1]: raise ValueError('bad on-diagonal submatrix') if submat[0,1] != -submat[1,0]: @@ -1767,42 +1869,10 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): z = submat[0,0] + submat[0,1]*i elements.append(z) - return matrix(F, n/2, elements) - - - @classmethod - def trace_inner_product(cls,X,Y): - """ - Compute a matrix inner product in this algebra directly from - its real embedding. - - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA - - TESTS: - - This gives the same answer as the slow, default method implemented - in :class:`MatrixEuclideanJordanAlgebra`:: - - sage: set_random_seed() - sage: J = ComplexHermitianEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: Xe = x.to_matrix() - sage: Ye = y.to_matrix() - sage: X = ComplexHermitianEJA.real_unembed(Xe) - sage: Y = ComplexHermitianEJA.real_unembed(Ye) - sage: expected = (X*Y).trace().real() - sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye) - sage: actual == expected - True - - """ - return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2 + return matrix(F, n/d, elements) -class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, - ComplexMatrixEuclideanJordanAlgebra): +class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -1879,7 +1949,6 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: field = QuadraticField(2, 'sqrt2') sage: B = ComplexHermitianEJA._denormalized_basis(n) sage: all( M.is_symmetric() for M in B) True @@ -1897,18 +1966,27 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, # * The diagonal will (as a result) be real. # S = [] + Eij = matrix.zero(F,n) for i in range(n): for j in range(i+1): - Eij = matrix(F, n, lambda k,l: k==i and l==j) + # "build" E_ij + Eij[i,j] = 1 if i == j: Sij = cls.real_embed(Eij) S.append(Sij) else: # The second one has a minus because it's conjugated. - Sij_real = cls.real_embed(Eij + Eij.transpose()) + Eij[j,i] = 1 # Eij = Eij + Eij.transpose() + Sij_real = cls.real_embed(Eij) S.append(Sij_real) - Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose()) + # Eij = I*Eij - I*Eij.transpose() + Eij[i,j] = I + Eij[j,i] = -I + Sij_imag = cls.real_embed(Eij) S.append(Sij_imag) + Eij[j,i] = 0 + # "erase" E_ij + Eij[i,j] = 0 # Since we embedded these, we can drop back to the "field" that we # started with instead of the complex extension "F". @@ -1924,8 +2002,12 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, self.jordan_product, self.trace_inner_product, **kwargs) + # TODO: this could be factored out somehow, but is left here + # because the MatrixEJA is not presently a subclass of the + # FDEJA class that defines rank() and one(). self.rank.set_cache(n) - # TODO: pre-cache the identity! + idV = matrix.identity(ZZ, self.dimension_over_reals()*n) + self.one.set_cache(self(idV)) @staticmethod def _max_random_instance_size(): @@ -1939,9 +2021,32 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) -class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): +class QuaternionMatrixEJA(MatrixEJA): + + # A manual dictionary-cache for the quaternion_extension() method, + # since apparently @classmethods can't also be @cached_methods. + _quaternion_extension = {} + + @classmethod + def quaternion_extension(cls,field): + r""" + The quaternion field that we embed/unembed, as an extension + of the given ``field``. + """ + if field in cls._quaternion_extension: + return cls._quaternion_extension[field] + + Q = QuaternionAlgebra(field,-1,-1) + + cls._quaternion_extension[field] = Q + return Q + @staticmethod - def real_embed(M): + def dimension_over_reals(): + return 4 + + @classmethod + def real_embed(cls,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 @@ -1951,8 +2056,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): SETUP:: - sage: from mjo.eja.eja_algebra import \ - ....: QuaternionMatrixEuclideanJordanAlgebra + sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA EXAMPLES:: @@ -1960,7 +2064,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): 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) + sage: QuaternionMatrixEJA.real_embed(M) [ 1 2 3 4] [-2 1 -4 3] [-3 4 1 -2] @@ -1973,17 +2077,16 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): 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 = QuaternionMatrixEJA.real_embed(X) + sage: Ye = QuaternionMatrixEJA.real_embed(Y) + sage: XYe = QuaternionMatrixEJA.real_embed(X*Y) sage: Xe*Ye == XYe True """ + super(QuaternionMatrixEJA,cls).real_embed(M) 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() @@ -1997,7 +2100,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): d = t[3] cplxM = matrix(F, 2, [[ a + b*i, c + d*i], [-c + d*i, a - b*i]]) - realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM) + realM = ComplexMatrixEJA.real_embed(cplxM) blocks.append(realM) # We should have real entries by now, so use the realest field @@ -2006,15 +2109,14 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): - @staticmethod - def real_unembed(M): + @classmethod + def real_unembed(cls,M): """ The inverse of _embed_quaternion_matrix(). SETUP:: - sage: from mjo.eja.eja_algebra import \ - ....: QuaternionMatrixEuclideanJordanAlgebra + sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA EXAMPLES:: @@ -2022,7 +2124,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): ....: [-2, 1, -4, 3], ....: [-3, 4, 1, -2], ....: [-4, -3, 2, 1]]) - sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M) + sage: QuaternionMatrixEJA.real_unembed(M) [1 + 2*i + 3*j + 4*k] TESTS: @@ -2032,31 +2134,28 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): 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 + sage: Me = QuaternionMatrixEJA.real_embed(M) + sage: QuaternionMatrixEJA.real_unembed(Me) == M True """ + super(QuaternionMatrixEJA,cls).real_unembed(M) n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(4).is_zero(): - raise ValueError("the matrix 'M' must be a quaternion embedding") + d = cls.dimension_over_reals() # Use the base ring of the matrix to ensure that its entries can be # multiplied by elements of the quaternion algebra. - field = M.base_ring() - Q = QuaternionAlgebra(field,-1,-1) + Q = cls.quaternion_extension(M.base_ring()) i,j,k = Q.gens() # Go top-left to bottom-right (reading order), converting every # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1 # quaternion block. elements = [] - for l in range(n/4): - for m in range(n/4): - submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed( - M[4*l:4*l+4,4*m:4*m+4] ) + for l in range(n/d): + for m in range(n/d): + submat = ComplexMatrixEJA.real_unembed( + M[d*l:d*l+d,d*m:d*m+d] ) if submat[0,0] != submat[1,1].conjugate(): raise ValueError('bad on-diagonal submatrix') if submat[0,1] != -submat[1,0].conjugate(): @@ -2067,42 +2166,10 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): z += submat[0,1].imag()*k elements.append(z) - return matrix(Q, n/4, elements) - - - @classmethod - def trace_inner_product(cls,X,Y): - """ - Compute a matrix 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.to_matrix() - sage: Ye = y.to_matrix() - sage: X = QuaternionHermitianEJA.real_unembed(Xe) - sage: Y = QuaternionHermitianEJA.real_unembed(Ye) - sage: expected = (X*Y).trace().coefficient_tuple()[0] - sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye) - sage: actual == expected - True - - """ - return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4 + return matrix(Q, n/d, elements) -class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, - QuaternionMatrixEuclideanJordanAlgebra): +class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the @@ -2193,23 +2260,39 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, # * The diagonal will (as a result) be real. # S = [] + Eij = matrix.zero(Q,n) for i in range(n): for j in range(i+1): - Eij = matrix(Q, n, lambda k,l: k==i and l==j) + # "build" E_ij + Eij[i,j] = 1 if i == j: Sij = cls.real_embed(Eij) S.append(Sij) else: # The second, third, and fourth ones have a minus # because they're conjugated. - Sij_real = cls.real_embed(Eij + Eij.transpose()) + # Eij = Eij + Eij.transpose() + Eij[j,i] = 1 + Sij_real = cls.real_embed(Eij) S.append(Sij_real) - Sij_I = cls.real_embed(I*Eij - I*Eij.transpose()) + # Eij = I*(Eij - Eij.transpose()) + Eij[i,j] = I + Eij[j,i] = -I + Sij_I = cls.real_embed(Eij) S.append(Sij_I) - Sij_J = cls.real_embed(J*Eij - J*Eij.transpose()) + # Eij = J*(Eij - Eij.transpose()) + Eij[i,j] = J + Eij[j,i] = -J + Sij_J = cls.real_embed(Eij) S.append(Sij_J) - Sij_K = cls.real_embed(K*Eij - K*Eij.transpose()) + # Eij = K*(Eij - Eij.transpose()) + Eij[i,j] = K + Eij[j,i] = -K + Sij_K = cls.real_embed(Eij) S.append(Sij_K) + Eij[j,i] = 0 + # "erase" E_ij + Eij[i,j] = 0 # Since we embedded these, we can drop back to the "field" that we # started with instead of the quaternion algebra "Q". @@ -2225,8 +2308,13 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, self.jordan_product, self.trace_inner_product, **kwargs) + # TODO: this could be factored out somehow, but is left here + # because the MatrixEJA is not presently a subclass of the + # FDEJA class that defines rank() and one(). self.rank.set_cache(n) - # TODO: cache one()! + idV = matrix.identity(ZZ, self.dimension_over_reals()*n) + self.one.set_cache(self(idV)) + @staticmethod def _max_random_instance_size(): @@ -2244,7 +2332,7 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, return cls(n, **kwargs) -class HadamardEJA(ConcreteEuclideanJordanAlgebra): +class HadamardEJA(ConcreteEJA): """ Return the Euclidean Jordan Algebra corresponding to the set `R^n` under the Hadamard product. @@ -2285,11 +2373,16 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra): """ def __init__(self, n, **kwargs): - def jordan_product(x,y): - P = x.parent() - return P(tuple( xi*yi for (xi,yi) in zip(x,y) )) - def inner_product(x,y): - return x.inner_product(y) + if n == 0: + jordan_product = lambda x,y: x + inner_product = lambda x,y: x + else: + def jordan_product(x,y): + P = x.parent() + return P( xi*yi for (xi,yi) in zip(x,y) ) + + def inner_product(x,y): + return (x.T*y)[0,0] # New defaults for keyword arguments. Don't orthonormalize # because our basis is already orthonormal with respect to our @@ -2300,12 +2393,12 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra): if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - - standard_basis = FreeModule(ZZ, n).basis() - super(HadamardEJA, self).__init__(standard_basis, - jordan_product, - inner_product, - **kwargs) + column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() ) + super().__init__(column_basis, + jordan_product, + inner_product, + associative=True, + **kwargs) self.rank.set_cache(n) if n == 0: @@ -2329,7 +2422,7 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra): return cls(n, **kwargs) -class BilinearFormEJA(ConcreteEuclideanJordanAlgebra): +class BilinearFormEJA(ConcreteEJA): r""" The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` with the half-trace inner product and jordan product ``x*y = @@ -2409,31 +2502,38 @@ class BilinearFormEJA(ConcreteEuclideanJordanAlgebra): ....: for j in range(n-1) ] sage: actual == expected True + """ def __init__(self, B, **kwargs): - if not B.is_positive_definite(): - raise ValueError("bilinear form is not positive-definite") + # The matrix "B" is supplied by the user in most cases, + # so it makes sense to check whether or not its positive- + # definite unless we are specifically asked not to... + if ("check_axioms" not in kwargs) or kwargs["check_axioms"]: + if not B.is_positive_definite(): + raise ValueError("bilinear form is not positive-definite") + + # However, all of the other data for this EJA is computed + # by us in manner that guarantees the axioms are + # satisfied. So, again, unless we are specifically asked to + # verify things, we'll skip the rest of the checks. + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False def inner_product(x,y): - return (B*x).inner_product(y) + return (y.T*B*x)[0,0] def jordan_product(x,y): P = x.parent() - x0 = x[0] - xbar = x[1:] - y0 = y[0] - ybar = y[1:] - z0 = inner_product(x,y) + x0 = x[0,0] + xbar = x[1:,0] + y0 = y[0,0] + ybar = y[1:,0] + z0 = inner_product(y,x) zbar = y0*xbar + x0*ybar - return P((z0,) + tuple(zbar)) - - # We know this is a valid EJA, but will double-check - # if the user passes check_axioms=True. - if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + return P([z0] + zbar.list()) n = B.nrows() - standard_basis = FreeModule(ZZ, n).basis() - super(BilinearFormEJA, self).__init__(standard_basis, + column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() ) + super(BilinearFormEJA, self).__init__(column_basis, jordan_product, inner_product, **kwargs) @@ -2561,7 +2661,7 @@ class JordanSpinEJA(BilinearFormEJA): return cls(n, **kwargs) -class TrivialEJA(ConcreteEuclideanJordanAlgebra): +class TrivialEJA(ConcreteEJA): """ The trivial Euclidean Jordan algebra consisting of only a zero element. @@ -2614,155 +2714,335 @@ class TrivialEJA(ConcreteEuclideanJordanAlgebra): # inappropriate for us. return cls(**kwargs) -class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): + +class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct, + FiniteDimensionalEJA): r""" - The external (orthogonal) direct sum of two other Euclidean Jordan - algebras. Essentially the Cartesian product of its two factors. - Every Euclidean Jordan algebra decomposes into an orthogonal - direct sum of simple Euclidean Jordan algebras, so no generality - is lost by providing only this construction. + The external (orthogonal) direct sum of two or more Euclidean + Jordan algebras. Every Euclidean Jordan algebra decomposes into an + orthogonal direct sum of simple Euclidean Jordan algebras which is + then isometric to a Cartesian product, so no generality is lost by + providing only this construction. SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, + ....: CartesianProductEJA, ....: HadamardEJA, - ....: RealSymmetricEJA, - ....: DirectSumEJA) + ....: JordanSpinEJA, + ....: RealSymmetricEJA) - EXAMPLES:: + EXAMPLES: + + The Jordan product is inherited from our factors and implemented by + our CombinatorialFreeModule Cartesian product superclass:: + sage: set_random_seed() sage: J1 = HadamardEJA(2) - sage: J2 = RealSymmetricEJA(3) - sage: J = DirectSumEJA(J1,J2) - sage: J.dimension() - 8 + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: x,y = J.random_elements(2) + sage: x*y in J + True + + The ability to retrieve the original factors is implemented by our + CombinatorialFreeModule Cartesian product superclass:: + + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = JordanSpinEJA(3, field=QQ) + sage: J = cartesian_product([J1,J2]) + sage: J.cartesian_factors() + (Euclidean Jordan algebra of dimension 2 over Rational Field, + Euclidean Jordan algebra of dimension 3 over Rational Field) + + You can provide more than two factors:: + + sage: J1 = HadamardEJA(2) + sage: J2 = JordanSpinEJA(3) + sage: J3 = RealSymmetricEJA(3) + sage: cartesian_product([J1,J2,J3]) + Euclidean Jordan algebra of dimension 2 over Algebraic Real + Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic + Real Field (+) Euclidean Jordan algebra of dimension 6 over + Algebraic Real Field + + Rank is additive on a Cartesian product:: + + sage: J1 = HadamardEJA(1) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J1.rank.clear_cache() + sage: J2.rank.clear_cache() + sage: J.rank.clear_cache() + sage: J.rank() + 3 + sage: J.rank() == J1.rank() + J2.rank() + True + + The same rank computation works over the rationals, with whatever + basis you like:: + + sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False) + sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False) + sage: J = cartesian_product([J1,J2]) + sage: J1.rank.clear_cache() + sage: J2.rank.clear_cache() + sage: J.rank.clear_cache() sage: J.rank() - 5 + 3 + sage: J.rank() == J1.rank() + J2.rank() + True + + The product algebra will be associative if and only if all of its + components are associative:: + + sage: J1 = HadamardEJA(2) + sage: J1.is_associative() + True + sage: J2 = HadamardEJA(3) + sage: J2.is_associative() + True + sage: J3 = RealSymmetricEJA(3) + sage: J3.is_associative() + False + sage: CP1 = cartesian_product([J1,J2]) + sage: CP1.is_associative() + True + sage: CP2 = cartesian_product([J1,J3]) + sage: CP2.is_associative() + False TESTS: - The external direct sum construction is only valid when the two factors - have the same base ring; an error is raised otherwise:: + All factors must share the same base field:: - sage: set_random_seed() - sage: J1 = random_eja(field=AA) - sage: J2 = random_eja(field=QQ,orthonormalize=False) - sage: J = DirectSumEJA(J1,J2) + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = RealSymmetricEJA(2) + sage: CartesianProductEJA((J1,J2)) Traceback (most recent call last): ... - ValueError: algebras must share the same base field + ValueError: all factors must share the same base field + + The cached unit element is the same one that would be computed:: + + sage: set_random_seed() # long time + sage: J1 = random_eja() # long time + sage: J2 = random_eja() # long time + sage: J = cartesian_product([J1,J2]) # long time + sage: actual = J.one() # long time + sage: J.one.clear_cache() # long time + sage: expected = J.one() # long time + sage: actual == expected # long time + True """ - def __init__(self, J1, J2, **kwargs): - if J1.base_ring() != J2.base_ring(): - raise ValueError("algebras must share the same base field") - field = J1.base_ring() - - self._factors = (J1, J2) - n1 = J1.dimension() - n2 = J2.dimension() - n = n1+n2 - V = VectorSpace(field, n) - mult_table = [ [ V.zero() for j in range(i+1) ] - for i in range(n) ] - for i in range(n1): - for j in range(i+1): - p = (J1.monomial(i)*J1.monomial(j)).to_vector() - mult_table[i][j] = V(p.list() + [field.zero()]*n2) + def __init__(self, algebras, **kwargs): + CombinatorialFreeModule_CartesianProduct.__init__(self, + algebras, + **kwargs) + field = algebras[0].base_ring() + if not all( J.base_ring() == field for J in algebras ): + raise ValueError("all factors must share the same base field") + + associative = all( m.is_associative() for m in algebras ) + + # The definition of matrix_space() and self.basis() relies + # only on the stuff in the CFM_CartesianProduct class, which + # we've already initialized. + Js = self.cartesian_factors() + m = len(Js) + MS = self.matrix_space() + basis = tuple( + MS(tuple( self.cartesian_projection(i)(b).to_matrix() + for i in range(m) )) + for b in self.basis() + ) - for i in range(n2): - for j in range(i+1): - p = (J2.monomial(i)*J2.monomial(j)).to_vector() - mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list()) - - # TODO: build the IP table here from the two constituent IP - # matrices (it'll be block diagonal, I think). - ip_table = [ [ field.zero() for j in range(i+1) ] - for i in range(n) ] - super(DirectSumEJA, self).__init__(field, - mult_table, - ip_table, - check_axioms=False, - **kwargs) - self.rank.set_cache(J1.rank() + J2.rank()) - - - def factors(self): + # Define jordan/inner products that operate on that matrix_basis. + def jordan_product(x,y): + return MS(tuple( + (Js[i](x[i])*Js[i](y[i])).to_matrix() for i in range(m) + )) + + def inner_product(x, y): + return sum( + Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m) + ) + + # There's no need to check the field since it already came + # from an EJA. Likewise the axioms are guaranteed to be + # satisfied, unless the guy writing this class sucks. + # + # If you want the basis to be orthonormalized, orthonormalize + # the factors. + FiniteDimensionalEJA.__init__(self, + basis, + jordan_product, + inner_product, + field=field, + orthonormalize=False, + associative=associative, + cartesian_product=True, + check_field=False, + check_axioms=False) + + ones = tuple(J.one() for J in algebras) + self.one.set_cache(self._cartesian_product_of_elements(ones)) + self.rank.set_cache(sum(J.rank() for J in algebras)) + + def matrix_space(self): r""" - Return the pair of this algebra's factors. + Return the space that our matrix basis lives in as a Cartesian + product. SETUP:: sage: from mjo.eja.eja_algebra import (HadamardEJA, - ....: JordanSpinEJA, - ....: DirectSumEJA) + ....: RealSymmetricEJA) EXAMPLES:: - sage: J1 = HadamardEJA(2, field=QQ) - sage: J2 = JordanSpinEJA(3, field=QQ) - sage: J = DirectSumEJA(J1,J2) - sage: J.factors() - (Euclidean Jordan algebra of dimension 2 over Rational Field, - Euclidean Jordan algebra of dimension 3 over Rational Field) + sage: J1 = HadamardEJA(1) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.matrix_space() + The Cartesian product of (Full MatrixSpace of 1 by 1 dense + matrices over Algebraic Real Field, Full MatrixSpace of 2 + by 2 dense matrices over Algebraic Real Field) """ - return self._factors + from sage.categories.cartesian_product import cartesian_product + return cartesian_product( [J.matrix_space() + for J in self.cartesian_factors()] ) - def projections(self): + @cached_method + def cartesian_projection(self, i): r""" - Return a pair of projections onto this algebra's factors. - SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: ComplexHermitianEJA, - ....: DirectSumEJA) + sage: from mjo.eja.eja_algebra import (random_eja, + ....: JordanSpinEJA, + ....: HadamardEJA, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA) - EXAMPLES:: + EXAMPLES: + + The projection morphisms are Euclidean Jordan algebra + operators:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.cartesian_projection(0) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [1 0 0 0 0] + [0 1 0 0 0] + Domain: Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field (+) Euclidean Jordan algebra of dimension 3 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field + sage: J.cartesian_projection(1) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [0 0 1 0 0] + [0 0 0 1 0] + [0 0 0 0 1] + Domain: Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field (+) Euclidean Jordan algebra of dimension 3 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic + Real Field + + The projections work the way you'd expect on the vector + representation of an element:: sage: J1 = JordanSpinEJA(2) sage: J2 = ComplexHermitianEJA(2) - sage: J = DirectSumEJA(J1,J2) - sage: (pi_left, pi_right) = J.projections() - sage: J.one().to_vector() - (1, 0, 1, 0, 0, 1) + sage: J = cartesian_product([J1,J2]) + sage: pi_left = J.cartesian_projection(0) + sage: pi_right = J.cartesian_projection(1) sage: pi_left(J.one()).to_vector() (1, 0) sage: pi_right(J.one()).to_vector() (1, 0, 0, 1) + sage: J.one().to_vector() + (1, 0, 1, 0, 0, 1) + + TESTS: + + The answer never changes:: + + sage: set_random_seed() + sage: J1 = random_eja() + sage: J2 = random_eja() + sage: J = cartesian_product([J1,J2]) + sage: P0 = J.cartesian_projection(0) + sage: P1 = J.cartesian_projection(0) + sage: P0 == P1 + True """ - (J1,J2) = self.factors() - m = J1.dimension() - n = J2.dimension() - V_basis = self.vector_space().basis() - # Need to specify the dimensions explicitly so that we don't - # wind up with a zero-by-zero matrix when we want e.g. a - # zero-by-two matrix (important for composing things). - P1 = matrix(self.base_ring(), m, m+n, V_basis[:m]) - P2 = matrix(self.base_ring(), n, m+n, V_basis[m:]) - pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1) - pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2) - return (pi_left, pi_right) - - def inclusions(self): - r""" - Return the pair of inclusion maps from our factors into us. + Ji = self.cartesian_factors()[i] + # Requires the fix on Trac 31421/31422 to work! + Pi = super().cartesian_projection(i) + return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix()) + @cached_method + def cartesian_embedding(self, i): + r""" SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, ....: JordanSpinEJA, - ....: RealSymmetricEJA, - ....: DirectSumEJA) + ....: HadamardEJA, + ....: RealSymmetricEJA) - EXAMPLES:: + EXAMPLES: + + The embedding morphisms are Euclidean Jordan algebra + operators:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.cartesian_embedding(0) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [1 0] + [0 1] + [0 0] + [0 0] + [0 0] + Domain: Euclidean Jordan algebra of dimension 2 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 2 over + Algebraic Real Field (+) Euclidean Jordan algebra of + dimension 3 over Algebraic Real Field + sage: J.cartesian_embedding(1) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [0 0 0] + [0 0 0] + [1 0 0] + [0 1 0] + [0 0 1] + Domain: Euclidean Jordan algebra of dimension 3 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 2 over + Algebraic Real Field (+) Euclidean Jordan algebra of + dimension 3 over Algebraic Real Field + + The embeddings work the way you'd expect on the vector + representation of an element:: sage: J1 = JordanSpinEJA(3) sage: J2 = RealSymmetricEJA(2) - sage: J = DirectSumEJA(J1,J2) - sage: (iota_left, iota_right) = J.inclusions() + sage: J = cartesian_product([J1,J2]) + sage: iota_left = J.cartesian_embedding(0) + sage: iota_right = J.cartesian_embedding(1) sage: iota_left(J1.zero()) == J.zero() True sage: iota_right(J2.zero()) == J.zero() @@ -2780,6 +3060,17 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): TESTS: + The answer never changes:: + + sage: set_random_seed() + sage: J1 = random_eja() + sage: J2 = random_eja() + sage: J = cartesian_product([J1,J2]) + sage: E0 = J.cartesian_embedding(0) + sage: E1 = J.cartesian_embedding(0) + sage: E0 == E1 + True + Composing a projection with the corresponding inclusion should produce the identity map, and mismatching them should produce the zero map:: @@ -2787,9 +3078,11 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: set_random_seed() sage: J1 = random_eja() sage: J2 = random_eja() - sage: J = DirectSumEJA(J1,J2) - sage: (iota_left, iota_right) = J.inclusions() - sage: (pi_left, pi_right) = J.projections() + sage: J = cartesian_product([J1,J2]) + sage: iota_left = J.cartesian_embedding(0) + sage: iota_right = J.cartesian_embedding(1) + sage: pi_left = J.cartesian_projection(0) + sage: pi_right = J.cartesian_projection(1) sage: pi_left*iota_left == J1.one().operator() True sage: pi_right*iota_right == J2.one().operator() @@ -2800,58 +3093,49 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): True """ - (J1,J2) = self.factors() - m = J1.dimension() - n = J2.dimension() - V_basis = self.vector_space().basis() - # Need to specify the dimensions explicitly so that we don't - # wind up with a zero-by-zero matrix when we want e.g. a - # two-by-zero matrix (important for composing things). - I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m]) - I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:]) - iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1) - iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2) - return (iota_left, iota_right) + Ji = self.cartesian_factors()[i] + # Requires the fix on Trac 31421/31422 to work! + Ei = super().cartesian_embedding(i) + return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix()) - def inner_product(self, x, y): + + def _element_constructor_(self, elt): r""" - The standard Cartesian inner-product. + Construct an element of this algebra from an ordered tuple. - We project ``x`` and ``y`` onto our factors, and add up the - inner-products from the subalgebras. + We just apply the element constructor from each of our factors + to the corresponding component of the tuple, and package up + the result. SETUP:: - sage: from mjo.eja.eja_algebra import (HadamardEJA, - ....: QuaternionHermitianEJA, - ....: DirectSumEJA) - - EXAMPLE:: - - sage: J1 = HadamardEJA(3,field=QQ) - sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False) - sage: J = DirectSumEJA(J1,J2) - sage: x1 = J1.one() - sage: x2 = x1 - sage: y1 = J2.one() - sage: y2 = y1 - sage: x1.inner_product(x2) - 3 - sage: y1.inner_product(y2) - 2 - sage: J.one().inner_product(J.one()) - 5 + ....: RealSymmetricEJA) + EXAMPLES:: + + sage: J1 = HadamardEJA(3) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) ) + e(0, 1) + e(1, 2) """ - (pi_left, pi_right) = self.projections() - x1 = pi_left(x) - x2 = pi_right(x) - y1 = pi_left(y) - y2 = pi_right(y) + m = len(self.cartesian_factors()) + try: + z = tuple( self.cartesian_factors()[i](elt[i]) for i in range(m) ) + return self._cartesian_product_of_elements(z) + except: + raise ValueError("not an element of this algebra") - return (x1.inner_product(y1) + x2.inner_product(y2)) + Element = CartesianProductEJAElement +FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA -random_eja = ConcreteEuclideanJordanAlgebra.random_instance +random_eja = ConcreteEJA.random_instance +#def random_eja(*args, **kwargs): +# from sage.categories.cartesian_product import cartesian_product +# J1 = HadamardEJA(1, **kwargs) +# J2 = RealSymmetricEJA(2, **kwargs) +# J = cartesian_product([J1,J2]) +# return J