# Root isolation of real-rooted integer polynomials
**URL:** https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421
**Category:** General Usage
**Tags:** question, roots, polynomials, nemo
**Created:** [December 7, 2020, 6:02pm UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421 "2020-12-07T18:02:46Z")
**Posts on this page:** 13
**Page:** 1
### Author:  [@malacroi](https://discourse.julialang.org/u/malacroi)
#### Post date: [December 7, 2020, 6:02pm UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/1 "2020-12-07T18:02:46Z")
I have a collection of polynomials with non-negative integer coefficients that I know to be real-rooted with no repeated roots by examining Sturm chains for each. I would like to show that they form an interlacing family (hence their convex combinations are real-rooted). Coming from `Maple`, I would use the `realroot` command (based on Descartes’ rule of signs) to find rationally bounded isolating intervals for the real roots of the polynomials. Is there a package that implements something equivalent in `Julia`? For a ballpark idea of scale, I’d like to be able to work with 10^5 polynomials with degree 80, and coefficients on the order of 2^480, but something working out-of-the-box on smaller instances would still be helpful.
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [December 7, 2020, 6:28pm UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/2 "2020-12-07T18:28:03Z")
You will likely find `IntervalRootFinding` to be helpful here. It will give rigorous bounds. (The issue with `roots` from Polynomials is the high degree will likely yield numeric issues). Here is a sample usage, though it sounds like you might use the `Polynomial(coeffs)` constructor:
```julia
julia> using Polynomials, IntervalRootFinding, IntervalArithmetic
julia> x = [2,3,4,5,7,8,9];
julia> p = fromroots(Polynomial, x)
Polynomial(-60480 + 100536*x - 68078*x^2 + 24433*x^3 - 5036*x^4 + 598*x^5 - 38*x^6 + x^7)
julia> IntervalRootFinding.roots(x->p(x), -20..20)
7-element Array{Root{Interval{Float64}},1}:
Root([2.99999, 3.00001], :unique)
Root([1.99999, 2.00001], :unique)
Root([7.99999, 8.00001], :unique)
Root([4.99999, 5.00001], :unique)
Root([6.99999, 7.00001], :unique)
Root([8.99999, 9.00001], :unique)
Root([3.99999, 4.00001], :unique)
```
The `:unique` shows there is a root and only 1 in the specified interval.
The choice of `20..20` can be widened without issue, and could be computed using some bounds on real roots based on coefficients.
To show the nesting, you might do something like this:
```julia
xs = IntervalRootFinding.roots(x->p(x), -20..20)
ys = IntervalRootFinding.roots(x->Polynomials.derivative(p)(x), -20..20)
xs = sort([i.interval for i in xs])# grab intervals only
ys = sort([i.interval for i in ys]);
all(xs[i] < ys[i] < xs[i+1] for i in 1:length(xs)-1) # true
```
---
### Author:  [@malacroi](https://discourse.julialang.org/u/malacroi)
#### Post date: [December 8, 2020, 1:24am UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/3 "2020-12-08T01:24:54Z")
`IntervalRootFinding.roots` doesn’t seem to be able to take advantage of the fact that the coefficients are integers. I’m finding it chokes on relatively low degree polynomials.
```julia
using IntervalRootFinding, IntervalArithmetic, Polynomials
bp = [
942438915208811912419937422298363203125
164182217245953398816894035758761902846875
4584900574568933770264468813466870772155175
48995332393110515074735075708247882042540865
266674183150777010544241114017741621823207005
852443280934837985352088128423887894438557515
1738546146302892245736990844587487000484756535
2381558158813900978436742173305983349418813145
2262003889258248241081177038309445610985409335
1516025051068561122302699855213604175083575145
720810987764796866354279087485114863858085005
241341213302325116160821160849326697595681275
55563538585119205063994483187179167406616375
8363912237256094118085070946083688310200625
740493466743082745510080711751444519503125
29215606371473169285018060091249259296875]
intervalroots(a) = IntervalRootFinding.roots(x->a(x), -big(4)..big(0))
poly = Polynomials.Polynomial(bp)
@btime intervalroots(poly)
```
this takes **minutes** at 100% CPU usage, though relative low memory. In contrast, in `Maple`, I can run
```nohighlight
cl := [ 942438915208811912419937422298363203125,
164182217245953398816894035758761902846875,
4584900574568933770264468813466870772155175,
48995332393110515074735075708247882042540865,
266674183150777010544241114017741621823207005,
852443280934837985352088128423887894438557515,
1738546146302892245736990844587487000484756535,
2381558158813900978436742173305983349418813145,
2262003889258248241081177038309445610985409335,
1516025051068561122302699855213604175083575145,
720810987764796866354279087485114863858085005,
241341213302325116160821160849326697595681275,
55563538585119205063994483187179167406616375,
8363912237256094118085070946083688310200625,
740493466743082745510080711751444519503125,
29215606371473169285018060091249259296875];
poly := add(cl[i]*x^(i-1), i=1..nops(cl));
realroot(poly,1/2^200);
```
to get interval that isolate to intervals of width 2^(-200) in on the order of 1/50th of a second.
I have a naive, hacked together, Sturm chain approach in Julia that gets me much wider intervals in a couple of seconds, but I can’t see how I’m going to make up that many orders of magnitude of difference.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [December 8, 2020, 1:38am UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/4 "2020-12-08T01:38:36Z")
I’m way out of my depth here, but have you looked at [Real-root isolation - Wikipedia](https://en.wikipedia.org/wiki/Real-root_isolation)?
---
### Author:  [@malacroi](https://discourse.julialang.org/u/malacroi)
#### Post date: [December 8, 2020, 2:17am UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/5 "2020-12-08T02:17:01Z")
That looks like about what I need. The _Recent Improvements_ section there, sure makes it look like there should be some interval arithmetic trick close by. I guess tomorrow’s project is to try and limit the number of intermediate values for all those change of variables clean when working with `BigInts`.
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [December 8, 2020, 6:40am UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/6 "2020-12-08T06:40:23Z")
I don’t see how the fact that the coefficients are integers will help? But it’s very interesting that Maple is so fast.
You might also be interested in HomotopyContinuation.jl
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [December 8, 2020, 6:53am UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/7 "2020-12-08T06:53:38Z")
> [@malacroi](#):
>
> ```julia
> intervalroots(a) = IntervalRootFinding.roots(x->a(x), -big(4)..big(0))
> poly = Polynomials.Polynomial(bp)
>
> ```
It seems that you can get some initial information quickly using Float64:
```julia
julia> bp2 = Interval{Float64}.(bp);
julia> poly2 = Polynomials.Polynomial(bp2);
julia> roots(x -> poly2(x), -Inf..Inf, Bisection, 1.0)
22-element Vector{Root{Interval{Float64}}}:
Root([-6.95784, -6.40872], :unknown)
Root([-9.129, -8.57124], :unknown)
[etc]
```
Then hopefully refining with BigFloat will be faster.
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [December 8, 2020, 6:57am UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/8 "2020-12-08T06:57:07Z")
For fast polynomial manipulations in Julia I believe you can use [https://github.com/Nemocas/AbstractAlgebra.jl](https://github.com/Nemocas/AbstractAlgebra.jl)
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [December 8, 2020, 7:03am UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/9 "2020-12-08T07:03:41Z")
If you have an approximate root, you can verify it, e.g.:
```julia
julia> using Roots
julia> rts = find_zeros(x -> poly(x), -10, 10); # approximate: Roots.jl
julia> big_roots = roots(x -> poly(x), big(rts[3] ± 1e-5), Newton, 1e-50) # rigorous: IntervalRootFinding.jl
1-element Vector{Root{Interval{BigFloat}}}:
Root([-3.39911, -3.3991]₂₅₆, :unique)
julia> diam(interval(big_roots[1]))
7.879417721966192043308957487195098830393286784342275257087148835580431618280549e-70
```
This finds an interval of with ~10^(-69) in which a unique root is proved to exist.
---
### Author:  [@thofma](https://discourse.julialang.org/u/thofma)
#### Post date: [December 8, 2020, 1:03pm UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/10 "2020-12-08T13:03:58Z")
You can also use `Hecke/Nemo`. Here are the timings on my very dated laptop:
```julia
julia> using Hecke
julia> Qx, x = PolynomialRing(QQ)
(Univariate Polynomial Ring in x over Rational Field, x)
julia> Qx, x = PolynomialRing(QQ);
julia> f = Qx(bp);
julia> Hecke._roots(f, 200) # the balls will have radious < 2^(-200)
15-element Array{acb,1}:
[-3.920616711153724253790204590529063281656171296624620351746200485310655120855543613013541271680907275607223979079 +/- 4.84e-112] + i*0
[-3.706111517064614641012473797293410340554769655675971950139715152216337896330719648922654909467811143497620428030 +/- 7.80e-112] + i*0
[-3.399101345956385184809590396663823245233910492969600172371668237761238821534096232299586114103416940774067560924 +/- 8.79e-112] + i*0
[-3.035412954945010421040655247164632911783923123071964813560272462311419306648072410827429032399236488318520657029 +/- 9.05e-112] + i*0
[-2.642241820837020072256321446404090852721123523348319017917035488525807038562259125380960186537327785108752839951 +/- 5.03e-112] + i*0
[-2.240386668347612606849921371192140625120311772030273906303162991172574142894086307045961801604703579352289007669 +/- 3.04e-112] + i*0
[-1.8459392275279762188118911328074682418106652190523122728768478452769447419809494076465189839407551707688825849549 +/- 9.17e-113] + i*0
[-1.4714364265485984347457103358546916040139602018877260092947033620076406186360124337515593788219523840023676498250 +/- 3.09e-113] + i*0
[-1.12667548938228648371659194465060333123610781848166932944102272185644329811383422404304887752307268162794786601481 +/- 6.41e-114] + i*0
[-0.819305999801997298040859319753994185201032318699276601025148883053951656453384749232848418410546346810549845454155 +/- 6.49e-115] + i*0
[-0.5552656322189970256982534594202473203239575869892228528269814900273624281956924496928014910474082527863381439566599 +/- 8.75e-116] + i*0
[-0.33910202510401947795464384494964958658273285483978651160525210708434910971533965735125598235860614612011137885597364 +/- 4.76e-117] + i*0
[-0.174209253846368310831930193848706110784111374105168275102590356950016921390554301061214369088829872621188752931599096 +/- 4.12e-118] + i*0
[-0.06299838108334910672357500092170790024222765409891867028284948240638078820689793485367480324463392976269126659634877413 +/- 9.09e-120] + i*0
[-0.007015398195279156379264993894770551409664575420748371105657755069724334773977198396431121866041220489066572460776894302 +/- 2.04e-121] + i*0
julia> @btime Hecke._roots(f, 200);
4.501 ms (3677 allocations: 204.58 KiB)
```
This is internally using [http://arblib.org/](http://arblib.org/).
P.S. To convert to `BigFloat` you can do `BigFloat(real(x))`, but you must handle the global precision yourself.
---
### Author:  [@malacroi](https://discourse.julialang.org/u/malacroi)
#### Post date: [December 8, 2020, 4:52pm UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/11 "2020-12-08T16:52:56Z")
> [@dpsanders](#):
>
> I don’t see how the fact that the coefficients are integers will help? But it’s very interesting that Maple is so fast.
I’m not really clear myself, I just know from browsing the `Maple` documentation, that they rely on integer (or rational) coefficients, presumably because it simplifies something. Page 2 of [https://arxiv.org/pdf/1605.00410.pdf](https://arxiv.org/pdf/1605.00410.pdf) seems to suggest that one of the problems is being able to guarantee that it’s possible to determine the sign of a polynomial near its roots. It also seems (page 8, paragraph 2), that there might be some implications to a heuristic that switches between MPFI and exact arithmetic depending on which is expected to be cheaper for a given subproblem.
Another guarantee from integer coefficients is that it’s possible to pre-filter and ensure no rational roots, so we don’t have to worry about corner cases where a root falls exactly at an end-point of an interval.
---
### Author:  [@malacroi](https://discourse.julialang.org/u/malacroi)
#### Post date: [December 8, 2020, 5:23pm UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/12 "2020-12-08T17:23:31Z")
Thanks. That’s definitely fast enough to fit into my pipeline.
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [December 9, 2020, 6:37pm UTC](https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/13 "2020-12-09T18:37:33Z")
I wanted to see if my initial suggestion to use `IntervalRootFinding` could be sped up paying more attention to the polynomial evaluation time and the underlying numbers used, so tried this:
```julia
using IntervalRootFinding, IntervalArithmetic, Polynomials
bp = [
942438915208811912419937422298363203125
164182217245953398816894035758761902846875
4584900574568933770264468813466870772155175
48995332393110515074735075708247882042540865
266674183150777010544241114017741621823207005
852443280934837985352088128423887894438557515
1738546146302892245736990844587487000484756535
2381558158813900978436742173305983349418813145
2262003889258248241081177038309445610985409335
1516025051068561122302699855213604175083575145
720810987764796866354279087485114863858085005
241341213302325116160821160849326697595681275
55563538585119205063994483187179167406616375
8363912237256094118085070946083688310200625
740493466743082745510080711751444519503125
29215606371473169285018060091249259296875]
bps = [convert(ArbFloat, string(b)) for b in bp]
p = ImmutablePolynomial(bps)
# error with show
# Error showing value of type Tuple{ArbFloat{116},ArbFloat{116}}:
# ERROR: MethodError: no method matching Array{String,1}(::Int64)
a,b = -convert(ArbFloat, "4"), convert(ArbFloat, "0");
IntervalRootFinding.roots(f,a..b) # error about a..b
#ERROR: MethodError: no method matching ArbFloat{116}(::String)
IntervalRootFinding.roots(f,Interval(a,b))
# ERROR: MethodError: *(::ArbFloat{116}, ::Interval{ArbFloat{116}}) is ambiguous. Candidates:
```
So a bunch of errors, the last one non-circumventable. I could file issues with `ArbFloats` for the first two, but the method ambiguity doesn’t have a natural place to resolve.