| Home | Trees | Indices | Help |
|---|
|
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 1)
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 http://www.sagemath.org:9001/days4schedule .
|
|||
| PrecisionError | |||
|
bernstein_polynomial_factory File: sage/rings/polynomial/real_roots.pyx (starting at line 2410) An abstract base class for bernstein_polynomial factories. |
|||
|
bernstein_polynomial_factory_ar File: sage/rings/polynomial/real_roots.pyx (starting at line 2601) 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. |
|||
|
bernstein_polynomial_factory_intlist File: sage/rings/polynomial/real_roots.pyx (starting at line 2447) 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. |
|||
|
bernstein_polynomial_factory_ratlist File: sage/rings/polynomial/real_roots.pyx (starting at line 2521) 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. |
|||
|
context File: sage/rings/polynomial/real_roots.pyx (starting at line 4138) A simple context class, which is passed through parts of the real root isolation algorithm to avoid global variables. |
|||
|
interval_bernstein_polynomial File: sage/rings/polynomial/real_roots.pyx (starting at line 170) An interval_bernstein_polynomial is an approximation to an exact polynomial. |
|||
|
interval_bernstein_polynomial_float File: sage/rings/polynomial/real_roots.pyx (starting at line 1299) This is the subclass of interval_bernstein_polynomial where polynomial coefficients are represented using floating-point numbers. |
|||
|
interval_bernstein_polynomial_integer File: sage/rings/polynomial/real_roots.pyx (starting at line 436) This is the subclass of interval_bernstein_polynomial where polynomial coefficients are represented using integers. |
|||
|
island File: sage/rings/polynomial/real_roots.pyx (starting at line 3112) This implements the island portion of my ocean-island root isolation algorithm. |
|||
|
linear_map File: sage/rings/polynomial/real_roots.pyx (starting at line 3635) A simple class to map linearly between original coordinates (ranging from [lower .. |
|||
|
ocean File: sage/rings/polynomial/real_roots.pyx (starting at line 2834) 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. |
|||
|
rr_gap File: sage/rings/polynomial/real_roots.pyx (starting at line 3617) A simple class representing the gaps between islands, in my ocean-island root isolation algorithm. |
|||
|
warp_map File: sage/rings/polynomial/real_roots.pyx (starting at line 3657) A class to map between original coordinates and ocean coordinates. |
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
|
|||
AA = |
|||
QQ = Rational Field
|
|||
RDF = Real Double Field
|
|||
RIF = Real Interval Field with 53 bits of precision
|
|||
RR = Real Field with 53 bits of precision
|
|||
ZZ = |
|||
dr_cache = |
|||
infinity = +Infinity
|
|||
lmap = <sage.rings.polynomial.real_roots.linear_map instance a
|
|||
realfield_rndu_cache = |
|||
|
|||
File: sage/rings/polynomial/real_roots.pyx (starting at line 1987)
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]
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4280)
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)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 2026)
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]
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 2217)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 2108)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 2183)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 1658)
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)
((0.7, 0.466666666667, 0.311111111111, 0.207407407407, 0.138271604938, 0.0921810699588), (0.0921810699588, 0.0, 0.0, 0.0, 0.0, 0.0), 15)
sage: de_casteljau_doublevec(c, 7/22)
((0.7, 0.477272727273, 0.32541322314, 0.221872652141, 0.151276808278, 0.103143278371), (0.103143278371, 0.0, 0.0, 0.0, 0.0, 0.0), 15)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 971)
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)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 1903)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4370)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4117)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 1239)
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
(bk + F1) * 2^d <= ak
and
ak + E <= (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.1102230246251565e-16, 0.37500000000000017, 3)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 1746)
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.76700000000000002
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 2086)
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
sage: maximum_root_first_lambda(x^2 - 1)
1.00000000000000
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 2161)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4417)
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)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4480)
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)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4451)
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)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4198)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 1632)
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: 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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 956)
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: mk_ibpi([50, 20, -90, -70, 200], error=5)
degree 4 IBP with 8-bit coefficients
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 1937)
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]
(3, [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]
[ -2/63 11/126 -2/63 -1/9 8/63 121/126], 2, ([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))
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 2301)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 3685)
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.00000000000000000 .. 0.00000000000000000], 1), ([1.0000000000000000 .. 1.0000000000000000], 1), ([2.0000000000000000 .. 2.0000000000000000], 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')
[([-1.5000000000000000 .. -1.0000000000000000], 1), ([1.0000000000000000 .. 1.5000000000000000], 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.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095530 .. -1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095519], 1), ([1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095519 .. 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095530], 1)]
sage: ar_rts = real_roots(x^2 - 2, retval='algebraic_real'); ar_rts
[([-1.4142135623730952 .. -1.4142135623730949], 1), ([1.4142135623730949 .. 1.4142135623730952], 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')
[([0.99999999999977251 .. 1.0000000000001137], 1), ([1.0000000000004545 .. 1.0000000000018190], 1)]
sage: v = 2^60
sage: real_roots((x-1) * (x-(v+1)/v), retval='interval')
[([0.999999999999999999941850329426260 .. 1.00000000000000000016070246737549], 1), ([1.00000000000000000037955460532465 .. 1.00000000000000000125496315712147], 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.0000000000000005 .. -2.9999999999999995], 1), ([-1.0000000000000003 .. -0.99999999999999988], 1), ([0.00000000000000000 .. 0.00000000000000000], 1), ([0.99999999999999988 .. 1.0000000000000003], 1), ([1.9999999999999997 .. 2.0000000000000005], 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
[([0.99999999999999988 .. 1.0000000000000003], 1), ([1.4142135623730949 .. 1.4142135623730952], 1), ([1.9999999999999997 .. 2.0000000000000005], 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.8284271247461917 .. -8.8284271247461898]*x^5 + [31.970562748477139 .. 31.970562748477143]*x^4 + [-60.779552621700475 .. -60.779552621700467]*x^3 + [63.985267632578008 .. 63.985267632578016]*x^2 + [-35.376134905855956 .. -35.376134905855948]*x + [8.0282842712474611 .. 8.0282842712474630]
sage: real_roots(p2, retval='interval')
[([0.99197568389057744 .. 1.0026279602750193], 1), ([1.0133947772657450 .. 1.0352795031055902], 1), ([1.3701989150090414 .. 1.3852957233848955], 1), ([1.4005860805860805 .. 1.4637593984962408], 1), ([1.9959314285714284 .. 2.0019352991697686], 1), ([2.0079632816982213 .. 2.0200921658986180], 1)]
sage: p = (x - 1) * (x - sqrt(AA(2)))^2 * (x - 2)^3 * sqrt(AA(3))
sage: real_roots(p, retval='interval')
[([0.99999999999999988 .. 1.0000000000000003], 1), ([1.4142135623730949 .. 1.4142135623730952], 2), ([1.9999999999999997 .. 2.0000000000000005], 3)]
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 1856)
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)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4095)
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)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 2237)
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, 0)
sage: root_bounds(x*(x+1))
(-1.00000000000000, 0)
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
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4040)
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)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 2709)
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: str(bps[0])
'<IBP: (999992, 999992, 999992) + [0 .. 15) over [8613397477114467984778830327/10633823966279326983230456482242756608 .. 591908168025934394813836527495938294787/730750818665451459101842416358141509827966271488]; level 2; slope_err [-1.9259259099049338e11 .. 1.9259259099049338e11]>'
sage: str(bps[1])
'<IBP: (-1562500, -1875001, -2222223, -2592593, -2969137, -3337450) + [0 .. 4) over [1/2 .. 2863311531/4294967296]>'
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4066)
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))
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4212)
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)
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 4261)
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]
|
File: sage/rings/polynomial/real_roots.pyx (starting at line 1769)
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
|
|
|||
AA
|
ZZ
|
lmap
|
| Home | Trees | Indices | Help |
|---|
| Generated by Epydoc 3.0beta1 on Thu Jul 17 04:23:28 2008 | http://epydoc.sourceforge.net |