Vincent-Collins-Akritas real root isolation algorithm.

gistfile1.txt

# Define a Poly as a vector of coefficients [c_0, c_1, ..., c_n]
# representing the polynomial c_0 + c_1*x + ... + c_n*x^n
const Poly = Vector{Rational{BigInt}}

# Evaluate polynomial at a given exact point using Horner's method
function evaluate(P::Poly, x::Rational{BigInt})
    res = zero(Rational{BigInt})
    for i in length(P):-1:1
        res = res * x + P[i]
    end
    return res
end

# Add two polynomials
function poly_add(A::Poly, B::Poly)
    n = max(length(A), length(B))
    C = zeros(Rational{BigInt}, n)
    for i in 1:length(A); C[i] += A[i]; end
    for i in 1:length(B); C[i] += B[i]; end
    return C
end

# Multiply two polynomials
function poly_mul(A::Poly, B::Poly)
    if isempty(A) || isempty(B) return Poly() end
    C = zeros(Rational{BigInt}, length(A) + length(B) - 1)
    for i in 1:length(A)
        for j in 1:length(B)
            C[i+j-1] += A[i] * B[j]
        end
    end
    return C
end

# Count the number of sign variations in the coefficient sequence
function sign_variations(P::Poly)
    vars = 0
    last_sign = 0
    for c in P
        if c > 0
            if last_sign < 0 vars += 1 end
            last_sign = 1
        elseif c < 0
            if last_sign > 0 vars += 1 end
            last_sign = -1
        end
    end
    return vars
end

# Computes Q(x) = (x+1)^n * P((a*x+b)/(x+1))
# This maps the interval (a, b) to (0, ∞)
function mobius_transform(P::Poly, a::Rational{BigInt}, b::Rational{BigInt})
    n = length(P) - 1
    Q = zeros(Rational{BigInt}, 1)
    
    num = Poly([b, a]) # Represents b + ax
    den = Poly([1, 1]) # Represents 1 + x
    
    # Precompute powers of the numerator and denominator to build the transformation
    num_pow = Vector{Poly}(undef, n + 1)
    den_pow = Vector{Poly}(undef, n + 1)
    num_pow[1] = Poly([1//1])
    den_pow[1] = Poly([1//1])
    
    for i in 1:n
        num_pow[i+1] = poly_mul(num_pow[i], num)
        den_pow[i+1] = poly_mul(den_pow[i], den)
    end
    
    for i in 0:n
        term = poly_mul(num_pow[i+1], den_pow[n-i+1])
        term = term .* P[i+1]
        Q = poly_add(Q, term)
    end
    
    return Q
end

# Calculate Cauchy's root bound to find the initial search interval
function root_bound(P::Poly)
    n = length(P)
    an = abs(P[n])
    max_val = zero(Rational{BigInt})
    for i in 1:(n-1)
        max_val = max(max_val, abs(P[i]))
    end
    return 1//1 + max_val / an
end

# Recursive bisection to isolate roots in (a, b)
function isolate_interval(P::Poly, a::Rational{BigInt}, b::Rational{BigInt}, roots_list)
    Q = mobius_transform(P, a, b)
    v = sign_variations(Q)
    
    if v == 0
        return # No roots in this interval
    elseif v == 1
        push!(roots_list, (a, b)) # Exactly one root isolated!
        return
    else
        # If v > 1, roots exist (or complex roots are mimicking real ones). Bisect!
        m = (a + b) / 2
        val_m = evaluate(P, m)
        
        if val_m == 0
            # m is an exact rational root
            isolate_interval(P, a, m, roots_list)
            push!(roots_list, (m, m))
            isolate_interval(P, m, b, roots_list)
        else
            isolate_interval(P, a, m, roots_list)
            isolate_interval(P, m, b, roots_list)
        end
    end
end

# 1. Add the refinement function
function refine_interval(P::Poly, a::Rational{BigInt}, b::Rational{BigInt}, tolerance::Rational{BigInt})
    # If it's already an exact root (a == b), no refinement needed
    if a == b return (a, b) end 

    while (b - a) > tolerance
        m = (a + b) / 2
        # Use evaluate to check signs and decide which half to keep
        if evaluate(P, a) * evaluate(P, m) <= 0
            b = m
        else
            a = m
        end
    end
    return (a, b)
end

# 2. Update the main wrapper to handle an optional tolerance
function collins_akritas(coeffs::Vector{<:Real}; tolerance::Union{Rational{BigInt}, Nothing} = nothing)
    P = Poly([Rational{BigInt}(c) for c in coeffs])
    
    zero_root = false
    while length(P) > 1 && P[1] == 0
        popfirst!(P)
        zero_root = true
    end
    
    zero_big = zero(Rational{BigInt}) 
    
    roots_list = Tuple{Rational{BigInt}, Rational{BigInt}}[]
    if length(P) <= 1
        if zero_root push!(roots_list, (zero_big, zero_big)) end
        return roots_list
    end
    
    R = root_bound(P)
    
    # Isolate Positive Roots
    pos_roots = Tuple{Rational{BigInt}, Rational{BigInt}}[]
    isolate_interval(P, zero_big, R, pos_roots)
    
    # Isolate Negative Roots
    P_neg = Poly([P[i] * (iseven(i-1) ? 1 : -1) for i in 1:length(P)])
    neg_roots_pos = Tuple{Rational{BigInt}, Rational{BigInt}}[]
    isolate_interval(P_neg, zero_big, R, neg_roots_pos)
    
    # Map back to negative intervals
    neg_roots = [(-b, -a) for (a, b) in reverse(neg_roots_pos)]
    
    # Combine results
    append!(roots_list, neg_roots)
    if zero_root push!(roots_list, (zero_big, zero_big)) end
    append!(roots_list, pos_roots)
    
    # ==========================================
    # THE REFINEMENT STEP
    # ==========================================
    if tolerance !== nothing
        # Use the deflated P, NOT P_original, to avoid evaluating to 0 
        # at the interval boundaries if 0 was an exact root!
        for i in 1:length(roots_list)
            a, b = roots_list[i]
            roots_list[i] = refine_interval(P, a, b, tolerance)
        end
    end
    
    return roots_list
end

# Helper to print the isolated intervals nicely
function print_isolation(name, coeffs, tolerance::Union{Rational{BigInt}, Nothing} = nothing)
    intervals = collins_akritas(coeffs, tolerance=tolerance)
    println("Polynomial: $name")
    println("Coefficients: $coeffs")
    println("Isolated Intervals (Exact Rational Boundaries):")
    for (a, b) in intervals
        if a == b
            println("  Exact root at x = $a")
        else
            println("  Root isolated in: ($a, $b)")
            println("     (Approx: $(Float64(a)) to $(Float64(b)))")
        end
    end
    println("-"^50)
end