X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=8005412e65da8b6192072420dffe794c7975b5d0;hb=f292189be8e5a79f9ae2b80ddaff76460e0d14c2;hp=ea62da8ca796031bd99382c5886bb1091d145b78;hpb=de451def1161cc9dfefcfc125523029881cb160a;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index ea62da8..8005412 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,7 +31,8 @@ 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 FiniteDimensionalEJAElement +from mjo.eja.eja_element import (CartesianProductEJAElement, + FiniteDimensionalEJAElement) from mjo.eja.eja_operator import FiniteDimensionalEJAOperator from mjo.eja.eja_utils import _mat2vec @@ -39,7 +42,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): INPUT: - - basis -- a tuple of basis elements in their matrix form. + - 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 @@ -60,6 +71,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): field=AA, orthonormalize=True, associative=False, + cartesian_product=False, check_field=True, check_axioms=True, prefix='e'): @@ -73,7 +85,12 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # If the basis given to us wasn't over the field that it's # supposed to be over, fix that. Or, you know, crash. - basis = tuple( b.change_ring(field) for b in basis ) + 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. @@ -96,6 +113,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): 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. @@ -113,10 +132,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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: - # Works on both column and square matrices... - degree = len(basis[0].list()) + degree = len(flatten(basis[0])) # Build an ambient space that fits our matrix basis when # written out as "long vectors." @@ -130,10 +156,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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(b.list()) for b in basis ) + deortho_vector_basis = tuple( V(flatten(b)) for b in basis ) from mjo.eja.eja_utils import gram_schmidt - basis = gram_schmidt(basis, inner_product) + basis = tuple(gram_schmidt(basis, inner_product)) # Save the (possibly orthonormalized) matrix basis for # later... @@ -142,7 +168,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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(b.list()) for b in basis ) + vector_basis = tuple( V(flatten(b)) for b in basis ) W = V.span_of_basis( vector_basis, check=check_axioms) if orthonormalize: @@ -158,7 +184,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # Now we actually compute the multiplication and inner-product # tables/matrices using the possibly-orthonormalized basis. - self._inner_product_matrix = matrix.zero(field, n) + self._inner_product_matrix = matrix.identity(field, n) self._multiplication_table = [ [0 for j in range(i+1)] for i in range(n) ] @@ -171,15 +197,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): q_i = basis[i] q_j = basis[j] - elt = jordan_product(q_i, q_j) - ip = inner_product(q_i, q_j) - # The jordan product returns a matrixy answer, so we # have to convert it to the algebra coordinates. - elt = W.coordinate_vector(V(elt.list())) + elt = jordan_product(q_i, q_j) + elt = W.coordinate_vector(V(flatten(elt))) self._multiplication_table[i][j] = self.from_vector(elt) - self._inner_product_matrix[i,j] = ip - self._inner_product_matrix[j,i] = ip + + 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() @@ -274,22 +305,32 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: y = J.random_element() sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) True + """ B = self._inner_product_matrix return (B*x.to_vector()).inner_product(y.to_vector()) - def _is_commutative(self): + def is_associative(self): r""" - Whether or not this algebra's multiplication table is commutative. + Return whether or not this algebra's Jordan product is associative. + + SETUP:: + + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + + EXAMPLES:: + + 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 - 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()) ) + return "Associative" in self.category().axioms() def _is_jordanian(self): r""" @@ -298,13 +339,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): 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 + :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)) + return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j]) == - (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j)) + (self.gens()[i])*((self.gens()[i]**2)*self.gens()[j]) for i in range(self.dimension()) for j in range(self.dimension()) ) @@ -315,7 +356,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): This method should of course always return ``True``, unless this algebra was constructed with ``check_axioms=False`` and - passed an invalid multiplication table. + passed an invalid Jordan or inner-product. """ # Used to check whether or not something is zero in an inexact @@ -326,9 +367,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): 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) + x = self.gens()[i] + y = self.gens()[j] + z = self.gens()[k] diff = (x*y).inner_product(z) - x.inner_product(y*z) if self.base_ring().is_exact(): @@ -651,8 +692,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # And to each subsequent row, prepend an entry that belongs to # the left-side "header column." - M += [ [self.monomial(i)] + [ self.product_on_basis(i,j) - for j in range(n) ] + M += [ [self.gens()[i]] + [ self.product_on_basis(i,j) + for j in range(n) ] for i in range(n) ] return table(M, header_row=True, header_column=True, frame=True) @@ -682,7 +723,7 @@ class FiniteDimensionalEJA(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. @@ -732,7 +773,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): if self.is_trivial(): return MatrixSpace(self.base_ring(), 0) else: - return self._matrix_basis[0].matrix_space() + return self.matrix_basis()[0].parent() @cached_method @@ -745,23 +786,57 @@ class FiniteDimensionalEJA(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 + + :: + + 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:: + 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() @@ -769,6 +844,27 @@ class FiniteDimensionalEJA(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:: @@ -780,6 +876,15 @@ class FiniteDimensionalEJA(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. @@ -924,14 +1029,12 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): if not c.is_idempotent(): raise ValueError("element is not idempotent: %s" % c) - 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 = FiniteDimensionalEJASubalgebra(self, ()) + trivial = self.subalgebra(()) J0 = trivial # eigenvalue zero J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half J1 = trivial # eigenvalue one @@ -941,9 +1044,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): J5 = eigspace else: gens = tuple( self.from_vector(b) for b in eigspace.basis() ) - subalg = FiniteDimensionalEJASubalgebra(self, - gens, - check_axioms=False) + subalg = self.subalgebra(gens, check_axioms=False) if eigval == 0: J0 = subalg elif eigval == 1: @@ -1032,6 +1133,21 @@ class FiniteDimensionalEJA(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() @@ -1041,7 +1157,7 @@ class FiniteDimensionalEJA(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) @@ -1067,10 +1183,17 @@ class FiniteDimensionalEJA(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): @@ -1131,7 +1254,7 @@ class FiniteDimensionalEJA(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 @@ -1140,6 +1263,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): return len(self._charpoly_coefficients()) + def subalgebra(self, basis, **kwargs): + r""" + Create a subalgebra of this algebra from the given basis. + + This is a simple wrapper around a subalgebra class constructor + that can be overridden in subclasses. + """ + from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra + return FiniteDimensionalEJASubalgebra(self, basis, **kwargs) + + def vector_space(self): """ Return the vector space that underlies this algebra. @@ -1187,9 +1321,7 @@ class RationalBasisEJA(FiniteDimensionalEJA): jordan_product, inner_product, field=AA, - orthonormalize=True, check_field=True, - check_axioms=True, **kwargs): if check_field: @@ -1198,6 +1330,7 @@ class RationalBasisEJA(FiniteDimensionalEJA): if not all( all(b_i in QQ for b_i in b.list()) for b in basis ): raise TypeError("basis not rational") + self._rational_algebra = None if field is not QQ: # There's no point in constructing the extra algebra if this # one is already rational. @@ -1212,15 +1345,13 @@ class RationalBasisEJA(FiniteDimensionalEJA): field=QQ, orthonormalize=False, check_field=False, - check_axioms=False, - **kwargs) + check_axioms=False) super().__init__(basis, jordan_product, inner_product, field=field, check_field=check_field, - check_axioms=check_axioms, **kwargs) @cached_method @@ -1260,7 +1391,14 @@ class RationalBasisEJA(FiniteDimensionalEJA): 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 @@ -1594,6 +1732,38 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): 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 dimension_over_reals(): return 2 @@ -1648,9 +1818,10 @@ class ComplexMatrixEJA(MatrixEJA): 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) @@ -1689,26 +1860,7 @@ class ComplexMatrixEJA(MatrixEJA): super(ComplexMatrixEJA,cls).real_unembed(M) n = ZZ(M.nrows()) d = cls.dimension_over_reals() - - # 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() - - # 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. - F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) + F = cls.complex_extension(M.base_ring()) i = F.gen() # Go top-left to bottom-right (reading order), converting every @@ -1804,7 +1956,6 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): 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 @@ -1822,18 +1973,27 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): # * 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". @@ -1869,6 +2029,25 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): return cls(n, **kwargs) 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 dimension_over_reals(): return 4 @@ -1973,8 +2152,7 @@ class QuaternionMatrixEJA(MatrixEJA): # 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 @@ -2089,23 +2267,39 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): # * 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". @@ -2207,7 +2401,11 @@ class HadamardEJA(ConcreteEJA): if "check_axioms" not in kwargs: kwargs["check_axioms"] = False column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() ) - super().__init__(column_basis, jordan_product, inner_product, **kwargs) + super().__init__(column_basis, + jordan_product, + inner_product, + associative=True, + **kwargs) self.rank.set_cache(n) if n == 0: @@ -2311,31 +2509,38 @@ class BilinearFormEJA(ConcreteEJA): ....: 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) @@ -2516,155 +2721,335 @@ class TrivialEJA(ConcreteEJA): # inappropriate for us. return cls(**kwargs) -class DirectSumEJA(ConcreteEJA): + +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 = FiniteDimensionalEJAOperator(self,J1,P1) - pi_right = FiniteDimensionalEJAOperator(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() @@ -2682,6 +3067,17 @@ class DirectSumEJA(ConcreteEJA): 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:: @@ -2689,9 +3085,11 @@ class DirectSumEJA(ConcreteEJA): 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() @@ -2702,58 +3100,49 @@ class DirectSumEJA(ConcreteEJA): 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 = FiniteDimensionalEJAOperator(J1,self,I1) - iota_right = FiniteDimensionalEJAOperator(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 = 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