Thiele's interpolation formula

Statement

Let $n \ge 1$, let $x_1, \dots, x_{n+1} \in \R$ be pairwise distinct, and let $f$ be a real function defined at these points. The reciprocal differences $\rho_k$ of $f$ are defined by $\rho_0(x_i) = f(x_i)$, $$\rho_1(x_i, x_{i+1}) = \frac{x_i - x_{i+1}}{f(x_i) - f(x_{i+1})},$$ and, for $k \ge 2$, $$\rho_k(x_i, \dots, x_{i+k}) = \frac{x_i - x_{i+k}}{\rho_{k-1}(x_i, \dots, x_{i+k-1}) - \rho_{k-1}(x_{i+1}, \dots, x_{i+k})} + \rho_{k-2}(x_{i+1}, \dots, x_{i+k-1}).$$ Set $a_0 = f(x_1)$, $a_1 = \rho_1(x_1, x_2)$ and $a_k = \rho_k(x_1, \dots, x_{k+1}) - \rho_{k-2}(x_1, \dots, x_{k-1})$ for $2 \le k \le n$. Assume that no denominator above is zero and that the finite continued fraction $$R(x) = a_0 + \cfrac{x - x_1}{a_1 + \cfrac{x - x_2}{a_2 + \cfrac{x - x_3}{\ddots + \cfrac{x - x_n}{a_n}}}}$$ can be evaluated at each node without division by zero. Then $R$ is a rational function and $R(x_i) = f(x_i)$ for $i = 1, \dots, n + 1$.