Coverage for local/lib/python2.7/site-packages/sage/rings/polynomial/real_roots.pyx : 92%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
""" Isolate Real Roots of Real Polynomials
AUTHOR:
- Carl Witty (2007-09-19): initial version
This is an implementation of real root isolation. That is, given a polynomial with exact real coefficients, we compute isolating intervals for the real roots of the polynomial. (Polynomials with integer, rational, or algebraic real coefficients are supported.)
We convert the polynomials into the Bernstein basis, and then use de Casteljau's algorithm and Descartes' rule of signs on the Bernstein basis polynomial (using interval arithmetic) to locate the roots. The algorithm is similar to that in "A Descartes Algorithm for Polynomials with Bit-Stream Coefficients", by Eigenwillig, Kettner, Krandick, Mehlhorn, Schmitt, and Wolpert, but has three crucial optimizations over the algorithm in that paper:
- Precision reduction: at certain points in the computation, we discard the low-order bits of the coefficients, widening the intervals.
- Degree reduction: at certain points in the computation, we find lower-degree polynomials that are approximately equal to our high-degree polynomial over the region of interest.
- When the intervals are too wide to continue (either because of a too-low initial precision, or because of precision or degree reduction), and we need to restart with higher precision, we recall which regions have already been proven not to have any roots and do not examine them again.
The best description of the algorithms used (other than this source code itself) is in the slides for my Sage Days 4 talk, currently available from https://wiki.sagemath.org/days4schedule . """
################################################################################ # Copyright (C) 2007 Carl Witty <cwitty@newtonlabs.com> # # Distributed under the terms of the GNU General Public License (GPL) # # http://www.gnu.org/licenses/ ################################################################################
# TODO: # These things would almost certainly improve the speed: # * Use Anatole Ruslanov's register tiling versions of de Casteljau's # algorithms when doing de Casteljau splitting at 1/2 in the integer case. # * Use a register tiling version of de Casteljau's algorithm for the # floating-point case...if you have n free FP registers (after allocating # registers to hold r and 1-r) you should be able to do about n rows at # once. # * Use vectorized FP instructions for de Casteljau's algorithm. # Using SSE2 to do 2 FP operations at once should be twice as fast. # * Faster degree elevation # Currently, when bernstein_expand expands from a degree-d1 to a # degree-d2 polynomial, it does O(d2^2 - d1^2) operations. I # have thought of (but not implemented) an approach that takes # O(d1*d2) operations, which should be much faster when d1 is 30 # and d2 is 1000. # Basically, the idea is to compute the d1 derivatives of the polynomial # at x=0. This can be done exactly, in O(d1^2) operations. Then # lift the d1'th derivative (a constant polynomial) to a polynomial of # formal degree d2-d1 (trivial; the coefficients are just d2-d1+1 copies # of the derivative). Then, compute the integral of that polynomial # d1 times, using the computed derivatives to give the constant offset. # (Computing the integral in Bernstein form involves divisions; these can # be done ahead of time, to the derivatives at x=0.) # The only tricky bit is tracking the error bounds and making sure you # use appropriate precisions at various parts of the computation.
# These things are interesting ideas that may or may not help: # * Partial degree reduction # Build a new subclass of interval_bernstein_polynomial, which # represents a polynomial as the sum of several # interval_bernstein_polynomial_integer (of different degrees). # Suppose you have a degree-1000 polynomial with 2000-bit coefficients. # When you try degree reduction, you get a degree-30 polynomial with # 2000-bit coefficients; but you get 300 bits of error. # Instead of just using the degree-30 polynomial and accepting the # loss of 300 bits of precision, you could represent the reduced polynomial # as the sum of a degree-30 polynomial with 2000-bit coefficients and # a degree-1000 polynomial with 300-bit coefficients. The combined # polynomial should still be much cheaper to work with than the original. # * Faster degree reduction # Perhaps somebody clever can figure out a way to do degree reduction # that doesn't involve matrix inversion, so that degree reduction can # target degrees larger than 30. # * Better heuristics # ** Better degree reduction heuristics # When should we degree reduce? What degrees should be targeted? # ** Better precision reduction heuristics # ** Better split point selection # Currently we attempt to split at 1/2, and then at a random offset. # We do not look at the current coefficients at all to try to select # a good split point. Would it be worthwhile to do so? (For instance, # it's a standard result that the polynomial in the region is bounded # by the convex hull of the coefficients.) # * Better initial transformation into Bernstein form # If a polynomial has roots which are very close to 0, or very large, # it may be better to rescale the roots (using a transformation like # p -> p(1000*x)). This can be done as part of the initial transformation # into Bernstein form. # If a polynomial has some roots which are very close to 0, and also # some roots which are very large, it might be better to isolate the roots # in two groups. # * More kinds of interval_bernstein_polynomial # Currently we have arbitrary-precision GMP coefficients or double-precision # floats. Would choices between these extremes help? (double-double? # quad-double? fixed-precision multi-word integer arithmetic?) # * Currently, some of our competitors are faster than this algorithm # for polynomials that are "easy" in some way. (This algorithm # seems to be vastly faster than everything else for "hard" polynomials.) # Maybe there is a way to start out using another algorithm # and then switch over to this algorithm once the polynomial # is determined to be hard. # * standard optimization # There's probably some more to be gained by just running the profiler # and optimizing individual functions. (To use more Cython features, for # instance.)
# Extra features: # * Specify a minimal island width: once islands are this narrow, # stop even if roots are not isolated. # * Do some sort of inexact version, for inexact input polynomials. # (For example, given a polynomial with (non-point) interval coefficients, # give a set of roots such that each root is guaranteed to be a root # of some polynomial in the set represented by the interval polynomial. # This may be vastly faster than the exact calculations carried out # by this algorithm! Is it enough faster to be faster than, say, # Pari's floating-point algorithms?) from __future__ import print_function, absolute_import
from sage.matrix.matrix_integer_dense cimport Matrix_integer_dense from sage.modules.vector_integer_dense cimport Vector_integer_dense from sage.modules.vector_real_double_dense cimport Vector_real_double_dense from sage.rings.integer cimport Integer from sage.rings.real_mpfr cimport RealNumber
cimport numpy
from libc.math cimport fabs, sqrt, ldexp, frexp
from sage.libs.gmp.mpz cimport * from sage.libs.gmp.mpq cimport * from sage.libs.mpfr cimport *
cdef class interval_bernstein_polynomial: r""" An interval_bernstein_polynomial is an approximation to an exact polynomial. This approximation is in the form of a Bernstein polynomial (a polynomial given as coefficients over a Bernstein basis) with interval coefficients.
The Bernstein basis of degree n over the region [a .. b] is the set of polynomials
.. MATH::
\binom{n}{k} (x-a)^k (b-x)^{n-k} / (b-a)^n
for `0 \le k \le n`.
A degree-n interval Bernstein polynomial P with its region [a .. b] can represent an exact polynomial p in two different ways: it can "contain" the polynomial or it can "bound" the polynomial.
We say that P contains p if, when p is represented as a degree-n Bernstein polynomial over [a .. b], its coefficients are contained in the corresponding interval coefficients of P. For instance, [0.9 .. 1.1]*x^2 (which is a degree-2 interval Bernstein polynomial over [0 .. 1]) contains x^2.
We say that P bounds p if, for all a <= x <= b, there exists a polynomial p' contained in P such that p(x) == p'(x). For instance, [0 .. 1]*x is a degree-1 interval Bernstein polynomial which bounds x^2 over [0 .. 1].
If P contains p, then P bounds p; but the converse is not necessarily true. In particular, if n < m, it is possible for a degree-n interval Bernstein polynomial to bound a degree-m polynomial; but it cannot contain the polynomial.
In the case where P bounds p, we maintain extra information, the "slope error". We say that P (over [a .. b]) bounds p with a slope error of E (where E is an interval) if there is a polynomial p' contained in P such that the derivative of (p - p') is bounded by E in the range [a .. b]. If P bounds p with a slope error of 0 then P contains p.
(Note that "contains" and "bounds" are not standard terminology; I just made them up.)
Interval Bernstein polynomials are useful in finding real roots because of the following properties:
- Given an exact real polynomial p, we can compute an interval Bernstein polynomial over an arbitrary region containing p.
- Given an interval Bernstein polynomial P over [a .. c], where a < b < c, we can compute interval Bernstein polynomials P1 over [a .. b] and P2 over [b .. c], where P1 and P2 contain (or bound) all polynomials that P contains (or bounds).
- Given a degree-n interval Bernstein polynomial P over [a .. b], and m < n, we can compute a degree-m interval Bernstein polynomial P' over [a .. b] that bounds all polynomials that P bounds.
- It is sometimes possible to prove that no polynomial bounded by P over [a .. b] has any roots in [a .. b]. (Roughly, this is possible when no polynomial contained by P has any complex roots near the line segment [a .. b], where "near" is defined relative to the length b-a.)
- It is sometimes possible to prove that every polynomial bounded by P over [a .. b] with slope error E has exactly one root in [a .. b]. (Roughly, this is possible when every polynomial contained by P over [a .. b] has exactly one root in [a .. b], there are no other complex roots near the line segment [a .. b], and every polynomial contained in P has a derivative which is bounded away from zero over [a .. b] by an amount which is large relative to E.)
- Starting from a sufficiently precise interval Bernstein polynomial, it is always possible to split it into polynomials which provably have 0 or 1 roots (as long as your original polynomial has no multiple real roots).
So a rough outline of a family of algorithms would be:
- Given a polynomial p, compute a region [a .. b] in which any real roots must lie. - Compute an interval Bernstein polynomial P containing p over [a .. b]. - Keep splitting P until you have isolated all the roots. Optionally, reduce the degree or the precision of the interval Bernstein polynomials at intermediate stages (to reduce computation time). If this seems not to be working, go back and try again with higher precision.
Obviously, there are many details to be worked out to turn this into a full algorithm, like:
- What initial precision is selected for computing P? - How do you decide when to reduce the degree of intermediate polynomials? - How do you decide when to reduce the precision of intermediate polynomials? - How do you decide where to split the interval Bernstein polynomial regions? - How do you decide when to give up and start over with higher precision?
Each set of answers to these questions gives a different algorithm (potentially with very different performance characteristics), but all of them can use this ``interval_bernstein_polynomial`` class as their basic building block.
To save computation time, all coefficients in an ``interval_bernstein_polynomial`` share the same interval width. (There is one exception: when creating an ``interval_bernstein_polynomial``, the first and last coefficients can be marked as "known positive" or "known negative". This has some of the same effect as having a (potentially) smaller interval width for these two coefficients, although it does not affect de Casteljau splitting.) To allow for widely varying coefficient magnitudes, all coefficients in an interval_bernstein_polynomial are scaled by `2^n` (where `n` may be positive, negative, or zero).
There are two representations for interval_bernstein_polynomials, integer and floating-point. These are the two subclasses of this class; ``interval_bernstein_polynomial`` itself is an abstract class.
``interval_bernstein_polynomial`` and its subclasses are not expected to be used outside this file. """
def variations(self): """ Consider a polynomial (written in either the normal power basis or the Bernstein basis). Take its list of coefficients, omitting zeroes. Count the number of positions in the list where the sign of one coefficient is opposite the sign of the next coefficient.
This count is the number of sign variations of the polynomial. According to Descartes' rule of signs, the number of real roots of the polynomial (counted with multiplicity) in a certain interval is always less than or equal to the number of sign variations, and the difference is always even. (If the polynomial is written in the power basis, the region is the positive reals; if the polynomial is written in the Bernstein basis over a particular region, then we count roots in that region.)
In particular, a polynomial with no sign variations has no real roots in the region, and a polynomial with one sign variation has one real root in the region.
In an interval Bernstein polynomial, we do not necessarily know the signs of the coefficients (if some of the coefficient intervals contain zero), so the polynomials contained by this interval polynomial may not all have the same number of sign variations. However, we can compute a range of possible numbers of sign variations.
This function returns the range, as a 2-tuple of integers. """
cdef void update_variations(self, interval_bernstein_polynomial bp1, interval_bernstein_polynomial bp2): """ Update the max_variations of bp1 and bp2 (which are assumed to be the result of splitting this polynomial).
If we knew the number of variations of self, bp1, and bp2 exactly, we would have self.variations == bp1.variations + bp2.variations + 2*n for some nonnegative integer n. Thus, we can use our information on min and max variations on self and bp1 (or bp2) to refine the range on bp2 (or bp1). """
def try_split(self, context ctx, logging_note): """ Try doing a de Casteljau split of this polynomial at 1/2, resulting in polynomials p1 and p2. If we see that the sign of this polynomial is determined at 1/2, then return (p1, p2, 1/2); otherwise, return None.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([50, 20, -90, -70, 200], error=5) sage: bp1, bp2, _ = bp.try_split(mk_context(), None) sage: bp1 <IBP: (50, 35, 0, -29, -31) + [0 .. 6) over [0 .. 1/2]> sage: bp2 <IBP: (-31, -33, -8, 65, 200) + [0 .. 6) over [1/2 .. 1]> sage: bp = mk_ibpf([0.5, 0.2, -0.9, -0.7, 0.99], neg_err=-0.1, pos_err=0.01) sage: bp1, bp2, _ = bp.try_split(mk_context(), None) sage: bp1 <IBP: (0.5, 0.35, 0.0, -0.2875, -0.369375) + [-0.1 .. 0.01] over [0 .. 1/2]> sage: bp2 <IBP: (-0.369375, -0.45125, -0.3275, 0.14500000000000002, 0.99) + [-0.1 .. 0.01] over [1/2 .. 1]> """
def try_rand_split(self, context ctx, logging_note): """ Compute a random split point r (using the random number generator embedded in ctx). We require 1/4 <= r < 3/4 (to ensure that recursive algorithms make progress).
Then, try doing a de Casteljau split of this polynomial at r, resulting in polynomials p1 and p2. If we see that the sign of this polynomial is determined at r, then return (p1, p2, r); otherwise, return None.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([50, 20, -90, -70, 200], error=5) sage: bp1, bp2, _ = bp.try_rand_split(mk_context(), None) sage: bp1 <IBP: (50, 29, -27, -56, -11) + [0 .. 6) over [0 .. 43/64]> sage: bp2 <IBP: (-11, 10, 49, 111, 200) + [0 .. 6) over [43/64 .. 1]> sage: bp1, bp2, _ = bp.try_rand_split(mk_context(seed=42), None) sage: bp1 <IBP: (50, 32, -11, -41, -29) + [0 .. 6) over [0 .. 583/1024]> sage: bp2 <IBP: (-29, -20, 13, 83, 200) + [0 .. 6) over [583/1024 .. 1]> sage: bp = mk_ibpf([0.5, 0.2, -0.9, -0.7, 0.99], neg_err=-0.1, pos_err=0.01) sage: bp1, bp2, _ = bp.try_rand_split(mk_context(), None) sage: bp1 # rel tol <IBP: (0.5, 0.2984375, -0.2642578125, -0.5511661529541015, -0.3145806974172592) + [-0.1 .. 0.01] over [0 .. 43/64]> sage: bp2 # rel tol <IBP: (-0.3145806974172592, -0.19903896331787108, 0.04135986328125002, 0.43546875, 0.99) + [-0.1 .. 0.01] over [43/64 .. 1]> """
# We want a split point which is a dyadic rational (denominator # is a power of 2), because that speeds up de_casteljau. We want # a small denominator, because that helps us find simpler isolating # intervals. However, if we require that our denominator be too # small, then that might keep the algorithm from terminating; # if the current polynomial has roots at (4/16, 5/16, ..., 12/16), # and we pick our split points from that same set, then we will # never find a split point with a determined sign. We avoid this # problem by making sure we have more possible split points # to choose from than our polynomial has roots.
# A different algorithm here might be more efficient.
div = div * 2
cdef int degree(self): raise NotImplementedError()
def region(self):
def region_width(self):
cdef class interval_bernstein_polynomial_integer(interval_bernstein_polynomial): """ This is the subclass of interval_bernstein_polynomial where polynomial coefficients are represented using integers.
In this integer representation, each coefficient is represented by a GMP arbitrary-precision integer A, and a (shared) interval width E (which is a machine integer). These represent the coefficients A*2^n <= c < (A+E)*2^n.
(Note that mk_ibpi is a simple helper function for creating elements of interval_bernstein_polynomial_integer in doctests.)
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([1, 2, 3], error=5); print(bp) degree 2 IBP with 2-bit coefficients sage: bp <IBP: (1, 2, 3) + [0 .. 5)> sage: bp.variations() (0, 0) sage: bp = mk_ibpi([-3, -1, 1, -1, -3, -1], lower=1, upper=5/4, usign=1, error=2, scale_log2=-3, level=2, slope_err=RIF(pi)); print(bp) degree 5 IBP with 2-bit coefficients sage: bp <IBP: ((-3, -1, 1, -1, -3, -1) + [0 .. 2)) * 2^-3 over [1 .. 5/4]; usign 1; level 2; slope_err 3.141592653589794?> sage: bp.variations() (3, 3) """
def __init__(self, Vector_integer_dense coeffs, Rational lower, Rational upper, int lsign, int usign, int error, int scale_log2, int level, RealIntervalFieldElement slope_err): """ Initialize an interval_bernstein_polynomial_integer.
INPUT:
- ``coeffs`` -- a coefficient vector for a polynomial in Bernstein form - ``lower`` -- the lower bound of the region over which the Bernstein basis is defined - ``upper`` -- the upper bound of the region over which the Bernstein basis is defined - ``lsign`` -- the sign of the polynomial at lower, if known - ``usign`` -- the sign of the polynomial at upper, if known - ``error`` -- the maximum error in the Bernstein coefficients - ``scale_log2`` -- the log2 of the scaling factor for the Bernstein coefficients - ``level`` -- the number of times we have performed degree reduction to get this polynomial - ``slope_err`` -- the maximum extra error in the derivative of this polynomial from degree reduction
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = interval_bernstein_polynomial_integer(vector(ZZ, [50, -30, -10]), -3/7, 4/7, 0, -1, 17, 3, 2, RIF(10^-30)) sage: print(bp) degree 2 IBP with 6-bit coefficients sage: bp <IBP: ((50, -30, -10) + [0 .. 17)) * 2^3 over [-3/7 .. 4/7]; usign -1; level 2; slope_err 1.0000000000000000?e-30> """
def __repr__(self): """ Reveal all the internals of this interval Bernstein polynomial.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([-11, 22, -33], upper=1/9, error=20, lsign=1) sage: repr(bp) '<IBP: (-11, 22, -33) + [0 .. 20) over [0 .. 1/9]; lsign 1>' """
def __str__(self): """ Return a short summary of this interval Bernstein polynomial.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([-11, 22, -33], upper=1/9, error=20, lsign=1) sage: print(bp) degree 2 IBP with 6-bit coefficients sage: str(bp) 'degree 2 IBP with 6-bit coefficients' """
def _type_code(self): """ Classifies this as either an integer representation ('i') or a floating-point representation ('f'). """
cdef void _set_bitsize(self): """ A private function that computes the maximum coefficient size of this Bernstein polynomial (in bits).
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: print(mk_ibpi([2^12345])) degree 0 IBP with 12346-bit coefficients sage: print(mk_ibpi([2^12345 - 1])) degree 0 IBP with 12345-bit coefficients """
cdef void _count_variations(self): """ A private function that counts the number of sign variations in this Bernstein polynomial. Since the coefficients are represented with intervals, not exactly, we cannot necessarily compute the exact number of sign variations; instead, we compute lower and upper bounds on this number.
TESTS::
sage: from sage.rings.polynomial.real_roots import * sage: mk_ibpi([-1, -2, -3], error=1).variations() (0, 0) sage: mk_ibpi([-1, -2, -3], error=2).variations() (0, 2) sage: mk_ibpi([-1, -2, -3], error=2, lsign=-1).variations() (0, 2) sage: mk_ibpi([-1, -2, -3], error=2, lsign=0).variations() (0, 2) sage: mk_ibpi([-1, -2, -3], error=2, lsign=1).variations() (1, 1) sage: mk_ibpi([-1, -2, -3], error=3).variations() (0, 2) sage: mk_ibpi([-1, -2, -3], error=3, lsign=-1).variations() (0, 2) sage: mk_ibpi([-1, -2, -3], error=3, lsign=1).variations() (1, 1) sage: mk_ibpi([-1, -2, -3], error=4).variations() (0, 2) sage: mk_ibpi([-1, -2, -3], error=4, lsign=-1, usign=1).variations() (1, 1) sage: mk_ibpi([-1, -2, -3], error=4, lsign=1, usign=-1).variations() (1, 1) """
cdef int count_maybe_pos, count_maybe_neg cdef int sign
cdef int new_count_maybe_pos, new_count_maybe_neg
cdef int lower_sgn, upper_sgn
cdef int i
else:
else:
cdef int degree(self): """ Returns the (formal) degree of this polynomial. """
def de_casteljau(self, context ctx, mid, msign=0): """ Uses de Casteljau's algorithm to compute the representation of this polynomial in a Bernstein basis over new regions.
INPUT:
- ``mid`` -- where to split the Bernstein basis region; 0 < mid < 1 - ``msign`` -- default 0 (unknown); the sign of this polynomial at mid
OUTPUT:
- ``bp1``, ``bp2`` -- the new interval Bernstein polynomials - ``ok`` -- a boolean; True if the sign of the original polynomial at mid is known
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([50, 20, -90, -70, 200], error=5) sage: ctx = mk_context() sage: bp1, bp2, ok = bp.de_casteljau(ctx, 1/2) sage: bp1 <IBP: (50, 35, 0, -29, -31) + [0 .. 6) over [0 .. 1/2]> sage: bp2 <IBP: (-31, -33, -8, 65, 200) + [0 .. 6) over [1/2 .. 1]> sage: bp1, bp2, ok = bp.de_casteljau(ctx, 2/3) sage: bp1 <IBP: (50, 30, -26, -55, -13) + [0 .. 6) over [0 .. 2/3]> sage: bp2 <IBP: (-13, 8, 47, 110, 200) + [0 .. 6) over [2/3 .. 1]> sage: bp1, bp2, ok = bp.de_casteljau(ctx, 7/39) sage: bp1 <IBP: (50, 44, 36, 27, 17) + [0 .. 6) over [0 .. 7/39]> sage: bp2 <IBP: (17, -26, -75, -22, 200) + [0 .. 6) over [7/39 .. 1]> """
cdef interval_bernstein_polynomial_integer bp1, bp2
(a, b, c, d) = self.lft bp1.lft = (a * mid, b, c * mid, d) bp2.lft = (a * (1-mid), b + a*mid, c * (1-mid), d + c*mid)
def as_float(self): """ Compute an interval_bernstein_polynomial_float which contains (or bounds) all the polynomials this interval polynomial contains (or bounds).
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([50, 20, -90, -70, 200], error=5) sage: print(bp.as_float()) degree 4 IBP with floating-point coefficients sage: bp.as_float() <IBP: ((0.1953125, 0.078125, -0.3515625, -0.2734375, 0.78125) + [-1.12757025938e-16 .. 0.01953125]) * 2^8> """
def get_msb_bit(self): """ Returns an approximation of the log2 of the maximum of the absolute values of the coefficients, as an integer. """
def down_degree(self, context ctx, max_err, exp_err_shift): """ Compute an interval_bernstein_polynomial_integer which bounds all the polynomials this interval polynomial bounds, but is of lesser degree.
During the computation, we find an "expected error" expected_err, which is the error inherent in our approach (this depends on the degrees involved, and is proportional to the error of the current polynomial).
We require that the error of the new interval polynomial be bounded both by max_err, and by expected_err << exp_err_shift. If we find such a polynomial p, then we return a pair of p and some debugging/logging information. Otherwise, we return the pair (None, None).
If the resulting polynomial would have error more than 2^17, then it is downscaled before returning.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([0, 100, 400, 903], error=2) sage: ctx = mk_context() sage: bp <IBP: (0, 100, 400, 903) + [0 .. 2)> sage: dbp, _ = bp.down_degree(ctx, 10, 32) sage: dbp <IBP: (-1, 148, 901) + [0 .. 4); level 1; slope_err 0.?e2> """
# return (None, None)
# v0 = dprod_matrow_vec(downmat, self.coeffs, 0) # vn = dprod_matrow_vec(downmat, self.coeffs, next) return (None, ('>', self.scale_log2 + bitsize(vn_err)))
else: # slope_err could be smaller by actually computing # the derivative of the error polynomial else: # warning: lsign and usign might have changed, invalidate them to 0
else:
def down_degree_iter(self, context ctx, max_scale): """ Compute a degree-reduced version of this interval polynomial, by iterating down_degree.
We stop when degree reduction would give a polynomial which is too inaccurate, meaning that either we think the current polynomial may have more roots in its region than the degree of the reduced polynomial, or that the least significant accurate bit in the result (on the absolute scale) would be larger than 1 << max_scale.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([0, 100, 400, 903, 1600, 2500], error=2) sage: ctx = mk_context() sage: bp <IBP: (0, 100, 400, 903, 1600, 2500) + [0 .. 2)> sage: rbp = bp.down_degree_iter(ctx, 6) sage: rbp <IBP: (-4, 249, 2497) + [0 .. 9); level 2; slope_err 0.?e3> """
# global dd_count_no # dd_count_no += 1 else: # global dd_count_yes # dd_count_yes += 1
def downscale(self, bits): """ Compute an interval_bernstein_polynomial_integer which contains (or bounds) all the polynomials this interval polynomial contains (or bounds), but uses "bits" fewer bits.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([0, 100, 400, 903], error=2) sage: bp.downscale(5) <IBP: ((0, 3, 12, 28) + [0 .. 1)) * 2^5> """
def slope_range(self): """ Compute a bound on the derivative of this polynomial, over its region.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([0, 100, 400, 903], error=2) sage: bp.slope_range().str(style='brackets') '[294.00000000000000 .. 1515.0000000000000]' """ else:
""" A simple wrapper for creating interval_bernstein_polynomial_integer objects with coercions, defaults, etc.
For use in doctests.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: print(mk_ibpi([50, 20, -90, -70, 200], error=5)) degree 4 IBP with 8-bit coefficients """
""" Given a polynomial in Bernstein form with integer coefficients over the region [0 .. 1], and a split point x, use de Casteljau's algorithm to give polynomials in Bernstein form over [0 .. x] and [x .. 1].
This function will work for an arbitrary rational split point x, as long as 0 < x < 1; but it has specialized code paths that make some values of x faster than others. If x == a/(a + b), there are special efficient cases for a==1, b==1, a+b fits in a machine word, a+b is a power of 2, a fits in a machine word, b fits in a machine word. The most efficient case is x==1/2.
Given split points x == a/(a + b) and y == c/(c + d), where min(a, b) and min(c, d) fit in the same number of machine words and a+b and c+d are both powers of two, then x and y should be equally fast split points.
If use_ints is nonzero, then instead of checking whether numerators and denominators fit in machine words, we check whether they fit in ints (32 bits, even on 64-bit machines). This slows things down, but allows for identical results across machines.
INPUT:
- ``c`` -- vector of coefficients of polynomial in Bernstein form - ``c_bitsize`` -- approximate size of coefficients in c (in bits) - ``x`` -- rational splitting point; 0 < x < 1
OUTPUT:
- ``c1`` -- coefficients of polynomial over range [0 .. x] - ``c2`` -- coefficients of polynomial over range [x .. 1] - ``err_inc`` -- amount by which error intervals widened
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: c = vector(ZZ, [1048576, 0, 0, 0, 0, 0]) sage: de_casteljau_intvec(c, 20, 1/2, 1) ((1048576, 524288, 262144, 131072, 65536, 32768), (32768, 0, 0, 0, 0, 0), 1) sage: de_casteljau_intvec(c, 20, 1/3, 1) ((1048576, 699050, 466033, 310689, 207126, 138084), (138084, 0, 0, 0, 0, 0), 1) sage: de_casteljau_intvec(c, 20, 7/22, 1) ((1048576, 714938, 487457, 332357, 226607, 154505), (154505, 0, 0, 0, 0, 0), 1) """
cdef Vector_integer_dense c1, c2, den_powers
cdef int i, j
cdef mpz_t num, den, diff, tmp, tmp2
cdef unsigned long num_ui, den_ui, diff_ui cdef int num_fits_ui, den_fits_ui, diff_fits_ui cdef int den_is_pow2, den_log2 cdef int num_less_diff
cdef int max_den_power
# We want to compute (diff*a + num*b)/den. This costs # 2*mul, 1*div, 1*add/sub. # We see below a way to trade a multiplication for a subtraction, # for a cost of 1*mul, 1*div, 2*add/sub. # Another possibility is to postpone divisions, for a cost of # 2*mul, 1/6*div, 1*add/sub. Clearly, this is worthwhile if # 5/6*div + 1*add/sub > 1*mul. This is very likely true whenever # the denominator is not a power of 2 (so that you have to do # real divisions, instead of shifts), and might be true in other cases.
# If one of the multiplications is by 1, then the costs become: # straightforward: 1*mul, 1*div, 1*add/sub # trade mul->sub: 1*div, 2*add/sub # postpone division: 1*mul, 1/6*div, 1*add/sub # for the same result.
# If both multiplications are by 1, then the costs are: # straightforward: 1*div, 1*add/sub # trade mul->sub: 1*div, 2*add/sub # postpone division: 1/6*div, 1*add/sub # with postpone division being the obvious winner.
# For now, we use "postpone division" whenever the denominator fits # in an unsigned long. This is not going to be drastically bad, # even when dividing by a power of two. For full generality, # it would also be important to use "postpone division" for # non-power-of-two denominators that do not fit in an unsigned long; # however, the current calling code essentially never does that, # so we'll stick with simpler code here.
else:
# These settings are slower than the above on laguerre(1000), but that's # the only experiment I've done so far... more testing is needed. # cdef int max_den_bits = 3 * c_bitsize / 2 # if max_den_bits < 100: max_den_bits = 300
else:
# x == 1/2 else: else: c1._entries[i-cur_den_steps+j+1], j*den_log2) c2._entries[n-i+cur_den_steps-j-2], j*den_log2) cur_den_steps*den_log2) else: c1._entries[i-cur_den_steps+j+1], den_powers._entries[j]) c2._entries[n-i+cur_den_steps-j-2], den_powers._entries[j]) mpz_fdiv_q(c2._entries[j], c2._entries[j], den_powers._entries[cur_den_steps]) else: # We want to compute (diff*a + num*b)/den, where # a == c2._entries[j] and b == c2._entries[j+1]. # The result goes in c2._entries[j]. Since den == diff + num, # this is equal to a + num*(b-a)/den or diff*(a-b)/den + b. # If num<diff, we compute the former, otherwise the latter.
elif den_fits_ui: mpz_fdiv_q_ui(tmp2, tmp, den_ui) else: mpz_fdiv_q(tmp2, tmp, den) else: else: mpz_mul(tmp2, tmp, num) elif den_fits_ui: mpz_fdiv_q_ui(tmp, tmp2, den_ui) else: mpz_fdiv_q(tmp, tmp2, den) else: else: else: mpz_mul(tmp, c2._entries[j], diff) elif den_fits_ui: mpz_fdiv_q_ui(c2._entries[j], tmp, den_ui) else: mpz_fdiv_q(c2._entries[j], tmp, den)
# An ULP is a "unit in the last place"; it is the (varying) unit for # how much adjacent floating-point numbers differ from each other. # A half-ULP is half this amount; it is the maximum rounding error # in the basic operations (+,-,*,/) in a correctly-operating IEEE # floating-point unit. # (Note that by default, the x86 does not use IEEE double precision; # instead, it uses extra precision, which can (counterintuitively) # actually give worse double-precision results in some rare cases. # To avoid this, we change the x86 floating-point unit into true # double-precision mode in places where it matters; that is, in # functions that add, subtract, multiply, or divide floating-point numbers. # Functions whose floating-point operations are limited to negation # and comparison do not require special treatment, since those operations # give the same results in double-precision or extended precision.)
# This is the value of a half-ULP for numbers in the range 0.5 <= x < 1. # This is actually slightly more than a half-ULP because of possible # double-rounding on x86 PCs.
""" Given a vector of integers A = [a1, ..., an], and an integer error bound E, returns a vector of floating-point numbers B = [b1, ..., bn], lower and upper error bounds F1 and F2, and a scaling factor d, such that
.. MATH::
(bk + F1) * 2^d \le ak
and
.. MATH::
ak + E \le (bk + F2) * 2^d
If bj is the element of B with largest absolute value, then 0.5 <= abs(bj) < 1.0 .
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: intvec_to_doublevec(vector(ZZ, [1, 2, 3, 4, 5]), 3) ((0.125, 0.25, 0.375, 0.5, 0.625), -1.1275702593849246e-16, 0.37500000000000017, 3) """ cdef unsigned int cwf # fpu_fix_start(&cwf)
cdef long cur_exp
cdef int i
cdef double d cdef int new_exp
# 0.5 <= d < 1; b._entries[i] ~= d*2^cur_exp # db[i] = ldexp(d, new_exp)
# The true value can be off by an ulp because mpz_get_d_2exp truncates. # (If we created a version of mpz_get_d_2exp that rounded instead, # we could do a little better.) # The third half_ulp in the positive direction is to compensate for # possible bad rounding when adding ldexp(err, delta).
# fpu_fix_end(&cwf)
cdef class interval_bernstein_polynomial_float(interval_bernstein_polynomial): """ This is the subclass of interval_bernstein_polynomial where polynomial coefficients are represented using floating-point numbers.
In the floating-point representation, each coefficient is represented as an IEEE double-precision float A, and the (shared) lower and upper interval widths E1 and E2. These represent the coefficients (A+E1)*2^n <= c <= (A+E2)*2^n.
Note that we always have E1 <= 0 <= E2. Also, each floating-point coefficient has absolute value less than one.
(Note that mk_ibpf is a simple helper function for creating elements of interval_bernstein_polynomial_float in doctests.)
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpf([0.1, 0.2, 0.3], pos_err=0.5); print(bp) degree 2 IBP with floating-point coefficients sage: bp <IBP: (0.1, 0.2, 0.3) + [0.0 .. 0.5]> sage: bp.variations() (0, 0) sage: bp = mk_ibpf([-0.3, -0.1, 0.1, -0.1, -0.3, -0.1], lower=1, upper=5/4, usign=1, pos_err=0.2, scale_log2=-3, level=2, slope_err=RIF(pi)); print(bp) degree 5 IBP with floating-point coefficients sage: bp <IBP: ((-0.3, -0.1, 0.1, -0.1, -0.3, -0.1) + [0.0 .. 0.2]) * 2^-3 over [1 .. 5/4]; usign 1; level 2; slope_err 3.141592653589794?> sage: bp.variations() (3, 3) """
def __init__(self, Vector_real_double_dense coeffs, Rational lower, Rational upper, int lsign, int usign, double neg_err, double pos_err, int scale_log2, int level, RealIntervalFieldElement slope_err): """ Initialize an interval_bernstein_polynomial_float.
INPUT:
- ``coeffs`` -- a coefficient vector for a polynomial in Bernstein form (all coefficients must have absolute value less than one) - ``lower`` -- the lower bound of the region over which the Bernstein basis is defined - ``upper`` -- the upper bound of the region over which the Bernstein basis is defined - ``lsign`` -- the sign of the polynomial at lower, if known - ``usign`` -- the sign of the polynomial at upper, if known - ``neg_err`` -- the minimum error in the Bernstein coefficients - ``pos_err`` -- the maximum error in the Bernstein coefficients (so a Bernstein coefficient x really represents the range [neg_err+x .. pos_err+x] - ``scale_log2`` -- the log2 of the scaling factor for the Bernstein coefficients - ``level`` -- the number of times we have performed degree reduction to get this polynomial - ``slope_err`` -- the maximum extra error in the derivative of this polynomial from degree reduction
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = interval_bernstein_polynomial_float(vector(RDF, [0.50, -0.30, -0.10]), -3/7, 4/7, 0, -1, -0.02, 0.17, 3, 2, RIF(10^-30)) sage: print(bp) degree 2 IBP with floating-point coefficients sage: bp <IBP: ((0.5, -0.3, -0.1) + [-0.02 .. 0.17]) * 2^3 over [-3/7 .. 4/7]; usign -1; level 2; slope_err 1.0000000000000000?e-30> """ cdef int exp
def __repr__(self): """ Reveal all the internals of this interval Bernstein polynomial.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpf([-0.11, 0.22, -0.33], upper=1/9, neg_err=-0.3, pos_err=0.1, lsign=1) sage: repr(bp) '<IBP: (-0.11, 0.22, -0.33) + [-0.3 .. 0.1] over [0 .. 1/9]>' """ s += "; lsign %d" % self.lsign
def __str__(self): """ Return a short summary of this interval Bernstein polynomial.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpf([-0.11, 0.22, -0.33], upper=1/9, neg_err=-0.1, pos_err=0.2, lsign=1) sage: print(bp) degree 2 IBP with floating-point coefficients sage: str(bp) 'degree 2 IBP with floating-point coefficients' """
def _type_code(self): """ Classifies this as either an integer representation ('i') or a floating-point representation ('f'). """
cdef void _count_variations(self): """ A private function that counts the number of sign variations in this Bernstein polynomial. Since the coefficients are represented with intervals, not exactly, we cannot necessarily compute the exact number of sign variations; instead, we compute lower and upper bounds on this number.
"""
cdef int count_maybe_pos cdef int count_maybe_neg cdef int sign
cdef int new_count_maybe_pos, new_count_maybe_neg
cdef int i
cdef double val
else:
else:
cdef int degree(self): """ Returns the (formal) degree of this polynomial. """
def de_casteljau(self, context ctx, mid, msign=0): """ Uses de Casteljau's algorithm to compute the representation of this polynomial in a Bernstein basis over new regions.
INPUT:
- ``mid`` -- where to split the Bernstein basis region; 0 < mid < 1 - ``msign`` -- default 0 (unknown); the sign of this polynomial at mid
OUTPUT:
- ``bp1``, ``bp2`` -- the new interval Bernstein polynomials - ``ok`` -- a boolean; True if the sign of the original polynomial at mid is known
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: ctx = mk_context() sage: bp = mk_ibpf([0.5, 0.2, -0.9, -0.7, 0.99], neg_err=-0.1, pos_err=0.01) sage: bp1, bp2, ok = bp.de_casteljau(ctx, 1/2) sage: bp1 <IBP: (0.5, 0.35, 0.0, -0.2875, -0.369375) + [-0.1 .. 0.01] over [0 .. 1/2]> sage: bp2 <IBP: (-0.369375, -0.45125, -0.3275, 0.14500000000000002, 0.99) + [-0.1 .. 0.01] over [1/2 .. 1]> sage: bp1, bp2, ok = bp.de_casteljau(ctx, 2/3) sage: bp1 # rel tol 2e-16 <IBP: (0.5, 0.30000000000000004, -0.2555555555555555, -0.5444444444444444, -0.32172839506172846) + [-0.1 .. 0.01] over [0 .. 2/3]> sage: bp2 # rel tol 3e-15 <IBP: (-0.32172839506172846, -0.21037037037037046, 0.028888888888888797, 0.4266666666666666, 0.99) + [-0.1 .. 0.01] over [2/3 .. 1]> sage: bp1, bp2, ok = bp.de_casteljau(ctx, 7/39) sage: bp1 # rel tol <IBP: (0.5, 0.4461538461538461, 0.36653517422748183, 0.27328680523946786, 0.1765692706232836) + [-0.1 .. 0.01] over [0 .. 7/39]> sage: bp2 # rel tol <IBP: (0.1765692706232836, -0.26556803047927313, -0.7802038132807364, -0.3966666666666666, 0.99) + [-0.1 .. 0.01] over [7/39 .. 1]> """
elif sign != 0: assert(msign == sign)
# As long as new_neg and new_pos have # magnitudes less than 0.5, these computations # are exact. This will be the case for any sensible # usage of this class.
# Give up on this computation...it's horribly inaccurate anyway.
cdef interval_bernstein_polynomial_float bp1, bp2
(a, b, c, d) = self.lft bp1.lft = (a * mid, b, c * mid, d) bp2.lft = (a * (1-mid), b + a*mid, c * (1-mid), d + c*mid)
def as_float(self): return self
def get_msb_bit(self): """ Returns an approximation of the log2 of the maximum of the absolute values of the coefficients, as an integer. """
def slope_range(self): """ Compute a bound on the derivative of this polynomial, over its region.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpf([0.5, 0.2, -0.9, -0.7, 0.99], neg_err=-0.1, pos_err=0.01) sage: bp.slope_range().str(style='brackets') '[-4.8400000000000017 .. 7.2000000000000011]' """ cdef unsigned int cwf # fpu_fix_start(&cwf)
# 2 half_ulp's because subtracting two numbers with absolute values # in (-1 .. 1) can give a number in (-2 .. 2), and the subtraction # can have an error of up to half an ulp in that range, which # is 2 half ulps for (-1 .. 1). else:
# fpu_fix_end(&cwf)
""" A simple wrapper for creating interval_bernstein_polynomial_float objects with coercions, defaults, etc.
For use in doctests.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: print(mk_ibpf([0.5, 0.2, -0.9, -0.7, 0.99], pos_err=0.1, neg_err=-0.01)) degree 4 IBP with floating-point coefficients """
""" Given a polynomial in Bernstein form with floating-point coefficients over the region [0 .. 1], and a split point x, use de Casteljau's algorithm to give polynomials in Bernstein form over [0 .. x] and [x .. 1].
This function will work for an arbitrary rational split point x, as long as 0 < x < 1; but it has a specialized code path for x==1/2.
INPUT:
- ``c`` -- vector of coefficients of polynomial in Bernstein form - ``x`` -- rational splitting point; 0 < x < 1
OUTPUT:
- ``c1`` -- coefficients of polynomial over range [0 .. x] - ``c2`` -- coefficients of polynomial over range [x .. 1] - ``err_inc`` -- number of half-ulps by which error intervals widened
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: c = vector(RDF, [0.7, 0, 0, 0, 0, 0]) sage: de_casteljau_doublevec(c, 1/2) ((0.7, 0.35, 0.175, 0.0875, 0.04375, 0.021875), (0.021875, 0.0, 0.0, 0.0, 0.0, 0.0), 5) sage: de_casteljau_doublevec(c, 1/3) # rel tol ((0.7, 0.4666666666666667, 0.31111111111111117, 0.20740740740740746, 0.13827160493827165, 0.09218106995884777), (0.09218106995884777, 0.0, 0.0, 0.0, 0.0, 0.0), 15) sage: de_casteljau_doublevec(c, 7/22) # rel tol ((0.7, 0.4772727272727273, 0.3254132231404959, 0.22187265214124724, 0.15127680827812312, 0.10314327837144759), (0.10314327837144759, 0.0, 0.0, 0.0, 0.0, 0.0), 15) """
cdef Vector_real_double_dense c1, c2
cdef unsigned int cwf # fpu_fix_start(&cwf)
cdef int i, j
cdef double rx, rnx
cdef int extra_err
# The following code lets us avoid O(n^2) floating-point multiplications # in favor of O(n) calls to ldexp(). In one experiment, though, it seems # to actually slow down the code by about 10%. (On an Intel Core 2 Duo # in 32-bit mode, on the test chebyt2(200).) # cur_den_steps = cur_den_steps + 1 # if cur_den_steps == 1020 or i == n-1: # for j from 1 <= j < cur_den_steps: # c1d[i-cur_den_steps+j+1] = ldexp(c1d[i-cur_den_steps+j+1], -j) # c2d[n-i+cur_den_steps-j-2] = ldexp(c2d[n-i+cur_den_steps-j-2], -j) # for j from 0 <= j < n-i-1: # c2d[j] = ldexp(c2d[j], -cur_den_steps)
else:
# fpu_fix_end(&cwf)
""" Given a floating-point vector, return the maximum of the absolute values of its elements.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: max_abs_doublevec(vector(RDF, [0.1, -0.767, 0.3, 0.693])) 0.767 """
cdef double a
""" Given rationals a and b, selects a de Casteljau split point r between a and b. An attempt is made to select an efficient split point (according to the criteria mentioned in the documentation for de_casteljau_intvec), with a bias towards split points near a.
In full detail:
Takes as input two rationals, a and b, such that 0<=a<=1, 0<=b<=1, and a!=b. Returns rational r, such that a<=r<=b or b<=r<=a. The denominator of r is a power of 2. Let m be min(r, 1-r), nm be numerator(m), and dml be log2(denominator(m)). The return value r is taken from the first of the following classes to have any members between a and b (except that if a <= 1/8, or 7/8 <= a, then class 2 is preferred to class 1).
1. dml < wordsize 2. bitsize(nm) <= wordsize 3. bitsize(nm) <= 2*wordsize 4. bitsize(nm) <= 3*wordsize
...
k. bitsize(nm) <= (k-1)*wordsize
From the first class to have members between a and b, r is chosen as the element of the class which is closest to a.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: wordsize_rational(1/5, 1/7, 32) 429496729/2147483648 sage: wordsize_rational(1/7, 1/5, 32) 306783379/2147483648 sage: wordsize_rational(1/5, 1/7, 64) 1844674407370955161/9223372036854775808 sage: wordsize_rational(1/7, 1/5, 64) 658812288346769701/4611686018427387904 sage: wordsize_rational(1/17, 1/19, 32) 252645135/4294967296 sage: wordsize_rational(1/17, 1/19, 64) 1085102592571150095/18446744073709551616 sage: wordsize_rational(1/1234567890, 1/1234567891, 32) 933866427/1152921504606846976 sage: wordsize_rational(1/1234567890, 1/1234567891, 64) 4010925763784056541/4951760157141521099596496896 """
# assert 0 <= a <= 1
# fld = RealField(cur_size, rnd='RNDU') cdef RealNumber rf # rf2 = RealField(cur_size + exp - 1, rnd='RNDU')(a) if rf <= -(fld(-b)): break cur_size = cur_size + wordsize fld = RealField(cur_size, rnd='RNDU')
# assert 0 <= r <= 1
""" INPUT:
- ``(al, ah)`` -- pair of rationals - ``(bl, bh)`` -- pair of rationals
OUTPUT:
- ``(cl, ch)`` -- pair of rationals
Computes the linear transformation that maps (al, ah) to (0, 1); then applies this transformation to (bl, bh) and returns the result.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: relative_bounds((1/7, 1/4), (1/6, 1/5)) (2/9, 8/15) """
cdef int bitsize(Integer a): """ Compute the number of bits required to write a given integer (the sign is ignored).
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bitsize_doctest(0) 1 sage: bitsize_doctest(1) 1 sage: bitsize_doctest(2) 2 sage: bitsize_doctest(-2) 2 sage: bitsize_doctest(65535) 16 sage: bitsize_doctest(65536) 17 """
""" Given n (a polynomial degree), returns either a smaller integer or None. This defines the sequence of degrees followed by our degree reduction implementation.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: degree_reduction_next_size(1000) 30 sage: degree_reduction_next_size(20) 15 sage: degree_reduction_next_size(3) 2 sage: degree_reduction_next_size(2) is None True """
# Virtually any descending sequence would be "correct" here; no # testing has been done to see if another sequence would be better.
# Constraints: must return None for n <= 2; degree reduction to # degrees > 30 may be very, very slow. (Part of reducing from # degree n to degree k is computing the exact inverse of a k by k # rational matrix. Fortunately, this matrix depends only on # n and k, so the inverted matrix can be cached; but still, # computing the exact inverse of a k by k matrix seems infeasible # for k much larger than 30.)
""" Compute and cache the matrices used for degree reduction, starting from degree n.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: precompute_degree_reduction_cache(5) sage: dr_cache[5] ( [121/126 8/63 -1/9 -2/63 11/126 -2/63] [ -3/7 37/42 16/21 1/21 -3/7 1/6] [ 1/6 -3/7 1/21 16/21 37/42 -3/7] 3, [ -2/63 11/126 -2/63 -1/9 8/63 121/126], 2, <BLANKLINE> ([121 16 -14 -4 11 -4] [-54 111 96 6 -54 21] [ 21 -54 6 96 111 -54] [ -4 11 -4 -14 16 121], 126) ) """
# We can compute a degree n -> degree k reduction by looking at # as few as k+1 of the coefficients of the degree n polynomial. # Using more samples reduces the error involved in degree # reduction, but is slower. More testing might reveal a better # way to select the number of samples here.
# For each coefficient in the reduced polynomial, we see how # it varies as the (sampled) coefficients in the original # polynomial change by +/- 1. Then we take the reduced # coefficient which changes the most, and call that the expected # error. We can multiply this by the error in the original # polynomial, and be fairly certain (absolutely certain?) that # the error in the reduced polynomial will be no better # than this product.
# bdd = bd.denominator() # bdi = MatrixSpace(ZZ, next+1, samps, sparse=False)(bd * bdd)
""" Given polynomial degrees d1 and d2 (where d1 < d2), and a number of samples s, computes a matrix bd.
If you have a Bernstein polynomial of formal degree d2, and select s of its coefficients (according to subsample_vec), and multiply the resulting vector by bd, then you get the coefficients of a Bernstein polynomial of formal degree d1, where this second polynomial is a good approximation to the first polynomial over the region of the Bernstein basis.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bernstein_down(3, 8, 5) [ 612/245 -348/245 -37/49 338/245 -172/245] [-724/441 132/49 395/441 -290/147 452/441] [ 452/441 -290/147 395/441 132/49 -724/441] [-172/245 338/245 -37/49 -348/245 612/245] """
# The use of the pseudoinverse means that the coefficients of the # reduced polynomial are selected according to a least-squares fit. # We would prefer to minimize the maximum error, rather than the # sum of the squares of the errors; but I don't know how to do that.
# Also, this pseudoinverse is very slow to compute if d1 is large. # Since the coefficients of bernstein_up(...) can be represented # with a fairly simple formula (see the implementation of # bernstein_up()), it seems possible that there is some sort of # formula for the coefficients of the pseudoinverse that would be # faster than the computation here.
""" Given polynomial degrees d1 and d2, where d1 < d2, compute a matrix bu.
If you have a Bernstein polynomial of formal degree d1, and multiply its coefficient vector by bu, then the result is the coefficient vector of the same polynomial represented as a Bernstein polynomial of formal degree d2.
If s is not None, then it represents a number of samples; then the product only gives s of the coefficients of the new Bernstein polynomial, selected according to subsample_vec.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bernstein_down(3, 7, 4) [ 12/5 -4 3 -2/5] [-13/15 16/3 -4 8/15] [ 8/15 -4 16/3 -13/15] [ -2/5 3 -4 12/5] """
cdef int subsample_vec(int a, int slen, int llen): """ Given a vector of length llen, and slen < llen, we want to select slen of the elements of the vector, evenly spaced.
Given an index 0 <= a < slen, this function computes the index in the original vector of length llen corresponding to a.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: [subsample_vec_doctest(a, 10, 100) for a in range(10)] [5, 15, 25, 35, 45, 54, 64, 74, 84, 94] sage: [subsample_vec_doctest(a, 3, 4) for a in range(3)] [1, 2, 3] """
# round((a + 0.5) * (llen - 1) / slen) # round((2*a + 1) * (llen - 1) / (2 * slen) # floor(((2*a + 1) * (llen - 1) + slen) / (2 * slen))
""" Given a polynomial with real coefficients, computes an upper bound on its largest real root, using the first-\lambda algorithm from "Implementations of a New Theorem for Computing Bounds for Positive Roots of Polynomials", by Akritas, Strzebo\'nski, and Vigklas.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(ZZ) sage: maximum_root_first_lambda((x-1)*(x-2)*(x-3)) 6.00000000000001 sage: maximum_root_first_lambda((x+1)*(x+2)*(x+3)) 0.000000000000000 sage: maximum_root_first_lambda(x^2 - 1) 1.00000000000000 """
r""" Given a polynomial represented by a list of its coefficients (as RealIntervalFieldElements), compute an upper bound on its largest real root.
Uses the first-\lambda algorithm from "Implementations of a New Theorem for Computing Bounds for Positive Roots of Polynomials", by Akritas, Strzebo\'nski, and Vigklas.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: cl_maximum_root_first_lambda([RIF(-1), RIF(0), RIF(1)]) 1.00000000000000
TESTS::
sage: bnd = cl_maximum_root_first_lambda(list(map(RIF, [0, 0, 0, 14, 1]))) sage: bnd, bnd.parent() (0.000000000000000, Real Field with 53 bits of precision and rounding RNDU) """ else: else:
r""" Given a polynomial with real coefficients, computes an upper bound on its largest real root, using the local-max algorithm from "Implementations of a New Theorem for Computing Bounds for Positive Roots of Polynomials", by Akritas, Strzebo\'nski, and Vigklas.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(ZZ) sage: maximum_root_local_max((x-1)*(x-2)*(x-3)) 12.0000000000001 sage: maximum_root_local_max((x+1)*(x+2)*(x+3)) 0.000000000000000 sage: maximum_root_local_max(x^2 - 1) 1.41421356237310 """
r""" Given a polynomial represented by a list of its coefficients (as RealIntervalFieldElements), compute an upper bound on its largest real root.
Uses the local-max algorithm from "Implementations of a New Theorem for Computing Bounds for Positive Roots of Polynomials", by Akritas, Strzebo\'nski, and Vigklas.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: cl_maximum_root_local_max([RIF(-1), RIF(0), RIF(1)]) 1.41421356237310 """
r""" Given a polynomial represented by a list of its coefficients (as RealIntervalFieldElements), compute an upper bound on its largest real root.
Uses two algorithms of Akritas, Strzebo\'nski, and Vigklas, and picks the better result.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: cl_maximum_root([RIF(-1), RIF(0), RIF(1)]) 1.00000000000000 """
r""" Given a polynomial with real coefficients, computes a lower and upper bound on its real roots. Uses algorithms of Akritas, Strzebo\'nski, and Vigklas.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(ZZ) sage: root_bounds((x-1)*(x-2)*(x-3)) (0.545454545454545, 6.00000000000001) sage: root_bounds(x^2) (0.000000000000000, 0.000000000000000) sage: root_bounds(x*(x+1)) (-1.00000000000000, 0.000000000000000) sage: root_bounds((x+2)*(x-3)) (-2.44948974278317, 3.46410161513776) sage: root_bounds(x^995 * (x^2 - 9999) - 1) (-99.9949998749937, 141.414284992713) sage: root_bounds(x^995 * (x^2 - 9999) + 1) (-141.414284992712, 99.9949998749938)
If we can see that the polynomial has no real roots, return None. ::
sage: root_bounds(x^2 + 1) is None True """
# not RIF.zero().endpoints() because of MPFI's convention that the # upper bound is -0.
else: return (lb, ub)
""" Given a polynomial p with real coefficients, computes rationals a and b, such that for every real root r of p, a < r < b. We try to find rationals which bound the roots somewhat tightly, yet are simple (have small numerators and denominators).
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(ZZ) sage: rational_root_bounds((x-1)*(x-2)*(x-3)) (0, 7) sage: rational_root_bounds(x^2) (-1/2, 1/2) sage: rational_root_bounds(x*(x+1)) (-3/2, 1/2) sage: rational_root_bounds((x+2)*(x-3)) (-3, 6) sage: rational_root_bounds(x^995 * (x^2 - 9999) - 1) (-100, 1000/7) sage: rational_root_bounds(x^995 * (x^2 - 9999) + 1) (-142, 213/2)
If we can see that the polynomial has no real roots, return None. sage: rational_root_bounds(x^2 + 7) is None True """
# There are two stages to the root isolation process. First, # we convert the polynomial to the Bernstein basis given by # the root bounds, exactly; then we isolate the roots from # that Bernstein basis, using interval arithmetic.
# We want the root bounds to be as small as possible, because that # speeds up the second phase; but we also want them to be as # simple as possible, because that speeds up the first phase.
# As a heuristic, given a polynomial of degree d with floating-point # root bounds lb and ub, we compute simple rational root bounds # rlb and rub such that lb - (ub - lb)/sqrt(d) <= rlb <= lb, # and ub <= rub <= ub + (ub - lb)/sqrt(d). Very little testing # was done to come up with this heuristic, and it can probably # be improved.
return (QQ(-1), QQ(1))
# The code below would give overly precise bounds in this case, # giving very non-simple isolating intervals in the result. # We don't need tight root bounds to quickly find the roots # of a linear polynomial, so go for simple root bounds.
# XXX A gross hack. Sometimes, root_bounds() gives a lower or # upper bound which is a root. We want bounds which are # not roots, so we just spread out the bounds a tiny bit. # (It might be more efficient to check the bounds from root_bounds() # and, if they are roots, use that information to factor out a linear # polynomial. However, as far as I can tell, root_bounds() only # gives roots as bounds on toy examples; so this is not too inefficient.)
# Given the current implementation of to_bernstein, we want rlb # and (rub/rlb - 1) to be simple rationals -- we don't really care # about the simplicity of rub by itself. (Or else we want rlb to # be zero, and rub to be a simple rational.)
# We want to compute an interval for the upper bound, # iub = [ub .. ub + (ub - lb)/sqrtd], # and then compute (iub/ilb - 1).simplest_rational(). # However, we want to compute this # interval with inward instead of outward rounding. Instead, # we compute the lower and upper bounds of the interval separately.
else:
pass
""" An abstract base class for bernstein_polynomial factories. That is, elements of subclasses represent Bernstein polynomials (exactly), and are responsible for creating interval_bernstein_polynomial_integer approximations at arbitrary precision.
Supports four methods, coeffs_bitsize(), bernstein_polynomial(), lsign(), and usign(). The coeffs_bitsize() method gives an integer approximation to the log2 of the max of the absolute values of the Bernstein coefficients. The bernstein_polynomial(scale_log2) method gives an approximation where the maximum coefficient has approximately coeffs_bitsize() - scale_log2 bits. The lsign() and usign() methods give the (exact) sign of the first and last coefficient, respectively. """
return 0
""" Returns the sign of the first coefficient of this Bernstein polynomial. """
""" Returns the sign of the last coefficient of this Bernstein polynomial. """
""" This class holds an exact Bernstein polynomial (represented as a list of integer coefficients), and returns arbitrarily-precise interval approximations of this polynomial on demand. """
""" Initializes a bernstein_polynomial_factory_intlist, given a list of integer coefficients.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bernstein_polynomial_factory_intlist([1, 2, 3]) degree 2 Bernstein factory with 2-bit integer coefficients """
""" Return a short summary of this Bernstein polynomial factory.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bernstein_polynomial_factory_intlist([1, 2, 3, 4000]) degree 3 Bernstein factory with 12-bit integer coefficients """
""" Computes the approximate log2 of the maximum of the absolute values of the coefficients.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bernstein_polynomial_factory_intlist([1, 2, 3, -60000]).coeffs_bitsize() 16 """ # return intlist_size(self.coeffs)
""" Compute an interval_bernstein_polynomial_integer that approximates this polynomial, using the given scale_log2. (Smaller scale_log2 values give more accurate approximations.)
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bpf = bernstein_polynomial_factory_intlist([10, -20, 30, -40]) sage: print(bpf.bernstein_polynomial(0)) degree 3 IBP with 6-bit coefficients sage: bpf.bernstein_polynomial(20) <IBP: ((0, -1, 0, -1) + [0 .. 1)) * 2^20; lsign 1> sage: bpf.bernstein_polynomial(0) <IBP: (10, -20, 30, -40) + [0 .. 1)> sage: bpf.bernstein_polynomial(-20) <IBP: ((10485760, -20971520, 31457280, -41943040) + [0 .. 1)) * 2^-20> """ else:
# return bp_of_intlist(self.coeffs, scale_log2)
""" This class holds an exact Bernstein polynomial (represented as a list of rational coefficients), and returns arbitrarily-precise interval approximations of this polynomial on demand. """
""" Initializes a bernstein_polynomial_factory_intlist, given a list of rational coefficients.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bernstein_polynomial_factory_ratlist([1, 1/2, 1/3]) degree 2 Bernstein factory with 0-bit rational coefficients """
""" Return a short summary of this Bernstein polynomial factory.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bernstein_polynomial_factory_ratlist([1, 2, 3, 4000/17]) degree 3 Bernstein factory with 7-bit rational coefficients """
""" Computes the approximate log2 of the maximum of the absolute values of the coefficients.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bernstein_polynomial_factory_ratlist([1, 2, 3, -60000]).coeffs_bitsize() 15 sage: bernstein_polynomial_factory_ratlist([65535/65536]).coeffs_bitsize() -1 sage: bernstein_polynomial_factory_ratlist([65536/65535]).coeffs_bitsize() 1 """
# This returns an estimate of the log base 2 of the rational in the # list with largest absolute value. The estimate is an integer, # and within +/- 1 of the correct answer.
""" Compute an interval_bernstein_polynomial_integer that approximates this polynomial, using the given scale_log2. (Smaller scale_log2 values give more accurate approximations.)
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bpf = bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]) sage: print(bpf.bernstein_polynomial(0)) degree 3 IBP with 3-bit coefficients sage: bpf.bernstein_polynomial(20) <IBP: ((0, -1, 0, -1) + [0 .. 1)) * 2^20; lsign 1> sage: bpf.bernstein_polynomial(0) <IBP: (0, -4, 2, -2) + [0 .. 1); lsign 1> sage: bpf.bernstein_polynomial(-20) <IBP: ((349525, -3295525, 2850354, -1482835) + [0 .. 1)) * 2^-20> """
else:
# return bp_of_ratlist(self.coeffs, scale_log2)
""" This class holds an exact Bernstein polynomial (represented as a list of algebraic real coefficients), and returns arbitrarily-precise interval approximations of this polynomial on demand. """
""" Initializes a bernstein_polynomial_factory_ar, given a polynomial with algebraic real coefficients. If neg is True, then gives the Bernstein polynomial for the negative half-line; if neg is False, the positive.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(AA) sage: p = (x - 1) * (x - sqrt(AA(2))) * (x - 2) sage: bernstein_polynomial_factory_ar(p, False) degree 3 Bernstein factory with algebraic real coefficients """ else:
""" Return a short summary of this Bernstein polynomial factory.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(AA) sage: p = (x - 1) * (x - sqrt(AA(2))) * (x - 2) sage: bernstein_polynomial_factory_ar(p, False) degree 3 Bernstein factory with algebraic real coefficients """
""" Computes the approximate log2 of the maximum of the absolute values of the coefficients.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(AA) sage: p = (x - 1) * (x - sqrt(AA(2))) * (x - 2) sage: bernstein_polynomial_factory_ar(p, False).coeffs_bitsize() 1 """
""" Compute an interval_bernstein_polynomial_integer that approximates this polynomial, using the given scale_log2. (Smaller scale_log2 values give more accurate approximations.)
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(AA) sage: p = (x - 1) * (x - sqrt(AA(2))) * (x - 2) sage: bpf = bernstein_polynomial_factory_ar(p, False) sage: print(bpf.bernstein_polynomial(0)) degree 3 IBP with 2-bit coefficients sage: bpf.bernstein_polynomial(-20) <IBP: ((-2965821, 2181961, -1542880, 1048576) + [0 .. 1)) * 2^-20> sage: bpf = bernstein_polynomial_factory_ar(p, True) sage: bpf.bernstein_polynomial(-20) <IBP: ((-2965821, -2181962, -1542880, -1048576) + [0 .. 1)) * 2^-20> sage: p = x^2 - 1 sage: bpf = bernstein_polynomial_factory_ar(p, False) sage: bpf.bernstein_polynomial(-10) <IBP: ((-1024, 0, 1024) + [0 .. 1)) * 2^-10> """
else:
# Compute lower and upper such that lower <= intv_c < upper
""" Given an interval Bernstein polynomial over a particular region (assumed to be a (not necessarily proper) subregion of [0 .. 1]), and a list of targets, uses de Casteljau's method to compute representations of the Bernstein polynomial over each target. Uses degree reduction as often as possible while maintaining the requested precision.
Each target is of the form (lgap, ugap, b). Suppose lgap.region() is (l1, l2), and ugap.region() is (u1, u2). Then we will compute an interval Bernstein polynomial over a region [l .. u], where l1 <= l <= l2 and u1 <= u <= u2. (split_for_targets() is free to select arbitrary region endpoints within these bounds; it picks endpoints which make the computation easier.) The third component of the target, b, is the maximum allowed scale_log2 of the result; this is used to decide when degree reduction is allowed.
The pair (l1, l2) can be replaced by None, meaning [-infinity .. 0]; or, (u1, u2) can be replaced by None, meaning [1 .. infinity].
There is another constraint on the region endpoints selected by split_for_targets() for a target ((l1, l2), (u1, u2), b). We set a size goal g, such that (u - l) <= g * (u1 - l2). Normally g is 256/255, but if precise is True, then g is 65536/65535.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: bp = mk_ibpi([1000000, -2000000, 3000000, -4000000, -5000000, -6000000]) sage: ctx = mk_context() sage: bps = split_for_targets(ctx, bp, [(rr_gap(1/1234567893, 1/1234567892, 1), rr_gap(1/1234567891, 1/1234567890, 1), 12), (rr_gap(1/3, 1/2, -1), rr_gap(2/3, 3/4, -1), 6)]) sage: bps[0] <IBP: (999992, 999992, 999992) + [0 .. 15) over [8613397477114467984778830327/10633823966279326983230456482242756608 .. 591908168025934394813836527495938294787/730750818665451459101842416358141509827966271488]; level 2; slope_err 0.?e12> sage: bps[1] <IBP: (-1562500, -1875001, -2222223, -2592593, -2969137, -3337450) + [0 .. 4) over [1/2 .. 2863311531/4294967296]> """ return []
cdef rr_gap l cdef rr_gap r
split_targets += [(QQ(0), None, 0)] else: else:
goal = Integer(65535)/65536 else:
(a, b, c, d) = bp.lft left = b/d if c+d != 0: right = (a+b)/(c+d) width = right-left err = (a/2 + b)/(c/2 + d) - (left+right)/2 else: width = RR(infinity) err = RR(infinity) slope = (a*d - b*c)/(d*d) ctx.dc_log_append(('lft', a, b, c, d, RR(left), RR(width), RR(slope/width), RR(err/width/width)))
if True: # p1.region_width() / bp.region_width() < tiny: else: if True: # p2.region_width() / bp.region_width() < tiny: else:
cdef class ocean: """ Given the tools we've defined so far, there are many possible root isolation algorithms that differ on where to select split points, what precision to work at when, and when to attempt degree reduction.
Here we implement one particular algorithm, which I call the ocean-island algorithm. We start with an interval Bernstein polynomial defined over the region [0 .. 1]. This region is the "ocean". Using de Casteljau's algorithm and Descartes' rule of signs, we divide this region into subregions which may contain roots, and subregions which are guaranteed not to contain roots. Subregions which may contain roots are "islands"; subregions known not to contain roots are "gaps".
All the real root isolation work happens in class island. See the documentation of that class for more information.
An island can be told to refine itself until it contains only a single root. This may not succeed, if the island's interval Bernstein polynomial does not have enough precision. The ocean basically loops, refining each of its islands, then increasing the precision of islands which did not succeed in isolating a single root; until all islands are done.
Increasing the precision of unsuccessful islands is done in a single pass using split_for_target(); this means it is possible to share work among multiple islands. """
def __init__(self, ctx, bpf, mapping): """ Initialize an ocean from a context and a Bernstein polynomial factory.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: ocean(mk_context(), bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]), lmap) ocean with precision 120 and 1 island(s) """
# Islands and gaps are maintained in a doubly-linked (circular) # list. Islands point to the gaps on either side, and gaps # point to the islands on either side.
# The list starts with self.lgap, which is the gap that starts # at 0; it ends at self.endpoint, which is a fictional island # after the gap that ends at 1.
# The initial gaps are unique in being 0-width; all gaps that # are created during execution of the algorithm have positive # width.
# Note that the constructor for islands side-effects its argument # gaps, so that they point to the island as their neighbor.
def __repr__(self): """ Return a short summary of this root isolator.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: ocean(mk_context(), bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]), lmap) ocean with precision 120 and 1 island(s) """
def _islands(self): """ Return a list of the islands in this ocean (for debugging purposes).
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]), lmap) sage: oc._islands() [island of width 1.00000000000000] """
def approx_bp(self, scale_log2): """ Returns an approximation to our Bernstein polynomial with the given scale_log2.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]), lmap) sage: oc.approx_bp(0) <IBP: (0, -4, 2, -2) + [0 .. 1); lsign 1> sage: oc.approx_bp(-20) <IBP: ((349525, -3295525, 2850354, -1482835) + [0 .. 1)) * 2^-20> """
def find_roots(self): """ Isolate all roots in this ocean.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]), lmap) sage: oc ocean with precision 120 and 1 island(s) sage: oc.find_roots() sage: oc ocean with precision 120 and 3 island(s) sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([1, 0, -1111/2, 0, 11108889/14, 0, 0, 0, 0, -1]), lmap) sage: oc.find_roots() sage: oc ocean with precision 240 and 3 island(s) """
def roots(self): """ Return the locations of all islands in this ocean. (If run after find_roots(), this is the location of all roots in the ocean.)
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]), lmap) sage: oc.find_roots() sage: oc.roots() [(1/32, 1/16), (1/2, 5/8), (3/4, 7/8)] sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([1, 0, -1111/2, 0, 11108889/14, 0, 0, 0, 0, -1]), lmap) sage: oc.find_roots() sage: oc.roots() [(95761241267509487747625/9671406556917033397649408, 191522482605387719863145/19342813113834066795298816), (1496269395904347376805/151115727451828646838272, 374067366568272936175/37778931862957161709568), (31/32, 63/64)] """
def refine_all(self): """ Refine all islands which are not done (which are not known to contain exactly one root).
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]), lmap) sage: oc ocean with precision 120 and 1 island(s) sage: oc.refine_all() sage: oc ocean with precision 120 and 3 island(s) """
def all_done(self): """ Returns true iff all islands are known to contain exactly one root.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]), lmap) sage: oc.all_done() False sage: oc.find_roots() sage: oc.all_done() True """
def reset_root_width(self, int isle_num, target_width): """ Require that the isle_num island have a width at most target_width.
If this is followed by a call to find_roots(), then the corresponding root will be refined to the specified width.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([-1, -1, 1]), lmap) sage: oc.find_roots() sage: oc.roots() [(1/2, 3/4)] sage: oc.reset_root_width(0, 1/2^200) sage: oc.find_roots() sage: oc ocean with precision 240 and 1 island(s) sage: RR(RealIntervalField(300)(oc.roots()[0]).absolute_diameter()).log2() -232.668979560890 """
def increase_precision(self): """ Increase the precision of the interval Bernstein polynomial held by any islands which are not done. (In normal use, calls to this function are separated by calls to self.refine_all().)
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([1/3, -22/7, 193/71, -140/99]), lmap) sage: oc ocean with precision 120 and 1 island(s) sage: oc.increase_precision() sage: oc.increase_precision() sage: oc.increase_precision() sage: oc ocean with precision 960 and 1 island(s) """
cdef class island: """ This implements the island portion of my ocean-island root isolation algorithm. See the documentation for class ocean, for more information on the overall algorithm.
Island root refinement starts with a Bernstein polynomial whose region is the whole island (or perhaps slightly more than the island in certain cases). There are two subalgorithms; one when looking at a Bernstein polynomial covering a whole island (so we know that there are gaps on the left and right), and one when looking at a Bernstein polynomial covering the left segment of an island (so we know that there is a gap on the left, but the right is in the middle of an island). An important invariant of the left-segment subalgorithm over the region [l .. r] is that it always finds a gap [r0 .. r] ending at its right endpoint.
Ignoring degree reduction, downscaling (precision reduction), and failures to split, the algorithm is roughly:
Whole island:
1. If the island definitely has exactly one root, then return. 2. Split the island in (approximately) half. 3. If both halves definitely have no roots, then remove this island from its doubly-linked list (merging its left and right gaps) and return. 4. If either half definitely has no roots, then discard that half and call the whole-island algorithm with the other half, then return. 5. If both halves may have roots, then call the left-segment algorithm on the left half. 6. We now know that there is a gap immediately to the left of the right half, so call the whole-island algorithm on the right half, then return.
Left segment:
1. Split the left segment in (approximately) half. 2. If both halves definitely have no roots, then extend the left gap over the segment and return. 3. If the left half definitely has no roots, then extend the left gap over this half and call the left-segment algorithm on the right half, then return. 4. If the right half definitely has no roots, then split the island in two, creating a new gap. Call the whole-island algorithm on the left half, then return. 5. Both halves may have roots. Call the left-segment algorithm on the left half. 6. We now know that there is a gap immediately to the left of the right half, so call the left-segment algorithm on the right half, then return.
Degree reduction complicates this picture only slightly. Basically, we use heuristics to decide when degree reduction might be likely to succeed and be helpful; whenever this is the case, we attempt degree reduction.
Precision reduction and split failure add more complications. The algorithm maintains a stack of different-precision representations of the interval Bernstein polynomial. The base of the stack is at the highest (currently known) precision; each stack entry has approximately half the precision of the entry below it. When we do a split, we pop off the top of the stack, split it, then push whichever half we're interested in back on the stack (so the different Bernstein polynomials may be over different regions). When we push a polynomial onto the stack, we may heuristically decide to push further lower-precision versions of the same polynomial onto the stack.
In the algorithm above, whenever we say "split in (approximately) half", we attempt to split the top-of-stack polynomial using try_split() and try_rand_split(). However, these will fail if the sign of the polynomial at the chosen split point is unknown (if the polynomial is not known to high enough precision, or if the chosen split point actually happens to be a root of the polynomial). If this fails, then we discard the top-of-stack polynomial, and try again with the next polynomial down (which has approximately twice the precision). This next polynomial may not be over the same region; if not, we split it using de Casteljau's algorithm to get a polynomial over (approximately) the same region first.
If we run out of higher-precision polynomials (if we empty out the entire stack), then we give up on root refinement for this island. The ocean class will notice this, provide the island with a higher-precision polynomial, and restart root refinement. Basically the only information kept in that case is the lower and upper bounds on the island. Since these are updated whenever we discover a "half" (of an island or a segment) that definitely contains no roots, we never need to re-examine these gaps. (We could keep more information. For example, we could keep a record of split points that succeeded and failed. However, a split point that failed at lower precision is likely to succeed at higher precision, so it's not worth avoiding. It could be useful to select split points that are known to succeed, but starting from a new Bernstein polynomial over a slightly different region, hitting such split points would require de Casteljau splits with non-power-of-two denominators, which are much much slower.) """
def __init__(self, interval_bernstein_polynomial bp, rr_gap lgap, rr_gap rgap): """ Initialize an island from a Bernstein polynomial, and the gaps to the left and right of the island. """
def __repr__(self): """ Return a short summary of this island. """
def _info(self): # A python accessor for the information in this island. # For debugging purposes. return (self.bp, self.ancestors, self.target_width, self.lgap, self.rgap, self.known_done)
def shrink_bp(self, context ctx): """ If the island's Bernstein polynomial covers a region much larger than the island itself (in particular, if either the island's left gap or right gap are totally contained in the polynomial's region) then shrink the polynomial down to cover the island more tightly. """
else:
def refine(self, context ctx): """ Attempts to shrink and/or split this island into sub-island that each definitely contain exactly one root. """ pass
def refine_recurse(self, context ctx, interval_bernstein_polynomial bp, ancestors, history, rightmost): """ This implements the root isolation algorithm described in the class documentation for class island. This is the implementation of both the whole-island and the left-segment algorithms; if the flag rightmost is True, then it is the whole-island algorithm, otherwise the left-segment algorithm.
The precision-reduction stack is (ancestors + [bp]); that is, the top-of-stack is maintained separately. """ cdef interval_bernstein_polynomial p1, p2 cdef rr_gap mgap cdef island lisland
# If you examine the description of this algorithm in the # class island class documentation, you see that several of # the calls involved are tail recursions (that is, they are # of the form "call (this algorithm) recursively, then return"). # We perform tail-recursion optimization by hand: such calls # assign new values to the method parameters, and then fall off # the end of the following "while True" loop, and return to # the top of the method here.
# This optimization is VITAL; without it, several test polynomials # (in particular, polynomials with roots that are very, very # close together) give stack overflow errors.
# No roots! Make the island disappear. self.lgap.risland = self.rgap.risland self.rgap.risland.lgap = self.lgap self.lgap.upper = self.rgap.upper else:
# This is our heuristic for deciding when to do degree # reduction.
# In general, given a degree-d Bernstein polynomial, # if you perform a de Casteljau split at 1/2, the # coefficients of the resulting polynomials are about d bits # smaller than the coefficients of the original polynomial.
# Conversely, if you do a split and the coefficients are # k bits smaller than the coefficients of the original # polynomial, this indicates that the original polynomial # may be very close to some degree-k polynomial.
# Here, we look at the amount that the coefficients got smaller # over the last 3 splits. If they got more than 100 bits # smaller (that is, an average of more than 33 bits per split), # then we expect that the original polynomial is not # close to any polynomial of degree 30 or less. Since degree # reduction works by finding a polynomial of (currently) # degree 30 that is close to the original polynomial, # then this drastic reduction in coefficient size means that # degree reduction is likely to fail, so we don't bother # attempting it.
# Note that this does not take into account the effects of # random splitting; maybe it should? For example, if the # last three splits were random splits with split points near # 1/4, and we're in the left-hand branch of all these splits, # then we would expect twice as much coefficient reduction; # so a coefficient size drop of 100 bits would actually # mean that we are near a degree-17 polynomial.
# Also, we are willing to throw away up to 70 + msb_delta # bits of precision during degree reduction. I don't remember # how I selected this heuristic...
# Heuristically push more lower-precision polynomials # on the "ancestors" stack
# Currently, we try a random split only once before giving up # and trying for higher precision; and then we try a 1/2 split # at the higher precision. Maybe we should try more random # splits before going to higher precision? Maybe we should # try less 1/2 splits?
# ctx.dc_log_append(('split_results', p1.variations(), p2.variations()))
# No roots! Make the island disappear. else: # return to top of function (tail recursion optimization) else: # return to top of function else: # Split the island! pass else: # return to top of function (tail recursion optimization)
def less_bits(self, ancestors, interval_bernstein_polynomial bp): """ Heuristically pushes lower-precision polynomials on the polynomial stack. See the class documentation for class island for more information. """ # The current heuristic is to always push one lower-precision # polynomial, unless the current polynomial is floating-point. # If the current polynomial has bitsize < 130, then the new # polynomial is floating-point, otherwise it is integer, # with half the precision of the original.
else: else:
def more_bits(self, context ctx, ancestors, interval_bernstein_polynomial bp, rightmost): """ Find a Bernstein polynomial on the "ancestors" stack with more precision than bp; if it is over a different region, then shrink its region to (approximately) match that of bp. (If this is rightmost -- if bp covers the whole island -- then we only require that the new region cover the whole island fairly tightly; if this is not rightmost, then the new region will have exactly the same right boundary as bp, although the left boundary may vary slightly.) """
cdef interval_bernstein_polynomial anc cdef interval_bernstein_polynomial ancestor_val
return (ancestors + [ancestor_val], ancestor_val.as_float()) else:
first_msb, first_lsb, rel_width_rr, new_lsb, target_lsb_h, target_lsb, target_lsb_l))
else: ancestor_val = ancestor_val.down_degree_iter(ctx, target_lsb_h)
# if rel_lbounds[1] > 0: # left_split = -exact_rational(simple_wordsize_float(-rel_lbounds[1], -rel_lbounds[0])) # (_, ancestor_val, _) = ancestor_val.de_casteljau(ctx, left_split) # ctx.dc_log_append(('pull_left', left_split))
def reset_root_width(self, target_width): """ Modify the criteria for this island to require that it is not "done" until its width is less than or equal to target_width. """
def bp_done(self, interval_bernstein_polynomial bp): """ Examine the given Bernstein polynomial to see if it is known to have exactly one root in its region. (In addition, we require that the polynomial region not include 0 or 1. This makes things work if the user gives explicit bounds to real_roots(), where the lower or upper bound is a root of the polynomial. real_roots() deals with this by explicitly detecting it, dividing out the appropriate linear polynomial, and adding the root to the returned list of roots; but then if the island considers itself "done" with a region including 0 or 1, the returned root regions can overlap with each other.) """
def done(self, context ctx): """ Check to see if the island is known to contain zero roots or is known to contain one root. """
else:
def has_root(self): """ Assuming that the island is done (has either 0 or 1 roots), reports whether the island has a root. """
cdef class rr_gap: """ A simple class representing the gaps between islands, in my ocean-island root isolation algorithm. Named "rr_gap" for "real roots gap", because "gap" seemed too short and generic. """
def __init__(self, lower, upper, sign): """ Initialize an rr_gap element. """
def region(self):
""" A simple class to map linearly between original coordinates (ranging from [lower .. upper]) and ocean coordinates (ranging from [0 .. 1]). """
""" A class to map between original coordinates and ocean coordinates. If neg is False, then the original->ocean transform is x -> x/(x+1), and the ocean->original transform is x/(1-x); this maps between [0 .. infinity] and [0 .. 1]. If neg is True, then the original->ocean transform is x -> -x/(1-x), and the ocean->original transform is the same thing: -x/(1-x). This maps between [0 .. -infinity] and [0 .. 1]. """
else:
return (-u/(1-u), -l/(1-l)) else:
""" Compute the real roots of a given polynomial with exact coefficients (integer, rational, and algebraic real coefficients are supported). Returns a list of pairs of a root and its multiplicity.
The root itself can be returned in one of three different ways. If retval=='rational', then it is returned as a pair of rationals that define a region that includes exactly one root. If retval=='interval', then it is returned as a RealIntervalFieldElement that includes exactly one root. If retval=='algebraic_real', then it is returned as an AlgebraicReal. In the former two cases, all the intervals are disjoint.
An alternate high-level algorithm can be used by selecting strategy='warp'. This affects the conversion into Bernstein polynomial form, but still uses the same ocean-island algorithm as the default algorithm. The 'warp' algorithm performs the conversion into Bernstein polynomial form much more quickly, but performs the rest of the computation slightly slower in some benchmarks. The 'warp' algorithm is particularly likely to be helpful for low-degree polynomials.
Part of the algorithm is randomized; the seed parameter gives a seed for the random number generator. (By default, the same seed is used for every call, so that results are repeatable.) The random seed may affect the running time, or the exact intervals returned, but the results are correct regardless of the seed used.
The bounds parameter lets you find roots in some proper subinterval of the reals; it takes a pair of a rational lower and upper bound and only roots within this bound will be found. Currently, specifying bounds does not work if you select strategy='warp', or if you use a polynomial with algebraic real coefficients.
By default, the algorithm will do a squarefree decomposition to get squarefree polynomials. The skip_squarefree parameter lets you skip this step. (If this step is skipped, and the polynomial has a repeated real root, then the algorithm will loop forever! However, repeated non-real roots are not a problem.)
For integer and rational coefficients, the squarefree decomposition is very fast, but it may be slow for algebraic reals. (It may trigger exact computation, so it might be arbitrarily slow. The only other way that this algorithm might trigger exact computation on algebraic real coefficients is that it checks the constant term of the input polynomial for equality with zero.)
Part of the algorithm works (approximately) by splitting numbers into word-size pieces (that is, pieces that fit into a machine word). For portability, this defaults to always selecting pieces suitable for a 32-bit machine; the wordsize parameter lets you make choices suitable for a 64-bit machine instead. (This affects the running time, and the exact intervals returned, but the results are correct on both 32- and 64-bit machines even if the wordsize is chosen "wrong".)
The precision of the results can be improved (at the expense of time, of course) by specifying the max_diameter parameter. If specified, this sets the maximum diameter() of the intervals returned. (Sage defines diameter() to be the relative diameter for intervals that do not contain 0, and the absolute diameter for intervals containing 0.) This directly affects the results in rational or interval return mode; in algebraic_real mode, it increases the precision of the intervals passed to the algebraic number package, which may speed up some operations on that algebraic real.
Some logging can be enabled with do_logging=True. If logging is enabled, then the normal values are not returned; instead, a pair of the internal context object and a list of all the roots in their internal form is returned.
ALGORITHM: We convert the polynomial into the Bernstein basis, and then use de Casteljau's algorithm and Descartes' rule of signs (using interval arithmetic) to locate the roots.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(ZZ) sage: real_roots(x^3 - x^2 - x - 1) [((7/4, 19/8), 1)] sage: real_roots((x-1)*(x-2)*(x-3)*(x-5)*(x-8)*(x-13)*(x-21)*(x-34)) [((11/16, 33/32), 1), ((11/8, 33/16), 1), ((11/4, 55/16), 1), ((77/16, 165/32), 1), ((11/2, 33/4), 1), ((11, 55/4), 1), ((165/8, 341/16), 1), ((22, 44), 1)] sage: real_roots(x^5 * (x^2 - 9999)^2 - 1) [((-29274496381311/9007199254740992, 419601125186091/2251799813685248), 1), ((2126658450145849453951061654415153249597/21267647932558653966460912964485513216, 4253316902721330018853696359533061621799/42535295865117307932921825928971026432), 1), ((1063329226287740282451317352558954186101/10633823966279326983230456482242756608, 531664614358685696701445201630854654353/5316911983139663491615228241121378304), 1)] sage: real_roots(x^5 * (x^2 - 9999)^2 - 1, seed=42) [((-123196838480289/18014398509481984, 293964743458749/9007199254740992), 1), ((8307259573979551907841696381986376143/83076749736557242056487941267521536, 16614519150981033789137940378745325503/166153499473114484112975882535043072), 1), ((519203723562592617581015249797434335/5192296858534827628530496329220096, 60443268924081068060312183/604462909807314587353088), 1)] sage: real_roots(x^5 * (x^2 - 9999)^2 - 1, wordsize=64) [((-62866503803202151050003/19342813113834066795298816, 901086554512564177624143/4835703278458516698824704), 1), ((544424563237337315214990987922809050101157/5444517870735015415413993718908291383296, 1088849127096660194637118845654929064385439/10889035741470030830827987437816582766592), 1), ((272212281929661439711063928866060007142141/2722258935367507707706996859454145691648, 136106141275823501959100399337685485662633/1361129467683753853853498429727072845824), 1)] sage: real_roots(x) [((-47/256, 81/512), 1)] sage: real_roots(x * (x-1)) [((-47/256, 81/512), 1), ((1/2, 1201/1024), 1)] sage: real_roots(x-1) [((209/256, 593/512), 1)] sage: real_roots(x*(x-1)*(x-2), bounds=(0, 2)) [((0, 0), 1), ((81/128, 337/256), 1), ((2, 2), 1)] sage: real_roots(x*(x-1)*(x-2), bounds=(0, 2), retval='algebraic_real') [(0, 1), (1, 1), (2, 1)] sage: v = 2^40 sage: real_roots((x^2-1)^2 * (x^2 - (v+1)/v)) [((-12855504354077768210885019021174120740504020581912910106032833/12855504354071922204335696738729300820177623950262342682411008, -6427752177038884105442509510587059395588605840418680645585479/6427752177035961102167848369364650410088811975131171341205504), 1), ((-1125899906842725/1125899906842624, -562949953421275/562949953421312), 2), ((62165404551223330269422781018352603934643403586760330761772204409982940218804935733653/62165404551223330269422781018352605012557018849668464680057997111644937126566671941632, 3885337784451458141838923813647037871787041539340705594199885610069035709862106085785/3885337784451458141838923813647037813284813678104279042503624819477808570410416996352), 2), ((509258994083853105745586001837045839749063767798922046787130823804169826426726965449697819/509258994083621521567111422102344540262867098416484062659035112338595324940834176545849344, 25711008708155536421770038042348240136257704305733983563630791/25711008708143844408671393477458601640355247900524685364822016), 1)] sage: real_roots(x^2 - 2) [((-3/2, -1), 1), ((1, 3/2), 1)] sage: real_roots(x^2 - 2, retval='interval') [(-2.?, 1), (2.?, 1)] sage: real_roots(x^2 - 2, max_diameter=1/2^30) [((-22506280506048041472675379598886543645348790970912519198456805737131269246430553365310109/15914343565113172548972231940698266883214596825515126958094847260581103904401068017057792, -45012561012096082945350759197773086524448972309421182031053065599548946985601579935498343/31828687130226345097944463881396533766429193651030253916189694521162207808802136034115584), 1), ((45012561012096082945350759197773086524448972309421182031053065599548946985601579935498343/31828687130226345097944463881396533766429193651030253916189694521162207808802136034115584, 22506280506048041472675379598886543645348790970912519198456805737131269246430553365310109/15914343565113172548972231940698266883214596825515126958094847260581103904401068017057792), 1)] sage: real_roots(x^2 - 2, retval='interval', max_diameter=1/2^500) [(-1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571470109559971605970274534596862014728517418640889198609552?, 1), (1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571470109559971605970274534596862014728517418640889198609552?, 1)] sage: ar_rts = real_roots(x^2 - 2, retval='algebraic_real'); ar_rts [(-1.414213562373095?, 1), (1.414213562373095?, 1)] sage: ar_rts[0][0]^2 - 2 == 0 True sage: v = 2^40 sage: real_roots((x-1) * (x-(v+1)/v), retval='interval') [(1.000000000000?, 1), (1.000000000001?, 1)] sage: v = 2^60 sage: real_roots((x-1) * (x-(v+1)/v), retval='interval') [(1.000000000000000000?, 1), (1.000000000000000001?, 1)] sage: real_roots((x-1) * (x-2), strategy='warp') [((499/525, 1173/875), 1), ((337/175, 849/175), 1)] sage: real_roots((x+3)*(x+1)*x*(x-1)*(x-2), strategy='warp') [((-1713/335, -689/335), 1), ((-2067/2029, -689/1359), 1), ((0, 0), 1), ((499/525, 1173/875), 1), ((337/175, 849/175), 1)] sage: real_roots((x+3)*(x+1)*x*(x-1)*(x-2), strategy='warp', retval='algebraic_real') [(-3.000000000000000?, 1), (-1.000000000000000?, 1), (0, 1), (1.000000000000000?, 1), (2.000000000000000?, 1)] sage: ar_rts = real_roots(x-1, retval='algebraic_real') sage: ar_rts[0][0] == 1 True
If the polynomial has no real roots, we get an empty list.
::
sage: (x^2 + 1).real_root_intervals() []
We can compute Conway's constant (see http://mathworld.wolfram.com/ConwaysConstant.html) to arbitrary precision. ::
sage: p = x^71 - x^69 - 2*x^68 - x^67 + 2*x^66 + 2*x^65 + x^64 - x^63 - x^62 - x^61 - x^60 - x^59 + 2*x^58 + 5*x^57 + 3*x^56 - 2*x^55 - 10*x^54 - 3*x^53 - 2*x^52 + 6*x^51 + 6*x^50 + x^49 + 9*x^48 - 3*x^47 - 7*x^46 - 8*x^45 - 8*x^44 + 10*x^43 + 6*x^42 + 8*x^41 - 5*x^40 - 12*x^39 + 7*x^38 - 7*x^37 + 7*x^36 + x^35 - 3*x^34 + 10*x^33 + x^32 - 6*x^31 - 2*x^30 - 10*x^29 - 3*x^28 + 2*x^27 + 9*x^26 - 3*x^25 + 14*x^24 - 8*x^23 - 7*x^21 + 9*x^20 + 3*x^19 - 4*x^18 - 10*x^17 - 7*x^16 + 12*x^15 + 7*x^14 + 2*x^13 - 12*x^12 - 4*x^11 - 2*x^10 + 5*x^9 + x^7 - 7*x^6 + 7*x^5 - 4*x^4 + 12*x^3 - 6*x^2 + 3*x - 6 sage: cc = real_roots(p, retval='algebraic_real')[2][0] # long time sage: RealField(180)(cc) # long time 1.3035772690342963912570991121525518907307025046594049
Now we play with algebraic real coefficients. ::
sage: x = polygen(AA) sage: p = (x - 1) * (x - sqrt(AA(2))) * (x - 2) sage: real_roots(p) [((499/525, 2171/1925), 1), ((1173/875, 2521/1575), 1), ((337/175, 849/175), 1)] sage: ar_rts = real_roots(p, retval='algebraic_real'); ar_rts [(1.000000000000000?, 1), (1.414213562373095?, 1), (2.000000000000000?, 1)] sage: ar_rts[1][0]^2 == 2 True sage: ar_rts = real_roots(x*(x-1), retval='algebraic_real') sage: ar_rts[0][0] == 0 True sage: p2 = p * (p - 1/100); p2 x^6 - 8.82842712474619?*x^5 + 31.97056274847714?*x^4 - 60.77955262170047?*x^3 + 63.98526763257801?*x^2 - 35.37613490585595?*x + 8.028284271247462? sage: real_roots(p2, retval='interval') [(1.00?, 1), (1.1?, 1), (1.38?, 1), (1.5?, 1), (2.00?, 1), (2.1?, 1)] sage: p = (x - 1) * (x - sqrt(AA(2)))^2 * (x - 2)^3 * sqrt(AA(3)) sage: real_roots(p, retval='interval') [(1.000000000000000?, 1), (1.414213562373095?, 2), (2.000000000000000?, 3)]
TESTS:
Check that :trac:`20269` is fixed::
sage: real_roots(polygen(AA))[0][0][0].parent() Rational Field
Check that :trac:`10803` is fixed::
sage: f = 2503841067*x^13 - 15465014877*x^12 + 37514382885*x^11 - 44333754994*x^10 + 24138665092*x^9 - 2059014842*x^8 - 3197810701*x^7 + 803983752*x^6 + 123767204*x^5 - 26596986*x^4 - 2327140*x^3 + 75923*x^2 + 7174*x + 102 sage: len(real_roots(f,seed=1)) 13 """
else: raise ValueError("Don't know how to isolate roots for " + str(p.parent()))
raise NotImplementedError("Cannot set your own bounds with algebraic real input")
raise NotImplementedError("Cannot set your own bounds with strategy=warp")
factors = [(p, 1)] else:
cdef ocean oc
else:
else: else: else: else: # Bad things happen if the bounds are roots themselves. # Avoid this by dividing out linear polynomials if # the bounds are roots.
# Check to make sure that no intervals are too wide.
# We use half_diameter, because if we ended up with a rational # interval that was exactly max_diameter, we might not # be able to coerce it into an interval small enough.
else:
# Check to be sure that all intervals are disjoint.
return ctx, all_roots
# The following line should work, but does not due to a Cython # bug (it calls PyObject_Cmp() for comparison operators, # instead of PyObject_RichCompare()).
# if not (intv_roots[j][0] < intv_roots[j+1][0]):
raise ValueError("Illegal retval parameter " + retval)
""" Given a vector of integers c of length n+1, and a rational k == kn / kd, multiplies each element c[i] by (kd^i)*(kn^(n-i)).
Modifies the input vector; has no return value.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: v = vector(ZZ, [1, 1, 1, 1]) sage: scale_intvec_var(v, 3/4) sage: v (64, 48, 36, 27) """
# XXX could use direct mpz calls for more speed
cdef int i
""" Given a vector of integers c of length d+1, representing the coefficients of a degree-d polynomial p, modify the vector to perform a Taylor shift by 1 (that is, p becomes p(x+1)).
This is the straightforward algorithm, which is not asymptotically optimal.
Modifies the input vector; has no return value.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(ZZ) sage: p = (x-1)*(x-2)*(x-3) sage: v = vector(ZZ, p.list()) sage: p, v (x^3 - 6*x^2 + 11*x - 6, (-6, 11, -6, 1)) sage: taylor_shift1_intvec(v) sage: p(x+1), v (x^3 - 3*x^2 + 2*x, (0, 2, -3, 1)) """
cdef int i, k
""" Given a vector of integers, reverse the vector (like the reverse() method on lists).
Modifies the input vector; has no return value.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: v = vector(ZZ, [1, 2, 3, 4]); v (1, 2, 3, 4) sage: reverse_intvec(v) sage: v (4, 3, 2, 1) """ cdef int i
""" A simple cache for RealField fields (with rounding set to round-to-positive-infinity).
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: get_realfield_rndu(20) Real Field with 20 bits of precision and rounding RNDU sage: get_realfield_rndu(53) Real Field with 53 bits of precision and rounding RNDU sage: get_realfield_rndu(20) Real Field with 20 bits of precision and rounding RNDU """
cdef class context: """ A simple context class, which is passed through parts of the real root isolation algorithm to avoid global variables.
Holds logging information, a random number generator, and the target machine wordsize. """
def __init__(self, do_logging, seed, wordsize): """ Initialize a context class. """
def __repr__(self): """ Return a short summary of this context.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: mk_context() root isolation context: seed=0 sage: mk_context(do_logging=True, seed=37, wordsize=64) root isolation context: seed=37; do_logging=True; wordsize=64 """
cdef void dc_log_append(self, x): """ Optional logging for the root isolation algorithm. """ self.dc_log.append(x)
cdef void be_log_append(self, x): """ Optional logging for degree reduction in the root isolation algorithm. """ self.be_log.append(x)
def get_dc_log(self): return self.dc_log
def get_be_log(self): return self.be_log
""" A simple wrapper for creating context objects with coercions, defaults, etc.
For use in doctests.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: mk_context(do_logging=True, seed=3, wordsize=64) root isolation context: seed=3; do_logging=True; wordsize=64 """
""" Given a polynomial p with integer coefficients, and rational bounds low and high, compute the exact rational Bernstein coefficients of p over the region [low .. high]. The optional parameter degree can be used to give a formal degree higher than the actual degree.
The return value is a pair (c, scale); c represents the same polynomial as p*scale. (If you only care about the roots of the polynomial, then of course scale can be ignored.)
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(ZZ) sage: to_bernstein(x) ([0, 1], 1) sage: to_bernstein(x, degree=5) ([0, 1/5, 2/5, 3/5, 4/5, 1], 1) sage: to_bernstein(x^3 + x^2 - x - 1, low=-3, high=3) ([-16, 24, -32, 32], 1) sage: to_bernstein(x^3 + x^2 - x - 1, low=3, high=22/7) ([296352, 310464, 325206, 340605], 9261) """
raise ValueError('Bernstein degree must be at least polynomial degree') else:
""" Given a polynomial p with rational coefficients, compute the exact rational Bernstein coefficients of p(x/(x+1)).
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: x = polygen(ZZ) sage: to_bernstein_warp(1 + x + x^2 + x^3 + x^4 + x^5) [1, 1/5, 1/10, 1/10, 1/5, 1] """
""" Given an integer vector representing a Bernstein polynomial p, and a degree d2, compute the representation of p as a Bernstein polynomial of formal degree d2.
This is similar to multiplying by the result of bernstein_up, but should be faster for large d2 (this has about the same number of multiplies, but in this version all the multiplies are by single machine words).
Returns a pair consisting of the expanded polynomial, and the maximum error E. (So if an element of the returned polynomial is a, and the true value of that coefficient is b, then a <= b < a + E.)
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: c = vector(ZZ, [1000, 2000, -3000]) sage: bernstein_expand(c, 3) ((1000, 1666, 333, -3000), 1) sage: bernstein_expand(c, 4) ((1000, 1500, 1000, -500, -3000), 1) sage: bernstein_expand(c, 20) ((1000, 1100, 1168, 1205, 1210, 1184, 1126, 1036, 915, 763, 578, 363, 115, -164, -474, -816, -1190, -1595, -2032, -2500, -3000), 1) """
cdef int i, j
cdef mpz_t tmp cdef mpz_t divisor
# XXX do experimentation here on how to decide when to divide
cdef int max_bitsize_intvec(Vector_integer_dense b): """ Given an integer vector, find the approximate log2 of the maximum of the absolute values of the elements.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: max_bitsize_intvec_doctest(vector(ZZ, [1, 2, 3, 1024])) 11 """ cdef int i cdef int size
""" Computes the dot product of row k of the matrix m with the vector v (that is, compute one element of the product m*v).
If v has more elements than m has columns, then elements of v are selected using subsample_vec.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: m = matrix(3, range(9)) sage: dprod_imatrow_vec(m, vector(ZZ, [1, 0, 0, 0]), 1) 0 sage: dprod_imatrow_vec(m, vector(ZZ, [0, 1, 0, 0]), 1) 3 sage: dprod_imatrow_vec(m, vector(ZZ, [0, 0, 1, 0]), 1) 4 sage: dprod_imatrow_vec(m, vector(ZZ, [0, 0, 0, 1]), 1) 5 sage: dprod_imatrow_vec(m, vector(ZZ, [1, 0, 0]), 1) 3 sage: dprod_imatrow_vec(m, vector(ZZ, [0, 1, 0]), 1) 4 sage: dprod_imatrow_vec(m, vector(ZZ, [0, 0, 1]), 1) 5 sage: dprod_imatrow_vec(m, vector(ZZ, [1, 2, 3]), 1) 26 """
cdef mpz_t tmp
cdef int ra cdef int a
""" Given two integer vectors a and b (of equal, nonzero length), return a pair of the minimum and maximum values taken on by a[i] - b[i].
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: a = vector(ZZ, [10, -30]) sage: b = vector(ZZ, [15, -60]) sage: min_max_delta_intvec(a, b) (30, -5) """
cdef mpz_t tmp
cdef int i
""" Given an integer vector b = (b0, ..., bn), compute the minimum and maximum values of b_{j+1} - b_j.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: min_max_diff_intvec(vector(ZZ, [1, 7, -2])) (-9, 6) """
""" Given a floating-point vector b = (b0, ..., bn), compute the minimum and maximum values of b_{j+1} - b_j.
EXAMPLES::
sage: from sage.rings.polynomial.real_roots import * sage: min_max_diff_doublevec(vector(RDF, [1, 7, -2])) (-9.0, 6.0) """
cdef double diff
|