Interval times Int64 and Interval times BigInt discrepancy

General Usage
,
1

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.

2

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.

3

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.

4

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

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.

4 Likes
5

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

6

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

7

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?

8

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

2 Likes