X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=9eba6699819956048747f648fd18bc80847fbbb6;hb=c9db5fa5d2f30833e4da85d32418048d34483e43;hp=44d83147bbeceb4acae52d573749545dec3c4b99;hpb=9708043704809263d2bd543de38c46b458b873cb;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 44d8314..9eba669 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -1040,7 +1040,7 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge n = len(basis) vector_basis = basis - from sage.matrix.matrix import is_Matrix + from sage.structure.element import is_Matrix basis_is_matrices = False degree = 0 @@ -1054,30 +1054,56 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge V = VectorSpace(field, degree) - self._deorthonormalization_matrix = matrix.identity(field,n) + # Compute this from "Q" (obtained from Gram-Schmidt) below as + # R = Q.solve_right(A), where the rows of "Q" are the + # orthonormalized vector_basis and and the rows of "A" are the + # original vector_basis. + self._deorthonormalization_matrix = None + if orthonormalize: + from mjo.eja.eja_utils import gram_schmidt A = matrix(field, vector_basis) - # uh oh, this is only the "usual" inner product - Q,R = A.gram_schmidt(orthonormal=True) - self._deorthonormalization_matrix = R.inverse().transpose() - vector_basis = Q.rows() + vector_basis = gram_schmidt(vector_basis, inner_product) W = V.span_of_basis( vector_basis ) + Q = matrix(field, vector_basis) + # A = QR <==> A.T == R.T*Q.T + # So, Q.solve_right() is equivalent to the Q.T.solve_left() + # that we want. + self._deorthonormalization_matrix = Q.solve_right(A) + if basis_is_matrices: from mjo.eja.eja_utils import _vec2mat basis = tuple( map(_vec2mat,vector_basis) ) + W = V.span_of_basis( vector_basis ) + 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) ] + # 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): - # do another mat2vec because the multiplication - # table is in terms of vectors - elt = _mat2vec(jordan_product(basis[i],basis[j])) + # ortho basis w.r.t. ambient coords + q_i = vector_basis[i] + q_j = vector_basis[j] + + if basis_is_matrices: + q_i = _vec2mat(q_i) + q_j = _vec2mat(q_j) + + elt = jordan_product(q_i, q_j) + ip = inner_product(q_i, q_j) + + if basis_is_matrices: + # do another mat2vec because the multiplication + # table is in terms of vectors + elt = _mat2vec(elt) + elt = W.coordinate_vector(elt) mult_table[i][j] = elt mult_table[j][i] = elt - ip = inner_product(basis[i],basis[j]) ip_table[i][j] = ip ip_table[j][i] = ip @@ -1086,12 +1112,14 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge if basis_is_matrices: for m in basis: m.set_immutable() + else: + basis = tuple( x.column() for x in basis ) super().__init__(field, mult_table, prefix, category, - basis, + basis, # matrix basis check_field, check_axioms) @@ -2109,7 +2137,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, return cls(n, field, **kwargs) -class HadamardEJA(RationalBasisEuclideanJordanAlgebra, +class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg, ConcreteEuclideanJordanAlgebra): """ Return the Euclidean Jordan Algebra corresponding to the set @@ -2152,22 +2180,17 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra, """ def __init__(self, n, field=AA, **kwargs): V = VectorSpace(field, n) - mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ] - for i in range(n) ] + basis = V.basis() - # Inner products are real numbers and not algebra - # elements, so once we turn the algebra element - # into a vector in inner_product(), we never go - # back. As a result -- contrary to what we do with - # self._multiplication_table -- we store the inner - # product table as a plain old matrix and not as - # an algebra operator. - ip_table = matrix.identity(field,n) - self._inner_product_matrix = ip_table + def jordan_product(x,y): + return V([ xi*yi for (xi,yi) in zip(x,y) ]) + def inner_product(x,y): + return x.inner_product(y) super(HadamardEJA, self).__init__(field, - mult_table, - check_axioms=False, + basis, + jordan_product, + inner_product, **kwargs) self.rank.set_cache(n) @@ -2192,7 +2215,7 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra, return cls(n, field, **kwargs) -class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, +class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg, ConcreteEuclideanJordanAlgebra): r""" The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` @@ -2273,39 +2296,28 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, True """ def __init__(self, B, field=AA, **kwargs): - n = B.nrows() - if not B.is_positive_definite(): raise ValueError("bilinear form is not positive-definite") + n = B.nrows() V = VectorSpace(field, n) - mult_table = [[V.zero() for j in range(n)] for i in range(n)] - for i in range(n): - for j in range(n): - x = V.gen(i) - y = V.gen(j) - x0 = x[0] - xbar = x[1:] - y0 = y[0] - ybar = y[1:] - z0 = (B*x).inner_product(y) - zbar = y0*xbar + x0*ybar - z = V([z0] + zbar.list()) - mult_table[i][j] = z - - # Inner products are real numbers and not algebra - # elements, so once we turn the algebra element - # into a vector in inner_product(), we never go - # back. As a result -- contrary to what we do with - # self._multiplication_table -- we store the inner - # product table as a plain old matrix and not as - # an algebra operator. - ip_table = B - self._inner_product_matrix = ip_table + + def inner_product(x,y): + return (B*x).inner_product(y) + + def jordan_product(x,y): + 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()) super(BilinearFormEJA, self).__init__(field, - mult_table, - check_axioms=False, + V.basis(), + jordan_product, + inner_product, **kwargs) # The rank of this algebra is two, unless we're in a