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):
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()
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
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
sage: J = random_eja()
sage: x = J.random_element()
sage: y = J.random_element()
+ sage: Lx = x.operator_matrix()
+ sage: Lxx = (x*x).operator_matrix()
+ sage: Qx = x.quadratic_representation()
+ sage: Qy = y.quadratic_representation()
+ sage: Qxy = x.quadratic_representation(y)
+ sage: Qex = J.one().quadratic_representation(x)
+ sage: n = ZZ.random_element(10)
+ sage: Qxn = (x^n).quadratic_representation()
Property 1:
- sage: actual = x.quadratic_representation(y)
- sage: expected = ( (x+y).quadratic_representation()
- ....: -x.quadratic_representation()
- ....: -y.quadratic_representation() ) / 2
- sage: actual == expected
+ sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
True
Property 2:
sage: alpha = QQ.random_element()
- sage: actual = (alpha*x).quadratic_representation()
- sage: expected = (alpha^2)*x.quadratic_representation()
- sage: actual == expected
+ sage: (alpha*x).quadratic_representation() == (alpha^2)*Qx
+ True
+
+ Property 3:
+
+ sage: not x.is_invertible() or (
+ ....: Qx*x.inverse().vector() == x.vector() )
+ True
+
+ sage: not x.is_invertible() or (
+ ....: Qx.inverse()
+ ....: ==
+ ....: x.inverse().quadratic_representation() )
+ True
+
+ sage: Qxy*(J.one().vector()) == (x*y).vector()
+ True
+
+ Property 4:
+
+ sage: not x.is_invertible() or (
+ ....: x.quadratic_representation(x.inverse())*Qx
+ ....: == Qx*x.quadratic_representation(x.inverse()) )
+ True
+
+ sage: not x.is_invertible() or (
+ ....: x.quadratic_representation(x.inverse())*Qx
+ ....: ==
+ ....: 2*x.operator_matrix()*Qex - Qx )
+ True
+
+ sage: 2*x.operator_matrix()*Qex - Qx == Lxx
True
Property 5:
- sage: Qy = y.quadratic_representation()
- sage: actual = J(Qy*x.vector()).quadratic_representation()
- sage: expected = Qy*x.quadratic_representation()*Qy
- sage: actual == expected
+ sage: J(Qy*x.vector()).quadratic_representation() == Qy*Qx*Qy
True
Property 6:
- sage: k = ZZ.random_element(1,10)
- sage: actual = (x^k).quadratic_representation()
- sage: expected = (x.quadratic_representation())^k
- sage: actual == expected
+ sage: Qxn == (Qx)^n
+ True
+
+ Property 7:
+
+ sage: not x.is_invertible() or (
+ ....: Qx*x.inverse().operator_matrix() == Lx )
+ True
+
+ Property 8:
+
+ sage: not x.operator_commutes_with(y) or (
+ ....: J(Qx*y.vector())^n == J(Qxn*(y^n).vector()) )
True
"""
TESTS::
sage: set_random_seed()
- sage: J = RealCartesianProductEJA(5)
- sage: c = J.random_element().subalgebra_idempotent()
- sage: c^2 == c
- True
- sage: J = JordanSpinEJA(5)
- sage: c = J.random_element().subalgebra_idempotent()
+ sage: J = random_eja()
+ sage: x = J.random_element()
+ sage: while x.is_nilpotent():
+ ....: x = J.random_element()
+ sage: c = x.subalgebra_idempotent()
sage: c^2 == c
True
projects / sage.d.git / blobdiff