X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=c35bf256453cfa9a3758c034a47a043745eb0c2b;hb=7dc61e259969f86d4f31706fcfe411e1677389f2;hp=1250fbda7981aba5474af0a3d2615c969bc9d1e1;hpb=4c8f9aac69d1cb4097b60b10e5b198b6372ec55e;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 1250fbd..c35bf25 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -113,60 +113,36 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # we see in things like x = 1*e1 + 2*e2. vector_basis = basis - from sage.structure.element import is_Matrix - basis_is_matrices = False - degree = 0 if n > 0: - if is_Matrix(basis[0]): - if basis[0].is_square(): - # TODO: this ugly is_square() hack works around the problem - # of passing to_matrix()ed vectors in as the basis from a - # subalgebra. They aren't REALLY matrices, at least not of - # the type that we assume here... Ugh. - basis_is_matrices = True - from mjo.eja.eja_utils import _vec2mat - vector_basis = tuple( map(_mat2vec,basis) ) - degree = basis[0].nrows()**2 - else: - # convert from column matrices to vectors, yuck - basis = tuple( map(_mat2vec,basis) ) - vector_basis = basis - degree = basis[0].degree() - else: - degree = basis[0].degree() + # Works on both column and square matrices... + degree = len(basis[0].list()) - # Build an ambient space that fits... + # Build an ambient space that fits our matrix basis when + # written out as "long vectors." V = VectorSpace(field, degree) - # 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 + # The matrix that will hole the orthonormal -> unorthonormal + # coordinate transformation. 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 + # 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 ) - if basis_is_matrices: - basis = tuple( map(_vec2mat,basis) ) + from mjo.eja.eja_utils import gram_schmidt + basis = tuple(gram_schmidt(basis, inner_product)) - # Save the matrix "basis" for later... this is the last time we'll - # reference it in this constructor. - if basis_is_matrices: - self._matrix_basis = basis - else: - MS = MatrixSpace(self.base_ring(), degree, 1) - self._matrix_basis = tuple( MS(b) for b in basis ) + # Save the (possibly orthonormalized) matrix basis for + # later... + self._matrix_basis = basis - # Now create the vector space for the algebra... + # 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 ) W = V.span_of_basis( vector_basis, check=check_axioms) if orthonormalize: @@ -182,40 +158,33 @@ 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._multiplication_table = [ [0 for j in range(i+1)] for i in range(n) ] + self._inner_product_matrix = matrix.identity(field, n) + self._multiplication_table = [ [0 for j in range(i+1)] + for i in range(n) ] - print("vector_basis:") - print(vector_basis) # 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) + 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) - 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) - - # TODO: the jordan product turns things back into - # matrices here even if they're supposed to be - # vectors. ugh. Can we get rid of vectors all together - # please? - elt = W.coordinate_vector(elt) + elt = W.coordinate_vector(V(elt.list())) 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() @@ -351,7 +320,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 @@ -768,7 +737,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 @@ -781,23 +750,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() @@ -805,6 +808,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:: @@ -816,6 +840,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. @@ -1068,6 +1101,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() @@ -1103,10 +1151,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): @@ -1167,7 +1222,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 @@ -1223,9 +1278,7 @@ class RationalBasisEJA(FiniteDimensionalEJA): jordan_product, inner_product, field=AA, - orthonormalize=True, check_field=True, - check_axioms=True, **kwargs): if check_field: @@ -1234,6 +1287,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. @@ -1248,15 +1302,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 @@ -1296,7 +1348,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 @@ -1630,6 +1689,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 @@ -1684,9 +1775,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) @@ -1725,26 +1817,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 @@ -1840,7 +1913,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 @@ -1858,18 +1930,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". @@ -1905,6 +1986,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 @@ -2009,8 +2109,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 @@ -2125,23 +2224,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". @@ -2222,11 +2337,16 @@ class HadamardEJA(ConcreteEJA): """ 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 @@ -2237,12 +2357,8 @@ class HadamardEJA(ConcreteEJA): 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, **kwargs) self.rank.set_cache(n) if n == 0: @@ -2346,31 +2462,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) @@ -2551,7 +2674,8 @@ class TrivialEJA(ConcreteEJA): # inappropriate for us. return cls(**kwargs) -class DirectSumEJA(ConcreteEJA): + +class DirectSumEJA(FiniteDimensionalEJA): r""" The external (orthogonal) direct sum of two other Euclidean Jordan algebras. Essentially the Cartesian product of its two factors. @@ -2575,6 +2699,10 @@ class DirectSumEJA(ConcreteEJA): 8 sage: J.rank() 5 + sage: J.matrix_space() + The Cartesian product of (Full MatrixSpace of 2 by 1 dense matrices + over Algebraic Real Field, Full MatrixSpace of 3 by 3 dense matrices + over Algebraic Real Field) TESTS: @@ -2595,199 +2723,213 @@ class DirectSumEJA(ConcreteEJA): 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) - - 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): - r""" - Return the pair of this algebra's factors. - - SETUP:: + M = J1.matrix_space().cartesian_product(J2.matrix_space()) + self._cartprod_algebra = J1.cartesian_product(J2) - sage: from mjo.eja.eja_algebra import (HadamardEJA, - ....: JordanSpinEJA, - ....: DirectSumEJA) + self._matrix_basis = tuple( [M((a,0)) for a in J1.matrix_basis()] + + [M((0,b)) for b in J2.matrix_basis()] ) - EXAMPLES:: + n = len(self._matrix_basis) + self._sets = None + CombinatorialFreeModule.__init__( + self, + field, + range(n), + category=self._cartprod_algebra.category(), + bracket=False, + **kwargs) + self.rank.set_cache(J1.rank() + J2.rank()) - 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) - """ - return self._factors - def projections(self): + def product(self,x,y): r""" - Return a pair of projections onto this algebra's factors. - SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: ComplexHermitianEJA, ....: DirectSumEJA) - EXAMPLES:: - - 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: pi_left(J.one()).to_vector() - (1, 0) - sage: pi_right(J.one()).to_vector() - (1, 0, 0, 1) - - """ - (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. - - SETUP:: - - sage: from mjo.eja.eja_algebra import (random_eja, - ....: JordanSpinEJA, - ....: RealSymmetricEJA, - ....: DirectSumEJA) - - EXAMPLES:: - - sage: J1 = JordanSpinEJA(3) - sage: J2 = RealSymmetricEJA(2) - sage: J = DirectSumEJA(J1,J2) - sage: (iota_left, iota_right) = J.inclusions() - sage: iota_left(J1.zero()) == J.zero() - True - sage: iota_right(J2.zero()) == J.zero() - True - sage: J1.one().to_vector() - (1, 0, 0) - sage: iota_left(J1.one()).to_vector() - (1, 0, 0, 0, 0, 0) - sage: J2.one().to_vector() - (1, 0, 1) - sage: iota_right(J2.one()).to_vector() - (0, 0, 0, 1, 0, 1) - sage: J.one().to_vector() - (1, 0, 0, 1, 0, 1) - - TESTS: - - Composing a projection with the corresponding inclusion should - produce the identity map, and mismatching them should produce - the zero map:: + TESTS:: sage: set_random_seed() - sage: J1 = random_eja() - sage: J2 = random_eja() + sage: J1 = JordanSpinEJA(3, field=QQ) + sage: J2 = ComplexHermitianEJA(2, field=QQ, orthonormalize=False) sage: J = DirectSumEJA(J1,J2) - sage: (iota_left, iota_right) = J.inclusions() - sage: (pi_left, pi_right) = J.projections() - sage: pi_left*iota_left == J1.one().operator() - True - sage: pi_right*iota_right == J2.one().operator() - True - sage: (pi_left*iota_right).is_zero() - True - sage: (pi_right*iota_left).is_zero() + sage: J.random_element()*J.random_element() in J 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) + xv = self._cartprod_algebra.from_vector(x.to_vector()) + yv = self._cartprod_algebra.from_vector(y.to_vector()) + return self.from_vector((xv*yv).to_vector()) - def inner_product(self, x, y): - r""" - The standard Cartesian inner-product. - We project ``x`` and ``y`` onto our factors, and add up the - inner-products from the subalgebras. + def cartesian_factors(self): + r""" + Return the pair of this algebra's factors. SETUP:: - sage: from mjo.eja.eja_algebra import (HadamardEJA, - ....: QuaternionHermitianEJA, + ....: JordanSpinEJA, ....: DirectSumEJA) - EXAMPLE:: + EXAMPLES:: - sage: J1 = HadamardEJA(3,field=QQ) - sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = JordanSpinEJA(3, field=QQ) 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 + sage: J.cartesian_factors() + (Euclidean Jordan algebra of dimension 2 over Rational Field, + Euclidean Jordan algebra of dimension 3 over Rational Field) """ - (pi_left, pi_right) = self.projections() - x1 = pi_left(x) - x2 = pi_right(x) - y1 = pi_left(y) - y2 = pi_right(y) - - return (x1.inner_product(y1) + x2.inner_product(y2)) + return self._cartprod_algebra.cartesian_factors() + + +# def projections(self): +# r""" +# Return a pair of projections onto this algebra's factors. + +# SETUP:: + +# sage: from mjo.eja.eja_algebra import (JordanSpinEJA, +# ....: ComplexHermitianEJA, +# ....: DirectSumEJA) + +# EXAMPLES:: + +# 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: pi_left(J.one()).to_vector() +# (1, 0) +# sage: pi_right(J.one()).to_vector() +# (1, 0, 0, 1) + +# """ +# (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. + +# SETUP:: + +# sage: from mjo.eja.eja_algebra import (random_eja, +# ....: JordanSpinEJA, +# ....: RealSymmetricEJA, +# ....: DirectSumEJA) + +# EXAMPLES:: + +# sage: J1 = JordanSpinEJA(3) +# sage: J2 = RealSymmetricEJA(2) +# sage: J = DirectSumEJA(J1,J2) +# sage: (iota_left, iota_right) = J.inclusions() +# sage: iota_left(J1.zero()) == J.zero() +# True +# sage: iota_right(J2.zero()) == J.zero() +# True +# sage: J1.one().to_vector() +# (1, 0, 0) +# sage: iota_left(J1.one()).to_vector() +# (1, 0, 0, 0, 0, 0) +# sage: J2.one().to_vector() +# (1, 0, 1) +# sage: iota_right(J2.one()).to_vector() +# (0, 0, 0, 1, 0, 1) +# sage: J.one().to_vector() +# (1, 0, 0, 1, 0, 1) + +# TESTS: + +# Composing a projection with the corresponding inclusion should +# produce the identity map, and mismatching them should produce +# the zero map:: + +# 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: pi_left*iota_left == J1.one().operator() +# True +# sage: pi_right*iota_right == J2.one().operator() +# True +# sage: (pi_left*iota_right).is_zero() +# True +# sage: (pi_right*iota_left).is_zero() +# 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) + +# def inner_product(self, x, y): +# r""" +# The standard Cartesian inner-product. + +# We project ``x`` and ``y`` onto our factors, and add up the +# inner-products from the subalgebras. + +# 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 + +# """ +# (pi_left, pi_right) = self.projections() +# x1 = pi_left(x) +# x2 = pi_right(x) +# y1 = pi_left(y) +# y2 = pi_right(y) + +# return (x1.inner_product(y1) + x2.inner_product(y2))