Barycentric interpolation formula with interval arithmetic

Barycentric interpolation is a numerically stable algorithm to evaluate Lagrange interpolating polynomials. An implementation in Julia is available at BarycentricInterpolation.jl.

The procedure takes as input a set of nodes X_k = \{ x_1, x_2, \ldots x_k\} and initially computes a set of weights W_k = \{ w_1, w_2, \ldots w_k \}.

Then, evaluating the k polynomials at the point x involves:

(This is pseudocode adopted from the above mentioned package just for description purposes)

for i = 1:k
      vals[i] = Wn[i] / (x - Xn[i])

(The complete algorithm checks for the case x == Xn[i] and then sets vals[i] = 1.0 and all other entries of vals to 0.0. Also, there is a final normalization step.)

I was wondering if there is a “Barycentric algorithm” that deals with the case when x::Interval (from IntervalArithmetic.jl).

Specifically, if x contains Xn[i] then the output bounds become arbitrarily loose [0,+/-∞]. Any inputs in dealing with this will be very helpful!

(It’s a little amusingly ironic that the barycentric algorithm is used because it is numerically robust and avoids cancellation errors, yet seems to not work immediately with interval arithmetic.)