X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feuclidean_jordan_algebra.py;h=a469bb0d7ba1302ccf30289c9b1c0c5518ed527f;hb=8f7ecef93d07a3f969d1c944c80b75a7052a4975;hp=b31322f9b1a2ee23a4451e249efd8734455d4595;hpb=644eb599855ddd51e8731104cc337eb9a4a33abe;p=sage.d.git diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py index b31322f..a469bb0 100644 --- a/mjo/eja/euclidean_jordan_algebra.py +++ b/mjo/eja/euclidean_jordan_algebra.py @@ -87,6 +87,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): return fmt.format(self.degree(), self.base_ring()) + def _a_regular_element(self): + """ + Guess a regular element. Needed to compute the basis for our + characteristic polynomial coefficients. + """ + gs = self.gens() + z = self.sum( (i+1)*gs[i] for i in range(len(gs)) ) + if not z.is_regular(): + raise ValueError("don't know a regular element") + return z + + + @cached_method + def _charpoly_basis_space(self): + """ + Return the vector space spanned by the basis used in our + characteristic polynomial coefficients. This is used not only to + compute those coefficients, but also any time we need to + evaluate the coefficients (like when we compute the trace or + determinant). + """ + z = self._a_regular_element() + V = z.vector().parent().ambient_vector_space() + V1 = V.span_of_basis( (z**k).vector() for k in range(self.rank()) ) + b = (V1.basis() + V1.complement().basis()) + return V.span_of_basis(b) + @cached_method def _charpoly_coeff(self, i): @@ -98,25 +125,38 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): store the trace/determinant (a_{r-1} and a_{0} respectively) separate from the entire characteristic polynomial. """ - (A_of_x, x) = self._charpoly_matrix() + (A_of_x, x, xr, detA) = self._charpoly_matrix_system() R = A_of_x.base_ring() - A_cols = A_of_x.columns() - A_cols[i] = (x**self.rank()).vector() - numerator = column_matrix(A_of_x.base_ring(), A_cols).det() - denominator = A_of_x.det() + if i >= self.rank(): + # Guaranteed by theory + return R.zero() + + # Danger: the in-place modification is done for performance + # reasons (reconstructing a matrix with huge polynomial + # entries is slow), but I don't know how cached_method works, + # so it's highly possible that we're modifying some global + # list variable by reference, here. In other words, you + # probably shouldn't call this method twice on the same + # algebra, at the same time, in two threads + Ai_orig = A_of_x.column(i) + A_of_x.set_column(i,xr) + numerator = A_of_x.det() + A_of_x.set_column(i,Ai_orig) # We're relying on the theory here to ensure that each a_i is # indeed back in R, and the added negative signs are to make # the whole charpoly expression sum to zero. - return R(-numerator/denominator) + return R(-numerator/detA) @cached_method - def _charpoly_matrix(self): + def _charpoly_matrix_system(self): """ Compute the matrix whose entries A_ij are polynomials in - X1,...,XN. This same matrix is used in more than one method and - it's not so fast to construct. + X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector + corresponding to `x^r` and the determinent of the matrix A = + [A_ij]. In other words, all of the fixed (cachable) data needed + to compute the coefficients of the characteristic polynomial. """ r = self.rank() n = self.dimension() @@ -130,11 +170,30 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): idmat = identity_matrix(J.base_ring(), n) + W = self._charpoly_basis_space() + W = W.change_ring(R.fraction_field()) + + # Starting with the standard coordinates x = (X1,X2,...,Xn) + # and then converting the entries to W-coordinates allows us + # to pass in the standard coordinates to the charpoly and get + # back the right answer. Specifically, with x = (X1,X2,...,Xn), + # we have + # + # W.coordinates(x^2) eval'd at (standard z-coords) + # = + # W-coords of (z^2) + # = + # W-coords of (standard coords of x^2 eval'd at std-coords of z) + # + # We want the middle equivalent thing in our matrix, but use + # the first equivalent thing instead so that we can pass in + # standard coordinates. x = J(vector(R, R.gens())) - l1 = [column_matrix((x**k).vector()) for k in range(r)] + l1 = [column_matrix(W.coordinates((x**k).vector())) for k in range(r)] l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)] A_of_x = block_matrix(R, 1, n, (l1 + l2)) - return (A_of_x, x) + xr = W.coordinates((x**r).vector()) + return (A_of_x, x, xr, A_of_x.det()) @cached_method @@ -640,14 +699,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): True """ + if not self.is_invertible(): + raise ValueError("element not invertible") + if self.parent().is_associative(): elt = FiniteDimensionalAlgebraElement(self.parent(), self) - return elt.inverse() - - # TODO: we can do better once the call to is_invertible() - # doesn't crash on irregular elements. - #if not self.is_invertible(): - # raise ValueError('element is not invertible') + # elt is in the right coordinates, but has the wrong class. + return self.parent()(elt.inverse().vector()) # We do this a little different than the usual recursive # call to a finite-dimensional algebra element, because we