X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=e32bc24cbb2741fca9de2285de703b7adf1497af;hb=3c7644ecfe369b6f83aa707b87d7a1f9aa246e27;hp=14666c2f2a4d47b2f1a5cde1448dcb3b778aa72a;hpb=8207ab8350d947b557bf6682da9e3c3ffe638523;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 14666c2..e32bc24 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -119,10 +119,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): The ``field`` we're given must be real with ``check_field=True``:: - sage: JordanSpinEJA(2,QQbar) + 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 The multiplication table must be square with ``check_axioms=True``:: @@ -173,16 +175,27 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): 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") + 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") + self._matrix_basis = matrix_basis if category is None: @@ -214,14 +227,18 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): 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. - self._inner_product_matrix = matrix(field, inner_product_table) - self._inner_product_matrix._cache = {'hermitian': False} + 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 check_axioms: if not self._is_jordanian(): @@ -276,6 +293,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: x = J.random_element() sage: J(x.to_vector().column()) == x True + """ msg = "not an element of this algebra" if elt == 0: @@ -298,8 +316,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # element's ring because the basis space might be an algebraic # closure whereas the base ring of the 3-by-3 identity matrix # could be QQ instead of QQbar. + # + # We pass check=False because the matrix basis is "guaranteed" + # to be linearly independent... right? Ha ha. V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols()) - W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() ) + W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()), + check=False) try: coords = W.coordinate_vector(_mat2vec(elt)) @@ -378,7 +400,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # Used to check whether or not something is zero in an inexact # ring. This number is sufficient to allow the construction of - # QuaternionHermitianEJA(2, RDF) with check_axioms=True. + # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True. epsilon = 1e-16 for i in range(self.dimension()): @@ -475,7 +497,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = HadamardEJA(2) sage: J.coordinate_polynomial_ring() Multivariate Polynomial Ring in X1, X2... - sage: J = RealSymmetricEJA(3,QQ) + sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False) sage: J.coordinate_polynomial_ring() Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6... @@ -1117,7 +1139,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Element = FiniteDimensionalEuclideanJordanAlgebraElement -class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra): +class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): r""" New class for algebras whose supplied basis elements have all rational entries. @@ -1140,16 +1162,26 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge """ def __init__(self, - field, basis, jordan_product, inner_product, + field=AA, orthonormalize=True, prefix='e', category=None, check_field=True, check_axioms=True): + if check_field: + # Abuse the check_field parameter to check that the entries of + # out basis (in ambient coordinates) are in the field QQ. + 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 @@ -1160,6 +1192,7 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge 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: @@ -1173,27 +1206,97 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge # 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._deortho_multiplication_table = None - self._deortho_inner_product_table = None + 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 - vector_basis = gram_schmidt(vector_basis, inner_product) - W = V.span_of_basis( vector_basis ) + 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: - from mjo.eja.eja_utils import _vec2mat basis = tuple( map(_vec2mat,basis) ) - W = V.span_of_basis( vector_basis ) + W = V.span_of_basis( vector_basis, check=check_axioms) - mult_table = [ [0 for i in range(n)] for j in range(n) ] - ip_table = [ [0 for i in range(n)] for j in range(n) ] + # 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 ) - # Note: the Jordan and inner- products are defined in terms + # 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): @@ -1216,9 +1319,12 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge elt = W.coordinate_vector(elt) mult_table[i][j] = elt - mult_table[j][i] = elt ip_table[i][j] = ip - ip_table[j][i] = 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: @@ -1235,26 +1341,13 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge check_field, check_axioms) -class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): - r""" - Algebras whose basis consists of vectors with rational - entries. Equivalently, algebras whose multiplication tables - contain only rational coefficients. - - When an EJA has a basis that can be made rational, we can speed up - the computation of its characteristic polynomial by doing it over - ``QQ``. All of the named EJA constructors that we provide fall - into this category. - """ @cached_method def _charpoly_coefficients(self): r""" - Override the parent method with something that tries to compute - over a faster (non-extension) field. - SETUP:: - sage: from mjo.eja.eja_algebra import JordanSpinEJA + sage: from mjo.eja.eja_algebra import (BilinearFormEJA, + ....: JordanSpinEJA) EXAMPLES: @@ -1272,29 +1365,31 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr Algebraic Real Field """ - if self.base_ring() is QQ: + if self.base_ring() is QQ or self._rational_algebra is None: # There's no need to construct *another* algebra over the - # rationals if this one is already over the rationals. + # 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() - mult_table = tuple( - tuple(map(lambda x: x.to_vector(), ls)) - for ls in self._multiplication_table - ) - # Do the computation over the rationals. The answer will be - # the same, because our basis coordinates are (essentially) - # rational. - J = FiniteDimensionalEuclideanJordanAlgebra(QQ, - mult_table, - check_field=False, - check_axioms=False) - a = J._charpoly_coefficients() - return tuple(map(lambda x: x.change_ring(self.base_ring()), a)) + # the same, because all we've done is a change of basis. + # Then, change back from QQ to our real base ring + 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 + # deorthonormalized ones. + R = self.coordinate_polynomial_ring() + from sage.modules.free_module_element import vector + X = vector(R, R.gens()) + BX = self._deortho_matrix*X -class ConcreteEuclideanJordanAlgebra: + 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): r""" A class for the Euclidean Jordan algebras that we know by name. @@ -1347,7 +1442,7 @@ class ConcreteEuclideanJordanAlgebra: raise NotImplementedError @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, *args, **kwargs): """ Return a random instance of this type of algebra. @@ -1355,127 +1450,14 @@ class ConcreteEuclideanJordanAlgebra: """ from sage.misc.prandom import choice eja_class = choice(cls.__subclasses__()) - return eja_class.random_instance(field) - - -class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): - - def __init__(self, field, basis, normalize_basis=True, **kwargs): - """ - Compared to the superclass constructor, we take a basis instead of - a multiplication table because the latter can be computed in terms - of the former when the product is known (like it is here). - """ - # Used in this class's fast _charpoly_coefficients() override. - self._basis_normalizers = None - - # We're going to loop through this a few times, so now's a good - # time to ensure that it isn't a generator expression. - basis = tuple(basis) - - algebra_dim = len(basis) - degree = 0 # size of the matrices - if algebra_dim > 0: - degree = basis[0].nrows() - - if algebra_dim > 1 and normalize_basis: - # We'll need sqrt(2) to normalize the basis, and this - # winds up in the multiplication table, so the whole - # algebra needs to be over the field extension. - R = PolynomialRing(field, 'z') - z = R.gen() - p = z**2 - 2 - if p.is_irreducible(): - field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt()) - basis = tuple( s.change_ring(field) for s in basis ) - self._basis_normalizers = tuple( - ~(self.matrix_inner_product(s,s).sqrt()) for s in basis ) - basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers)) - - # Now compute the multiplication and inner product tables. - # We have to do this *after* normalizing the basis, because - # scaling affects the answers. - V = VectorSpace(field, degree**2) - W = V.span_of_basis( _mat2vec(s) for s in basis ) - mult_table = [[W.zero() for j in range(algebra_dim)] - for i in range(algebra_dim)] - ip_table = [[field.zero() for j in range(algebra_dim)] - for i in range(algebra_dim)] - for i in range(algebra_dim): - for j in range(algebra_dim): - mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 - mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) - - try: - # HACK: ignore the error here if we don't need the - # inner product (as is the case when we construct - # a dummy QQ-algebra for fast charpoly coefficients. - ip_table[i][j] = self.matrix_inner_product(basis[i], - basis[j]) - except: - pass - - super(MatrixEuclideanJordanAlgebra, self).__init__(field, - mult_table, - ip_table, - matrix_basis=basis, - **kwargs) - - if algebra_dim == 0: - self.one.set_cache(self.zero()) - else: - n = basis[0].nrows() - # The identity wrt (A,B) -> (AB + BA)/2 is independent of the - # details of this algebra. - self.one.set_cache(self(matrix.identity(field,n))) - - @cached_method - def _charpoly_coefficients(self): - r""" - Override the parent method with something that tries to compute - over a faster (non-extension) field. - """ - if self._basis_normalizers is None or self.base_ring() is QQ: - # We didn't normalize, or the basis we started with had - # entries in a nice field already. Just compute the thing. - return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients() - - basis = ( (b/n) for (b,n) in zip(self.matrix_basis(), - self._basis_normalizers) ) - - # Do this over the rationals and convert back at the end. - # Only works because we know the entries of the basis are - # integers. The argument ``check_axioms=False`` is required - # because the trace inner-product method for this - # class is a stub and can't actually be checked. - J = MatrixEuclideanJordanAlgebra(QQ, - basis, - normalize_basis=False, - check_field=False, - check_axioms=False) - a = J._charpoly_coefficients() - - # Unfortunately, changing the basis does change the - # coefficients of the characteristic polynomial, but since - # these are really the coefficients of the "characteristic - # polynomial of" function, everything is still nice and - # unevaluated. It's therefore "obvious" how scaling the - # basis affects the coordinate variables X1, X2, et - # cetera. Scaling the first basis vector up by "n" adds a - # factor of 1/n into every "X1" term, for example. So here - # we simply undo the basis_normalizer scaling that we - # performed earlier. - # - # The a[0] access here is safe because trivial algebras - # won't have any basis normalizers and therefore won't - # make it to this "else" branch. - XS = a[0].parent().gens() - subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i] - for i in range(len(XS)) } - return tuple( a_i.subs(subs_dict) for a_i in a ) + # These all bubble up to the RationalBasisEuclideanJordanAlgebra + # superclass constructor, so any (kw)args valid there are also + # valid here. + return eja_class.random_instance(*args, **kwargs) +class MatrixEuclideanJordanAlgebra: @staticmethod def real_embed(M): """ @@ -1500,8 +1482,12 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): """ raise NotImplementedError + @staticmethod + def jordan_product(X,Y): + return (X*Y + Y*X)/2 + @classmethod - def matrix_inner_product(cls,X,Y): + def trace_inner_product(cls,X,Y): Xu = cls.real_unembed(X) Yu = cls.real_unembed(Y) tr = (Xu*Yu).trace() @@ -1534,8 +1520,8 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return M -class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, + RealMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1558,9 +1544,9 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, In theory, our "field" can be any subfield of the reals:: - sage: RealSymmetricEJA(2, RDF) + sage: RealSymmetricEJA(2, field=RDF) Euclidean Jordan algebra of dimension 3 over Real Double Field - sage: RealSymmetricEJA(2, RR) + sage: RealSymmetricEJA(2, field=RR) Euclidean Jordan algebra of dimension 3 over Real Field with 53 bits of precision @@ -1601,7 +1587,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, """ @classmethod - def _denormalized_basis(cls, n, field): + def _denormalized_basis(cls, n): """ Return a basis for the space of real symmetric n-by-n matrices. @@ -1613,7 +1599,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: B = RealSymmetricEJA._denormalized_basis(n,QQ) + sage: B = RealSymmetricEJA._denormalized_basis(n) sage: all( M.is_symmetric() for M in B) True @@ -1623,13 +1609,13 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, S = [] for i in range(n): for j in range(i+1): - Eij = matrix(field, n, lambda k,l: k==i and l==j) + Eij = matrix(ZZ, n, lambda k,l: k==i and l==j) if i == j: Sij = Eij else: Sij = Eij + Eij.transpose() S.append(Sij) - return S + return tuple(S) @staticmethod @@ -1637,20 +1623,20 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, return 4 # Dimension 10 @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) + return cls(n, **kwargs) - def __init__(self, n, field=AA, **kwargs): - basis = self._denormalized_basis(n, field) - super(RealSymmetricEJA, self).__init__(field, - basis, - check_axioms=False, + def __init__(self, n, **kwargs): + super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n), + self.jordan_product, + self.trace_inner_product, **kwargs) self.rank.set_cache(n) + self.one.set_cache(self(matrix.identity(ZZ,n))) class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @@ -1781,7 +1767,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @classmethod - def matrix_inner_product(cls,X,Y): + def trace_inner_product(cls,X,Y): """ Compute a matrix inner product in this algebra directly from its real embedding. @@ -1803,16 +1789,16 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: X = ComplexHermitianEJA.real_unembed(Xe) sage: Y = ComplexHermitianEJA.real_unembed(Ye) sage: expected = (X*Y).trace().real() - sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye) + sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye) sage: actual == expected True """ - return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2 + return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2 -class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, + ComplexMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -1827,9 +1813,9 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, In theory, our "field" can be any subfield of the reals:: - sage: ComplexHermitianEJA(2, RDF) + sage: ComplexHermitianEJA(2, field=RDF) Euclidean Jordan algebra of dimension 4 over Real Double Field - sage: ComplexHermitianEJA(2, RR) + sage: ComplexHermitianEJA(2, field=RR) Euclidean Jordan algebra of dimension 4 over Real Field with 53 bits of precision @@ -1871,7 +1857,7 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, """ @classmethod - def _denormalized_basis(cls, n, field): + def _denormalized_basis(cls, n): """ Returns a basis for the space of complex Hermitian n-by-n matrices. @@ -1890,15 +1876,16 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: field = QuadraticField(2, 'sqrt2') - sage: B = ComplexHermitianEJA._denormalized_basis(n, field) + sage: B = ComplexHermitianEJA._denormalized_basis(n) sage: all( M.is_symmetric() for M in B) True """ + field = ZZ R = PolynomialRing(field, 'z') z = R.gen() F = field.extension(z**2 + 1, 'I') - I = F.gen() + I = F.gen(1) # This is like the symmetric case, but we need to be careful: # @@ -1921,28 +1908,28 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, # Since we embedded these, we can drop back to the "field" that we # started with instead of the complex extension "F". - return ( s.change_ring(field) for s in S ) + return tuple( s.change_ring(field) for s in S ) - def __init__(self, n, field=AA, **kwargs): - basis = self._denormalized_basis(n,field) - super(ComplexHermitianEJA,self).__init__(field, - basis, - check_axioms=False, - **kwargs) + def __init__(self, n, **kwargs): + super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n), + self.jordan_product, + self.trace_inner_product, + **kwargs) self.rank.set_cache(n) + # TODO: pre-cache the identity! @staticmethod def _max_random_instance_size(): return 3 # Dimension 9 @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) + return cls(n, **kwargs) class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @staticmethod @@ -2076,7 +2063,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @classmethod - def matrix_inner_product(cls,X,Y): + def trace_inner_product(cls,X,Y): """ Compute a matrix inner product in this algebra directly from its real embedding. @@ -2098,16 +2085,16 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: X = QuaternionHermitianEJA.real_unembed(Xe) sage: Y = QuaternionHermitianEJA.real_unembed(Ye) sage: expected = (X*Y).trace().coefficient_tuple()[0] - sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye) + sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye) sage: actual == expected True """ - return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4 + return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4 -class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, + QuaternionMatrixEuclideanJordanAlgebra): r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the @@ -2122,9 +2109,9 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, In theory, our "field" can be any subfield of the reals:: - sage: QuaternionHermitianEJA(2, RDF) + sage: QuaternionHermitianEJA(2, field=RDF) Euclidean Jordan algebra of dimension 6 over Real Double Field - sage: QuaternionHermitianEJA(2, RR) + sage: QuaternionHermitianEJA(2, field=RR) Euclidean Jordan algebra of dimension 6 over Real Field with 53 bits of precision @@ -2165,7 +2152,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, """ @classmethod - def _denormalized_basis(cls, n, field): + def _denormalized_basis(cls, n): """ Returns a basis for the space of quaternion Hermitian n-by-n matrices. @@ -2183,11 +2170,12 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, sage: set_random_seed() sage: n = ZZ.random_element(1,5) - sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ) + sage: B = QuaternionHermitianEJA._denormalized_basis(n) sage: all( M.is_symmetric() for M in B ) True """ + field = ZZ Q = QuaternionAlgebra(QQ,-1,-1) I,J,K = Q.gens() @@ -2217,16 +2205,16 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, # Since we embedded these, we can drop back to the "field" that we # started with instead of the quaternion algebra "Q". - return ( s.change_ring(field) for s in S ) + return tuple( s.change_ring(field) for s in S ) - def __init__(self, n, field=AA, **kwargs): - basis = self._denormalized_basis(n,field) - super(QuaternionHermitianEJA,self).__init__(field, - basis, - check_axioms=False, - **kwargs) + def __init__(self, n, **kwargs): + super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n), + self.jordan_product, + self.trace_inner_product, + **kwargs) self.rank.set_cache(n) + # TODO: cache one()! @staticmethod def _max_random_instance_size(): @@ -2236,16 +2224,15 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, return 2 # Dimension 6 @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) + return cls(n, **kwargs) -class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg, - ConcreteEuclideanJordanAlgebra): +class HadamardEJA(ConcreteEuclideanJordanAlgebra): """ Return the Euclidean Jordan Algebra corresponding to the set `R^n` under the Hadamard product. @@ -2285,19 +2272,25 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg, (r0, r1, r2) """ - def __init__(self, n, field=AA, **kwargs): - V = VectorSpace(field, n) - basis = V.basis() - + def __init__(self, n, **kwargs): def jordan_product(x,y): - return V([ xi*yi for (xi,yi) in zip(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) - super(HadamardEJA, self).__init__(field, - basis, + # Don't orthonormalize because our basis is already + # orthonormal with respect to our inner-product. + if not 'orthonormalize' in kwargs: + kwargs['orthonormalize'] = False + + # But also don't pass check_field=False here, because the user + # can pass in a field! + standard_basis = FreeModule(ZZ, n).basis() + super(HadamardEJA, self).__init__(standard_basis, jordan_product, inner_product, + check_axioms=False, **kwargs) self.rank.set_cache(n) @@ -2314,16 +2307,15 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg, return 5 @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) + return cls(n, **kwargs) -class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg, - ConcreteEuclideanJordanAlgebra): +class BilinearFormEJA(ConcreteEuclideanJordanAlgebra): 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 = @@ -2404,27 +2396,26 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg, sage: actual == expected True """ - def __init__(self, B, field=AA, **kwargs): + def __init__(self, B, **kwargs): if not B.is_positive_definite(): raise ValueError("bilinear form is not positive-definite") - n = B.nrows() - V = VectorSpace(field, n) - def inner_product(x,y): return (B*x).inner_product(y) 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) zbar = y0*xbar + x0*ybar - return V([z0] + zbar.list()) + return P((z0,) + tuple(zbar)) - super(BilinearFormEJA, self).__init__(field, - V.basis(), + n = B.nrows() + standard_basis = FreeModule(ZZ, n).basis() + super(BilinearFormEJA, self).__init__(standard_basis, jordan_product, inner_product, **kwargs) @@ -2447,28 +2438,28 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg, return 5 @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, **kwargs): """ Return a random instance of this algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) if n.is_zero(): - B = matrix.identity(field, n) - return cls(B, field, **kwargs) + B = matrix.identity(ZZ, n) + return cls(B, **kwargs) - B11 = matrix.identity(field,1) - M = matrix.random(field, n-1) - I = matrix.identity(field, n-1) - alpha = field.zero() + B11 = matrix.identity(ZZ, 1) + M = matrix.random(ZZ, n-1) + I = matrix.identity(ZZ, n-1) + alpha = ZZ.zero() while alpha.is_zero(): - alpha = field.random_element().abs() + alpha = ZZ.random_element().abs() B22 = M.transpose()*M + alpha*I from sage.matrix.special import block_matrix B = block_matrix(2,2, [ [B11, ZZ(0) ], [ZZ(0), B22 ] ]) - return cls(B, field, **kwargs) + return cls(B, **kwargs) class JordanSpinEJA(BilinearFormEJA): @@ -2521,11 +2512,21 @@ class JordanSpinEJA(BilinearFormEJA): True """ - def __init__(self, n, field=AA, **kwargs): - # This is a special case of the BilinearFormEJA with the identity - # matrix as its bilinear form. - B = matrix.identity(field, n) - super(JordanSpinEJA, self).__init__(B, field, **kwargs) + def __init__(self, n, **kwargs): + # This is a special case of the BilinearFormEJA with the + # identity matrix as its bilinear form. + B = matrix.identity(ZZ, n) + + # Don't orthonormalize because our basis is already + # orthonormal with respect to our inner-product. + if not 'orthonormalize' in kwargs: + kwargs['orthonormalize'] = False + + # But also don't pass check_field=False here, because the user + # can pass in a field! + super(JordanSpinEJA, self).__init__(B, + check_axioms=False, + **kwargs) @staticmethod def _max_random_instance_size(): @@ -2535,18 +2536,17 @@ class JordanSpinEJA(BilinearFormEJA): return 5 @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. Needed here to override the implementation for ``BilinearFormEJA``. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) + return cls(n, **kwargs) -class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, - ConcreteEuclideanJordanAlgebra): +class TrivialEJA(ConcreteEuclideanJordanAlgebra): """ The trivial Euclidean Jordan algebra consisting of only a zero element. @@ -2575,13 +2575,13 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, 0 """ - def __init__(self, field=AA, **kwargs): - mult_table = [] - ip_table = [] - super(TrivialEJA, self).__init__(field, - mult_table, - ip_table, - check_axioms=False, + def __init__(self, **kwargs): + jordan_product = lambda x,y: x + inner_product = lambda x,y: 0 + basis = () + super(TrivialEJA, self).__init__(basis, + jordan_product, + inner_product, **kwargs) # The rank is zero using my definition, namely the dimension of the # largest subalgebra generated by any element. @@ -2589,10 +2589,10 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, self.one.set_cache( self.zero() ) @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, **kwargs): # We don't take a "size" argument so the superclass method is # inappropriate for us. - return cls(field, **kwargs) + return cls(**kwargs) class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): r""" @@ -2625,8 +2625,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): have the same base ring; an error is raised otherwise:: sage: set_random_seed() - sage: J1 = random_eja(AA) - sage: J2 = random_eja(QQ) + sage: J1 = random_eja(field=AA) + sage: J2 = random_eja(field=QQ,orthonormalize=False) sage: J = DirectSumEJA(J1,J2) Traceback (most recent call last): ... @@ -2643,21 +2643,22 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): n2 = J2.dimension() n = n1+n2 V = VectorSpace(field, n) - mult_table = [ [ V.zero() for j in range(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(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(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 = None + ip_table = [ [ field.zero() for j in range(i+1) ] + for i in range(n) ] super(DirectSumEJA, self).__init__(field, mult_table, ip_table, @@ -2678,8 +2679,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): EXAMPLES:: - sage: J1 = HadamardEJA(2,QQ) - sage: J2 = JordanSpinEJA(3,QQ) + 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, @@ -2808,8 +2809,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): EXAMPLE:: - sage: J1 = HadamardEJA(3,QQ) - sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False) + 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