- (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
- R = A_of_x.base_ring()
-
- if i == self.rank():
- return R.one()
- 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/detA)
-
-
- @cached_method
- def _charpoly_matrix_system(self):
- """
- Compute the matrix whose entries A_ij are polynomials in
- 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()
-
- # Turn my vector space into a module so that "vectors" can
- # have multivatiate polynomial entries.
- names = tuple('X' + str(i) for i in range(1,n+1))
- R = PolynomialRing(self.base_ring(), names)
-
- # Using change_ring() on the parent's vector space doesn't work
- # here because, in a subalgebra, that vector space has a basis
- # and change_ring() tries to bring the basis along with it. And
- # that doesn't work unless the new ring is a PID, which it usually
- # won't be.
- V = FreeModule(R,n)
-
- # Now let x = (X1,X2,...,Xn) be the vector whose entries are
- # indeterminates...
- x = V(names)
-
- # And figure out the "left multiplication by x" matrix in
- # that setting.
- lmbx_cols = []
- monomial_matrices = [ self.monomial(i).operator().matrix()
- for i in range(n) ] # don't recompute these!
- for k in range(n):
- ek = self.monomial(k).to_vector()
- lmbx_cols.append(
- sum( x[i]*(monomial_matrices[i]*ek)
- for i in range(n) ) )
- Lx = matrix.column(R, lmbx_cols)
-
- # Now we can compute powers of x "symbolically"
- x_powers = [self.one().to_vector(), x]
- for d in range(2, r+1):
- x_powers.append( Lx*(x_powers[-1]) )
-
- idmat = matrix.identity(R, 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_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
- l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
- A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
- return (A_of_x, x, x_powers[r], A_of_x.det())
projects / sage.d.git / blobdiff