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 .

exception sage.rings.polynomial.real_roots.PrecisionError

Bases: ValueError

sage.rings.polynomial.real_roots.bernstein_down(d1, d2, s)

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

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]

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bernstein_down(Integer(3), Integer(8), Integer(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]

sage.rings.polynomial.real_roots.bernstein_expand(c, d2)

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).

This 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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> c = vector(ZZ, [Integer(1000), Integer(2000), -Integer(3000)])
>>> bernstein_expand(c, Integer(3))
((1000, 1666, 333, -3000), 1)
>>> bernstein_expand(c, Integer(4))
((1000, 1500, 1000, -500, -3000), 1)
>>> bernstein_expand(c, Integer(20))
((1000, 1100, 1168, 1205, 1210, 1184, 1126, 1036, 915, 763, 578, 363, 115, -164, -474, -816, -1190, -1595, -2032, -2500, -3000), 1)

class sage.rings.polynomial.real_roots.bernstein_polynomial_factory

Bases: object

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.

lsign()

Return the sign of the first coefficient of this Bernstein polynomial.

usign()

Return the sign of the last coefficient of this Bernstein polynomial.

class sage.rings.polynomial.real_roots.bernstein_polynomial_factory_ar(poly, neg)

Bases: bernstein_polynomial_factory

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(scale_log2)

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

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>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(AA)
>>> p = (x - Integer(1)) * (x - sqrt(AA(Integer(2)))) * (x - Integer(2))
>>> bpf = bernstein_polynomial_factory_ar(p, False)
>>> print(bpf.bernstein_polynomial(Integer(0)))
degree 3 IBP with 2-bit coefficients
>>> bpf.bernstein_polynomial(-Integer(20))
<IBP: ((-2965821, 2181961, -1542880, 1048576) + [0 .. 1)) * 2^-20>
>>> bpf = bernstein_polynomial_factory_ar(p, True)
>>> bpf.bernstein_polynomial(-Integer(20))
<IBP: ((-2965821, -2181962, -1542880, -1048576) + [0 .. 1)) * 2^-20>
>>> p = x**Integer(2) - Integer(1)
>>> bpf = bernstein_polynomial_factory_ar(p, False)
>>> bpf.bernstein_polynomial(-Integer(10))
<IBP: ((-1024, 0, 1024) + [0 .. 1)) * 2^-10>

coeffs_bitsize()

Compute the approximate log2 of the maximum of the absolute values of the coefficients.

EXAMPLES:

Sage

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(AA)
>>> p = (x - Integer(1)) * (x - sqrt(AA(Integer(2)))) * (x - Integer(2))
>>> bernstein_polynomial_factory_ar(p, False).coeffs_bitsize()
1

class sage.rings.polynomial.real_roots.bernstein_polynomial_factory_intlist(coeffs)

Bases: bernstein_polynomial_factory

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(scale_log2)

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

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>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bpf = bernstein_polynomial_factory_intlist([Integer(10), -Integer(20), Integer(30), -Integer(40)])
>>> print(bpf.bernstein_polynomial(Integer(0)))
degree 3 IBP with 6-bit coefficients
>>> bpf.bernstein_polynomial(Integer(20))
<IBP: ((0, -1, 0, -1) + [0 .. 1)) * 2^20; lsign 1>
>>> bpf.bernstein_polynomial(Integer(0))
<IBP: (10, -20, 30, -40) + [0 .. 1)>
>>> bpf.bernstein_polynomial(-Integer(20))
<IBP: ((10485760, -20971520, 31457280, -41943040) + [0 .. 1)) * 2^-20>

coeffs_bitsize()

Compute the approximate log2 of the maximum of the absolute values of the coefficients.

EXAMPLES:

Sage

sage: from sage.rings.polynomial.real_roots import *
sage: bernstein_polynomial_factory_intlist([1, 2, 3, -60000]).coeffs_bitsize()
16

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bernstein_polynomial_factory_intlist([Integer(1), Integer(2), Integer(3), -Integer(60000)]).coeffs_bitsize()
16

class sage.rings.polynomial.real_roots.bernstein_polynomial_factory_ratlist(coeffs)

Bases: bernstein_polynomial_factory

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.

bernstein_polynomial(scale_log2)

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

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>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bpf = bernstein_polynomial_factory_ratlist([Integer(1)/Integer(3), -Integer(22)/Integer(7), Integer(193)/Integer(71), -Integer(140)/Integer(99)])
>>> print(bpf.bernstein_polynomial(Integer(0)))
degree 3 IBP with 3-bit coefficients
>>> bpf.bernstein_polynomial(Integer(20))
<IBP: ((0, -1, 0, -1) + [0 .. 1)) * 2^20; lsign 1>
>>> bpf.bernstein_polynomial(Integer(0))
<IBP: (0, -4, 2, -2) + [0 .. 1); lsign 1>
>>> bpf.bernstein_polynomial(-Integer(20))
<IBP: ((349525, -3295525, 2850354, -1482835) + [0 .. 1)) * 2^-20>

coeffs_bitsize()

Compute the approximate log2 of the maximum of the absolute values of the coefficients.

EXAMPLES:

Sage

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bernstein_polynomial_factory_ratlist([Integer(1), Integer(2), Integer(3), -Integer(60000)]).coeffs_bitsize()
15
>>> bernstein_polynomial_factory_ratlist([Integer(65535)/Integer(65536)]).coeffs_bitsize()
-1
>>> bernstein_polynomial_factory_ratlist([Integer(65536)/Integer(65535)]).coeffs_bitsize()
1

sage.rings.polynomial.real_roots.bernstein_up(d1, d2, s=None)

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

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]

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bernstein_down(Integer(3), Integer(7), Integer(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]

sage.rings.polynomial.real_roots.bitsize_doctest(n)

sage.rings.polynomial.real_roots.cl_maximum_root(cl)

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

sage: from sage.rings.polynomial.real_roots import *
sage: cl_maximum_root([RIF(-1), RIF(0), RIF(1)])
1.00000000000000

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> cl_maximum_root([RIF(-Integer(1)), RIF(Integer(0)), RIF(Integer(1))])
1.00000000000000

sage.rings.polynomial.real_roots.cl_maximum_root_first_lambda(cl)

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

sage: from sage.rings.polynomial.real_roots import *
sage: cl_maximum_root_first_lambda([RIF(-1), RIF(0), RIF(1)])
1.00000000000000

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> cl_maximum_root_first_lambda([RIF(-Integer(1)), RIF(Integer(0)), RIF(Integer(1))])
1.00000000000000

sage.rings.polynomial.real_roots.cl_maximum_root_local_max(cl)

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

sage: from sage.rings.polynomial.real_roots import *
sage: cl_maximum_root_local_max([RIF(-1), RIF(0), RIF(1)])
1.41421356237310

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> cl_maximum_root_local_max([RIF(-Integer(1)), RIF(Integer(0)), RIF(Integer(1))])
1.41421356237310

class sage.rings.polynomial.real_roots.context

Bases: object

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.

get_be_log()

get_dc_log()

sage.rings.polynomial.real_roots.de_casteljau_doublevec(c, x)

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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> c = vector(RDF, [RealNumber('0.7'), Integer(0), Integer(0), Integer(0), Integer(0), Integer(0)])
>>> de_casteljau_doublevec(c, Integer(1)/Integer(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)
>>> de_casteljau_doublevec(c, Integer(1)/Integer(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)
>>> de_casteljau_doublevec(c, Integer(7)/Integer(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)

sage.rings.polynomial.real_roots.de_casteljau_intvec(c, c_bitsize, x, use_ints)

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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> c = vector(ZZ, [Integer(1048576), Integer(0), Integer(0), Integer(0), Integer(0), Integer(0)])
>>> de_casteljau_intvec(c, Integer(20), Integer(1)/Integer(2), Integer(1))
((1048576, 524288, 262144, 131072, 65536, 32768), (32768, 0, 0, 0, 0, 0), 1)
>>> de_casteljau_intvec(c, Integer(20), Integer(1)/Integer(3), Integer(1))
((1048576, 699050, 466033, 310689, 207126, 138084), (138084, 0, 0, 0, 0, 0), 1)
>>> de_casteljau_intvec(c, Integer(20), Integer(7)/Integer(22), Integer(1))
((1048576, 714938, 487457, 332357, 226607, 154505), (154505, 0, 0, 0, 0, 0), 1)

sage.rings.polynomial.real_roots.degree_reduction_next_size(n)

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

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> degree_reduction_next_size(Integer(1000))
30
>>> degree_reduction_next_size(Integer(20))
15
>>> degree_reduction_next_size(Integer(3))
2
>>> degree_reduction_next_size(Integer(2)) is None
True

sage.rings.polynomial.real_roots.dprod_imatrow_vec(m, v, k)

Compute 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

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> m = matrix(Integer(3), range(Integer(9)))
>>> dprod_imatrow_vec(m, vector(ZZ, [Integer(1), Integer(0), Integer(0), Integer(0)]), Integer(1))
0
>>> dprod_imatrow_vec(m, vector(ZZ, [Integer(0), Integer(1), Integer(0), Integer(0)]), Integer(1))
3
>>> dprod_imatrow_vec(m, vector(ZZ, [Integer(0), Integer(0), Integer(1), Integer(0)]), Integer(1))
4
>>> dprod_imatrow_vec(m, vector(ZZ, [Integer(0), Integer(0), Integer(0), Integer(1)]), Integer(1))
5
>>> dprod_imatrow_vec(m, vector(ZZ, [Integer(1), Integer(0), Integer(0)]), Integer(1))
3
>>> dprod_imatrow_vec(m, vector(ZZ, [Integer(0), Integer(1), Integer(0)]), Integer(1))
4
>>> dprod_imatrow_vec(m, vector(ZZ, [Integer(0), Integer(0), Integer(1)]), Integer(1))
5
>>> dprod_imatrow_vec(m, vector(ZZ, [Integer(1), Integer(2), Integer(3)]), Integer(1))
26

sage.rings.polynomial.real_roots.get_realfield_rndu(n)

A simple cache for RealField fields (with rounding set to round-to-positive-infinity).

EXAMPLES:

Sage

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> get_realfield_rndu(Integer(20))
Real Field with 20 bits of precision and rounding RNDU
>>> get_realfield_rndu(Integer(53))
Real Field with 53 bits of precision and rounding RNDU
>>> get_realfield_rndu(Integer(20))
Real Field with 20 bits of precision and rounding RNDU

class sage.rings.polynomial.real_roots.interval_bernstein_polynomial

Bases: object

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

\[\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.

region()

region_width()

try_rand_split(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

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.10000000000000069 .. 0.010000000000000677] over [0 .. 43/64]>
sage: bp2  # rel tol
<IBP: (-0.3145806974172592, -0.19903896331787108, 0.04135986328125002, 0.43546875, 0.99) + [-0.10000000000000069 .. 0.010000000000000677] over [43/64 .. 1]>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(50), Integer(20), -Integer(90), -Integer(70), Integer(200)], error=Integer(5))
>>> bp1, bp2, _ = bp.try_rand_split(mk_context(), None)
>>> bp1
<IBP: (50, 29, -27, -56, -11) + [0 .. 6) over [0 .. 43/64]>
>>> bp2
<IBP: (-11, 10, 49, 111, 200) + [0 .. 6) over [43/64 .. 1]>
>>> bp1, bp2, _ = bp.try_rand_split(mk_context(seed=Integer(42)), None)
>>> bp1
<IBP: (50, 32, -11, -41, -29) + [0 .. 6) over [0 .. 583/1024]>
>>> bp2
<IBP: (-29, -20, 13, 83, 200) + [0 .. 6) over [583/1024 .. 1]>
>>> bp = mk_ibpf([RealNumber('0.5'), RealNumber('0.2'), -RealNumber('0.9'), -RealNumber('0.7'), RealNumber('0.99')], neg_err=-RealNumber('0.1'), pos_err=RealNumber('0.01'))
>>> bp1, bp2, _ = bp.try_rand_split(mk_context(), None)
>>> bp1  # rel tol
<IBP: (0.5, 0.2984375, -0.2642578125, -0.5511661529541015, -0.3145806974172592) + [-0.10000000000000069 .. 0.010000000000000677] over [0 .. 43/64]>
>>> bp2  # rel tol
<IBP: (-0.3145806974172592, -0.19903896331787108, 0.04135986328125002, 0.43546875, 0.99) + [-0.10000000000000069 .. 0.010000000000000677] over [43/64 .. 1]>

try_split(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

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.10000000000000023 .. 0.010000000000000226] over [0 .. 1/2]>
sage: bp2
<IBP: (-0.369375, -0.45125, -0.3275, 0.14500000000000002, 0.99) + [-0.10000000000000023 .. 0.010000000000000226] over [1/2 .. 1]>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(50), Integer(20), -Integer(90), -Integer(70), Integer(200)], error=Integer(5))
>>> bp1, bp2, _ = bp.try_split(mk_context(), None)
>>> bp1
<IBP: (50, 35, 0, -29, -31) + [0 .. 6) over [0 .. 1/2]>
>>> bp2
<IBP: (-31, -33, -8, 65, 200) + [0 .. 6) over [1/2 .. 1]>
>>> bp = mk_ibpf([RealNumber('0.5'), RealNumber('0.2'), -RealNumber('0.9'), -RealNumber('0.7'), RealNumber('0.99')], neg_err=-RealNumber('0.1'), pos_err=RealNumber('0.01'))
>>> bp1, bp2, _ = bp.try_split(mk_context(), None)
>>> bp1
<IBP: (0.5, 0.35, 0.0, -0.2875, -0.369375) + [-0.10000000000000023 .. 0.010000000000000226] over [0 .. 1/2]>
>>> bp2
<IBP: (-0.369375, -0.45125, -0.3275, 0.14500000000000002, 0.99) + [-0.10000000000000023 .. 0.010000000000000226] over [1/2 .. 1]>

variations()

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.

class sage.rings.polynomial.real_roots.interval_bernstein_polynomial_float

Bases: 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

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],                         # needs sage.symbolic
....:              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                                                                        # needs sage.symbolic
<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()                                                           # needs sage.symbolic
(3, 3)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpf([RealNumber('0.1'), RealNumber('0.2'), RealNumber('0.3')], pos_err=RealNumber('0.5')); print(bp)
degree 2 IBP with floating-point coefficients
>>> bp
<IBP: (0.1, 0.2, 0.3) + [0.0 .. 0.5]>
>>> bp.variations()
(0, 0)
>>> bp = mk_ibpf([-RealNumber('0.3'), -RealNumber('0.1'), RealNumber('0.1'), -RealNumber('0.1'), -RealNumber('0.3'), -RealNumber('0.1')],                         # needs sage.symbolic
...              lower=Integer(1), upper=Integer(5)/Integer(4), usign=Integer(1), pos_err=RealNumber('0.2'),
...              scale_log2=-Integer(3), level=Integer(2), slope_err=RIF(pi)); print(bp)
degree 5 IBP with floating-point coefficients
>>> bp                                                                        # needs sage.symbolic
<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?>
>>> bp.variations()                                                           # needs sage.symbolic
(3, 3)

as_float()

de_casteljau(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 – boolean; True if the sign of the original polynomial at mid is known

EXAMPLES:

Sage

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.10000000000000023 .. 0.010000000000000226] over [0 .. 1/2]>
sage: bp2
<IBP: (-0.369375, -0.45125, -0.3275, 0.14500000000000002, 0.99) + [-0.10000000000000023 .. 0.010000000000000226] 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.10000000000000069 .. 0.010000000000000677] over [0 .. 2/3]>
sage: bp2  # rel tol 3e-15
<IBP: (-0.32172839506172846, -0.21037037037037046, 0.028888888888888797, 0.4266666666666666, 0.99) + [-0.10000000000000069 .. 0.010000000000000677] 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.10000000000000069 .. 0.010000000000000677] over [0 .. 7/39]>
sage: bp2  # rel tol
<IBP: (0.1765692706232836, -0.26556803047927313, -0.7802038132807364, -0.3966666666666666, 0.99) + [-0.10000000000000069 .. 0.010000000000000677] over [7/39 .. 1]>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> ctx = mk_context()
>>> bp = mk_ibpf([RealNumber('0.5'), RealNumber('0.2'), -RealNumber('0.9'), -RealNumber('0.7'), RealNumber('0.99')], neg_err=-RealNumber('0.1'), pos_err=RealNumber('0.01'))
>>> bp1, bp2, ok = bp.de_casteljau(ctx, Integer(1)/Integer(2))
>>> bp1
<IBP: (0.5, 0.35, 0.0, -0.2875, -0.369375) + [-0.10000000000000023 .. 0.010000000000000226] over [0 .. 1/2]>
>>> bp2
<IBP: (-0.369375, -0.45125, -0.3275, 0.14500000000000002, 0.99) + [-0.10000000000000023 .. 0.010000000000000226] over [1/2 .. 1]>
>>> bp1, bp2, ok = bp.de_casteljau(ctx, Integer(2)/Integer(3))
>>> bp1 # rel tol 2e-16
<IBP: (0.5, 0.30000000000000004, -0.2555555555555555, -0.5444444444444444, -0.32172839506172846) + [-0.10000000000000069 .. 0.010000000000000677] over [0 .. 2/3]>
>>> bp2  # rel tol 3e-15
<IBP: (-0.32172839506172846, -0.21037037037037046, 0.028888888888888797, 0.4266666666666666, 0.99) + [-0.10000000000000069 .. 0.010000000000000677] over [2/3 .. 1]>
>>> bp1, bp2, ok = bp.de_casteljau(ctx, Integer(7)/Integer(39))
>>> bp1  # rel tol
<IBP: (0.5, 0.4461538461538461, 0.36653517422748183, 0.27328680523946786, 0.1765692706232836) + [-0.10000000000000069 .. 0.010000000000000677] over [0 .. 7/39]>
>>> bp2  # rel tol
<IBP: (0.1765692706232836, -0.26556803047927313, -0.7802038132807364, -0.3966666666666666, 0.99) + [-0.10000000000000069 .. 0.010000000000000677] over [7/39 .. 1]>

get_msb_bit()

Return an approximation of the log2 of the maximum of the absolute values of the coefficients, as an integer.

slope_range()

Compute a bound on the derivative of this polynomial, over its region.

EXAMPLES:

Sage

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]'

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpf([RealNumber('0.5'), RealNumber('0.2'), -RealNumber('0.9'), -RealNumber('0.7'), RealNumber('0.99')], neg_err=-RealNumber('0.1'), pos_err=RealNumber('0.01'))
>>> bp.slope_range().str(style='brackets')
'[-4.8400000000000017 .. 7.2000000000000011]'

class sage.rings.polynomial.real_roots.interval_bernstein_polynomial_integer

Bases: 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

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,        # needs sage.symbolic
....:              error=2, scale_log2=-3, level=2, slope_err=RIF(pi)); print(bp)
degree 5 IBP with 2-bit coefficients
sage: bp                                                                        # needs sage.symbolic
<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()                                                           # needs sage.symbolic
(3, 3)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(1), Integer(2), Integer(3)], error=Integer(5)); print(bp)
degree 2 IBP with 2-bit coefficients
>>> bp
<IBP: (1, 2, 3) + [0 .. 5)>
>>> bp.variations()
(0, 0)
>>> bp = mk_ibpi([-Integer(3), -Integer(1), Integer(1), -Integer(1), -Integer(3), -Integer(1)], lower=Integer(1), upper=Integer(5)/Integer(4), usign=Integer(1),        # needs sage.symbolic
...              error=Integer(2), scale_log2=-Integer(3), level=Integer(2), slope_err=RIF(pi)); print(bp)
degree 5 IBP with 2-bit coefficients
>>> bp                                                                        # needs sage.symbolic
<IBP: ((-3, -1, 1, -1, -3, -1) + [0 .. 2)) * 2^-3 over [1 .. 5/4]; usign 1;
      level 2; slope_err 3.141592653589794?>
>>> bp.variations()                                                           # needs sage.symbolic
(3, 3)

as_float()

Compute an interval_bernstein_polynomial_float which contains (or bounds) all the polynomials this interval polynomial contains (or bounds).

EXAMPLES:

Sage

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.1275702593849246e-16 .. 0.01953125000000017]) * 2^8>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(50), Integer(20), -Integer(90), -Integer(70), Integer(200)], error=Integer(5))
>>> print(bp.as_float())
degree 4 IBP with floating-point coefficients
>>> bp.as_float()
<IBP: ((0.1953125, 0.078125, -0.3515625, -0.2734375, 0.78125) + [-1.1275702593849246e-16 .. 0.01953125000000017]) * 2^8>

de_casteljau(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 – boolean; True if the sign of the original polynomial at mid is known

EXAMPLES:

Sage

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]>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(50), Integer(20), -Integer(90), -Integer(70), Integer(200)], error=Integer(5))
>>> ctx = mk_context()
>>> bp1, bp2, ok = bp.de_casteljau(ctx, Integer(1)/Integer(2))
>>> bp1
<IBP: (50, 35, 0, -29, -31) + [0 .. 6) over [0 .. 1/2]>
>>> bp2
<IBP: (-31, -33, -8, 65, 200) + [0 .. 6) over [1/2 .. 1]>
>>> bp1, bp2, ok = bp.de_casteljau(ctx, Integer(2)/Integer(3))
>>> bp1
<IBP: (50, 30, -26, -55, -13) + [0 .. 6) over [0 .. 2/3]>
>>> bp2
<IBP: (-13, 8, 47, 110, 200) + [0 .. 6) over [2/3 .. 1]>
>>> bp1, bp2, ok = bp.de_casteljau(ctx, Integer(7)/Integer(39))
>>> bp1
<IBP: (50, 44, 36, 27, 17) + [0 .. 6) over [0 .. 7/39]>
>>> bp2
<IBP: (17, -26, -75, -22, 200) + [0 .. 6) over [7/39 .. 1]>

down_degree(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

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>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(0), Integer(100), Integer(400), Integer(903)], error=Integer(2))
>>> ctx = mk_context()
>>> bp
<IBP: (0, 100, 400, 903) + [0 .. 2)>
>>> dbp, _ = bp.down_degree(ctx, Integer(10), Integer(32))
>>> dbp
<IBP: (-1, 148, 901) + [0 .. 4); level 1; slope_err 0.?e2>

down_degree_iter(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

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>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(0), Integer(100), Integer(400), Integer(903), Integer(1600), Integer(2500)], error=Integer(2))
>>> ctx = mk_context()
>>> bp
<IBP: (0, 100, 400, 903, 1600, 2500) + [0 .. 2)>
>>> rbp = bp.down_degree_iter(ctx, Integer(6))
>>> rbp
<IBP: (-4, 249, 2497) + [0 .. 9); level 2; slope_err 0.?e3>

downscale(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

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>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(0), Integer(100), Integer(400), Integer(903)], error=Integer(2))
>>> bp.downscale(Integer(5))
<IBP: ((0, 3, 12, 28) + [0 .. 1)) * 2^5>

get_msb_bit()

Return an approximation of the log2 of the maximum of the absolute values of the coefficients, as an integer.

slope_range()

Compute a bound on the derivative of this polynomial, over its region.

EXAMPLES:

Sage

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]'

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(0), Integer(100), Integer(400), Integer(903)], error=Integer(2))
>>> bp.slope_range().str(style='brackets')
'[294.00000000000000 .. 1515.0000000000000]'

sage.rings.polynomial.real_roots.intvec_to_doublevec(b, err)

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 \le ak\]

and

\[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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> intvec_to_doublevec(vector(ZZ, [Integer(1), Integer(2), Integer(3), Integer(4), Integer(5)]), Integer(3))
((0.125, 0.25, 0.375, 0.5, 0.625), -1.1275702593849246e-16, 0.37500000000000017, 3)

class sage.rings.polynomial.real_roots.island

Bases: object

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.)

bp_done(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.)

done(ctx)

Check to see if the island is known to contain zero roots or is known to contain one root.

has_root()

Assuming that the island is done (has either 0 or 1 roots), reports whether the island has a root.

less_bits(ancestors, bp)

Heuristically push lower-precision polynomials on the polynomial stack. See the class documentation for class island for more information.

more_bits(ctx, ancestors, 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.)

refine(ctx)

Attempts to shrink and/or split this island into sub-island that each definitely contain exactly one root.

refine_recurse(ctx, 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.

reset_root_width(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.

shrink_bp(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.

class sage.rings.polynomial.real_roots.linear_map(lower, upper)

Bases: object

A simple class to map linearly between original coordinates (ranging from [lower .. upper]) and ocean coordinates (ranging from [0 .. 1]).

from_ocean(region)

to_ocean(region)

sage.rings.polynomial.real_roots.max_abs_doublevec(c)

Given a floating-point vector, return the maximum of the absolute values of its elements.

EXAMPLES:

Sage

sage: from sage.rings.polynomial.real_roots import *
sage: max_abs_doublevec(vector(RDF, [0.1, -0.767, 0.3, 0.693]))
0.767

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> max_abs_doublevec(vector(RDF, [RealNumber('0.1'), -RealNumber('0.767'), RealNumber('0.3'), RealNumber('0.693')]))
0.767

sage.rings.polynomial.real_roots.max_bitsize_intvec_doctest(b)

sage.rings.polynomial.real_roots.maximum_root_first_lambda(p)

Given a polynomial with real coefficients, computes an upper bound on its largest real root.

This is 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

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(ZZ)
>>> maximum_root_first_lambda((x-Integer(1))*(x-Integer(2))*(x-Integer(3)))
6.00000000000001
>>> maximum_root_first_lambda((x+Integer(1))*(x+Integer(2))*(x+Integer(3)))
0.000000000000000
>>> maximum_root_first_lambda(x**Integer(2) - Integer(1))
1.00000000000000

sage.rings.polynomial.real_roots.maximum_root_local_max(p)

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

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(ZZ)
>>> maximum_root_local_max((x-Integer(1))*(x-Integer(2))*(x-Integer(3)))
12.0000000000001
>>> maximum_root_local_max((x+Integer(1))*(x+Integer(2))*(x+Integer(3)))
0.000000000000000
>>> maximum_root_local_max(x**Integer(2) - Integer(1))
1.41421356237310

sage.rings.polynomial.real_roots.min_max_delta_intvec(a, b)

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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> a = vector(ZZ, [Integer(10), -Integer(30)])
>>> b = vector(ZZ, [Integer(15), -Integer(60)])
>>> min_max_delta_intvec(a, b)
(30, -5)

sage.rings.polynomial.real_roots.min_max_diff_doublevec(c)

Given a floating-point vector b = (b0, …, bn), compute the minimum and maximum values of b_{j+1} - b_j.

EXAMPLES:

Sage

sage: from sage.rings.polynomial.real_roots import *
sage: min_max_diff_doublevec(vector(RDF, [1, 7, -2]))
(-9.0, 6.0)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> min_max_diff_doublevec(vector(RDF, [Integer(1), Integer(7), -Integer(2)]))
(-9.0, 6.0)

sage.rings.polynomial.real_roots.min_max_diff_intvec(b)

Given an integer vector b = (b0, …, bn), compute the minimum and maximum values of b_{j+1} - b_j.

EXAMPLES:

Sage

sage: from sage.rings.polynomial.real_roots import *
sage: min_max_diff_intvec(vector(ZZ, [1, 7, -2]))
(-9, 6)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> min_max_diff_intvec(vector(ZZ, [Integer(1), Integer(7), -Integer(2)]))
(-9, 6)

sage.rings.polynomial.real_roots.mk_context(do_logging=False, seed=0, wordsize=32)

A simple wrapper for creating context objects with coercions, defaults, etc.

For use in doctests.

EXAMPLES:

Sage

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> mk_context(do_logging=True, seed=Integer(3), wordsize=Integer(64))
root isolation context: seed=3; do_logging=True; wordsize=64

sage.rings.polynomial.real_roots.mk_ibpf(coeffs, lower=0, upper=1, lsign=0, usign=0, neg_err=0, pos_err=0, scale_log2=0, level=0, slope_err=None)

A simple wrapper for creating interval_bernstein_polynomial_float objects with coercions, defaults, etc.

For use in doctests.

EXAMPLES:

Sage

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> print(mk_ibpf([RealNumber('0.5'), RealNumber('0.2'), -RealNumber('0.9'), -RealNumber('0.7'), RealNumber('0.99')], pos_err=RealNumber('0.1'), neg_err=-RealNumber('0.01')))
degree 4 IBP with floating-point coefficients

sage.rings.polynomial.real_roots.mk_ibpi(coeffs, lower=0, upper=1, lsign=0, usign=0, error=1, scale_log2=0, level=0, slope_err=None)

A simple wrapper for creating interval_bernstein_polynomial_integer objects with coercions, defaults, etc.

For use in doctests.

EXAMPLES:

Sage

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> print(mk_ibpi([Integer(50), Integer(20), -Integer(90), -Integer(70), Integer(200)], error=Integer(5)))
degree 4 IBP with 8-bit coefficients

class sage.rings.polynomial.real_roots.ocean

Bases: object

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.

all_done()

Return True iff all islands are known to contain exactly one root.

EXAMPLES:

Sage

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([Integer(1)/Integer(3), -Integer(22)/Integer(7), Integer(193)/Integer(71), -Integer(140)/Integer(99)]), lmap)
>>> oc.all_done()
False
>>> oc.find_roots()
>>> oc.all_done()
True

approx_bp(scale_log2)

Return an approximation to our Bernstein polynomial with the given scale_log2.

EXAMPLES:

Sage

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>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([Integer(1)/Integer(3), -Integer(22)/Integer(7), Integer(193)/Integer(71), -Integer(140)/Integer(99)]), lmap)
>>> oc.approx_bp(Integer(0))
<IBP: (0, -4, 2, -2) + [0 .. 1); lsign 1>
>>> oc.approx_bp(-Integer(20))
<IBP: ((349525, -3295525, 2850354, -1482835) + [0 .. 1)) * 2^-20>

find_roots()

Isolate all roots in this ocean.

EXAMPLES:

Sage

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([Integer(1)/Integer(3), -Integer(22)/Integer(7), Integer(193)/Integer(71), -Integer(140)/Integer(99)]), lmap)
>>> oc
ocean with precision 120 and 1 island(s)
>>> oc.find_roots()
>>> oc
ocean with precision 120 and 3 island(s)
>>> oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([Integer(1), Integer(0), -Integer(1111)/Integer(2), Integer(0), Integer(11108889)/Integer(14), Integer(0), Integer(0), Integer(0), Integer(0), -Integer(1)]), lmap)
>>> oc.find_roots()
>>> oc
ocean with precision 240 and 3 island(s)

increase_precision()

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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([Integer(1)/Integer(3), -Integer(22)/Integer(7), Integer(193)/Integer(71), -Integer(140)/Integer(99)]), lmap)
>>> oc
ocean with precision 120 and 1 island(s)
>>> oc.increase_precision()
>>> oc.increase_precision()
>>> oc.increase_precision()
>>> oc
ocean with precision 960 and 1 island(s)

refine_all()

Refine all islands which are not done (which are not known to contain exactly one root).

EXAMPLES:

Sage

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([Integer(1)/Integer(3), -Integer(22)/Integer(7), Integer(193)/Integer(71), -Integer(140)/Integer(99)]), lmap)
>>> oc
ocean with precision 120 and 1 island(s)
>>> oc.refine_all()
>>> oc
ocean with precision 120 and 3 island(s)

reset_root_width(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

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([-Integer(1), -Integer(1), Integer(1)]), lmap)
>>> oc.find_roots()
>>> oc.roots()
[(1/2, 3/4)]
>>> oc.reset_root_width(Integer(0), Integer(1)/Integer(2)**Integer(200))
>>> oc.find_roots()
>>> oc
ocean with precision 240 and 1 island(s)
>>> RR(RealIntervalField(Integer(300))(oc.roots()[Integer(0)]).absolute_diameter()).log2()
-232.668979560890

roots()

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

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)]

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([Integer(1)/Integer(3), -Integer(22)/Integer(7), Integer(193)/Integer(71), -Integer(140)/Integer(99)]), lmap)
>>> oc.find_roots()
>>> oc.roots()
[(1/32, 1/16), (1/2, 5/8), (3/4, 7/8)]
>>> oc = ocean(mk_context(), bernstein_polynomial_factory_ratlist([Integer(1), Integer(0), -Integer(1111)/Integer(2), Integer(0), Integer(11108889)/Integer(14), Integer(0), Integer(0), Integer(0), Integer(0), -Integer(1)]), lmap)
>>> oc.find_roots()
>>> oc.roots()
[(95761241267509487747625/9671406556917033397649408, 191522482605387719863145/19342813113834066795298816), (1496269395904347376805/151115727451828646838272, 374067366568272936175/37778931862957161709568), (31/32, 63/64)]

sage.rings.polynomial.real_roots.precompute_degree_reduction_cache(n)

Compute and cache the matrices used for degree reduction, starting from degree n.

EXAMPLES:

Sage

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,

([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)
)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> precompute_degree_reduction_cache(Integer(5))
>>> dr_cache[Integer(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)
)

sage.rings.polynomial.real_roots.pseudoinverse(m)

sage.rings.polynomial.real_roots.rational_root_bounds(p)

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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(ZZ)
>>> rational_root_bounds((x-Integer(1))*(x-Integer(2))*(x-Integer(3)))
(0, 7)
>>> rational_root_bounds(x**Integer(2))
(-1/2, 1/2)
>>> rational_root_bounds(x*(x+Integer(1)))
(-3/2, 1/2)
>>> rational_root_bounds((x+Integer(2))*(x-Integer(3)))
(-3, 6)
>>> rational_root_bounds(x**Integer(995) * (x**Integer(2) - Integer(9999)) - Integer(1))
(-100, 1000/7)
>>> rational_root_bounds(x**Integer(995) * (x**Integer(2) - Integer(9999)) + Integer(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

sage.rings.polynomial.real_roots.real_roots(p, bounds=None, seed=None, skip_squarefree=False, do_logging=False, wordsize=32, retval='rational', strategy=None, max_diameter=None)

Compute the real roots of a given polynomial with exact coefficients (integer, rational, and algebraic real coefficients are supported).

This 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

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(ZZ)
>>> real_roots(x**Integer(3) - x**Integer(2) - x - Integer(1))
[((7/4, 19/8), 1)]
>>> real_roots((x-Integer(1))*(x-Integer(2))*(x-Integer(3))*(x-Integer(5))*(x-Integer(8))*(x-Integer(13))*(x-Integer(21))*(x-Integer(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)]
>>> real_roots(x**Integer(5) * (x**Integer(2) - Integer(9999))**Integer(2) - Integer(1))
[((-29274496381311/9007199254740992, 419601125186091/2251799813685248), 1), ((2126658450145849453951061654415153249597/21267647932558653966460912964485513216, 4253316902721330018853696359533061621799/42535295865117307932921825928971026432), 1), ((1063329226287740282451317352558954186101/10633823966279326983230456482242756608, 531664614358685696701445201630854654353/5316911983139663491615228241121378304), 1)]
>>> real_roots(x**Integer(5) * (x**Integer(2) - Integer(9999))**Integer(2) - Integer(1), seed=Integer(42))
[((-123196838480289/18014398509481984, 293964743458749/9007199254740992), 1), ((8307259573979551907841696381986376143/83076749736557242056487941267521536, 16614519150981033789137940378745325503/166153499473114484112975882535043072), 1), ((519203723562592617581015249797434335/5192296858534827628530496329220096, 60443268924081068060312183/604462909807314587353088), 1)]
>>> real_roots(x**Integer(5) * (x**Integer(2) - Integer(9999))**Integer(2) - Integer(1), wordsize=Integer(64))
[((-62866503803202151050003/19342813113834066795298816, 901086554512564177624143/4835703278458516698824704), 1), ((544424563237337315214990987922809050101157/5444517870735015415413993718908291383296, 1088849127096660194637118845654929064385439/10889035741470030830827987437816582766592), 1), ((272212281929661439711063928866060007142141/2722258935367507707706996859454145691648, 136106141275823501959100399337685485662633/1361129467683753853853498429727072845824), 1)]
>>> real_roots(x)
[((-47/256, 81/512), 1)]
>>> real_roots(x * (x-Integer(1)))
[((-47/256, 81/512), 1), ((1/2, 1201/1024), 1)]
>>> real_roots(x-Integer(1))
[((209/256, 593/512), 1)]
>>> real_roots(x*(x-Integer(1))*(x-Integer(2)), bounds=(Integer(0), Integer(2)))
[((0, 0), 1), ((81/128, 337/256), 1), ((2, 2), 1)]
>>> real_roots(x*(x-Integer(1))*(x-Integer(2)), bounds=(Integer(0), Integer(2)), retval='algebraic_real')
[(0, 1), (1, 1), (2, 1)]
>>> v = Integer(2)**Integer(40)
>>> real_roots((x**Integer(2)-Integer(1))**Integer(2) * (x**Integer(2) - (v+Integer(1))/v))
[((-12855504354077768210885019021174120740504020581912910106032833/12855504354071922204335696738729300820177623950262342682411008, -6427752177038884105442509510587059395588605840418680645585479/6427752177035961102167848369364650410088811975131171341205504), 1), ((-1125899906842725/1125899906842624, -562949953421275/562949953421312), 2), ((62165404551223330269422781018352603934643403586760330761772204409982940218804935733653/62165404551223330269422781018352605012557018849668464680057997111644937126566671941632, 3885337784451458141838923813647037871787041539340705594199885610069035709862106085785/3885337784451458141838923813647037813284813678104279042503624819477808570410416996352), 2), ((509258994083853105745586001837045839749063767798922046787130823804169826426726965449697819/509258994083621521567111422102344540262867098416484062659035112338595324940834176545849344, 25711008708155536421770038042348240136257704305733983563630791/25711008708143844408671393477458601640355247900524685364822016), 1)]
>>> real_roots(x**Integer(2) - Integer(2))
[((-3/2, -1), 1), ((1, 3/2), 1)]
>>> real_roots(x**Integer(2) - Integer(2), retval='interval')
[(-2.?, 1), (2.?, 1)]
>>> real_roots(x**Integer(2) - Integer(2), max_diameter=Integer(1)/Integer(2)**Integer(30))
[((-22506280506048041472675379598886543645348790970912519198456805737131269246430553365310109/15914343565113172548972231940698266883214596825515126958094847260581103904401068017057792, -45012561012096082945350759197773086524448972309421182031053065599548946985601579935498343/31828687130226345097944463881396533766429193651030253916189694521162207808802136034115584), 1), ((45012561012096082945350759197773086524448972309421182031053065599548946985601579935498343/31828687130226345097944463881396533766429193651030253916189694521162207808802136034115584, 22506280506048041472675379598886543645348790970912519198456805737131269246430553365310109/15914343565113172548972231940698266883214596825515126958094847260581103904401068017057792), 1)]
>>> real_roots(x**Integer(2) - Integer(2), retval='interval', max_diameter=Integer(1)/Integer(2)**Integer(500))
[(-1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571470109559971605970274534596862014728517418640889198609552?, 1), (1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571470109559971605970274534596862014728517418640889198609552?, 1)]
>>> ar_rts = real_roots(x**Integer(2) - Integer(2), retval='algebraic_real'); ar_rts
[(-1.414213562373095?, 1), (1.414213562373095?, 1)]
>>> ar_rts[Integer(0)][Integer(0)]**Integer(2) - Integer(2) == Integer(0)
True
>>> v = Integer(2)**Integer(40)
>>> real_roots((x-Integer(1)) * (x-(v+Integer(1))/v), retval='interval')
[(1.000000000000?, 1), (1.000000000001?, 1)]
>>> v = Integer(2)**Integer(60)
>>> real_roots((x-Integer(1)) * (x-(v+Integer(1))/v), retval='interval')
[(1.000000000000000000?, 1), (1.000000000000000001?, 1)]
>>> real_roots((x-Integer(1)) * (x-Integer(2)), strategy='warp')
[((499/525, 1173/875), 1), ((337/175, 849/175), 1)]
>>> real_roots((x+Integer(3))*(x+Integer(1))*x*(x-Integer(1))*(x-Integer(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)]
>>> real_roots((x+Integer(3))*(x+Integer(1))*x*(x-Integer(1))*(x-Integer(2)), strategy='warp', retval='algebraic_real')
[(-3.000000000000000?, 1), (-1.000000000000000?, 1), (0, 1), (1.000000000000000?, 1), (2.000000000000000?, 1)]
>>> ar_rts = real_roots(x-Integer(1), retval='algebraic_real')
>>> ar_rts[Integer(0)][Integer(0)] == Integer(1)
True

If the polynomial has no real roots, we get an empty list.

Sage

sage: (x^2 + 1).real_root_intervals()
[]

Python

>>> from sage.all import *
>>> (x**Integer(2) + Integer(1)).real_root_intervals()
[]

We can compute Conway’s constant (see http://mathworld.wolfram.com/ConwaysConstant.html) to arbitrary precision.

Sage

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

Python

>>> from sage.all import *
>>> p = x**Integer(71) - x**Integer(69) - Integer(2)*x**Integer(68) - x**Integer(67) + Integer(2)*x**Integer(66) + Integer(2)*x**Integer(65) + x**Integer(64) - x**Integer(63) - x**Integer(62) - x**Integer(61) - x**Integer(60) - x**Integer(59) + Integer(2)*x**Integer(58) + Integer(5)*x**Integer(57) + Integer(3)*x**Integer(56) - Integer(2)*x**Integer(55) - Integer(10)*x**Integer(54) - Integer(3)*x**Integer(53) - Integer(2)*x**Integer(52) + Integer(6)*x**Integer(51) + Integer(6)*x**Integer(50) + x**Integer(49) + Integer(9)*x**Integer(48) - Integer(3)*x**Integer(47) - Integer(7)*x**Integer(46) - Integer(8)*x**Integer(45) - Integer(8)*x**Integer(44) + Integer(10)*x**Integer(43) + Integer(6)*x**Integer(42) + Integer(8)*x**Integer(41) - Integer(5)*x**Integer(40) - Integer(12)*x**Integer(39) + Integer(7)*x**Integer(38) - Integer(7)*x**Integer(37) + Integer(7)*x**Integer(36) + x**Integer(35) - Integer(3)*x**Integer(34) + Integer(10)*x**Integer(33) + x**Integer(32) - Integer(6)*x**Integer(31) - Integer(2)*x**Integer(30) - Integer(10)*x**Integer(29) - Integer(3)*x**Integer(28) + Integer(2)*x**Integer(27) + Integer(9)*x**Integer(26) - Integer(3)*x**Integer(25) + Integer(14)*x**Integer(24) - Integer(8)*x**Integer(23) - Integer(7)*x**Integer(21) + Integer(9)*x**Integer(20) + Integer(3)*x**Integer(19) - Integer(4)*x**Integer(18) - Integer(10)*x**Integer(17) - Integer(7)*x**Integer(16) + Integer(12)*x**Integer(15) + Integer(7)*x**Integer(14) + Integer(2)*x**Integer(13) - Integer(12)*x**Integer(12) - Integer(4)*x**Integer(11) - Integer(2)*x**Integer(10) + Integer(5)*x**Integer(9) + x**Integer(7) - Integer(7)*x**Integer(6) + Integer(7)*x**Integer(5) - Integer(4)*x**Integer(4) + Integer(12)*x**Integer(3) - Integer(6)*x**Integer(2) + Integer(3)*x - Integer(6)
>>> cc = real_roots(p, retval='algebraic_real')[Integer(2)][Integer(0)] # long time
>>> RealField(Integer(180))(cc)                                # long time
1.3035772690342963912570991121525518907307025046594049

Now we play with algebraic real coefficients.

Sage

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)]

Python

>>> from sage.all import *
>>> x = polygen(AA)
>>> p = (x - Integer(1)) * (x - sqrt(AA(Integer(2)))) * (x - Integer(2))
>>> real_roots(p)
[((499/525, 2171/1925), 1), ((1173/875, 2521/1575), 1), ((337/175, 849/175), 1)]
>>> ar_rts = real_roots(p, retval='algebraic_real'); ar_rts
[(1.000000000000000?, 1), (1.414213562373095?, 1), (2.000000000000000?, 1)]
>>> ar_rts[Integer(1)][Integer(0)]**Integer(2) == Integer(2)
True
>>> ar_rts = real_roots(x*(x-Integer(1)), retval='algebraic_real')
>>> ar_rts[Integer(0)][Integer(0)] == Integer(0)
True
>>> p2 = p * (p - Integer(1)/Integer(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?
>>> real_roots(p2, retval='interval')
[(1.00?, 1), (1.1?, 1), (1.38?, 1), (1.5?, 1), (2.00?, 1), (2.1?, 1)]
>>> p = (x - Integer(1)) * (x - sqrt(AA(Integer(2))))**Integer(2) * (x - Integer(2))**Integer(3) * sqrt(AA(Integer(3)))
>>> real_roots(p, retval='interval')
[(1.000000000000000?, 1), (1.414213562373095?, 2), (2.000000000000000?, 3)]

sage.rings.polynomial.real_roots.relative_bounds(a, b)

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

sage: from sage.rings.polynomial.real_roots import *
sage: relative_bounds((1/7, 1/4), (1/6, 1/5))
(2/9, 8/15)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> relative_bounds((Integer(1)/Integer(7), Integer(1)/Integer(4)), (Integer(1)/Integer(6), Integer(1)/Integer(5)))
(2/9, 8/15)

sage.rings.polynomial.real_roots.reverse_intvec(c)

Given a vector of integers, reverse the vector (like the reverse() method on lists).

Modifies the input vector; has no return value.

EXAMPLES:

Sage

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> v = vector(ZZ, [Integer(1), Integer(2), Integer(3), Integer(4)]); v
(1, 2, 3, 4)
>>> reverse_intvec(v)
>>> v
(4, 3, 2, 1)

sage.rings.polynomial.real_roots.root_bounds(p)

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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(ZZ)
>>> root_bounds((x-Integer(1))*(x-Integer(2))*(x-Integer(3)))
(0.545454545454545, 6.00000000000001)
>>> root_bounds(x**Integer(2))
(0.000000000000000, 0.000000000000000)
>>> root_bounds(x*(x+Integer(1)))
(-1.00000000000000, 0.000000000000000)
>>> root_bounds((x+Integer(2))*(x-Integer(3)))
(-2.44948974278317, 3.46410161513776)
>>> root_bounds(x**Integer(995) * (x**Integer(2) - Integer(9999)) - Integer(1))
(-99.9949998749937, 141.414284992713)
>>> root_bounds(x**Integer(995) * (x**Integer(2) - Integer(9999)) + Integer(1))
(-141.414284992712, 99.9949998749938)

If we can see that the polynomial has no real roots, return None.

Sage

sage: root_bounds(x^2 + 1) is None
True

Python

>>> from sage.all import *
>>> root_bounds(x**Integer(2) + Integer(1)) is None
True

class sage.rings.polynomial.real_roots.rr_gap

Bases: object

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.

region()

sage.rings.polynomial.real_roots.scale_intvec_var(c, k)

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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> v = vector(ZZ, [Integer(1), Integer(1), Integer(1), Integer(1)])
>>> scale_intvec_var(v, Integer(3)/Integer(4))
>>> v
(64, 48, 36, 27)

sage.rings.polynomial.real_roots.split_for_targets(ctx, bp, target_list, precise=False)

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

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]>

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> bp = mk_ibpi([Integer(1000000), -Integer(2000000), Integer(3000000), -Integer(4000000), -Integer(5000000), -Integer(6000000)])
>>> ctx = mk_context()
>>> bps = split_for_targets(ctx, bp, [(rr_gap(Integer(1)/Integer(1234567893), Integer(1)/Integer(1234567892), Integer(1)), rr_gap(Integer(1)/Integer(1234567891), Integer(1)/Integer(1234567890), Integer(1)), Integer(12)), (rr_gap(Integer(1)/Integer(3), Integer(1)/Integer(2), -Integer(1)), rr_gap(Integer(2)/Integer(3), Integer(3)/Integer(4), -Integer(1)), Integer(6))])
>>> bps[Integer(0)]
<IBP: (999992, 999992, 999992) + [0 .. 15) over [8613397477114467984778830327/10633823966279326983230456482242756608 .. 591908168025934394813836527495938294787/730750818665451459101842416358141509827966271488]; level 2; slope_err 0.?e12>
>>> bps[Integer(1)]
<IBP: (-1562500, -1875001, -2222223, -2592593, -2969137, -3337450) + [0 .. 4) over [1/2 .. 2863311531/4294967296]>

sage.rings.polynomial.real_roots.subsample_vec_doctest(a, slen, llen)

sage.rings.polynomial.real_roots.taylor_shift1_intvec(c)

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

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))

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(ZZ)
>>> p = (x-Integer(1))*(x-Integer(2))*(x-Integer(3))
>>> v = vector(ZZ, p.list())
>>> p, v
(x^3 - 6*x^2 + 11*x - 6, (-6, 11, -6, 1))
>>> taylor_shift1_intvec(v)
>>> p(x+Integer(1)), v
(x^3 - 3*x^2 + 2*x, (0, 2, -3, 1))

sage.rings.polynomial.real_roots.to_bernstein(p, low=0, high=1, degree=None)

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

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)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(ZZ)
>>> to_bernstein(x)
([0, 1], 1)
>>> to_bernstein(x, degree=Integer(5))
([0, 1/5, 2/5, 3/5, 4/5, 1], 1)
>>> to_bernstein(x**Integer(3) + x**Integer(2) - x - Integer(1), low=-Integer(3), high=Integer(3))
([-16, 24, -32, 32], 1)
>>> to_bernstein(x**Integer(3) + x**Integer(2) - x - Integer(1), low=Integer(3), high=Integer(22)/Integer(7))
([296352, 310464, 325206, 340605], 9261)

sage.rings.polynomial.real_roots.to_bernstein_warp(p)

Given a polynomial p with rational coefficients, compute the exact rational Bernstein coefficients of p(x/(x+1)).

EXAMPLES:

Sage

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]

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> x = polygen(ZZ)
>>> to_bernstein_warp(Integer(1) + x + x**Integer(2) + x**Integer(3) + x**Integer(4) + x**Integer(5))
[1, 1/5, 1/10, 1/10, 1/5, 1]

class sage.rings.polynomial.real_roots.warp_map(neg)

Bases: object

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].

from_ocean(region)

to_ocean(region)

sage.rings.polynomial.real_roots.wordsize_rational(a, b, wordsize)

Given rationals a and b, select 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:

This takes as input two rationals, a and b, such that 0<=a<=1, 0<=b<=1, and a!=b. This 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

…

11. 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

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

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.real_roots import *
>>> wordsize_rational(Integer(1)/Integer(5), Integer(1)/Integer(7), Integer(32))
429496729/2147483648
>>> wordsize_rational(Integer(1)/Integer(7), Integer(1)/Integer(5), Integer(32))
306783379/2147483648
>>> wordsize_rational(Integer(1)/Integer(5), Integer(1)/Integer(7), Integer(64))
1844674407370955161/9223372036854775808
>>> wordsize_rational(Integer(1)/Integer(7), Integer(1)/Integer(5), Integer(64))
658812288346769701/4611686018427387904
>>> wordsize_rational(Integer(1)/Integer(17), Integer(1)/Integer(19), Integer(32))
252645135/4294967296
>>> wordsize_rational(Integer(1)/Integer(17), Integer(1)/Integer(19), Integer(64))
1085102592571150095/18446744073709551616
>>> wordsize_rational(Integer(1)/Integer(1234567890), Integer(1)/Integer(1234567891), Integer(32))
933866427/1152921504606846976
>>> wordsize_rational(Integer(1)/Integer(1234567890), Integer(1)/Integer(1234567891), Integer(64))
4010925763784056541/4951760157141521099596496896