On the numerical stability of linear barycentric rational interpolation

The barycentric forms of polynomial and rational interpolation have recently gained popularity, because they can be computed with simple, efficient, and numerically stable algorithms. In this paper, we show more generally that the evaluation of any function that can be expressed as r(x)=∑i=0nai(x)fi/∑j=0mbj(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r(x)=\sum _{i=0}^n a_i(x) f_i\big /\sum _{j=0}^m b_j(x)$$\end{document} in terms of data values fi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_i$$\end{document} and some functions ai\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_i$$\end{document} and bj\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j$$\end{document} for i=0,…,n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=0,\ldots ,n$$\end{document} and j=0,⋯,m\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=0,\dots ,m$$\end{document} with a simple algorithm that first sums up the terms in the numerator and the denominator, followed by a final division, is forward and backward stable under certain assumptions. This result includes the two barycentric forms of rational interpolation as special cases. Our analysis further reveals that the stability of the second barycentric form depends on the Lebesgue constant associated with the interpolation nodes, which typically grows with n, whereas the stability of the first barycentric form depends on a similar, but different quantity, that can be bounded in terms of the mesh ratio, regardless of n. We support our theoretical results with numerical experiments.


Introduction
Given n + 1 distinct interpolation nodes x 0 , . . . , x n with associated data f 0 , . . . , f n , the classical interpolation problem consists in finding a function r : R → R that interpolates f i at x i , that is, In the case of polynomial interpolation of degree at most n, this problem has a unique solution, which can be expressed in Lagrange form as While this form is advantageous for theoretical analysis, its evaluation requires O(n 2 ) operations and can be numerically unstable. It is advisable to consider instead the first barycentric form of r , with the Lagrange weights w i defined as The first barycentric form is more efficient than the Lagrange form, as it can be evaluated in O(n) operations, after computing the w i , which are independent of x, in O(n 2 ) operations in a preprocessing step. Higham [1] further shows that this evaluation is backward stable with respect to perturbations of the data f i . Another means of evaluating r is given by the second barycentric form of r , which can be derived from (2) by noticing that 1 = (x) n i=0 w i x−x i . Evaluating this formula also requires O(n) operations, but it comes with the advantage that the w i can be rescaled by a common factor to avoid underflow and overflow [2]. Moreover, the second barycentric form is forward stable, as long as the Lebesgue constant associated with the interpolation nodes x i is small [1], which is the case, for example, for Chebyshev nodes of the first and the second kind, but not for equidistant nodes.
In the case of rational interpolation, the interpolation problem (1) no longer has a unique solution, but Berrut and Mittelmann [3] show that every rational interpolant of degree at most n can be expressed in the second barycentric form (4) for a specific choice of weights w i . Vice versa, Schneider and Werner [4] note that for any set of non-zero weights w i , the function r in (4) is a rational interpolant of degree at most n. An important subset of these barycentric rational interpolants are those that do not have any poles in R. This is obviously true for the Lagrange weights in (3), but also for the Berrut weights [5] and, assuming the interpolation nodes to be in ascending order, that is, x 0 < · · · < x n , for the family of weights with parameter d ∈ {0, . . . , n}, proposed by Floater and Hormann [6]. Note that this family includes the Berrut and the Lagrange weights as special cases for d = 0 and d = n. For the weights in (6), Floater and Hormann further observe that the barycentric rational interpolant in (4) with the weights (6) can also be written as where As the formula in (7) simplifies to (2) for d = n, we refer to it as the first barycentric form of the corresponding rational interpolant. Note that the first and the second form are identical for Berrut's interpolant, that is, for d = 0. The result of Higham [1] can be extended to show that the evaluation of the second barycentric form (4) is forward stable, not only in the case of polynomial interpolation, but for general barycentric rational interpolants, provided that the weights w i can be computed with a forward stable algorithm and that the corresponding Lebesgue constant is small [7]. For Berrut's rational interpolant with weights in (5), this is the case for all well-spaced interpolation nodes [8], including equidistant and Chebyshev nodes. For the family of barycentric rational interpolants with weights in (6), the Lebesgue constant is known to grow logarithmically in n, for any constant d > 0 and equidistant [9] as well as quasi-equidistant [10] nodes, and the formula in (6) turns out to be a forward stable means of computing the weights w i [11].
In this paper, we further generalize the proof of Salazar Celis [7], such that it can also be used for proving the forward stability of the first barycentric form (7), with two important changes (Sect. 3). On the one hand, the result relies on the fact that not only the w i , but also the λ i can be evaluated with a forward stable algorithm. This, however is indeed the case, both for the formula in (8), which requires O(d(n − d)) operations for 0 < d < n, as well as a more efficient, but slightly less stable formula that gets by with O(n) operations (Sect. 4). On the other hand, the Lebesgue constant must be replaced by a similar quantity that depends on the functions λ i , for which we prove that it is at most on the order of O(μ d ), where μ is the mesh ratio of the interpolation nodes (Sect. 6). Moreover, we show that a more efficient formula [12] for computing the weights in (6) is forward stable, too (Sect. 4).
Regarding backward stability, Mascarenhas and Camargo [13] show that the second barycentric form is backward stable under the same assumptions, namely forward stable weights w i and small Lebesgue constant. Moreover, Camargo [11] proves that the first barycentric form is backward stable, as long as the denominator in (7) is computed with a special algorithm in O(nd) operations.
In this paper, we derive a substantially different approach that can be used to show the backward stability of both barycentric forms (Sect. 5). For the second barycentric form, our result provides an upper bound on the perturbation of the data f i that is smaller than the upper bound by Mascarenhas and Camargo [13]. For the first barycentric form, our upper bound is larger than the one found by Camargo [11], but it comes with the advantage of holding for a more efficient way of computing the denominator in (7) in O(n) operations, which is based on our new O(n) algorithm for evaluating the λ i (Sect. 4).
It is important to note that our results hold under the assumption that the input values to the algorithm, x i , f i , and x, are given as floating point numbers, so they do not introduce any additional error when we compute both forms (4) and (7). Our stability analysis does not cover errors that result from initially rounding the given values to floating point numbers.
Our numerical experiments (Sect. 7) confirm that this leads to an evaluation of the first barycentric form (7) of the rational interpolant with weights in (6), which is considerably faster than the algorithm proposed by Camargo [11], especially for larger d, at the price of only marginally larger forward errors. Evaluating the interpolant using the second barycentric form (4) is even faster and can be as stable, but may also result in significantly larger errors for certain choices of interpolation nodes. However, we also report a case in which the second form is stable and the first is not.

Preliminaries
Following Trefethen [14], a general problem g can be seen as a function g : U → V from a normed vector space (U , · U ) of data to a normed vector space (V , · V ) of solutions and a numerical algorithm for solving the problem on a computer as another functionĝ : U → V that approximates g.
For our analysis, we consider a computer that uses a set F of floating point numbers with the corresponding machine epsilon and let fl : R → F be the rounding function that maps each x ∈ R to the closest floating point approximation fl(x) ∈ F. Moreover, we denote by the floating point analogue of the arithmetic operation * ∈ {+, −, ×, ÷}, that is, for x, y ∈ F, and recall [14,Lecture 13] that fl and introduce a relative error of size at most . That is, for all x ∈ R there exists some δ ∈ R with |δ| < , such that and likewise for all x, y ∈ F, which follows immediately from (9). We further note that fl(−x) = −x for all x ∈ F, so that multiplying a floating point number by (−1) i or taking its absolute value does not entail any rounding error. In order to analyse and to compare the quality of numerical algorithms, the usual quantities to consider are the absolute forward error ĝ(u) − g(u) V and the rela- where the notation x = O( ) means that there exists some positive constant C, such that |x| ≤ C as → 0. In the case of forward stability, this constant usually depends on the relative condition number of the problem g at u and on the dimensions of U and V . In this paper, we analyse algorithms for evaluating barycentric rational interpolants, based on the formulas (4) and (7). Following Higham [1], we assume the evaluation point x ∈ F and the interpolation nodes x = (x 0 , . . . , x n ) ∈ F n+1 to be fixed and prove the forward and backward stability with respect to the data f = ( f 0 , . . . , f n ) ∈ F n+1 . More precisely, we solve the problem by a numerical algorithm and we show that for any choice of x, x, and f , the relative forward error is bounded from above (Sect. 3), and that there exists some perturbed dataf with where the additional arguments of r are used to emphasize the dependence of the interpolant on the interpolation nodes and the data (Sect. 5). In both cases, is assumed to be sufficiently small, since this is enough to guarantee that the errors are O( ). We will see that both constants C 1 and C 2 depend linearly on n and a quantity β(x), which is independent of f and different for the first and the second barycentric form. It turns out that this quantity is usually smaller when the first barycentric form is used to evaluate r (Sect. 7). As mentioned above, the constant C 1 further depends on the condition number of the problem, and more specifically on the componentwise relative condition number [15] κ(x; For the derivation of these upper bounds, we shall frequently use the following basic facts: 1. By Taylor expansion, Moreover, for any δ ∈ R with |δ| ≤ , there exists some δ ∈ R with |δ | ≤ + O( 2 ), such that 1 This observation is useful for "moving" the perturbation (10) caused by a floating point operation from the denominator to the numerator.

For any
We use this observation to gather the perturbations caused by computing the product of m terms into a single perturbation 1 This follows from the previous observation, and we use it to estimate the rounding error introduced by simple iterative summation of m + 1 floating point numbers.

Forward stability
For analysing the relative forward error of barycentric rational interpolation, we first observe that (4) and (7) can both be written in the common form where a i (x) = w i /(x − x i ) for both forms, while m = n and b j (x) = a j (x) for the second form and m = n − d and b j (x) = λ j (x) for the first form. Next, we define the functions and Assuming now that we have forward stable algorithms for computing a i (x) asâ i (x) and b j (x) asb j (x), we can derive a general bound on the relative forward error for the function r in (16).
for some constants A and B. Then, assuming f ∈ F n+1 , the relative forward error of r in (16) satisfies for small enough.
Proof We first notice thatr is given bŷ and δ ÷ are the relative errors introduced respectively by the product a i (x) f i , the sums in the numerator and the denominator, and the final division. It then follows from (10) such thatr Therefore,r Since, by the triangle inequality, we can use (12) to express this ratio aŝ and obtain the relative forward error of r as The upper bound in (19) then follows immediately by using again (22) and (23).
While Theorem 1 holds for any function r that can be expressed as in (16), we shall now focus on the special cases for the two different forms of the barycentric rational interpolant. For the second barycentric form, the only assumption we need is that the weights w i can be computed asŵ i with a forward stable algorithm. Moreover, we recall the definition of the Lebesgue function of the barycentric rational interpolant [9] as for some constant W . Then, assuming f ∈ F n+1 , the relative forward error of the second barycentric form in (4) satisfies for small enough.
Proof We first notice that a i (x) = w i /(x − x i ) can be computed with one subtraction and one division, so that, by (10), (13), (14), and (24), in case of the second barycentric form. From the latter, it also follows immediately that β(x) in (18) is equal to n (x; x) in this case, and it only remains to show that α(x; f ) in (17) is equal to κ(x; x, f ) in (11). To this end, we first use the triangle inequality to see that for any δ > 0 and any Dividing both sides of this inequality by H and n i=0 a i (x) f i , taking the supremum over all admissible h and the limit δ → 0, gives κ(x; For the first barycentric form, we additionally need to assume that the λ i (x) can be computed with a forward stable algorithm asλ i (x), and we note that β(x) in (18) is equal to . (26) in this case.

Corollary 2
Assume that the weights w 0 , . . . , w n can be computed as in (24) and that there exist for some constant C. Then, assuming f ∈ F n+1 , the relative forward error of the first barycentric form in (7) satisfies for small enough.
Proof As the numerator of the first and second barycentric form are identical, the only difference to the proof of Corollary 1 is that While upper bounds for the Lebesgue function n (x; x) can be found in the literature [8,9], we are unaware of any previous work bounding the function d (x; x), and we derive such an upper bound in Sect. 6. the error propagation that occurs in the implementation of different algorithms and further analyse them in terms of computational cost.
Regarding the w i , it was shown by Higham [1] that the Lagrange weights in (3) can be computed stably with W = 2n in (24), and the Berrut weights in (5) can be represented exactly in F, so that W = 0. The same holds for the weights in (6) if the interpolation nodes are equidistant, because they simplify to the integers in this special case [6]. For the general case, Camargo [11,Lemma 1] shows that W = 3d, if the w i are computed with a straightforward implementation of the formula in (6). While this construction requires O(nd 2 ) operations, Hormann and Schaefer [12] suggest an improved O(nd) pyramid algorithm, which turns out to have the same precision. Their algorithm starts from the values and iteratively computes for l = d − 1, d − 2, . . . , 0, tacitly assuming v l i = 0 for i < 0 and i > n − l. They show that the resulting values v 0 i are essentially the weights w i in (6), up to a factor of (−1) i−d .

and the intermediate value theorem further guarantees that
Let us now focus on the functions λ i that appear in the barycentric formula (7) and first study the error propagation when computing them straightforwardly, following the formula in (8).
Evaluating λ i in this way clearly has a computational cost of O(d) for d > 0 and O(1) for d = 0, so that computing all λ i (x) requires O(d(n − d)) operations for 0 < d < n and O(n) operations for d = 0 and d = n. However, for 0 < d < n this can be improved by exploiting the fact that λ i (x) and λ i+1 (x) have d common factors in the denominator, which in turn suggests to first compute the "central" λ m (x) for m = n−d 2 as above and then the remaining λ i (x) iteratively as Computing all λ i (x) this way requires only O(n) operations, but it comes at the price of a likely reduced precision.

Backward stability
Similar to how we established the forward stability of both barycentric forms in a unified way in Sect. 3, we can prove the backward stability in general for the function r in (16) and then derive upper bounds on the perturbation of the data for both forms as special cases.

. , m for some constants A and B. Then there exists for any
for small enough, such that the numerical evaluation of r in (16) satisfiesr (x; f ) = r (x;f ).
Proof Starting from (21), with the η i and μ j satisfying (20), we get, again with the help of (12), By (23), this means that there exist some ξ 0 , . . . , ξ n ∈ R with |ξ 0 |, . . . , and by (14) and (20) we can further assume the existence of some ϕ 0 , . . . , ϕ n with |ϕ 0 |, . . . , |ϕ n | ≤ (n , such that The statement then follows directly from The special cases of Theorem 2 for the two different forms of the barycentric rational interpolant then follow with the same reasoning as in Sect. 3.

Corollary 3
Assume that the weights w 0 , . . . , w n can be computed as in (24). Then there exists for any f ∈ F n+1 somef ∈ F n+1 with for small enough, such that the second barycentric form in (4) for small enough, such that the first barycentric form in (7) satisfiesr (x; x, f ) = r (x; x,f ).

Upper bound for 0 d
In case of the first barycentric form, the bounds for the forward and backward stability depend on the function d in (26), and we still need to show that this function is bounded from above. Note that 0 = n , because w i = (−1) i and λ i = , so that the bound for the Lebesgue constant of Berrut's interpolant [8] also holds for d in this case. In the following, we therefore assume d ≥ 1 and define . We further assume that x ∈ (x k , x k+1 ) for some k ∈ {0, . . . , n − 1}. It then follows from the definition in (8) that all λ i with index in where I = {0, . . . , n − d}, have the same sign as (−1) k+d and that the sign alternates for decreasing indices "to the left" and increasing indices "to the right" of I 3 . More precisely, the λ i with index in have the same sign as the ones with index in I 3 , while the sign is opposite, if the index is in Without loss of generality, we assume that the λ i are positive for i ∈ I 2 , I 3 , I 4 and negative for i ∈ I 1 , I 5 , since multiplying all λ i with a common constant does not change d . Letting we then find that To bound d , we need to bound the sum of the negative λ i with i ∈ I 1 , I 5 as well as the λ i with i ∈ I 3 , which in turn requires to first bound the terms x − x i . To this end, Proof The statement follows directly by noting that for the given t and because h min ≤ h i ≤ h max for any i ∈ {0, . . . , n − 1}.
For bounding the negative λ i from above, it turns out to be useful to consider them in pairs, with indices from I 1 and I 5 at the same "distance" from I 3 .

Lemma 4 For any j ∈ N and x
where we set λ i (x) = 0 for any i / ∈ I .
As g is clearly even and , it remains to show that g(t) ≤ g (1) for t ∈ [0, 1]. To this end, note that g(t) = p(t)/q(t) for and By the product rule, and For t ∈ [0, 1], we observe that p(t) > 0, q(t) > 0, p (t) ≥ 0, and q (t) ≤ 0, hence and it follows from the quotient rule that g is monotonically increasing over [0, 1].
Next, let us bound the λ i with indices in I 3 from below.

Lemma 5
For any i ∈ I 3 and x ∈ (x k , x k+1 ), Proof Since the denominator of λ i (x), for i ∈ I 3 , contains the factors x − x k−m for m = 0, . . . , k − i and the factors x k+1+l − x for l = 0, . . . , i + d − k − 1, we conclude from Proposition 1 that and the statement then follows, because a!b! ≤ (a + b − 1)! for any a, b ∈ N.
We are now ready to derive a general upper bound on the function d in (26), which turns out to depend on d and the mesh ratio for any set of ascending interpolation nodes x = (x 0 , . . . , x n ) ∈ R n+1 and any x ∈ [x 0 , x n ].
Proof If x = x k for some k ∈ {0, . . . , n}, then, after multiplying both N (x) and D(x) by n i=0 |x − x i |, we find that d (x; x) = 1, which is clearly smaller than the upper bound in (34). Otherwise, there exists some k ∈ {0, . . . , n − 1}, such that x ∈ (x k , x k+1 ), and it follows from Lemma 4 that where the sum of the series can be found in [17, p. 464]. Together with Lemma 5, this implies for any i ∈ I 3 . As λ i (x) ≤ S 3 ≤ D(x) and N (x) − D(x) = −2(S 1 + S 5 ), we conclude that and the statement then follows immediately.
In the case of equidistant nodes, when the mesh ratio is μ = 1, the upper bound in (34) is simply 1 + 1/(2d), so it becomes smaller as d grows. For other nodes, μ may depend on n, which may result in very large upper bounds. For example, in the case of Chebyshev nodes, one can show that μ grows asymptotically linear in n. However, our numerical experiments suggest that the function d is always small, and we believe that the upper bound in Theorem 3 can be improved significantly in future work.

Numerical experiments
We performed numerous experiments to verify the results proven in the previous sections numerically and report some representative results below. In particular, we analyze the various algorithms that implement the first barycentric form (7) both in terms of stability and computational cost (Sect. 7.1) and focus on the comparison between the first and the second form for an example where n d (Sect. 7.2). Our experimental platform is a Windows 10 laptop with 1.8 GHz Intel Core i7-10510U processor and 16 GB RAM, and we implemented all algorithms in C++. In what follows, the 'exact' values were computed in multiple-precision (1024 bit) floating point arithmetic using the MPFR library [18], while all other values were computed in standard double precision. Moreover, we took care of providing all input data (interpolation nodes, data values, and evaluation points) in double precision, so that they do not cause any additional error.

Comparison of algorithms for the first barycentric form
For the first example, we consider the case of n + 1 interpolation nodes x i = fl(y i ) ∈ [−1, 1] for i = 0, . . . , n, derived from the equidistant nodes y i = 2i/n−1 by rounding them to double precision, and associated (rounded) data f i = fl( f (y i )), sampled from the test function and we compare three ways of evaluating the first barycentric form in (7), which differ in the way the denominator is computed. The first algorithm simply evaluates the functions λ i as in (8), leading to the error mentioned in Lemma 2 and then sums up these values to get the denominator. The second algorithm by Camargo [11,Section 4.1] instead increases the stability of the summation by first computing sums of pairs of λ i 's such that all these sums have the same sign. The third algorithm implements our iterative strategy in (31) before taking the sum, which is more efficient than the first algorithm, but less precise (cf. Lemma 3). All three algorithms compute the numerator of the first barycentric form in the same way, first dividing the weights w i by x − x i , then multiplying the results by f i , and finally summing up these values. The w i themselves are precomputed with the pyramid algorithm in (29) and (30). Note that, although the weights for equidistant nodes are integer multiples of each other in theory [6], they do not have this property in this example, because the nodes x i are not exactly equidistant, because of the rounding.
To compare these three algorithms, we used them to evaluate the barycentric rational interpolant with weights in (6) for d ∈ {1, 5, 25} and an increasing number of interpolation nodes, n ∈ {39, 79, 159, 319, 639, 1279}, at 50,000 random points from [−1 + 10 3 , 1 − 10 3 ]\{x 0 , . . . , x n }. Figure 1 shows the corresponding running times and the distribution of the relative forward error of the computed values. For the latter, we chose a box plot, where the bottom and top of each box represent the interquartile range, from the first to the third quartile, and the line in the middle shows the median of the relative forward errors. The whiskers range from the minimum to the maximum, excluding those values that are more than 1.5 times the interquartile range greater than the third quartile, which are instead considered outliers and shown as isolated points.
On the one hand, we observe that Camargo's algorithm beats the standard algorithm in terms of running time, but that our efficient algorithm is the fastest, especially as d grows, because its time complexity does not depend on d. On the other hand, our efficient algorithm is less precise than the standard algorithm, as predicted by Lemma 3 and Camargo's algorithm gives the smallest errors, except for d = 5 and n ∈ {639, 1279}. Nevertheless, the computations confirm the forward stability for all three algorithms. For Camargo's algorithm, this follows from the backward stability, which is proven in [11], and for the other two algorithms it is implied by Corollary 2 and Lemmas 1-3.
The rather large errors of the outliers in the case d = 25 can be explained by the behaviour of the componentwise relative condition number κ, shown in Figure 2.
While κ is small for all x ∈ [−1, 1] in the case of d = 1, it starts to grow considerably close to the end points of the interval as d grows, up to 10 6 for n = 39 and 10 8 for n = 1279 in the case of d = 25, and so does the upper bound on the relative forward error in (28). While this upper bound also depends on d , Figure 2 shows that this function is always small and its maximum even decreases as d grows, independently of n. This is in agreement with Theorem 3, because μ ≈ 1 in this example. The fact that the maximum error still seems to decrease for d = 25 as n increases is simply due to the fact that the 50,000 sample points are not sufficiently many to "catch" the worst case, because the region near the endpoints where κ grows rapidly actually shrinks as n grows.
Of course, it is also possible to evaluate the rational interpolant using the second barycentric form (4), which is actually the best choice for this example, giving relative forward errors that are similar to the ones of the standard algorithm for the first barycentric form, but being roughly twice as fast as the efficient algorithm, both of which is not surprising. Regarding the efficiency, the second form is superior, because the denominator can be computed "on-the-fly" at almost no extra cost during the evaluation of the numerator. As for the error, we note that n is much smaller than κ for the nodes used in this example [10] and so the upper bound in (25) is dominated by κ, exactly as the upper bound in (28). However, for non-uniform nodes, the situation can be quite different, as the next example shows.

Worst-case comparison of first and second barycentric form
The aim of the second example is to compare the standard algorithm for the first barycentric form in (7), as described in Sect. 7.1, with a straightforward implementation of the second barycentric form, following the formula in (4). The weights w i are again precomputed with the pyramid algorithm.
We consider n = 29, d = 3, and interpolation nodes x i = fl(y i ) ∈ [0, 1] for i = 0, . . . , n, obtained by rounding to double precision the values y i = F(t i ), where and t i = i/n. We choose these nodes, because the functions n and d behave completely differently in this case, as shown in Fig. 3. While n reaches huge values, up to 10 17 , d is small (even though the latter is not guaranteed by Theorem 3). Hence, the upper bounds in Corollaries 1 and 2 suggest that we may see a big difference in the forward stability of the first and the second barycentric form, if we choose the data such that the condition number κ is small. One such choice, which is also presented by Higham for the case d = n with equidistant nodes [1], is to take the data which can be interpreted as having been sampled from the nth Lagrange basis polynomial l n (x) = n−1 i=0 x−x i x n −x i , that is, f i = fl(l n (y i )). For this data, we know that κ(x; x, f ) = 1 for all x ∈ [0, 1], so the upper bounds on the relative forward errors in (25) and (28) are dominated by n and d , respectively. Consequently, as shown in Fig. 4, the barycentric rational interpolant is reproduced faithfully by the first form, but not by the second, because the relative forward error for the first form is on the order of , while it can be on the order of 1 for the second form; see   (25) and (28) suggest that both forms can be unstable, even though the barycentric rational interpolant is simply r (x) = 1. Figure 5 (right) and Figure 6 confirm that this is indeed the case for the first barycentric form. However, the second barycentric form is perfectly stable, because the numerator and the denominator in (4) are identical and cancel out to always give the correct function value 1.

Conclusion
Barycentric interpolation offers a fast and stable means of evaluating the polynomial Lagrange interpolant using either the first or the second barycentric form. While the first form is always numerically stable, the second form is stable only for interpolation nodes with a small Lebesgue constant. Evaluating a rational interpolant via the second barycentric form comes with the same limitation, but for the special family of barycentric rational interpolants with weights in (6), a computationally more attractive first barycentric form is available. Instead of depending on the Lebesgue constant, both the forward and the backward stability of a straightforward implementation of this first barycentric form depend on the function d in (26).
Unlike the Lebesgue constant, Theorem 3 shows that the maximum of d is independent of n. Moreover, it is guaranteed to be very small for equidistant nodes, regardless of d, while the Lebesgue constant is known to grow logarithmically in n and exponentially in d in this case. Based on our numerical experiments, the maximum of d seems to be small for other distributions of interpolation nodes, too, even if the mesh ratio μ is big, as in the example in Sect. 7.2, and we believe that the upper bound in (34) can be improved considerably in future work. For example, if d = n, then d is just the constant one function, independent of μ.
Overall, we conclude that the most efficient way of stably evaluating a barycentric rational interpolant with weights in (6) is by determining the weights with the pyramid algorithm in (29) and (30) in a preprocessing step, and, for every evaluation point x, first computing the values λ i (x) with our iterative strategy in (31) and then the value of the interpolant using a straightforward implementation of the first barycentric form in (7). Alternatively computing the sum of the λ i (x) in the denominator with Camargo's algorithm [11] results in slightly smaller forward errors, but is significantly slower, especially for larger d.