I want to get the positive root of some polynomials with ridiculous precision. The polynomials I care I about are of the form \sum_{k=0}^{n-1} a_k x^k - x^n, with a_k's being natural numbers. I am trying to write a program that does that efficiently and I found a weird bug.

The following is a minimal example. Unfortunately, I couldn’t minimise it any more.

```
using IntervalArithmetic
struct CauchyPolynomial{T,S,I<:Union{Int,Int128,BigInt}}
c::NTuple{T,I}
d::NTuple{S,I}
end
function CauchyPolynomial(ar::NTuple{N,T}) where N where T <:Union{Int,Int128,BigInt}
d = zeros(T,length(ar)-1)
for i in 2:length(ar)
d[i-1] = (i-1)*ar[i]
end
CauchyPolynomial(ar, tuple(d...))
end
@generated function pval_g(c::NTuple{N,I}, x::R) where I<:Union{Int,Int128,BigInt} where R<:Real where N
r = :(-one(R))
for i in N:-1:1
r = :(muladd(x, $r, c[$i] ) )
end
return r
end
@generated function dval_g(d::NTuple{N,I}, x::R) where I<:Union{Int,Int128,BigInt} where R<:Real where N
r = :(-one(R)*(N+1))
for i in N:-1:1
r = :(muladd(x, $r, d[$i] ) )
end
return r
end
function pval(p::CauchyPolynomial{M,N}, x::R) where M where N where R<:Real
return pval_g(p.c, x)
end
function dval(p::CauchyPolynomial{M,N}, x::R) where M where N where R<:Real
return dval_g(p.d, x)
end
function root(p::CauchyPolynomial{M,N}, bts::Int) where M where N
println("root called: ",p," - ",bts)
m = max(p.c...)
x = @biginterval m
while -log2(diam(x)) > bts
for i in 1:100
pt = pval(p,x)
x = x - pt/dval(p,x)
if 0 in pt
return x
end
end
end
throw(ErrorException("Newton method timed out, p: "*string(p)))
end
setprecision(8384);
ts = (832152526037543613, 4191427560420539572, 1948805851146604192, 41780545826284041, 2952882476324826786);
tb = map(big, ts);
ps = CauchyPolynomial(ts);
pb = CauchyPolynomial(tb);
root(pb, 8303)
root(ps, 8303)
```

When I execute this, the line `root(pb, 8303)`

gives me the root but the line `root(ps, 8303)`

times out.

The only difference between the polynomials `ps`

and `pb`

is that the coefficients of `ps`

are `Int64`

and the ones of `pb`

are `BigInt`

. Since integers are represented exactly, I expected that the difference in types wouldn’t matter. What is happening here?

Also, separate to this issue, any advice on how to improve the code is more than welcomed.