_check_root_interval(f, a::Arf, b::Arf, fa_sign::Integer = Arblib.sgn_nonzero(f(Arb(a))), fb_sign::Integer = Arblib.sgn_nonzero(f(Arb(b))); check_unique = true, buffer1 = zero(Arb), buffer2 = zero(Arb))

Check if the function f has a zero on the interval [a, b].

Returns (maybe, unique). If maybe is false then f is proved to not have a zero on the interval, if it's true then f might or might not have a zero on the interval. If unique is true then f is proved to have exactly one zero on the interval, if it's false then no information is given.

The arguments fa_sign and fb_sign should be set to the sign of f on the endpoints a and b according to how Arblib.sgn_nonzero does it. They can be computed as

fa_sign = Arblib.sgn_nonzero(f(Arb(a)))
fb_sign = Arblib.sgn_nonzero(f(Arb(b)))

If check_unique = false then don't check for uniqueness, unique will always be false in this case.

The function f should satisfy the same properties as for isolate_roots.

The arguments buffer1 and buffer2 can optionally be given to be used as scratch space for the computations. This reduces the number of allocations required.

_check_signs(pa, pb, values = ())

Check the signs of pa, pb and values. Returns 1 if all of them definitely have the same sign, -1 if there are at least two who definitely have different signs and 0 otherwise.

_compute_endpoints(a::Arf, b::Arf, x::Arb)

Given a, b and x with x = Arb((a, b)) compute

c = a - midpoint(x)
d = b - midpoint(x)

rounded towards zero. It returns c, c_exact, d, d_exact where c and d are as above and c_exact is true if c was computed exactly and false if rounding was performed, similarly for d_exact.

This is used in extrema_series, minimum_series and maximum_series to compute the interval on which the extrema of the polynomial should be enclosed.

If c_exact is false then an enclosure of the endpoint can be computed with _enclose_inexact_endpoint.

_current_lower_upper_bound(minormax, low, upp, values_low, values_upp, values)

Returns a lower and upper bound of the extremum given that low and upp and lower and upper bounds of the extremum and so are all values in values_low and Values_upp respectively. It determines either the minimum or the maximum depending on if minormax = min or minormax = max.

This is mainly an internal function which implements behaviour which is common for extrema_enclosure, minimum_enclosure and maximum_enclosure.

_enclose_inexact_endpoint(endpoint::Arf)

Compute a ball centered at endpoint with a radius of 1 ulp of endpoint.

This is used to get an enclosure of the endpoint in case _compute_endpoints gives an inexact result. The semantics of Arb guarantees that the error of the endpoint is at most 1 ulp when using directed rounding and a ball with endpoint as midpoint and a radius set to the value of this 1 ulp will therefore always enclose the exact result.

_extrema_polynomial_low_degree(p::ArbPoly, a::Arf, b::Arf; abs_value = false)

Internal method used by extrema_polynomial for polynomials of degree at most 2.

It assumes that check_interval(a, b) holds but does not check it.

The method is split into separate cases depending on the degree of p. It is optimized for performance and sacrifices readability.

IMPROVE: Add support for polynomials of degree 3.

_maximum_polynomial_low_degree(p::ArbPoly, a::Arf, b::Arf; abs_value = false)

Internal method used by minimum_polynomial for polynomials of degree at most 2.

Similar to _extrema_polynomial_low_degree but only computes minimum.

_midpoint_interval(a::T, b::T) where {T<:Union{Arf,Arb}}

Return the midpoint of a and b.

_midpoint_interval_log!(buffer::Arb, a::Arf, b::Arf)

Return the logarithmic midpoint of a and b as described in bisect_interval.

The argument buffer is used as scratch space during the computations.

_midpoint_interval_log(a::Arf, b::Arf)

Return the logarithmic midpoint of a and b as described in bisect_interval.

_minimum_polynomial_low_degree(p::ArbPoly, a::Arf, b::Arf; abs_value = false)

Internal method used by minimum_polynomial for polynomials of degree at most 2.

Similar to _extrema_polynomial_low_degree but only computes minimum.

bisect_interval(a::Arf, b::Arf; log_midpoint::Bool = false)

Compute the midpoint mid of a and b and return two tuples, (a, mid) and (mid, b), which corresponds to splitting the interval in half.

If log_midpoint = true then the midpoint is computed in logarithmic scale. If a, b > 0 this sets the midpoint to exp((log(a) + log(b)) / 2), can also be written as sqrt(a * b). If a, b < 0 we instead get -sqrt(a * b). If zero is contained in the interval then currently a normal bisection is performed (i.e. the midpoint is (a + b) / 2).

Logarithmic bisection can be useful if the function in consideration has a logarithmic behaviour and the interval has numbers exponentially close to zero, for example [1e-10000, 1e-1]. In that case normal bisection could give very slow convergence.

TODO: Think about how to handle a logarithmic bisection when the interval contains zero.

The value of mid is aliased in the two tuples and care should therefore be taken if doing inplace operations on it.

bisect_interval_recursive(a::T, b::T, depth::Integer; log_midpoint::Bool = false) where {T<:Union{Arf,Arb}}

Recursively bisect the interval [a, b] a number of depth times. The bisection is done using bisect_interval and it returns a vector with 2^depth intervals.

It accepts the same arguments as bisect_interval.

bisect_intervals(intervals::Vector{NTuple{2,T}}, to_bisect::Union{BitVector,Vector{Bool}}; log_midpoint = false) where {T<:Union{Arf,Arb}}

Bisect all intervals in intervals for which to_bisect is true and return a vector with the bisected intervals.

See also bisect_interval.

bounded_by(f, a::Arf, b::Arf, C::Arf; degree, abs_value, log_bisection, depth_start, maxevals, depth, threaded, verbose)

Return true if the function f can be shown to be bounded by C on the interval [a, b], i.e. f(x) <= C for all x ∈ [a, b], otherwise return false.

This function is similar to first computing the maximum with maximum_enclosure and then check if the computed maximum satisfies the bound. However if the only thing needed is to check if the bound holds this method has a number of benefits

• it aborts early if the bound is shown to not hold
• it doesn't try to compute an accurate enclosure of the maximum, it only bisects as much as is needed for getting the bound.
• the implementation is simpler and easier to check for correctness.

The maximum of f is enclosed by the use of Taylor series through maximum_series. The degree of the expansion is kept constant and the interval is bisected until we can either conclude that the bound holds on the whole interval or there is subinterval where it doesn't hold.

The function f should support evaluation on both Arb and ArbSeries and should return an enclosure of the result in both cases. The degree used for ArbSeries can be set with the degree argument (defaults to 8). If degree is negative then it will fall back to direct evaluation with Arb and not make use of maximum_series, this is usually much slower but does not require the function to be implemented for ArbSeries.

If abs_value = true the instead consider the function abs(f(x)) on the interval [a, b].

If log_bisection = true then the intervals are bisected in a logarithmic scale, see bisect_interval for details.

The argument depth_start bisect the interval using bisect_interval_recursive before starting to compute the extrema. This can be useful if it is known beforehand that a certain number of bisections will be necessary before the enclosures get good enough. It defaults to 0 which corresponds to not bisecting the interval at all before starting.

The arguments maxevals and depth can be used to limit the number of function evaluations and the number of bisections of the interval respectively. Notice that depth takes depth_start into account, so the maximum number of iterations is depth - depth_start.

If threaded = true then evaluate f in parallel on the intervals using Threads.@threads.

If verbose = true then output information about the process.

check_interval(a::Arf, b::Arf)

Throw an error if a and b does not define a finite interval [a, b]. It checks that a and b are both finite and satisfy a <= b.

check_interval(Bool, a::Arf, b::Arf)

Return true if a and b defines a finite interval [a, b]. It checks that a and b are both finite and satisfy a <= b.

check_tolerance(x::Arb; atol = nothing, rtol = nothing)
check_tolerance(x::AbstractVector{Arb}; atol = nothing, rtol = nothing)

Return true if x satisfies the absolute tolerance atol or the relative tolerance rtol.

For x::AbstractVector{Arb} the tolerance is checked element-wise.

The absolute tolerance is satisfied if the diameter of x is less than or equal to atol.

The relative tolerance is satisfied if the diameter of x divided by the absolute value of x is less than or equal to rtol. If x contains zero this is only satisfied if x is exactly zero.

Both atol and rtol can be given as nothing, in which case that check is skipped. If both of them are nothing then it returns true if x is finite.

It always returns true if x is finite and the radius is zero.

compose_zero!(res::T, p::T, q::T) where {T<:Union{ArbPoly,AcbPoly,ArbSeries,AcbSeries}}

In-place version of compose_zero.

compose_zero(p::T, q::T) where {T<:Union{ArbPoly,AcbPoly,ArbSeries,AcbSeries}}

Compute compose(p, q), but with the constant term of q set to zero.

derivative_function(f, n = 1)

Return a function for computing the nth derivative of f.

The returned function accepts only Arb, Acb, ArbSeries and AcbSeries as input. For Arb and ArbSeries it calls f with ArbSeries. For Acb and AcbSeries it calls f with AcbSeries.

If f is a polynomial, then return the derivative of the polynomial.

IMPROVE: Avoid overflow in factorial function for large n.

enclosure_lbound(x::Arb)

Compute an enclosure of the lower and upper endpoints of x.

See enclosure_ubound.

enclosure_lbound(x::Arb)

Compute an enclosure of the lower endpoint of x.

See enclosure_ubound.

enclosure_series(f, x::Arb; degree = 0, abs_value = false, verbose = false)

Convenience function for computing an enclosure of f(x) using extrema_series.

The enclosure is computed by finding the minimum and maximum value using extrema_series and typically gives a tighter enclosures than just evaluating f(x) directly.

It is equivalent to

Arb(ArbExtras.extrema_series(f, getinterval(x)...; degree, abs_value, verbose)[1:2])

but has the default value degree = 0, which is different than extrema_series. The reason for the default degree = 0 is that it can pick up on monotonicity of f and is therefore enough in many cases.

enclosure_ubound(x::Arb)

Compute an enclosure of the upper endpoint of x.

The function ubound can be used to compute an upper bound of x but it doesn't give guarantee on how much larger than the upper endpoint of x this bound is.

One common use of this method is to compute enclosures of functions which have been proved to be monotone on x. Consider a function f and say we have checked that the enclosure of its derivative is positive on x. In that case f(enclosure_ubound(x)) gives an enclosure of the upper bound of f on x. In contrast f(ubound(Arb, x)) is not even guaranteed to give an upper bound of f on x. This is because ubound(Arb, x) can be larger than x and hence f has not been proved to be increasing there.

As an example the following code produces an x which satisfy x < π but for which ubound(x) > π.

x = Arb(3, prec = 16)
r = Mag(lbound(Arb(π, prec = 16) - x))
Arblib.set!(Arblib.radref(x), r)

x < Arb(π, prec = 256)

xᵤ = ubound(Arb, x)

Arb(π, prec = 256) < Arb(xᵤ)

extrema_enclosure(f, a::Arf, b::Arf; degree, atol, rtol, lbound_tol, ubound_tol, abs_value, log_bisection, point_value_min, point_value_max, depth_start, maxevals, depth, threaded, verbose)

Compute both the minimum and maximum of the function f on the interval [a, b] and return them as a 2-tuple.

The extrema are computed by enclosing the extrema of the function on an interval using Taylor series through extrema_series. The degree of the expansion is kept constant and the interval is bisected until the required tolerance is met.

The function f should support evaluation on both Arb and ArbSeries and should return an enclosure of the result in both cases. The degree used for ArbSeries can be set with the degree argument (defaults to 8). If degree is negative then it will fall back to direct evaluation with Arb and not make use of extrema_series, this is usually much slower but does not require the function to be implemented for ArbSeries.

The atol and rtol argument are used to set the tolerance for when a subinterval is considered done and not bisected further, see check_tolerance for details. The lbound_tol and ubound_tol arguments can optionally be set to give another criteria for when to stop bisection of subintervals. For the minimum a subinterval is considered done if the enclosure of its minimum is strictly greater than lbound_tol, for the maximum if its maximum is strictly less than ubound_tol. This can be useful if you mostly care about whether the minimum/maximum is less/greater than some specific value, then you don't need to bisect if the current enclosure is good enough. See also bounded_by if you only need to prove that f is bounded by some value.

If abs_value = true the compute the extrema of abs(f(x)) on the interval [a, b]. For the computation of the maximum this is mostly the same, only difference is that we have to take the absolute value of the evaluations. For the minimum we have to take into account that the polynomial might cross zero, in which case the minimum is zero.

If log_bisection = true then the intervals are bisected in a logarithmic scale, see bisect_interval for details.

The arguments point_value_min and point_value_max can optionally be set to a a priori upper or lower bounds of the min and max respectively. Typically this is computed by evaluating f on some (possibly well chosen) points before calling this method and taking the minimum and maximum of the evaluations respectively. This can allow the method to quicker discard subintervals where the extrema could not possibly be located.

The argument depth_start bisect the interval using bisect_interval_recursive before starting to compute the extrema. This can be useful if it is known beforehand that a certain number of bisections will be necessary before the enclosures get good enough. It defaults to 0 which corresponds to not bisecting the interval at all before starting.

The arguments maxevals and depth can be used to limit the number of function evaluations and the number of bisections of the interval respectively. Notice that depth takes depth_start into account, so the maximum number of iterations is depth - depth_start.

If threaded = true then evaluate f in parallel on the intervals using Threads.@threads.

If verbose = true then output information about the process.

TODO: Currently this always computes both minimum and maximum on all subintervals. Change it so that it takes into account if the minimum/maximum could be located in the subinterval or not.

extrema_polynomial(p::ArbPoly, a::Arf, b::Arf; abs_value = false, verbose = false)

Compute an enclosure of both the minimum and maximum of the polynomial p on the interval [a, b] and return them as a 2-tuple.

The extrema are computed by finding the zeros of the derivative of p using isolate_roots and refining them using refine_root, then evaluating p on the potential zeros as well as the endpoints a, b of the interval.

If p cannot be evaluated to a reasonable precision, because its coefficients are too wide, it does in general not make sense to try and isolate the roots. To check for this we compute min(radius(p(a)), radius(p(b))) and p(Arb((a, b))). If the quotient is larger than 0.99, meaning we don't even reduce the radius by one percent by evaluating at the endpoints instead of the whole interval, we don't try to isolate the roots but just evaluate the polynomial directly.

If abs_value = true the compute the extrema of abs(p(x)) on the interval [a, b]. For the computation of the maximum this is mostly the same, only difference is that we have to take the absolute value of the evaluations. For the minimum we have to take into account that the polynomial might cross zero, in which case the minimum is zero. To see if p crosses zero we compare the sign of p at the endpoints and all zeros of its derivative. If any of these have different signs we are guaranteed that p crosses zero and the minimum is zero, if all of them have the same sign we are instead guaranteed that the polynomial doesn't cross zero.

If verbose = true then output information about the process.

TODO: Allow passing arguments to isolate_roots?

extrema_series(f, a::Arf, b::Arf; degree, abs_value, verbose)

Compute both the minimum and maximum of the function f on the interval [a, b]. It returns a three tuple (fmin, fmax, fmid) where fmin and fmax are enclosures of the minimum and maximum respectively and fmid is an enclosure of f evaluated on the midpoint of the interval.

The extrema are computed by computing the extrema of the Taylor series of f on the midpoint of the interval and then adding an error term. The degree of the Taylor series used is set with the degree argument.

The reason that fmid is returned is that as part of computing the Taylor series we get this value for free. In some cases the Taylor series for at the midpoint is however not computed, in this case fmid will just be NaN instead.

If abs_value = true then compute the extrema of abs(f(x)) on the interval [a, b]. For the computation of the maximum this is mostly the same, only difference is that we have to compute the maximum of the absolute value of the Taylor series. For the minimum we have to take into account that the polynomial might cross zero, in which case the minimum is zero.

If verbose = true then output information about the process.

This method is mainly intended for internal use in the function extrema_enclosure but can be used directly as well.

TODO: Handle crossing zero when computing absolute minimum better. Currently we always add the rest term in the end.

format_interval(a, b; digits = 5)

Return a nicely formatted string representing the interval [a, b].

It prints the interval either on the form [a, b] or on the form [c +/- d] depending on the width of the width of the interval.

If a or b is non-finite it always prints an interval of the form [a, b]. Otherwise it checks if the relative error is less than 10^-digits, if it is then it prints it as a ball, otherwise as an interval. When printed in interval form it prints at most digits digits.

integrate(f, a::Arb, b::Arb; atol, rtol, depth_start, maxevals, depth, verbose)

Compute an enclosure of the integral of f on the interval[a, b]`.

It uses a Gauss-Legendre quadrature of order 2 through integrate_gauss_legendre together with bisection of the interval. On each subinterval the integral is computed using integrate_gauss_legendre and it checks if the result satisfies the required tolerance, if it doesn't the interval is bisected. For more details about the integration method see integrate_gauss_legendre.

Notice that the tolerance is checked on each subinterval and not on the interval as a whole. This means that even if it finishes before reaching the maximum depth or the maximum number of evaluations the result as a whole might not satisfy the given tolerance.

The argument depth_start bisect the interval using bisect_interval_recursive before starting to compute the integral. This can be useful if it is known beforehand that a certain number of bisections will be necessary before the enclosures get good enough. It defaults to 0 which corresponds to not bisecting the interval at all before starting.

The arguments maxevals and depth can be used to limit the number of function evaluations and the number of bisections of the interval respectively.

If verbose = true then output information about the process.

integrate_gauss_legendre(f, a::Arb, b::Arb)

Compute an enclosure of the integral of f on the interval [a, b] using a Gauss-Legendre quadrature of order 2.

The approximation is given by

(b - a) / 2 * (f(x₁) + f(x₂))

where

x₁ = (b - a) / 2 * sqrt(3) / 3 + (b + a) / 2
x₂ = -(b - a) / 2 * sqrt(3) / 3 + (b + a) / 2

are quadrature nodes. The remainder term is enclosed by

(b - a)^5 / 4320 * d4f

where d4f is an enclosure of the fourth derivative of f on the interval [a, b].

If it happens to be that d4f is not finite, for example because f is not differentiable on the interval, it instead computes a naive enclosure of the integral using the enclosure of f evaluated on the interval [a, b].

The function f should support evaluation on both Arb and ArbSeries and should return an enclosure of the result in both cases. For ArbSeries it is enough that it supports evaluation up to degree 4 and only the zeroth and fourth order term are ever used.

• TODO: Optimize for performance
• TODO: Add option for computing integral of absolute value?
• TODO: Allow giving an separate function for computing derivative?

iscpx(p::Union{ArbPoly,AcbPoly,ArbSeries,AcbSeries})

Return true if p is of the form c + x.

isolate_roots(f, a::Arf, b::Arf; depth = 10, check_unique = true, verbose = false)

Isolate the roots of f on the interval [a, b].

Returns a tuple (found, flags) where found is a vector of subintervals of [a, b] given by pairs of Arf and flags is a vector of booleans. The output has the following properties:

• The function has no roots on interval outside of the subintervals in found.
• Subintervals are sorted in increasing order (with no overlap except possibly starting and ending with the same point).
• Subintervals with a flag true has exactly one (simple) root.
• Subintervals with a flag false may or may not contain roots.

If all flags are true then all roots of the function on interval have been isolated. If there are output subintervals on which the existence or nonexistence of roots could not be determined, the user may attempt further searches on those subintervals (possibly with increased precision and/or increased bounds for the breaking criteria). Note that roots of multiplicity higher than one and roots located exactly at endpoints cannot be isolated by the algorithm.

This method closely mimics the function arb_calc_isolate_roots in Arb, both in behaviour and in implementation.

The function f should support evaluation on both Arb and ArbSeries and should return an enclosure of the result in both cases. The evaluation on ArbSeries is done to compute the derivative to check uniqueness, this can be skipped by setting check_unique = false, see below. Alternatively f can be an ArbPoly, in which case some optimizations are done.

If check_unique = false then don't check if roots are unique. In this case flags will always be false and it will simply continue to split the subintervals until it reaches the maximum depth. The only benefit to this is that the derivative of f doesn't have to be computed so it can be used when the function of interest doesn't support evaluating derivatives.

If verbose = true then output information after each iteration on the number of found roots and the number of remaining intervals that needs to be split further.

FUTURE WORK:

• Parallelize.
• Compute values on endpoints of subintervals once to not have to compute them multiple times.
• Should it use BitVector or Vector{Bool} for flags?
• The current implementation does a breadth-first-search. Consider rewriting it to do a depth-first-search instead to reduce memory usage in case of many splittings. This might however make it harder to parallelize.

maximum_enclosure(f, a::Arf, b::Arf; degree, atol, rtol, ubound_tol, abs_value, log_bisection, point_value_max, depth_start, maxevals, depth, threaded, verbose)

Compute the maximum of the function f on the interval [a, b].

Takes the same arguments as extrema_polynomial. The algorithm is also the same except that it only looks for the maximum.

maximum_polynomial(p::ArbPoly, a::Arf, b::Arf; abs_value = false, verbose = false)

Compute an enclosure of the maximum of the polynomial p on the interval [a, b].

Takes the same arguments as extrema_polynomial. The algorithm is also the same except that it only looks for the maximum.

maximum_series(f, a::Arf, b::Arf; degree, abs_value, verbose)

Compute the maximum of the function f on the interval [a, b].

Takes the same arguments as extrema_series. The algorithm is also the same except that it only looks for the maximum.

This method is mainly intended for internal use in the function maximum_enclosure but can be used directly as well.

minimum_enclosure(f, a::Arf, b::Arf; degree, atol, rtol, lbound_tol, abs_value, log_bisection, point_value_min, depth_start, maxevals, depth, threaded, verbose)

Compute the minimum of the function f on the interval [a, b].

Takes the same arguments as extrema_polynomial. The algorithm is also the same except that it only looks for the minimum.

minimum_polynomial(p::ArbPoly, a::Arf, b::Arf; abs_value = false, verbose = false)

Compute an enclosure of the minimum of the polynomial p on the interval [a, b].

Takes the same arguments as extrema_polynomial. The algorithm is also the same except that it only looks for the minimum.

minimum_series(f, a::Arf, b::Arf; degree, abs_value, verbose)

Compute the minimum of the function f on the interval [a, b].

Takes the same arguments as extrema_series. The algorithm is also the same except that it only looks for the minimum.

This method is mainly intended for internal use in the function minimum_enclosure but can be used directly as well.

refine_root(f, root::Arb; df, atol, rtol, min_iterations, max_iterations, strict)
refine_root(f::ArbPoly, root::Arb; df, atol, rtol, min_iterations, max_iterations, strict)

Refine the given root of the function f.

Given a ball root containing a unique root of the function f it refines the enclosure of the root using successive interval Newton iterations. The method can only handle simple roots contained in the interior of the enclosure and to work the derivative must be non-zero on the whole enclosure.

The function f can be either an ArbPoly or a regular function. By default the derivative is computed by differentiating the polynomial, in the case of ArbPoly, or using ArbSeries, in the case of a regular function. Alternatively the keyword argument df can be set to a function for computing the derivative.

At each iteration it checks if the required tolerance is met according to check_tolerance. The default tolerances are atol = 0 and rtol = 4eps(one(root)).

If the absolute error does not improve by at least a factor 1.5 between two iterations then it stops early since it's unlikely that further iterations would improve the result much. This can for example happen when the function is computed with too low precision. This check is only done if min_iterations have been performed, this is to avoid stopping early due to slow convergence in the beginning. To avoid this check entirely min_iterations can be put to the same as max_iterations.

If strict = true then only return an enclosure which has been proved to contain a unique root contained in the original enclosure, if it could not prove that this is the case it returns NaN. If strict = false it will return an enclosure in all cases, this enclosure is only valid if root indeed contains a unique root and the returned enclosure could just be the starting enclosure root itself. This is useful if you know that root contains a unique root and don't want to worry about this method failing.

If any iteration gives NaN as a result then it returns the enclosure from the previous iteration, possibly just the starting enclosure. This is unless strict = true and it has not been able to prove that there is a root, in which case NaN is returned.

If verbose = true then print the enclosure at each iteration and some more information in the end.

In rare cases when strict = false and the original enclosure does in fact not contain a root the method could converge to a potential root just outside the original enclosure. The reason this might happen is that when the intersection of two balls is computed the enclosure is typically slightly larger than the true intersection. This is not an issue if strict = false nor if the original enclosure does contain a root.

refine_root_bisection(f, root::Arb; atol, rtol, max_iterations, strict, verbose)

Refine the root of the function f on the interval [a, b] using bisection.

The function f is assumed to be continuous on the interval [a, b] and f(a) and f(b) should have different signs.

It iteratively bisects the interval and determines which subinterval to keep by checking the sign of f at the midpoint. It stops either when the required tolerance is met according to check_tolerance, it has reached the maximum number of iterations, it can't determine the sign of the midpoint or the midpoint is equal to one of the endpoints (meaning that we have bisected as much as possible at the given precision).

If the zero lies exactly at or very close to one of the bisection point it can fail to determine the sign there. We avoid this by in that case instead trying to bisect at a point between the left endpoint of the interval and the midpoint, if that also fails it stops.

When possible refine_root should be used instead since it converges much faster. This method is useful when the derivative is not available, i.e. evaluation with ArbSeries is not possible, or as a preprocessor for refine_root if the original enclosure of the root gives an enclosure on the derivative which contains zero.

If strict = true return NaN if the signs of f(a) and f(b) don't differ. In this case a finite result proves that there is a root of f on the interval, assuming that f is continuous there, but it says nothing about the uniqueness.

If verbose = true then print the enclosure at each iteration and some more information in the end.

• IMPROVE: Reduce the number of allocations.
• IMPROVE: Handle the case when the sign at the midpoint can't be determined better by using a perturbation of the midpoint.

taylor_remainder(p::ArbSeries, x::Arb)

Compute p[degree] * (x - midpoint(x))^degree, where degree is the degree of p

This corresponds to the Lagrange form of the remainder term in a Taylor expansion around the point midpoint(x) of degree degree - 1 valid on the full interval x.

For example, for a function f, a ball x, any degree >= 0 and

p = f(ArbSeries((x, 1), degree = degree + 1))
q = f(ArbSeries((midpoint(x), 1), degree = degree))

we have for any y ∈ x that

f(y) ∈ q(y) + taylor_remainder(p, x)