Interval times Int64 and Interval times BigInt discrepancy

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

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]
    CauchyPolynomial(ar, tuple(d...))

@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] ) )
    return r

@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] ) )
    return r

function pval(p::CauchyPolynomial{M,N}, x::R) where M where N where R<:Real
    return pval_g(p.c, x)

function dval(p::CauchyPolynomial{M,N}, x::R) where M where N where R<:Real
    return dval_g(p.d, x)

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
    throw(ErrorException("Newton method timed out, p: "*string(p)))


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.

For relatively short polynomials, the best way to evaluate this will almost certainly be
pval_g(c,x) = evalpoly(x, tuple(c... ,-1)).
That said, BigInt will always be way slower since Int is doing the computation modulo 2^64.

I think my function is faster since it is a generated function that takes into account that the highest order coefficient is always -1. If it is not, then I’m very curious to see why it is not.

(without looking at your code)
Ints (the default Integer type) are bounded by (typemin(Int), typemax(Int)) and they will wrap around:

julia> 0 > typemax(Int) + 2

So your evaluation logic may be causing the root finder to bounce back and forth ad infinitum when Ints are used with your specific calculation.

BigInts grow and grow, allowing the root finder to utilize whatever integer precision is required in this case.


The “point” to be evaluated is Interval{BigFloat}. So the result of any computation is also Interval{BigFloat}.

There was a recent thread about special methods to find roots of polynomials with integer coefficients.

I have seen the that thread, but it was not very useful to me.

Regardless, the integer type should not affect the result of the product. Am I missing something here? Is this not a bug?

I was an idiot, it overflows when I compute the coefficients of the derivative. :man_facepalming: