Rational Interpolation: Jacobi’s Approach Reminiscence

Abstract

We treat the interpolation problem $${\{f\left(x_j\right)=y_j\}}_{j=1}^N$$ for polynomial and rational functions. Developing the approach originated by C. Jacobi, we represent the interpolants by virtue of the Hankel polynomials generated by the sequences of special symmetric functions of the data set like $${\{{\sum }_{j=1}^N{x}_j^k{y}_j/W^{\prime }\left(x_j\right)\}}_{k\in \mathbb{N}}$$ and $${\{{\sum }_{j=1}^N{x}_j^k/\left(y_j{W}^{\prime }\left(x_j\right)\right)\}}_{k\in \mathbb{N}}$$; here, $$W\left(x\right)={\prod }_{j=1}^N(x-x_j)$$. We also review the results by Jacobi, Joachimsthal, Kronecker and Frobenius on the recursive procedure for computation of the sequence of Hankel polynomials. The problem of evaluation of the resultant of polynomials $$p\left(x\right)$$ and $$q\left(x\right)$$ given a set of values $${\{p\left(x_j\right)/q\left(x_j\right)\}}_{j=1}^N$$ is also tackled within the framework of this approach. An effective procedure is suggested for recomputation of rational interpolants in case of extension of the data set by an extra point.

1. Introduction

Given the data set for the variables x and y

(1)$\begin{array}{cccccc}x& & x_1& x_2& \dots & x_N\\ y& & y_1& y_2& \dots & y_N\end{array},$

with distinct nodes $\{x_1,\dots ,x_N\}\subset \mathbb{R}$, we treat the interpolation problem in a sense of finding a function $f\left(x\right)$ such that $f\left(x_j\right)=y_j\in \mathbb{R}$ for $j\in \{1,\dots ,N\}$. Systematic exploration of the problem was started since 17th century; for the historical overview and the review of numerous applications of the relevant theory, we refer to [1] and the references cited therein.

Depending on the concrete objectives, the solution to the problem can be looked for in various classes of interpolants like algebraic polynomial, trigonometric polynomial, sum of exponentials, rational functions, etc. The two interpolation problems dealt with in the present paper are polynomial and rational ones.

In comparison with the polynomial interpolation problem, the rational interpolation one has its beginning a century and a half later, with the first explicit formula due to Cauchy [2]. Further its development was made by Kronecker [3] and Netto [4] (we briefly discuss it in Section 5). The interest to the problem revives in the second part of the 20th century and is connected with its application in Approximation Theory, Control Theory (recovering of the transfer function from the frequency responses) and Error Correcting Codes; in the latter case, the problem is treated in finite fields. We refer to [5,6,7,8] for recent developments and further references on the rational interpolation problem as well as to its generalization, known as rational Hermite’s or osculatory rational interpolation problem, where the values of some derivatives for $f\left(x\right)$ are assumed to be specified at some nodes.

The present paper is devoted to an approach to the rational interpolation problem originated in 1846 by Carl Jacobi [9]. Within the framework of this approach, the numerator and the denominator of the rational interpolant

$f\left(x\right)=p\left(x\right)/q\left(x\right),degp\left(x\right)+degq\left(x\right)=N-1$

are constructed in the form of the so called Hankel polynomials, i.e., polynomials represented in the determinantal form as

(2)${\mathcal{H}}_k\left(x\right)=\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & c_2\hfill & \dots \hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & c_3\hfill & \dots \hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ c_{k-1}\hfill & c_k\hfill & c_{k+1}\hfill & \dots \hfill & c_{2k-1}\hfill \\ 1\hfill & x\hfill & x^2\hfill & \dots \hfill & x^k\hfill \end{array}\right|.$

with the entries $\{c_0,c_1,\cdots ,c_{2k-1}\}\subset \mathbb{C}$ called its generators. Numerator and denominator of the rational interpolant can be expressed as Hankel polynomials of appropriate orders with the generating sequences chosen in the form of suitable symmetric functions of the data set (1), namely

(3)${\tilde{\tau }}_k=\sum _{j=1}^N\frac{1}{y_j}\frac{x_j^k}{W^{\prime }\left(x_j\right)}\mathrm{and}{\tau }_k=\sum _{j=1}^N{y}_j\frac{x_j^k}{W^{\prime }\left(x_j\right)}\mathrm{for}k\in \{0,1,\dots \}$

correspondingly; here $W\left(x\right)={\prod }_{j=1}^N(x-x_j)$. Jacobi’s approach has not got much attention and further development in modern computational mathematics; perhaps the only exception is the paper [10]. This disregard is probably caused by an evident computational bottleneck, namely the computation of a parameter dependent determinant. In case of its large order, it is not a trivial task.

Fortunately a specific structure of the involved determinants helps a lot: there exists a recursive procedure for the Hankel polynomial computation. It is based on the identity linking the Hankel determinants of three consecutive orders:

(4)$A_k{\mathcal{H}}_{k-2}\left(x\right)+(B_k-x){\mathcal{H}}_{k-1}\left(x\right)+\frac{1}{A_k}{\mathcal{H}}_k\left(x\right)\equiv 0,$

here, $A_k$ and $B_k$ are some constants that can be expressed via the coefficients of ${\mathcal{H}}_{k-2}\left(x\right)$ and ${\mathcal{H}}_{k-1}\left(x\right)$. Formula (4) looks similar to the recurrent formula connecting orthogonal polynomials. It is indeed the case for the positive definite Hankel matrix, i.e., under this condition, the Hankel polynomials ${\mathcal{H}}_1\left(x\right),{\mathcal{H}}_2\left(x\right),\dots$ are orthogonal with respect to the inner product introduced via this matrix taken as the Gram matrix. In [11], this idea was extended to the Hankel matrices with the non-zero leading principal minors. Application of the recursive procedures for computation of the sequences of Hankel polynomials to the rational interpolation problem provides one with an opportunity to efficiently compute not only a single interpolant, but the whole family of fractions with all possible combinations of degrees for the numerator and the denominator satisfying the restriction $degp\left(x\right)+degq\left(x\right)=N-1$.

One issue of this scheme should be additionally discussed, namely constants’ evaluation in (4). For the classical orthogonal polynomials or for their counterparts from [11], this is done via the inner product computations involving polynomials ${\mathcal{H}}_{k-1}\left(x\right)$ and $x{\mathcal{H}}_{k-1}\left(x\right)$. However, there exists an alternative algorithm for the computation of the constants in (4). Surprisingly, this version of the identity (4) is also related to Jacobi with the original idea outlined in his paper [12] as of 1836, and further complete treatment executed by his disciple Ferdinand Joachimsthal [13]. This approach is free of the orthogonality concept and is based on purely determinantal formalism.

The present paper is organized as follows. We first discuss the properties of Hankel polynomials in Section 3 starting with the Jacobi–Joachimsthal result. The first two of the present authors have excavated this result and published it with the original Joachimsthal’s proof in the paper [14]. For the convenience of a reader of the present paper, we decided to include an English translation of the proof (original Jacobi’s paper was written in Latin, Joachimsthal’s—in German, the paper [14] is in Russian). Further investigation of the 19th century heredity, results in discovering the generalization of the Jacobi–Joachimsthal result (the cases of degeneracy of (4) like, for instance, $A_k=0$) in the works by Kronecker and Frobenius [15]; we include their results in Section 3. The necessity in detailed presentation of this classical algebra material (with the exception of Theorem 8, which is due to the present authors) is justified by its complete forgetfulness in the next two centuries. Indeed, we failed to find any references to them either in the monographs [16,17,18] or even in research papers. This state of affairs might seem unfair, but not an extraordinary one, if not for one circumstance. As a matter of fact, the referred results should be treated as the gist of the algorithm, which is nowadays known as the Berlekamp–Massey algorithm [19,20]; it was suggested for the decoding procedure in Bose–Chaudhury–Hocquenghem (BCH) codes and for finding the minimal polynomial of a linear recurrent sequence.

The deployment of this Hankel polynomial formalism to the main problem of the paper is started in Section 4 with the treatment of the polynomial interpolation. We discuss here some properties of the sums like (3), i.e., symmetric functions of the data set. The rational interpolant construction is dealt with in Section 5 where the efficiency of the two approaches (i.e., Jacobi–Joachimsthal and orthogonality-based) for computation of the Hankel polynomial via (4) is discussed.

The results of Section 6 and Section 7 should be regarded as completely original. In Section 6, we deal with the question which relates to the uniqueness problem for the rational interpolation. Having the data set (1) generated by some rational fraction $p\left(x\right)/q\left(x\right)$ with $degp\left(x\right)+degq\left(x\right)=N-1$, is it possible to establish that the fraction is irreducible, (or, equivalently, that polynomials $p\left(x\right)$ and $q\left(x\right)$ do not possess a common zero) avoiding, as an intermediate step, their explicit representation? The question is equivalent to the possibility of expressing the resultant of polynomials $p\left(x\right)$ and $q\left(x\right)$ in terms of the entries of the data set (1). We prove that the resultant can be expressed in the form of an appropriate Hankel determinant generated by any of the sequences (3).

In Section 7, we present an effective procedure for recomputation of Hankel polynomials for the data set (1) complemented by an additional point. It turns out that generically the kth order Hankel polynomial composed for the extended data set can be expressed as a linear combination of two of its counterparts (of the orders k and $k-1$) constructed for the initial set (1).

Relationship of the Jacobi’s approach to the rational interpolation problem with that one based on barycentric formula [6,21] is discussed in Section 8.

In Section 9, we briefly outline an application of the Hankel polynomial formalism to the problems of Coding Theory. We address here the Berlekamp–Welch algorithm [22] for the error detection and correction in a data sequence.

2. Algebraic Preliminaries

For the convenience of further references, we recall here some auxiliary results from the Matrix theory and Polynomial Field theory.

Theorem 1

(Kronecker [23,24]). If in a given matrix a certain $\mathfrak{r}$th order minor is not zero, and all the $(\mathfrak{r}+1)$th order minors containing that $\mathfrak{r}$th order minor are zero, then all the $(\mathfrak{r}+1)$th minors are zero, and, therefore, the rank of the matrix equals $\mathfrak{r}$.

Theorem 2

(Sylvester [16]). For a square matrix A of the order n denote by

$A\left(\begin{array}{cccc}{i}_1& i_2& \dots & i_n\\ j_1& j_2& \dots & j_n\end{array}\right)$

its minor standing in the rows $i_1,i_2,\dots ,i_k$ and in the columns $j_1,j_2,\dots ,j_k$. The following Sylvester’s identity is valid:

(5)$$\begin{gathered} A\left(\begin{array}{cccc}1& 2& \dots & n-2\\ 1& 2& \dots & n-2\end{array}\right)·detA \\ =A\left(\begin{array}{ccccc}1& 2& \dots & n-2& n-1\\ 1& 2& \dots & n-2& n-1\end{array}\right)·A\left(\begin{array}{ccccc}1& 2& \dots & n-2& n\\ 1& 2& \dots & n-2& n\end{array}\right) \\ -A\left(\begin{array}{ccccc}1& 2& \dots & n-2& n\\ 1& 2& \dots & n-2& n-1\end{array}\right)·A\left(\begin{array}{ccccc}1& 2& \dots & n-2& n-1\\ 1& 2& \dots & n-2& n\end{array}\right). \end{gathered}$$

We next recall the definition of the resultant of polynomials [8,23,25]. For the polynomials

(6)$p\left(x\right)=p_0{x}^n+p_1{x}^{n-1}+\cdots +p_n\mathrm{and}q\left(x\right)=q_0{x}^m+q_1{x}^{m-1}+\cdots +q_m$

with $p_0\ne 0,q_0\ne 0,n\ge 1,m\ge 1$ their resultant is formally defined as

(7)$\mathcal{R}(p\left(x\right),q\left(x\right))=p_0^m\prod _{j=1}^{n}q\left({\lambda }_j\right)$

where ${\lambda }_1,\dots ,{\lambda }_n$ denote the zeros of $p\left(x\right)$ (counted in accordance with their multiplicities). Equivalently, the resultant can be defined as

(8)$\mathcal{R}(p\left(x\right),q\left(x\right))={(-1)}^{mn}{q}_0^n\prod _{\ell =1}^{m}p\left({\mu }_{\ell }\right)$

where ${\mu }_1,\dots ,{\mu }_m$ denote the zeros of $q\left(x\right)$ (counted in accordance with their multiplicities). As for the constructive methods of computing the resultant as a polynomial function of the coefficients of $p\left(x\right)$ and $q\left(x\right)$, this can be done with the aid of several determinantal representations (like Sylvester’s, Bézout’s or Kronecker’s).

Example 1.

For the polynomials

$p\left(x\right)=p_0{x}^3+p_1{x}^2+p_{2}x+p_{3}andq\left(x\right)=q_0{x}^5+\cdots +q_5$

with $p_0\ne 0,q_0\ne 0$, their resultant in Sylvester’s form is the $(3+5)$-order determinant:

$\mathcal{R}(p\left(x\right),q\left(x\right))=-\left|\begin{array}{cccccccc}{p}_0& p_1& p_2& p_3& & & & \\ & p_0& p_1& p_2& p_3& & & \\ & & p_0& p_1& p_2& p_3& & \\ & & & p_0& p_1& p_2& p_3& \\ & & & & p_0& p_1& p_2& p_3\\ & & q_0& q_1& q_2& q_3& q_4& q_5\\ & q_0& q_1& q_2& q_3& q_4& q_5& \\ q_0& q_1& q_2& q_3& q_4& q_5& & \end{array}\right|\begin{array}{c}\left.\begin{array}{c}\\ \\ \\ \\ \end{array}\right\}5\hfill \\ \left.\begin{array}{c}\\ \\ \end{array}\right\}3\hfill \end{array}$

Not indicated entries of the determinant are assigned to zero.

Theorem 3.

Polynomials $p\left(x\right)$ and $q\left(x\right)$ possess a common zero if and only if their resultant vanishes: $\mathcal{R}\left(p\right(x),q(x\left)\right)=0$.

An important particular case related to the resultant of a polynomial and its derivative gives rise to a special notion: the expression

$\mathcal{D}\left(p\left(x\right)\right)=\frac{{(-1)}^{n(n-1)/2}}{p_0}\mathcal{R}(p\left(x\right),p^{\prime }\left(x\right))$

is known as the discriminant of the polynomial $p\left(x\right)$. It is a polynomial function in the coefficients $p_0,\dots ,p_n$.

Corollary 1.

Polynomial $p\left(x\right),degp\left(x\right)\ge 2$ possesses a multiple zero if and only if $\mathcal{D}\left(p\right(x\left)\right)=0$.

We conclude the list of algebraic preliminaries with Hermann Weyl’s result known as the Principle of the Irrelevance of Algebraic Inequalities ([26], p. 4, Lemma (1.1.A); in the following original statement, k stands for an infinite integral domain):

Theorem 4

(Weyl H.). A k-polynomial $F(x,y,\dots )$ vanishes identically if it vanishes numerically for all sets of rational values $x=\alpha ,y=\beta ,\dots$ subject to a number of algebraic inequalities

$R_1(\alpha ,\beta ,\dots )\ne 0,R_2(\alpha ,\beta ,\dots )\ne 0,\dots$

3. Hankel Determinants and Polynomials

For a (finite or infinite) sequence of complex numbers

(9)${\left\{c_j\right\}}_{j=0}^{\infty }=\left\{c_0,c_1,\dots \right\}$

a square matrix in which each ascending skew-diagonal from left to right is constant

(10)${\left[c_{i+j-2}\right]}_{i,j=1}^k={\left[\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill \\ c_1\hfill & c_2\hfill & c_3\hfill & \dots \hfill & c_k\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-3}\hfill \\ c_{k-1}\hfill & c_k\hfill & c_{k+1}\hfill & \dots \hfill & c_{2k-2}\hfill \end{array}\right]}_{k\times k}$

is called Hankel matrix of order k generated by the sequence (9); its determinant will be denoted by $H_k\left(\left\{c\right\}\right)$ or simply $H_k$ if it will not cause misunderstandings.

The following result is a particular case of Sylvester’s identity (5).

Theorem 5.

Hankel determinants of three consecutive orders are linked by the equality (sometimes also referred to as the Jacobi identity [17]):

(11) $$\begin{gathered} H_{k-2}{H}_k \\ =H_{k-1}\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots & c_{k-3}\hfill & \mid c_{k-1}\hfill \\ c_1\hfill & c_2\hfill & \dots & c_{k-2}\hfill & \mid c_k\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill \\ c_{k-3}\hfill & c_{k-2}\hfill & \dots & c_{2k-6}\hfill & \mid c_{2k-4}\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots & c_{2k-4}\hfill & \mid c_{2k-2}\hfill \end{array}\right|-{\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots & c_{k-3}\hfill & c_{k-1}\hfill \\ c_1\hfill & c_2\hfill & \dots & c_{k-2}\hfill & c_k\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill \\ c_{k-3}\hfill & c_{k-2}\hfill & \dots & c_{2k-6}\hfill & c_{2k-4}\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots & c_{2k-5}\hfill & c_{2k-3}\hfill \end{array}\right|}^2. \end{gathered}$$

If we replace the last row of the Hankel matrix of order $k+1$ with the row of powers of x, then the corresponding determinant

(12)${\mathcal{H}}_k(x;\left\{c\right\})={\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & c_2\hfill & \dots \hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & c_3\hfill & \dots \hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ c_{k-1}\hfill & c_k\hfill & c_{k+1}\hfill & \dots \hfill & c_{2k-1}\hfill \\ 1\hfill & x\hfill & x^2\hfill & \dots \hfill & x^k\hfill \end{array}\right|}_{(k+1)\times (k+1)}$

or simply ${\mathcal{H}}_k\left(x\right)$, is called the kth Hankel polynomial [17] generated by the sequence (9). Expansion of (12) by its last row yields

(13)${\mathcal{H}}_k\left(x\right)\equiv h_{k0}{x}^k+h_{k1}{x}^{k-1}+h_{k2}{x}^{k-2}+\dots \mathrm{with}{h}_{k0}=H_k.$

Thus, $deg{\mathcal{H}}_k\left(x\right)=k$ iff $H_k\ne 0$. We will also utilize an alternative representation for Hankel polynomial in the form of Hankel determinant resulting from linear transformation of the columns of the determinant (12):

(14)$\begin{array}{ccc}\hfill {\mathcal{H}}_k(x;\left\{c\right\})& =& {(-1)}^k{\left|\begin{array}{cccc}{c}_1-c_{0}x\hfill & c_2-c_{1}x\hfill & \dots \hfill & c_k-c_{k-1}x\hfill \\ c_2-c_{1}x\hfill & c_3-c_{2}x\hfill & \dots \hfill & c_{k+1}-c_{k}x\hfill \\ \cdots \hfill & \cdots \hfill & & \cdots \hfill \\ c_k-c_{k-1}x\hfill & c_{k+1}-c_{k}x\hfill & \dots \hfill & c_{2k-1}-c_{2k-2}x\hfill \end{array}\right|}_{k\times k}\hfill \\ & =& {(-1)}^k{H}_k\left({\{c_{j+1}-c_{j}x\}}_{j=0}^{\infty }\right)\hfill \end{array}$

It turns out that there exists an explicit linear relation between any three consecutive Hankel polynomials generated by arbitrary sequence (9).

Theorem 6

(Jacobi, Joachimsthal, refs. [12,13]). Any three consecutive Hankel polynomials

${\mathcal{H}}_{k-2}\left(x\right),{\mathcal{H}}_{k-1}\left(x\right),{\mathcal{H}}_k\left(x\right)$

are linked by the following identity

(15)$H_k^2{\mathcal{H}}_{k-2}\left(x\right)+\left(H_k{h}_{k-1,1}-H_{k-1}{h}_{k1}-H_k{H}_{k-1}x\right){\mathcal{H}}_{k-1}\left(x\right)+H_{k-1}^2{\mathcal{H}}_k\left(x\right)\equiv 0$

which hereinafter will be referred to as the JJ-identity.

Proof is a slightly modernized original proof from [13].

Lemma 1.

Let the Hankel polynomial ${\mathcal{H}}_k\left(x\right)$ be generated by the sequence

(16)$c_j=\sum _{\ell =1}^m{\lambda }_{\ell }^j,\mathrm{for}j\in \{0,\dots ,2k-1\}$

for arbitrary distinct ${\lambda }_1,\dots ,{\lambda }_m$ with $m>k$; hence, $c_0=m$. Then, the following equalities are valid:

(17)$\sum _{\ell =1}^m{\lambda }_{\ell }^j{\mathcal{H}}_k\left({\lambda }_{\ell }\right)=\left\{\begin{array}{cc}0\hfill & ifj\in \{0,\dots ,k-1\},\hfill \\ H_{k+1}\hfill & ifj=k.\hfill \end{array}\right.$

Proof. 

$$\begin{gathered} {\lambda }_1^j{\mathcal{H}}_k\left({\lambda }_1\right)+{\lambda }_2^j{\mathcal{H}}_k\left({\lambda }_2\right)+\cdots +{\lambda }_m^j{\mathcal{H}}_k\left({\lambda }_m\right) \\ ={\lambda }_1^j\left|\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-1}\hfill \\ 1\hfill & {\lambda }_1\hfill & \dots \hfill & {\lambda }_1^k\hfill \end{array}\right|+\cdots +{\lambda }_m^j\left|\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-1}\hfill \\ 1\hfill & {\lambda }_m\hfill & \dots \hfill & {\lambda }_m^k\hfill \end{array}\right|. \end{gathered}$$

Using the linear property of the determinant convert this linear combination of determinants into a single one:

$\left|\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-1}\hfill \\ c_j\hfill & c_{j+1}\hfill & \dots \hfill & c_{j+k}\hfill \end{array}\right|.$

If $j<k$, then the last determinant possesses two identical rows. Therefore, in this case, it is just zero. For $j=k$ the obtained determinant coincides with $H_{k+1}$.  □

Proof of Theorem 6. 

First consider the case where the generating sequence for the polynomials ${\mathcal{H}}_{k-2}\left(x\right),{\mathcal{H}}_{k-1}\left(x\right),{\mathcal{H}}_k\left(x\right)$ is given by (16).  □

Assuming that $H_{k-1}\ne 0$ (i.e., $deg{\mathcal{H}}_{k-1}\left(x\right)=k-1$) divide ${\mathcal{H}}_k\left(x\right)$ by ${\mathcal{H}}_{k-1}\left(x\right)$:

(18)${\mathcal{H}}_k\left(x\right)\equiv Q\left(x\right){\mathcal{H}}_{k-1}\left(x\right)+R\left(x\right).$

Here, the undetermined coefficients of the quotient $Q\left(x\right)$ can be determined from those of ${\mathcal{H}}_k\left(x\right)$ and ${\mathcal{H}}_{k-1}\left(x\right)$ via equating the coefficients of corresponding powers of x in the both sides of (18):

(19)$Q\left(x\right)=Q_0+Q_{1}x\mathrm{where}{Q}_1=\frac{H_k}{H_{k-1}},Q_0=\frac{H_{k-1}{h}_{k1}-H_k{h}_{k-1,1}}{H_{k-1}^2},$

or, equivalently,

(20)$Q\left(x\right)\equiv \frac{1}{H_{k-1}^2}\left|\begin{array}{ccc}0& H_{k-1}& H_k\\ H_{k-1}& h_{k-1,1}& h_{k1}\\ 1& x& 0\end{array}\right|.$

To find the coefficients of the remainder

(21)$R\left(x\right)=R_0+R_{1}x+\cdots +R_{k-2}{x}^{k-2}$

substitute $x={\lambda }_1,\dots ,x={\lambda }_m$ into (18):

(22)${\left\{{\mathcal{H}}_k\left({\lambda }_{\ell }\right)=\left(Q_1{\lambda }_{\ell }+Q_0\right){\mathcal{H}}_{k-1}\left({\lambda }_{\ell }\right)+\left(R_0+R_1{\lambda }_{\ell }+\cdots +R_{k-2}{\lambda }_{\ell }^{k-2}\right)\right\}}_{\ell =1}^m.$

Summation of these equalities yields

(23)$$\begin{gathered} \sum _{\ell =1}^m{\mathcal{H}}_k\left({\lambda }_{\ell }\right) \\ =Q_0\sum _{\ell =1}^m{\mathcal{H}}_{k-1}\left({\lambda }_{\ell }\right)+Q_1\sum _{\ell =1}^m{\lambda }_{\ell }{\mathcal{H}}_{k-1}\left({\lambda }_{\ell }\right)+(c_0{R}_0+c_1{R}_1+\cdots +c_{k-2}{R}_{k-2}). \end{gathered}$$

Due to (17), one gets

$0=c_0{R}_0+c_1{R}_1+\cdots +c_{k-2}{R}_{k-2}.$

Next multiply every equality (22) by the corresponding ${\lambda }_{\ell }^j$ and sum up the obtained equalities by ℓ. For $j\in \{1,\dots ,k-3\}$, the resulting equality looks similar to the previously deduced:

$0=c_j{R}_0+c_{j+1}{R}_1+\cdots +c_{j+k-2}{R}_{k-2}.$

Multiplication of equalities (22) by ${\lambda }_{\ell }^{k-2}$ yields something different:

$0=H_k{Q}_1+c_{k-2}{R}_0+c_{k-1}{R}_1+\cdots +c_{2k-4}{R}_{k-2}.$

Unifying the obtained relations for $R_0,\dots ,R_{k-2}$ with (21), one gets the linear system:

$\left\{\begin{array}{cccccc}\hfill c_0{R}_0& \hfill +c_1{R}_1& \hfill +\dots & \hfill +c_{k-2}{R}_{k-2}& & \hfill =0,\\ \hfill c_1{R}_0& \hfill +c_2{R}_1& \hfill +\dots & \hfill +c_{k-1}{R}_{k-2}& & \hfill =0,\\ \hfill \cdots & & & & & \hfill \cdots \\ \hfill c_{k-3}{R}_0& \hfill +c_{k-2}{R}_1& \hfill +\dots & \hfill +c_{2k-5}{R}_{k-2}& & \hfill =0,\\ \hfill c_{k-2}{R}_0& \hfill +c_{k-1}{R}_1& \hfill +\dots & \hfill +c_{2k-4}{R}_{k-2}& \hfill +H_k{Q}_1& \hfill =0,\\ \hfill R_0& \hfill +x{R}_1& \hfill +\dots & \hfill +x^{k-2}{R}_{k-2}& \hfill -R\left(x\right)& \hfill =0.\end{array}\right.$

Consider it as a system of homogeneous equations with respect to the variables $R_0$, $R_1,\dots ,R_{k-2}$, 1. Since it possesses a nontrivial solution, its determinant necessarily vanishes:

$\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill & 0\\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill & 0\\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\\ c_{k-3}\hfill & c_{k-2}\hfill & \dots \hfill & c_{2k-5}\hfill & 0\\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-4}\hfill & H_k{Q}_1\\ 1\hfill & x\hfill & \dots \hfill & x^{k-2}\hfill & -R\left(x\right)\end{array}\right|=0.$

Expansion of the determinant by its last column yields:

$H_{k-1}R\left(x\right)+H_k{Q}_1{\mathcal{H}}_{k-2}\left(x\right)\equiv 0.$

Together with the already obtained expression (19) for $Q_1$, this confirms the validity of (15) for the particular case of the generating sequence given by (16).

Consider now the case of arbitrary generating sequence (9). For any given sequence of complex numbers $c_1,\dots ,c_{2k-1}$ it is possible to find complex numbers ${\lambda }_1,\dots ,{\lambda }_m$ with $m>2k-1$ such that the equations (16) become consistent. These numbers can be chosen to be the zeros of a polynomial of the degree m whose first $2k-1$ Newton sums [25] coincide with ${\left\{c_j\right\}}_{j=1}^{2k-1}$.

To complete the proof of (15), one should fill one gap in the arguments of the previous paragraph. Whereas the numbers $c_1,\dots ,c_{2k-1}$ can be chosen arbitrarily, the number $c_0$ takes the positive integer value, namely m. Thus, the validity of (15) is proved only for any positive integer $c_0$. However, this equality is an algebraic one in $c_0$. Being valid for an infinite set of integers, it should be valid for any $c_0\in \mathbb{C}$.

Joachimsthal did not provide the proof of the result for the case $H_{k-1}=0$. Theorem 4 manages this case. The left-hand side of (15) and $H_{k-1}$ are polynomials in $c_0,\dots ,c_{2k-1}$.

The relationship (15) gives rise to a more symmetric form:

Corollary 1.

If $H_k\ne 0,H_{k-1}\ne 0$, then the JJ-identity can be written down as

(24) $\frac{H_k}{H_{k-1}}{\mathcal{H}}_{k-2}\left(x\right)-\left(x-\frac{h_{k-1,1}}{H_{k-1}}+\frac{h_{k1}}{H_k}\right){\mathcal{H}}_{k-1}\left(x\right)+\frac{H_{k-1}}{H_k}{\mathcal{H}}_k\left(x\right)\equiv 0.$

The next result will be used in one of the further proofs.

Corollary 2.

For the coefficient $h_{k-1,2}$, the following representation is valid:

(25) $h_{k-1,2}=\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-3}\hfill & \mid c_{k-1}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-2}\hfill & \mid c_k\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill & ⋮\hfill \\ c_{k-3}\hfill & c_{k-2}\hfill & \dots \hfill & c_{2k-6}\hfill & \mid c_{2k-4}\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-4}\hfill & \mid c_{2k-2}\hfill \end{array}\right|-\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-3}\hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-2}\hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-5}\hfill & c_{2k-2}\hfill \end{array}\right|.$

Proof. 

Take the coefficient of $x^{k-2}$ from the left-hand side of the JJ-identity:

$H_k^2{H}_{k-2}+\left(H_k{H}_{k-1,1}-H_{k-1}{h}_{k1}\right)h_{k-1,1}-H_k{H}_{k-1}{h}_{k-1,2}+H_{k-1}^2{h}_{k2}=0.$

Under assumption that $H_{k-1}\ne 0,H_k\ne 0$, it follows

$h_{k-1,2}=\frac{H_k{H}_{k-2}+h_{k-1,1}^2}{H_{k-1}}-\frac{h_{k1}{h}_{k-1,1}-H_{k-1}{h}_{k2}}{H_k}.$

Apply Sylvester’s identity (5) to both fractions in the right-hand side. For the first one, we take it in the version (11) while for the second one we apply it to the matrix

${\left(\begin{array}{cccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill & c_{k-1}\hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill & c_k\hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-4}\hfill & c_{2k-3}\hfill & c_{2k-2}\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-3}\hfill & c_{2k-2}\hfill & c_{2k-1}\hfill \\ 0\hfill & 0\hfill & \dots \hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right)}_{(k+1)\times (k+1)}$

with the determinant equal to $H_k$:

$\frac{h_{k1}{h}_{k-1,1}-H_{k-1}{h}_{k2}}{H_k}=\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-3}\hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-2}\hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-5}\hfill & c_{2k-2}\hfill \end{array}\right|.$

To cover the case where $H_{k-1}=0$ or $H_k=0$, one should apply the result of Theorem 4.  □

The JJ-identity permits one to generate the recursive procedure for computation of Hankel polynomials. Indeed, assume that the expressions for ${\mathcal{H}}_{k-2}\left(x\right)$ and ${\mathcal{H}}_{k-1}\left(x\right)$ are already computed and

(26)${\mathcal{H}}_{k-1}\left(x\right)\equiv h_{k-1,0}{x}^{k-1}+h_{k-1,1}{x}^{k-2}+\cdots +h_{k-1,k-1},h_{k-1,0}=H_{k-1}.$

Then in (15) all the constants are also evaluated except for $H_k$ and $h_{k1}$ for which one has only their determinantal representations:

$H_k=\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill & c_{k-1}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill & c_k\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-4}\hfill & c_{2k-3}\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-3}\hfill & c_{2k-2}\hfill \end{array}\right|,h_{k1}=-\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-4}\hfill & c_{2k-2}\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-3}\hfill & c_{2k-1}\hfill \end{array}\right|.$

These determinants differ from the transposed determinantal representation for ${\mathcal{H}}_{k-1}\left(x\right)$ only in their last columns. Expansions by the entries of these columns have the same values for corresponding cofactors, and, therefore, the following formulas:

(27)$\left\{\begin{array}{ccc}\hfill h_{k0}=H_k& =& c_{k-1}{h}_{k-1,k-1}+c_k{h}_{k-1,k-2}+\cdots +c_{2k-2}{h}_{k-1,0},\hfill \\ \hfill h_{k1}& =& -(c_k{h}_{k-1,k-1}+c_{k+1}{h}_{k-1,k-2}+\cdots +c_{2k-1}{h}_{k-1,0})\hfill \end{array}\right.$

allow one to evaluate $h_{k0}$ and $h_{k1}$ via the already computed coefficients of ${\mathcal{H}}_{k-1}\left(x\right)$.

However, the just-outlined algorithm for recursive computation of ${\mathcal{H}}_k\left(x\right)$ fails for the case where $H_{k-1}=0$.

Theorem 7.

Let $H_{k-2}\ne 0,H_{k-1}=0$. If $h_{k-1,1}=0$, then $H_k=0$ and

${\mathcal{H}}_{k-1}\left(x\right)\equiv 0,$

and

(28) ${\mathcal{H}}_k\left(x\right)\equiv \frac{h_{k2}}{H_{k-2}}{\mathcal{H}}_{k-2}\left(x\right).$

Otherwise,

(29) ${\mathcal{H}}_{k-1}\left(x\right)\equiv \frac{h_{k-1,1}}{H_{k-2}}{\mathcal{H}}_{k-2}\left(x\right)$

and

(30) ${\mathcal{H}}_k\left(x\right)\equiv \frac{1}{H_{k-2}^3}\left(H_k{H}_{k-2}{h}_{k-1,1}{\mathcal{H}}_{k-3}\left(x\right)+\left|\begin{array}{cccc}0& 0& H_{k-2}& H_k\\ 0& H_{k-2}& h_{k-2,1}& h_{k1}\\ H_{k-2}& h_{k-2,1}& h_{k-2,2}& h_{k2}\\ 1& x& x^2& 0\end{array}\right|{\mathcal{H}}_{k-2}\left(x\right)\right).$

Proof. 

If $H_{k-2}\ne 0$ and $H_{k-1}=0,h_{k-1,1}=0$ then all the $(k-1)$th order minors of the matrix

${\left(\begin{array}{cccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-3}\hfill & c_{k-2}\hfill & c_{k-1}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-2}\hfill & c_{k-1}\hfill & c_k\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ c_{k-3}\hfill & c_{k-2}\hfill & \dots \hfill & c_{2k-6}\hfill & c_{2k-5}\hfill & c_{2k-4}\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-5}\hfill & c_{2k-4}\hfill & c_{2k-3}\hfill \end{array}\right)}_{(k-1)\times k}$

are zero. This follows from Theorem 1. Therefore, all the coefficients of the polynomial ${\mathcal{H}}_{k-1}\left(x\right)$ are also zero since they are the $(k-1)$th minors of this matrix.

To establish the validity of (28), it is sufficient to prove that the remainder of the division of ${\mathcal{H}}_k\left(x\right)$ by ${\mathcal{H}}_{k-2}\left(x\right)$ is identically zero. This will be done later, while now we intend to prove (29). Since $H_{k-1}=0$, the identity (15) can be rewritten as

(31)$H_k(H_k{\mathcal{H}}_{k-2}\left(x\right)+h_{k-1,1}{\mathcal{H}}_{k-1}\left(x\right))\equiv 0.$

If $h_{k-1,1}\ne 0$, then the Sylvester identity (5) leads one to

$H_k{H}_{k-2}=-h_{k-1,1}^2,$

from which it follows that

$H_k\ne 0\mathrm{and}{H}_k=-h_{k-1,1}^2/H_{k-2}.$

Consequently, the identity (31) leads one to

${\mathcal{H}}_{k-1}\left(x\right)\equiv -\frac{H_k}{h_{k-1,1}}{\mathcal{H}}_{k-2}\left(x\right)\equiv \frac{h_{k-1,1}}{H_{k-2}}{\mathcal{H}}_{k-2}\left(x\right)$

which proves (29).

We now intend to prove (28) and (30). Divide ${\mathcal{H}}_k\left(x\right)$ by ${\mathcal{H}}_{k-2}\left(x\right)$:

(32)${\mathcal{H}}_k\left(x\right)\equiv Q\left(x\right){\mathcal{H}}_{k-2}\left(x\right)+R\left(x\right).$

Here, the quotient $Q\left(x\right)$ equals $h_{k2}/H_{k-2}$ if $h_{k-1,1}=0$, while for the case $h_{k-1,1}\ne 0$ one has

$Q\left(x\right)\equiv Q_0+Q_{1}x+Q_2{x}^2$

with coefficients determined by the equalities:

(33)$Q_2=\frac{H_k}{H_{k-2}},Q_1=\frac{H_{k-2}{h}_{k1}-H_k{h}_{k-2,1}}{H_{k-2}^2},Q_0=\frac{H_{k-2}^2{h}_{k2}-H_k{H}_{k-2}^2{h}_{k-2,2}-H_k{h}_{k-2,1}^2}{H_{k-2}^3}.$

To find the coefficients of the remainder

(34)$R\left(x\right)=R_0+R_{1}x+\cdots +R_{k-3}{x}^{k-3}$

use arguments similar to those used in the proof of Theorem 6. First, consider the case where the generating sequence (9) is given by (16). Substitute $x={\lambda }_{\ell }$ into (32)

(35)${\left\{{\mathcal{H}}_k\left({\lambda }_{\ell }\right)=(Q_0+Q_1{\lambda }_{\ell }+Q_2{\lambda }_{\ell }^2){\mathcal{H}}_{k-2}\left({\lambda }_{\ell }\right)+R_0+\cdots +R_{k-3}{\lambda }_{\ell }^{k-3}\right\}}_{\ell =1}^m,$

and sum up the obtained equalities. Due to (17) one gets

$0=c_0{R}_0+c_1{R}_1+\cdots +c_{k-3}{R}_{k-3}.$

Similar equalities result from multiplication of (35) by ${\lambda }_{\ell }^j$ for $j\in \{1,\dots ,k-5\}$ and further summation by ℓ:

$\begin{array}{ccccc}0=& \hfill c_1{R}_0& \hfill +c_2{R}_1& +\dots & \hfill +c_{k-2}{R}_{k-3},\\ \cdots & & \hfill \cdots & & \\ 0=& \hfill c_{k-5}{R}_0& \hfill +c_{k-4}{R}_1& +\dots & \hfill +c_{2k-8}{R}_{k-3}.\end{array}$

Multiplication of (35) by ${\lambda }_{\ell }^{k-4}$ and summation yields

$0=Q_2{H}_{k-1}+c_{k-4}{R}_0+c_{k-3}{R}_1+\cdots +c_{2k-7}{R}_{k-3},$

and, since by assumption $H_{k-1}=0$, the obtained equality looks similar to the previous ones. Multiplication of (35) by ${\lambda }_{\ell }^{k-3}$ and summation leads to

$0=Q_2\sum _{j=1}^m{\lambda }_j^{k-1}{\mathcal{H}}_{k-2}\left({\lambda }_j\right)+c_{k-3}{R}_0+c_{k-2}{R}_1+\cdots +c_{2k-6}{R}_{k-3}$

with

$$\begin{gathered} \sum _{\ell =1}^m{\lambda }_{\ell }^{k-1}{\mathcal{H}}_{k-2}\left({\lambda }_{\ell }\right) \\ =\left|\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ c_{k-3}\hfill & c_{k-2}\hfill & \dots \hfill & c_{2k-5}\hfill \\ {\displaystyle \sum _{\ell =1}^m{\lambda }_{\ell }^{k-1}}\hfill & {\displaystyle \sum _{\ell =1}^m{\lambda }_{\ell }^k}\hfill & \dots \hfill & {\displaystyle \sum _{\ell =1}^m{\lambda }_{\ell }^{2k-3}}\hfill \end{array}\right|=\left|\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ c_{k-3}\hfill & c_{k-2}\hfill & \dots \hfill & c_{2k-5}\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-3}\hfill \end{array}\right|=-h_{k-1,1}. \end{gathered}$$

If $h_{k-1,1}=0$ then the obtained linear system of equalities with respect to $R_0,\dots ,R_{k-3}$, namely

$c_i{R}_0+c_{i+1}{R}_1+\cdots +c_{k+i-3}{R}_{k-3}=0\mathrm{for}i\in \{0,\dots ,k-3\}$

implies that ${\{R_j=0\}}_{j=0}^{k-3}$ (since the determinant of the system equals $H_{k-2}\ne 0$). This proves (28). For the case $h_{k-1,1}\ne 0$, we unify all the obtained relationships with (34) and compose the linear system with respect to $R_0,\dots ,R_{k-3},1$. Since it is consistent, its determinant should vanish:

$\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill & 0\\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill & 0\\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\\ c_{k-3}\hfill & c_{k-2}\hfill & \dots \hfill & c_{2k-6}\hfill & -h_{k-1,1}{Q}_2\\ 1\hfill & x\hfill & \dots \hfill & x^{k-3}\hfill & -R\left(x\right)\end{array}\right|=0.$

Expansion of the determinant by its last column and the use of the Formula (33) lead to

$H_{k-2}R\left(x\right)\equiv \frac{h_{k-1,1}{H}_k}{H_{k-2}}{\mathcal{H}}_{k-3}\left(x\right)$

which completes the proof of (30).  □

Remark 1.

Formulas of Theorem 7 allow one to organize recursive computation for ${\mathcal{H}}_k\left(x\right)$ if the polynomials ${\mathcal{H}}_{k-2}\left(x\right)$ and ${\mathcal{H}}_{k-3}\left(x\right)$ are already computed. The involved factors, such as $h_{k-2,1}$, $h_{k-2,2}$, $h_{k-1,1}$ and $h_{k1}$, are either treated known as the coefficients of ${\mathcal{H}}_{k-2}\left(x\right),{\mathcal{H}}_{k-3}\left(x\right)$, or can be calculated by Formulas (27). The only exception is the value for $h_{k2}$. For its computation, we suggest the following original result.

Theorem 8.

Under the assumptions $H_{k-2}\ne 0,H_{k-1}=0$, one has:

(36) $$\begin{gathered} h_{k2}=-\frac{1}{H_{k-2}}{\left|\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ c_{k-3}\hfill & c_{k-2}\hfill & \dots \hfill & c_{2k-5}\hfill \\ c_k\hfill & c_{k+1}\hfill & \dots \hfill & c_{2k-2}\hfill \end{array}\right|}^2-\left|\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-1}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_k\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-3}\hfill \\ c_{k+1}\hfill & c_{k+2}\hfill & \dots \hfill & c_{2k}\hfill \end{array}\right| \\ =-\frac{1}{H_{k-2}}{\left(\sum _{j=0}^{k-2}{h}_{k-2,j}{c}_{2k-j-2}\right)}^2-\sum _{j=1}^{k-1}{h}_{k-1,j}{c}_{2k-j}. \end{gathered}$$

Proof. 

Is based on the equality (25). By replacing there $(k-1)\to k$, we start with

$$\begin{gathered} h_{k2}={\left|\begin{array}{cccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-3}\hfill & c_{k-1}\hfill & c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-2}\hfill & c_k\hfill & c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-5}\hfill & c_{2k-3}\hfill & c_{2k-2}\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-4}\hfill & c_{2k-2}\hfill & c_{2k-1}\hfill \end{array}\right|}_{k\times k} \\ =\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill & \mid c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill & \mid c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-4}\hfill & \mid c_{2k-2}\hfill \\ c_k\hfill & c_{k+1}\hfill & \dots \hfill & c_{2k-2}\hfill & \mid c_{2k}\hfill \end{array}\right|-\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill & c_{k-1}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill & c_k\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-4}\hfill & c_{2k-3}\hfill \\ c_{k+1}\hfill & c_{k+2}\hfill & \dots \hfill & c_{2k-1}\hfill & c_{2k}\hfill \end{array}\right|. \end{gathered}$$

We represent the first determinant with the aid of Sylvester’s identity (5):

(37)$\frac{1}{H_{k-2}}\left(H_{k-1}\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-3}\hfill & \mid c_k\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-2}\hfill & \mid c_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill & ⋮\hfill \\ c_{k-3}\hfill & c_{k-2}\hfill & \dots \hfill & c_{2k-6}\hfill & \mid c_{2k-4}\hfill \\ c_k\hfill & c_{k+1}\hfill & \dots \hfill & c_{2k-4}\hfill & \mid c_{2k-2}\hfill \end{array}\right|-{\left|\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-2}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_{k-1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ c_{k-3}\hfill & c_{k-2}\hfill & \dots \hfill & c_{2k-5}\hfill \\ c_k\hfill & c_{k+1}\hfill & \dots \hfill & c_{2k-2}\hfill \end{array}\right|}^2\right).$

If $H_{k-1}=h_{k-1,1}=0$, then the expansions of the remained determinants by their last rows lead to (8).  □

As for the extension of the results of the present section in the 19th century, we recall here one generalization due to Kronecker [3] and Frobenius [15]. It relates to the general degenerate case of the vanishing of several entries in the sequence of Hankel determinants $H_1,H_2,\dots$. The key point in their treatment is the number of consecutive zeros in this sequence.

Theorem 9

(Frobenius). Let

(38)$H_k\ne 0,H_{k+1}=0,\dots ,H_{k+\mathfrak{r}-1}=0,H_{k+\mathfrak{r}}\ne 0where\mathfrak{r}>2.$

Then

${\mathcal{H}}_{k+1}\left(x\right)\equiv 0,\dots ,{\mathcal{H}}_{k+\mathfrak{r}-2}\left(x\right)\equiv 0$

and

${\mathcal{H}}_{k+\mathfrak{r}-1}\left(x\right)\equiv {(-1)}^{\mathfrak{r}(\mathfrak{r}-1)/2}{\left[\frac{B_{0,\mathfrak{r}-1}}{H_k}\right]}^{\mathfrak{r}-1}{\mathcal{H}}_k\left(x\right),$

while

${\mathcal{H}}_{k+\mathfrak{r}}\left(x\right)\equiv \frac{1}{H_k^{\mathfrak{r}+1}}\left[Q\left(x\right){\mathcal{H}}_k\left(x\right)-{(-1)}^{\mathfrak{r}(\mathfrak{r}-1)/2}{B}_{0,\mathfrak{r}-1}^{\mathfrak{r}+1}{\mathcal{H}}_{k-1}\left(x\right)\right].$

Here

(39) $Q\left(x\right)=\left|\begin{array}{ccccc}{B}_{00}\hfill & B_{01}\hfill & \dots & B_{0,\mathfrak{r}-1}\hfill & \mid D_0\hfill \\ B_{10}\hfill & B_{11}\hfill & \dots & B_{1,\mathfrak{r}-1}\hfill & \mid D_{0}x+D_1\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill \\ B_{\mathfrak{r}-1,0}\hfill & B_{\mathfrak{r}-1,1}\hfill & \dots & B_{\mathfrak{r}-1,\mathfrak{r}-1}\hfill & \mid D_0{x}^{\mathfrak{r}-1}+D_1{x}^{\mathfrak{r}-2}+\cdots +D_{\mathfrak{r}-1}\hfill \\ B_{\mathfrak{r}0}\hfill & B_{\mathfrak{r}1}\hfill & \dots & B_{\mathfrak{r},\mathfrak{r}-1}\hfill & \mid D_0{x}^{\mathfrak{r}}+D_1{x}^{\mathfrak{r}-1}+\cdots +D_{\mathfrak{r}-1}x+D_{\mathfrak{r}}\hfill \end{array}\right|$

with

$B_{ij}=\left|\begin{array}{ccccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-1}\hfill & \mid c_{k+j}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_k\hfill & \mid c_{k+j+1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill & ⋮\hfill \\ c_{k-1}\hfill & c_k\hfill & \dots \hfill & c_{2k-2}\hfill & \mid c_{2k+j-1}\hfill \\ c_{k+i}\hfill & c_{k+i+1}\hfill & \dots \hfill & c_{2k+i-2}\hfill & \mid c_{2k+i+j}\hfill \end{array}\right|,D_i=\left|\begin{array}{cccc}{c}_0\hfill & c_1\hfill & \dots \hfill & c_{k-1}\hfill \\ c_1\hfill & c_2\hfill & \dots \hfill & c_k\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ c_{k-2}\hfill & c_{k-1}\hfill & \dots \hfill & c_{2k-3}\hfill \\ c_{k+i-1}\hfill & c_{k+i-2}\hfill & \dots \hfill & c_{2k+i-2}\hfill \end{array}\right|.$

The $\mathfrak{r}\times \mathfrak{r}$-submatrix of (39) standing in the upper left corner happens to possess the Hankel structure, i.e.,

$B_{i_1{j}_1}=B_{i_2{j}_2}if{i}_1+j_1=i_2+j_2,\{i_1,i_2,j_1,j_2\}\subset \{0,1,\cdots ,\mathfrak{r}-1\},$

with its first generating elements equal to zero: $B_{00}=0,B_{01}=0,\dots ,B_{0,\mathfrak{r}-2}=0$.

If the coefficients of ${\mathcal{H}}_{k-1}\left(x\right)$ and ${\mathcal{H}}_k\left(x\right)$ are known, then the expressions for $D_0$, $D_1,\dots$, $D_{\mathfrak{r}}$ and $B_{0,\mathfrak{r}-1},B_{\mathfrak{r},0}$ can be computed via the linear relationships similar to (27). As for the values $B_{ij}$ with $i+j\ge \mathfrak{r},i>0,j>0$, we are sure that there exist formulas for their evaluation similar to (8), i.e., their explicit representations as combinations of the coefficients of polynomials ${\mathcal{H}}_{k-1}\left(x\right)$ and ${\mathcal{H}}_k\left(x\right)$. We failed to find them in classical papers, we also have not yet succeeded in deducing them ourselves. Our confidence is based on the coincidence of the results of Theorems 7 and 9 with the algorithm known as the Berlekamp–Massey algorithm [19,20]. Massey suggested it for the problem of finding the minimal polynomial of a linear recurrent sequence. The relationship of this problem to that of finding some characteristics of an appropriate Hankel matrix is due to Kronecker and partially discussed in [18]; in the papers [27], the backgrounds of this algorithm in terms of systems of linear equations with Hankel matrices have been discussed, with some of deduced conclusions equivalent to those presented here in terms of Hankel polynomials. The paper [28] extends the orthogonality oriented approach to the rational interpolation problem initiated in [11] to the case of the vanishing of several leading principal minors of the Hankel matrix. The result is very similar to Theorem 8 (though it is obtained by a different technique).

As it is mentioned in the Introduction of the present paper, it looks like the results of Theorems 7 and 9 are lost in the past.

4. Polynomial Interpolation via Symmetric Functions of the Data Set

The classical polynomial interpolation problem consists in finding the polynomial

(40)$p\left(x\right)=p_0{x}^{N-1}+p_1{x}^{N-2}+\dots +p_{N-1}$

satisfying the data set (1), i.e.,

(41)$p\left(x_j\right)=y_j\mathrm{for}j\in \{1,\dots N\}.$

This problem always possesses a unique solution that can be found via resolving the system of linear relationships (41) with the Vandermonde matrix or, in other terms, via the determinant evaluation

(42)$p\left(x\right)\equiv \frac{1}{{\displaystyle \prod _{1\le j<k\le N}(x_k-x_j)}}\left|\begin{array}{cccccc}1& x_1& x_1^2& \dots & x_1^{N-1}& y_1\\ 1& x_2& x_2^2& \dots & x_2^{N-1}& y_2\\ ⋮& ⋮& ⋮& & ⋮& ⋮\\ 1& x_N& x_N^2& \dots & x_N^{N-1}& y_N\\ 1& x& x^2& \dots & x^{N-1}& 0\end{array}\right|.$

This computation is usually performed by virtue of some auxiliary constructions like representation of the interpolant in either the Newton or the Lagrange form. The last one stems from the expansion of the determinant from (42) by its last column:

(43)$p\left(x\right)\equiv \sum _{j=1}^N{y}_j\frac{W_j\left(x\right)}{W_j\left(x_j\right)}=\sum _{j=1}^N{y}_j\frac{W_j\left(x\right)}{W^{\prime }\left(x_j\right)}.$

where

(44)$W\left(x\right)=\prod _{j=1}^N(x-x_j)$

and

(45)$W_k\left(x\right)=\frac{W\left(x\right)}{x-x_k}\equiv \prod _{\genfrac{}{}{0pt}{}{j=1}{j\ne k}}^N(x-x_j)\mathrm{for}k\in \{1,\dots ,N\}.$

The following easily deduced equality will be used in some further proofs:

(46)$\prod _{j=1}^N{W}_j\left(x_j\right)=\prod _{j=1}^N{W}^{\prime }\left(x_j\right)={(-1)}^{N(N-1)/2}\prod _{1\le j<k\le N}{(x_k-x_j)}^2.$

It should be noted that the representation (43) does not immediately provide one with the canonical form for the interpolant, i.e., an explicit expression for its coefficients. In order to extract them from this formula, let us prove first a preliminary result

Theorem 10.

Let $G\left(x\right)\equiv G_0{x}^{N-1}+G_1{x}^{N-2}+\cdots +G_{N-1}$ be a polynomial of a degree at most $N-1$. Then the following equalities are valid

(47) $\sum _{j=1}^N\frac{G\left(x_j\right)}{W^{\prime }\left(x_j\right)}=\left\{\begin{array}{cc}0\hfill & ifdegG<N-1;\hfill \\ G_0\hfill & ifdegG=N-1.\hfill \end{array}\right.$

Proof. 

Use the following Lagrange formula:

$\frac{G\left(x\right)}{W\left(x\right)}\equiv \frac{G\left(x_1\right)}{W^{\prime }\left(x_1\right)(x-x_1)}+\frac{G\left(x_2\right)}{W^{\prime }\left(x_2\right)(x-x_2)}+\cdots +\frac{G\left(x_N\right)}{W^{\prime }\left(x_N\right)(x-x_N)}$

provided that ${\{G\left(x_j\right)\ne 0\}}_{j=1}^N$. From this the expansion of the fraction $G\left(x\right)/W\left(x\right)$ in the Laurent series in negative powers of x can be derived:

$\frac{G\left(x\right)}{W\left(x\right)}\equiv \left(\sum _{j=1}^N\frac{G\left(x_j\right)}{W^{\prime }\left(x_j\right)}\right)\frac{1}{x}+\left(\sum _{j=1}^N\frac{G\left(x_j\right)x_j}{W^{\prime }\left(x_j\right)}\right)\frac{1}{x^2}+\cdots +\left(\sum _{j=1}^N\frac{G\left(x_j\right)x_j^k}{W^{\prime }\left(x_j\right)}\right)\frac{1}{x^{k+1}}+\dots$

On the other hand, multiplying the formal expansion

$\frac{G\left(x\right)}{W\left(x\right)}\equiv \frac{c_0}{x}+\frac{c_1}{x^2}+\cdots +\frac{c_k}{x^{k+1}}+\dots$

by $W\left(x\right)\equiv x^N+w_1{x}^{N-1}+\cdots +w_N$, one obtains the following formulas for determining the coefficients $c_0,c_1,\dots$ recursively:

$G_0=c_0,G_1=c_0{w}_1+c_1,G_2=c_0{w}_2+c_1{w}_1+c_2,\dots$

Comparing the two obtained forms for the expansion of $G\left(x\right)/W\left(x\right)$ yields the claimed equalities (47).  □

Setting $G\left(x\right)\equiv x^k$ in (47) results in

Corollary 1.

The following Euler–Lagrange equalities are valid:

(48) $\sum _{j=1}^N\frac{x_j^k}{W^{\prime }\left(x_j\right)}=\left\{\begin{array}{cc}0\hfill & ifk\in \{0,\dots ,N-2\};\hfill \\ 1\hfill & ifk=N-1.\hfill \end{array}\right.$

Corollary 2.

Let

$G\left(x\right)\equiv G_0{x}^M+G_1{x}^{M-1}+\cdots +G_M,F\left(x\right)\equiv F_0{x}^n+F_1{x}^{n-1}+\cdots +F_n,$

$G_0\ne 0,F_0\ne 0$, and $F\left(x\right)$ possesses only simple zeros ${\lambda }_1,\dots ,{\lambda }_n$ not coinciding with $x_1,\dots ,x_N$. Then the following equalities are valid

(49) $\sum _{j=1}^N\frac{G\left(x_j\right)x_j^k}{F\left(x_j\right)W^{\prime }\left(x_j\right)}=\left\{\begin{array}{c}{\displaystyle -\sum _{\ell =1}^n\frac{G\left({\lambda }_{\ell }\right){\lambda }_{\ell }^k}{F^{\prime }\left({\lambda }_{\ell }\right)W\left({\lambda }_{\ell }\right)}ifk\in \{0,\dots ,N+n-M-2\};}\\ {\displaystyle \frac{G_0}{F_0}-\sum _{\ell =1}^n\frac{G\left({\lambda }_{\ell }\right){\lambda }_{\ell }^{N+n-M-1}}{F^{\prime }\left({\lambda }_{\ell }\right)W\left({\lambda }_{\ell }\right)}ifk=N+n-M-1.}\end{array}\right.$

Theorem 11.

Calculate two sequences of values:

(50) ${\sigma }_k=\sum _{j=1}^N\frac{x_j^{N+k-1}}{W^{\prime }\left(x_j\right)}fork\in \{1,\dots ,N\}$

and

(51) ${\tau }_k=\sum _{j=1}^N{y}_j\frac{x_j^k}{W^{\prime }\left(x_j\right)}fork\in \{0,\dots N-1\}.$

The following recursive formulas connect the values (50) and (51) with the coefficients of the interpolation polynomial:

${\tau }_0=p_0,{\tau }_k=p_0{\sigma }_k+p_1{\sigma }_{k-1}+\cdots +p_{k-1}{\sigma }_1+p_{k}fork\in \{1,\dots ,N-1\}.$

Proof. 

Due to equalities (48) one has

$\begin{array}{ccc}\hfill {\tau }_k& =& {\displaystyle \sum _{j=1}^N\frac{x_j^k{y}_j}{W^{\prime }\left(x_j\right)}=\sum _{j=1}^N\frac{(p_0{x}_j^{N-1}+p_1{x}_j^{N-2}+\cdots +p_k{x}_j^{N-k-1}+\cdots +p_{N-1})x_j^k}{W^{\prime }\left(x_j\right)}}\hfill \\ & =& p_0{\sigma }_k+p_1{\sigma }_{k-1}+\cdots +p_{k-1}{\sigma }_1+p_k.\hfill \end{array}$

 □

Remark 2.

In the 19th century publications, the expression (50) treated as a symmetric function of $x_1,\dots ,x_N$ was known as alephfunction ((Germ.) Alephfunktion ) and was denoted by ${\aleph }_k(x_1,\dots ,x_N)$ [29]. Despite of its rationally looking definition, it happens to be a complete homogeneous symmetric polynomial, i.e., it equals the sum of all possible monomials of the degree k in the variables $x_1,\dots ,x_N$. Thus, for instance,

${\sigma }_2(x_1,x_2,x_3)\equiv x_1^2+x_2^2+x_3^2+x_1{x}_2+x_1{x}_3+x_2{x}_3.$

The expression ${\tau }_k$ defined by (51) can also be treated as a symmetric function of pairs of variables $(x_1,y_1),(x_2,y_2),\dots ,(x_N,y_N)$, i.e., such a function whose value is unchanged when any of the pairs are interchanged. Although the present authors suppose that the result of Theorem 11 is contained in an old textbook in Algebra, this source remains still hidden from them.

We now turn to an alternative representation for the interpolant.

Theorem 12.

Assume that $y_j\ne 0$ for $j\in \{1,\dots N\}$. Calculate the sequence of values:

(52) ${\tilde{\tau }}_k=\sum _{j=0}^N\frac{1}{y_j}\frac{x_j^k}{W^{\prime }\left(x_j\right)}fork\in \{0,\dots ,2N-2\}.$

The interpolation polynomial can be represented in the Hankel polynomial form:

(53) $p\left(x\right)={(-1)}^{N(N-1)/2}\left(\prod _{j=1}^N{y}_j\right)\underset{{\mathcal{H}}_{N-1}(x;\left\{\tilde{\tau }\right\})}{\underbrace{\left|\begin{array}{ccccc}{\tilde{\tau }}_0\hfill & {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & \dots \hfill & {\tilde{\tau }}_{N-1}\hfill \\ {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & {\tilde{\tau }}_3\hfill & \dots \hfill & {\tilde{\tau }}_N\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ {\tilde{\tau }}_{N-2}\hfill & {\tilde{\tau }}_{N-1}\hfill & {\tilde{\tau }}_N\hfill & \dots \hfill & {\tilde{\tau }}_{2N-3}\hfill \\ 1\hfill & x\hfill & x^2\hfill & \dots \hfill & x^{N-1}\hfill \end{array}\right|}}.$

Proof. 

Assume that the interpolation polynomial (40) exists. Rewrite the equalities (41) in the form

(54)$p_{N-1}\frac{1}{y_j}+p_{N-2}\frac{x_j}{y_j}+\cdots +p_1\frac{x_j^{N-2}}{y_j}+p_0\frac{x_j^{N-1}}{y_j}=1\mathrm{for}j\in \{1,\dots ,N\}.$

Multiply each of these equalities by the corresponding multiplier $1/W^{\prime }\left(x_j\right)$ and sum the obtained results. Due to (48), one gets

$p_{N-1}{\tilde{\tau }}_0+p_{N-2}{\tilde{\tau }}_1+\cdots +p_1{\tilde{\tau }}_{N-2}+p_0{\tilde{\tau }}_{N-1}=0.$

Similar equalities result from multiplication of (54) by ${\{x_j^k/W^{\prime }\left(x_j\right)\}}_{k=1}^{N-2}$:

$\begin{array}{cccccc}{p}_{N-1}{\tilde{\tau }}_1\hfill & +p_{N-2}{\tilde{\tau }}_2\hfill & +\dots \hfill & +p_1{\tilde{\tau }}_{N-1}\hfill & +p_0{\tilde{\tau }}_N\hfill & =0,\\ p_{N-1}{\tilde{\tau }}_2\hfill & +p_{N-2}{\tilde{\tau }}_3\hfill & +\dots \hfill & +p_1{\tilde{\tau }}_N\hfill & +p_0{\tilde{\tau }}_{N+1}\hfill & =0,\\ \cdots \hfill & & & & & \cdots ,\\ p_{N-1}{\tilde{\tau }}_{N-2}\hfill & +p_{N-2}{\tilde{\tau }}_{N-1}\hfill & +\dots \hfill & +p_1{\tilde{\tau }}_{2N-4}\hfill & +p_0{\tilde{\tau }}_{2N-3}\hfill & =0.\end{array}$

Linear combination of (54) multiplied by $x_j^{N-1}/W^{\prime }\left(x_j\right)$ yields something different:

$p_{N-1}{\tilde{\tau }}_{N-1}+p_{N-2}{\tilde{\tau }}_N+\cdots +p_1{\tilde{\tau }}_{2N-3}+p_0{\tilde{\tau }}_{2N-2}=1.$

Unifying all the obtained equalities with (40), one obtains a system of linear equations with respect to $p_{N-1},p_{N-2},\dots ,p_0$. If the interpolant exists, then it should satisfy the following identity:

$\left|\begin{array}{cccccc}{\tilde{\tau }}_0\hfill & {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & \dots \hfill & {\tilde{\tau }}_{N-1}\hfill & 0\hfill \\ {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & {\tilde{\tau }}_3\hfill & \dots \hfill & {\tilde{\tau }}_N\hfill & 0\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & & ⋮\hfill & ⋮\hfill \\ {\tilde{\tau }}_{N-2}\hfill & {\tilde{\tau }}_{N-1}\hfill & {\tilde{\tau }}_N\hfill & \dots \hfill & {\tilde{\tau }}_{2N-3}\hfill & 0\hfill \\ {\tilde{\tau }}_{N-1}\hfill & {\tilde{\tau }}_N\hfill & {\tilde{\tau }}_{N+1}\hfill & \dots \hfill & {\tilde{\tau }}_{2N-2}\hfill & 1\hfill \\ 1\hfill & x\hfill & x^2\hfill & \dots \hfill & x^{N-1}\hfill & p\left(x\right)\hfill \end{array}\right|\equiv 0.$

Expansion of the determinant by its last column gives

(55)$p\left(x\right)\equiv \frac{{\mathcal{H}}_{N-1}(x,\left\{\tilde{\tau }\right\})}{H_N\left(\left\{\tilde{\tau }\right\}\right)}$

provided that $H_N\left(\left\{\tilde{\tau }\right\}\right)\ne 0$. To prove the latter, represent the determinant as the product:

(56)$\begin{array}{c}\hfill H_N\left(\left\{\tilde{\tau }\right\}\right)=\left|\begin{array}{cccc}1& 1& \dots & 1\\ x_1& x_2& \dots & x_N\\ x_1^2& x_2^2& \dots & x_N^2\\ ⋮& ⋮& & ⋮\\ x_1^{N-1}& x_2^{N-1}& \dots & x_N^{N-1}\end{array}\right|·\left|\begin{array}{cccc}\frac{1}{y_1{W}^{\prime }\left(x_1\right)}& 0& \dots & 0\\ 0& \frac{1}{y_2{W}^{\prime }\left(x_2\right)}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \frac{1}{y_N{W}^{\prime }\left(x_N\right)}\end{array}\right|\\ \hfill \times \left|\begin{array}{ccccc}1\hfill & x_1\hfill & x_1^2\hfill & \dots \hfill & x_1^{N-1}\hfill \\ 1\hfill & x_2\hfill & x_2^2\hfill & \dots \hfill & x_2^{N-1}\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ 1\hfill & x_N\hfill & x_N^2\hfill & \dots \hfill & x_N^{N-1}\hfill \end{array}\right|=\frac{{\displaystyle \prod _{1\le j<k\le N}{(x_k-x_j)}^2}}{{\displaystyle \prod _{j=1}^N{y}_j\prod _{j=1}^N{W}^{\prime }\left(x_j\right)}}\stackrel{\left(46\right)}{=}\frac{{(-1)}^{N(N-1)/2}}{{\displaystyle \prod _{j=1}^N{y}_j}}.\end{array}$

Therefore, the denominator in the right-hand side of (55) does not vanish.

One can utilize the last expression for the direct deduction of the validity of the Formula (53) as a solution to the polynomial interpolation problem. Indeed, let us substitute $x=x_1$ into the numerator of the fraction in the right-hand side of (53):

${\mathcal{H}}_{N-1}(x_1;\left\{\tilde{\tau }\right\})={\left|\begin{array}{ccccc}{\tilde{\tau }}_0\hfill & {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & \dots \hfill & {\tilde{\tau }}_{N-1}\hfill \\ {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & {\tilde{\tau }}_3\hfill & \dots \hfill & {\tilde{\tau }}_N\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ {\tilde{\tau }}_{N-2}\hfill & {\tilde{\tau }}_{N-1}\hfill & {\tilde{\tau }}_N\hfill & \dots \hfill & {\tilde{\tau }}_{2N-3}\hfill \\ 1\hfill & x_1\hfill & x_1^2\hfill & \dots \hfill & x_1^{N-1}\hfill \end{array}\right|}_{N\times N}$

$\stackrel{\left(14\right)}{\equiv }{(-1)}^{N+1}{\left|\begin{array}{cccc}{\tilde{\tau }}_1-x_1{\tilde{\tau }}_0\hfill & {\tilde{\tau }}_2-x_1{\tilde{\tau }}_1\hfill & \dots \hfill & {\tilde{\tau }}_{N-1}-x_1{\tilde{\tau }}_{N-2}\hfill \\ {\tilde{\tau }}_2-x_1{\tilde{\tau }}_1\hfill & {\tilde{\tau }}_3-x_1{\tilde{\tau }}_2\hfill & \dots \hfill & {\tilde{\tau }}_N-x_1{\tilde{\tau }}_{N-1}\hfill \\ \cdots \hfill & \cdots \hfill & & \cdots \hfill \\ {\tilde{\tau }}_{N-1}-x_1{\tilde{\tau }}_{N-2}\hfill & {\tilde{\tau }}_N-x_1{\tilde{\tau }}_{N-1}\hfill & \dots \hfill & {\tilde{\tau }}_{2N-3}-x_1{\tilde{\tau }}_{2N-4}\hfill \end{array}\right|}_{(N-1)\times (N-1)}.$

Consider an entry of the last determinant

$\begin{array}{ccc}\hfill {\tilde{\tau }}_k-x_1{\tilde{\tau }}_{k-1}& =& \frac{x_1^k}{y_1{W}^{\prime }\left(x_1\right)}+\frac{x_2^k}{y_2{W}^{\prime }\left(x_2\right)}+\cdots +\frac{x_N^k}{y_N{W}^{\prime }\left(x_N\right)}\hfill \\ & & -\frac{x_1^k}{y_1{W}^{\prime }\left(x_1\right)}-\frac{x_2^{k-1}{x}_1}{y_2{W}^{\prime }\left(x_2\right)}-\cdots -\frac{x_N^{k-1}{x}_1}{y_N{W}^{\prime }\left(x_N\right)}\hfill \\ & =& \frac{x_2^{k-1}(x_2-x_1)}{y_2{W}^{\prime }\left(x_2\right)}+\cdots +\frac{x_N^{k-1}(x_N-x_1)}{y_N{W}^{\prime }\left(x_N\right)}.\hfill \end{array}$

It can be easily verified that

$\frac{x_j-x_1}{W^{\prime }\left(x_j\right)}=\frac{1}{W_1^{\prime }\left(x_j\right)}\mathrm{for}j\in \{2,\dots ,N\}\mathrm{and}{W}_1\left(x\right)\stackrel{\left(45\right)}{\equiv }\prod _{k=2}^N(x-x_k).$

Thus, introducing the sequence

(57)${\tilde{T}}_k={\tilde{\tau }}_k-x_1{\tilde{\tau }}_{k-1}=\sum _{j=2}^N\frac{x_j^k}{y_j{W}_1^{\prime }\left(x_j\right)}\mathrm{for}k\in \{0,\dots ,2N-4\}$

one gets

(58)${\mathcal{H}}_{N-1}(x_1;\left\{\tilde{\tau }\right\})={(-1)}^{N-1}\left|\begin{array}{cccc}{\tilde{T}}_0\hfill & {\tilde{T}}_1\hfill & \dots \hfill & {\tilde{T}}_{N-2}\hfill \\ {\tilde{T}}_1\hfill & {\tilde{T}}_2\hfill & \dots \hfill & {\tilde{T}}_{N-1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ {\tilde{T}}_{N-2}\hfill & {\tilde{T}}_{N-1}\hfill & \dots \hfill & {\tilde{T}}_{2N-4}\hfill \end{array}\right|={(-1)}^{N-1}H\left(\tilde{T}\right).$

The last determinant can be evaluated by Formula (56):

${\mathcal{H}}_{N-1}(x_1;\left\{\tilde{\tau }\right\})={(-1)}^{N-1}{(-1)}^{(N-1)(N-2)/2}/\prod _{j=2}^N{y}_j.$

Combining this equality with (56), one gets in (55): $p\left(x_1\right)=y_1$. Other nodes $x_j$ are dealt with similarly.  □

At first glance, Formula (53) seems to have no advantage over the algorithm from Theorem 11. We dispel doubts in its practical utility in Section 9. The present section is concluded by a technical result.

Theorem 13.

Under the condition of Theorem 12, one has

(59) ${\mathcal{H}}_N(x;\left\{\tilde{\tau }\right\})\equiv H_N\left(\left\{\tilde{\tau }\right\}\right)\prod _{j=1}^N(x-x_j)\equiv \frac{{(-1)}^{N(N-1)/2}}{{\prod }_{j=1}^N{y}_j}\prod _{j=1}^N(x-x_j).$

Proof. 

Let us prove that

(60)${\mathcal{H}}_N(x_j;\left\{\tilde{\tau }\right\})=0,\mathrm{for}j\in \{1,\dots ,N\}.$

The linear combination

$\sum _{j=1}^N\frac{x_j^k}{W^{\prime }\left(x_j\right)}{\mathcal{H}}_N(x_j;\left\{\tilde{\tau }\right\})=\left|\begin{array}{cccc}{\tilde{\tau }}_0\hfill & {\tilde{\tau }}_1\hfill & \dots \hfill & {\tilde{\tau }}_N\hfill \\ {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & \dots \hfill & {\tilde{\tau }}_{N+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ {\tilde{\tau }}_{N-1}\hfill & {\tilde{\tau }}_N\hfill & \dots \hfill & {\tilde{\tau }}_{2N-1}\hfill \\ {\tilde{\tau }}_k\hfill & {\tilde{\tau }}_{k+1}\hfill & \dots \hfill & {\tilde{\tau }}_{N+k}\hfill \end{array}\right|.$

equals zero for $k\in \{0,\dots ,N-1\}$. Therefore we have N linear homogeneous equalities

$\sum _{j=1}^N\frac{x_j^{k-1}}{W^{\prime }\left(x_j\right)}{\mathcal{H}}_N(x_j;\left\{\tilde{\tau }\right\})=0,j\in \{1,\dots ,N\}$

valid for N values ${\left\{{\mathcal{H}}_N(x_j;\left\{\tilde{\tau }\right\})\right\}}_{j=1}^N$. Since

$det{\left[\frac{x_j^{k-1}}{W^{\prime }\left(x_j\right)}\right]}_{j,k=1}^N\stackrel{\left(46\right)}{=}{(-1)}^{N(N-1)/2}\prod _{1\le j<k\le N}(x_k-x_j)\ne 0,$

all the equalities (60) are valid. Thus, we know all the zeros of the polynomial ${\mathcal{H}}_N(x;\left\{\tilde{\tau }\right\})$, while its leading coefficient equals to $H_N\left(\left\{\tilde{\tau }\right\}\right)$. The latter has been already evaluated by (56). This completes the proof.  □

5. Rational Interpolation

We now search for polynomials (6) such that the rational function

(61)$r\left(x\right)=p\left(x\right)/q\left(x\right)$

satisfies the data set (1), i.e.,

(62)$r\left(x_j\right)=y_j\mathrm{for}j\in \{1,\dots ,N=n+m+1\}.$

Hereinafter, we will not distinguish the solutions to the problem with numerator and denominator multiplied by a common numerical factor.

The first solution to the problem was suggested in the form of generalization of the Lagrange representation (43) for the polynomial interpolant.

Theorem 14

(Cauchy, ref. [2]). Denote

$W_{\overline{j_1{j}_2\dots j_{\ell }}}\left(x\right)=\prod _{s=1}^{\ell }(x-x_{j_s})and{W}_{j_1{j}_2\dots j_{\ell }}\left(x\right)=\frac{W\left(x\right)}{W_{\overline{j_1{j}_2\dots j_{\ell }}}\left(x\right)}.$

Rational interpolant (61) can be found via the formulas

$p\left(x\right)=\sum _{(j_1,j_2,\dots ,j_{m+1})}\frac{{\prod }_{s=1}^{m+1}{y}_{j_s}}{{\displaystyle \prod _{s=1}^{m+1}{W}_{j_1{j}_2\dots j_{m+1}}\left(x_{j_s}\right)}}{W}_{j_1{j}_2\dots j_{m+1}}\left(x\right)$

and

$q\left(x\right)={(-1)}^{mn}\sum _{(j_1,j_2,\dots ,j_m)}\frac{{\prod }_{\ell =1}^m{y}_{j_{\ell }}}{{\displaystyle \prod _{j\in \{1,\dots ,N\}\setminus \{j_1,\dots ,j_m\}}{W}_{\overline{j_1{j}_2\dots j_m}}\left(x_j\right)}}{W}_{\overline{j_1{j}_2\dots j_m}}\left(x\right).$

Both sums are taken with respect to all combinations $m+1$ and, respectively, m integers from the set $\{1,\dots ,N\}$.

Being valid generically, Cauchy’s solution fails for some particular choices of the data set (1). Whereas the polynomial interpolation problem always has a solution, the rational interpolation one is not always resolvable. This defect was first discovered by Kronecker [3], and later discussed by Netto [4]. To exemplify this, we first generate from the condition (62) the system of equations

(63)$p\left(x_j\right)=y_{j}q\left(x_j\right)\mathrm{for}j\in \{1,\dots ,N\}$

or, equivalently,

(64)$p_n+p_{n-1}{x}_j+\cdots +p_0{x}_j^n=q_m{y}_j+q_{m-1}{x}_j{y}_j+\cdots +q_0{x}_j^m{y}_j\mathrm{for}j\in \{1,\dots ,N\}$

which is linear with respect to the $N+1$ coefficients of $p\left(x\right)$ and $q\left(x\right)$. The principal solvability of this system can be established with the aid of Linear Algebra methods, like, for instance via Gaussian elimination procedure.

Example 2.

Given the data set

$\begin{array}{ccccccc}x& & -1& 0& 1& 2& 3\\ y& & 1& 1& 1/3& 3& 1/13\end{array}$

find its rational interpolant in the form $r\left(x\right)=p\left(x\right)/q\left(x\right)$ with $degp\left(x\right)=1,degq\left(x\right)=3$.

Solution. Resolving the system (64), one gets the expressions:

$p\left(x\right)\equiv x-2,q\left(x\right)\equiv x^3-x^2-x-2.$

However, $p\left(2\right)=0$ and $q\left(2\right)=0$, and therefore the condition $r\left(2\right)=3$ is not satisfied. It is not satisfied even if we cancel the numerator and denominator by the common linear factor.

Explanation for this phenomenon of the unattainable points consists in the nonequivalence of the passage from (62) to (63) since for some node $x_j$ one might obtain a solution to the linear system (64) satisfying both conditions $p\left(x_j\right)=0$ and $q\left(x_j\right)=0$.

On the other hand, solution to the problem might be not unique.

Example 3.

For the data set

$\begin{array}{ccccccc}x& & -1& 0& 1& 2& 3\\ y& & 1& 1& 1/3& 1/7& 1/13\end{array}$

generated by the rational function $1/(x^2+x+1)$ there exists infinitely many rational interpolants with $degp\left(x\right)=1,degq\left(x\right)=3$ in the form $(x-\lambda )/\left((x-\lambda )(x^2+x+1)\right)$ where $\lambda \notin \{-1,0,1,2,3\}$.

We now pass to an alternative approach to the problem, due to Jacobi.

Theorem 15.

Let $y_j\ne 0$ for $j\in \{1,\dots ,N\}$. Compute the values

(65) ${\tau }_k=\sum _{j=1}^N{y}_j\frac{x_j^k}{W^{\prime }\left(x_j\right)}fork\in \{0,\dots ,2m\}$

and

(66) ${\tilde{\tau }}_k=\sum _{j=1}^N\frac{1}{y_j}\frac{x_j^k}{W^{\prime }\left(x_j\right)}fork\in \{0,\dots ,2n-1\},$

and generate the corresponding Hankel polynomials ${\mathcal{H}}_m(x;\left\{\tau \right\})$ and ${\mathcal{H}}_n(x;\left\{\tilde{\tau }\right\})$. If

(67) $H_n\left(\left\{\tilde{\tau }\right\}\right)\ne 0$

and

(68) ${\mathcal{H}}_m(x_j;\left\{\tau \right\})\ne 0forj\in \{1,\dots ,N\}$

then there exists a unique rational interpolant with $degp\left(x\right)=n,degq\left(x\right)\le m=N-n-1$ which can be expressed as:

(69) $$\begin{gathered} p\left(x\right)=H_{m+1}\left(\left\{\tau \right\}\right){\mathcal{H}}_n(x;\left\{\tilde{\tau }\right\}) \\ =\left|\begin{array}{cccc}{\tau }_0\hfill & {\tau }_1\hfill & \dots \hfill & {\tau }_m\hfill \\ {\tau }_1\hfill & {\tau }_2\hfill & \dots \hfill & {\tau }_{m+1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ {\tau }_{m-1}\hfill & {\tau }_m\hfill & \dots \hfill & {\tau }_{2m-1}\hfill \\ {\tau }_m\hfill & {\tau }_{m+1}\hfill & \dots \hfill & {\tau }_{2m}\hfill \end{array}\right|·\left|\begin{array}{cccc}{\tilde{\tau }}_0\hfill & {\tilde{\tau }}_1\hfill & \dots \hfill & {\tilde{\tau }}_n\hfill \\ {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & \dots \hfill & {\tilde{\tau }}_{n+1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ {\tilde{\tau }}_{n-1}\hfill & {\tilde{\tau }}_n\hfill & \dots \hfill & {\tilde{\tau }}_{2n-1}\hfill \\ 1\hfill & x\hfill & \dots \hfill & x^n\hfill \end{array}\right|, \end{gathered}$$

(70) $$\begin{gathered} q\left(x\right)=H_n\left(\left\{\tilde{\tau }\right\}\right){\mathcal{H}}_m(x;\left\{\tau \right\}) \\ =\left|\begin{array}{cccc}{\tilde{\tau }}_0\hfill & {\tilde{\tau }}_1\hfill & \dots \hfill & {\tilde{\tau }}_{n-1}\hfill \\ {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & \dots \hfill & {\tilde{\tau }}_n\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ {\tilde{\tau }}_{n-1}\hfill & {\tilde{\tau }}_n\hfill & \dots \hfill & {\tilde{\tau }}_{2n-2}\hfill \end{array}\right|·\left|\begin{array}{cccc}{\tau }_0\hfill & {\tau }_1\hfill & \dots \hfill & {\tau }_m\hfill \\ {\tau }_1\hfill & {\tau }_2\hfill & \dots \hfill & {\tau }_{m+1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ {\tau }_{m-1}\hfill & {\tau }_m\hfill & \dots \hfill & {\tau }_{2m-1}\hfill \\ 1\hfill & x\hfill & \dots \hfill & x^m\hfill \end{array}\right|. \end{gathered}$$

Proof. 

Can be found in [11], and we refer here only with the underlying idea leading to the Hankel polynomial appearance. If a desired interpolant exists, then the equalities (64) are valid. Multiply the jth equality by $x_j^k/W^{\prime }\left(x_j\right)$ for $k\in \{0,\dots ,m-1\}$ and sum the obtained equalities by j. Due to the Euler–Lagrange equalities (48), one arrives at a system of equations

$\left\{\begin{array}{ccccc}{q}_m{\tau }_0\hfill & +q_{m-1}{\tau }_1\hfill & +\dots \hfill & +q_0{\tau }_m\hfill & =0,\hfill \\ q_m{\tau }_1\hfill & +q_{m-1}{\tau }_2\hfill & +\dots \hfill & +q_0{\tau }_{m+1}\hfill & =0,\hfill \\ \dots \hfill & & & & \dots ,\hfill \\ q_m{\tau }_{m-1}\hfill & +q_{m-1}{\tau }_m\hfill & +\dots \hfill & +q_0{\tau }_{2m-1}\hfill & =0.\hfill \end{array}\right.$

Therefore, the denominator of the fraction should satisfy the relation

$Aq\left(x\right)\equiv \left|\begin{array}{cccc}{\tau }_0\hfill & {\tau }_1\hfill & \dots \hfill & {\tau }_m\hfill \\ {\tau }_1\hfill & {\tau }_2\hfill & \dots \hfill & {\tau }_{m+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ {\tau }_{m-1}\hfill & {\tau }_m\hfill & \dots \hfill & {\tau }_{2m-1}\hfill \\ 1\hfill & x\hfill & \dots \hfill & x^m\hfill \end{array}\right|$

for some constant factor A.

In a similar way, multiplying the equalities

$p_n\frac{1}{y_j}+p_{n-1}\frac{x_j}{y_j}+\cdots +p_0\frac{x_j^n}{y_j}=q_m+q_{m-1}{x}_j+\cdots +q_0{x}_j^m,j\in \{1,\dots ,N\}$

by $x_j^{\ell }/W^{\prime }\left(x_j\right)$ for $\ell \in \{0,\dots ,n-1\}$ and summarizing by j, one gets the equality for the numerator:

$Bp\left(x\right)\equiv \left|\begin{array}{cccc}{\tilde{\tau }}_0\hfill & {\tilde{\tau }}_1\hfill & \dots \hfill & {\tilde{\tau }}_n\hfill \\ {\tilde{\tau }}_1\hfill & {\tilde{\tau }}_2\hfill & \dots \hfill & {\tilde{\tau }}_{n+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ {\tilde{\tau }}_{n-1}\hfill & {\tilde{\tau }}_n\hfill & \dots \hfill & {\tilde{\tau }}_{2n-1}\hfill \\ 1\hfill & x\hfill & \dots \hfill & x^n\hfill \end{array}\right|\equiv H_n\left(\left\{\tilde{\tau }\right\}\right)x^n+\dots$

for some constant factor B. Due to assumption (67), $B\ne 0$ and $degp\left(x\right)=n$. This can be interpreted as a polynomial interpolation problem for the data set ${\left\{(x_j,1/y_j)\right\}}_{j=1}^N$. Thereby the rational interpolation problem is reduced to the polynomial interpolation problem discussed in Section 4.  □

Corollary 1.

The following relationship is valid:

(71) $H_n\left(\left\{\tilde{\tau }\right\}\right)H_m\left(\left\{\tau \right\}\right)=H_{n+1}\left(\left\{\tilde{\tau }\right\}\right)H_{m+1}\left(\left\{\tau \right\}\right).$

Example 4.

Given the data set

$\begin{array}{ccccccccc}x& & -3& -2& -1& 0& 1& 2& 3\\ y& & 3& 1& 4& 1& 5& 9& 3\end{array}$

find all rational functions $r\left(x\right)=p\left(x\right)/q\left(x\right)$ with $degp\left(x\right)+degq\left(x\right)\le 6$ satisfying it.

Solution. Since we do not know a priori the degrees of the numerator and the denominator of $r\left(x\right)$, we have to compute the sequences (65) and (66) for the maximal possible number of elements:

${\tau }_0=\frac{61}{720},{\tau }_1=-\frac{9}{80},{\tau }_2=-\frac{17}{240},\dots ,{\tau }_{12}=\frac{981007}{240};$

${\tilde{\tau }}_0=-\frac{77}{2880},{\tilde{\tau }}_1=\frac{119}{8640},{\tilde{\tau }}_2=-\frac{167}{8640},\dots ,{\tilde{\tau }}_{12}=\frac{3923929}{8640}.$

Now, compute the Hankel polynomials of the first and the second order:

${\mathcal{H}}_1(x;\left\{\tau \right\})=\frac{61}{720}x+\frac{9}{80},{\mathcal{H}}_2(x;\left\{\tau \right\})=\underset{h_{2,0}}{\underbrace{-\frac{403}{21600}}}{x}^2+\underset{h_{2,1}}{\underbrace{\frac{37}{720}}}x+\underset{h_{2,2}}{\underbrace{\frac{379}{7200}}}.$

Computation of ${\mathcal{H}}_3(x;\left\{\tau \right\})$ can be organized with the aid of the JJ-identity (24):

${\mathcal{H}}_3(x;\left\{\tau \right\})\equiv -{\left(\frac{h_{3,0}}{h_{2,0}}\right)}^2{\mathcal{H}}_1(x;\left\{\tau \right\})+\frac{h_{3,0}}{h_{2,0}}\left(x-\frac{h_{2,1}}{h_{2,0}}+\frac{h_{3,1}}{h_{3,0}}\right){\mathcal{H}}_2(x;\left\{\tau \right\})$

where all the constants are already known except for $h_{3,0}=H_3\left(\left\{\tau \right\}\right)$ and $h_{3,1}$. To find the latter, utilize the equalities (27)

$h_{3,0}=H_3\left(\left\{\tau \right\}\right)={\tau }_4{h}_{2,0}+{\tau }_3{h}_{2,1}+{\tau }_2{h}_{2,2}=-\frac{1379}{64800},$

$h_{3,1}=-({\tau }_5{h}_{2,0}+{\tau }_4{h}_{2,1}+{\tau }_3{h}_{2,2})=\frac{127}{10800}.$

Therefore,

${\mathcal{H}}_3(x;\left\{\tau \right\})\equiv -\frac{1379}{64800}{x}^3+\frac{127}{10800}{x}^2+\frac{5111}{64800}x-\frac{17}{1200}.$

Continuing the recursive utilization of the JJ-identity (24), we get further:

$\begin{array}{ccc}\hfill {\mathcal{H}}_4(x;\left\{\tau \right\})& \equiv & -\frac{1609}{21600}{x}^4-\frac{347}{10800}{x}^3+\frac{368}{675}{x}^2-\frac{2201}{5400}x-\frac{111}{160},\hfill \\ \hfill \\ \hfill {\mathcal{H}}_5(x;\left\{\tau \right\})& \equiv & -\frac{763}{1440}{x}^5+\frac{301}{144}{x}^4+\frac{1715}{288}{x}^3-\frac{2947}{144}{x}^2+\frac{301}{80}x+\frac{357}{16},\hfill \\ \hfill \\ \hfill {\mathcal{H}}_6(x;\left\{\tau \right\})& \equiv & \underset{h_{6,0}}{\underbrace{\frac{693}{16}}}{x}^6\underset{h_{6,1}}{\underbrace{-\frac{357}{16}}}{x}^5-\frac{9201}{16}{x}^4+\frac{3489}{16}{x}^3+\frac{7149}{4}{x}^2-\frac{621}{4}x\underset{h_{6,6}}{\underbrace{-1620}}.\hfill \end{array}$

We need one extra computation, namely

$H_7\left(\left\{\tau \right\}\right)={\tau }_{12}{h}_{6,0}+{\tau }_{11}{h}_{6,1}+\cdots +{\tau }_6{h}_{6,6}=-1620.$

Thus, all the potential denominators of the interpolation fractions are computed. In parallel, similar recursive procedure can be organized for the numerator’s computation:

${\mathcal{H}}_1(x;\left\{\tilde{\tau }\right\})=-\frac{77}{2880}x-\frac{119}{8640},{\mathcal{H}}_2(x;\left\{\tilde{\tau }\right\})=\underset{{\tilde{h}}_{2,0}=H_2\left(\left\{\tilde{\tau }\right\}\right)}{\underbrace{\frac{763}{2332800}}}{x}^2+\underset{{\tilde{h}}_{2,1}}{\underbrace{\frac{301}{233280}}}x+\underset{{\tilde{h}}_{2,2}}{\underbrace{\frac{37}{86400}}},$

${\tilde{h}}_{3,0}=H_3\left(\left\{\tilde{\tau }\right\}\right)={\tilde{\tau }}_4{\tilde{h}}_{2,0}+{\tilde{\tau }}_3{\tilde{h}}_{2,1}+{\tilde{\tau }}_2{\tilde{h}}_{2,2}=\frac{1609}{34992000},$

${\tilde{h}}_{3,1}=-({\tilde{\tau }}_5{\tilde{h}}_{2,0}+{\tilde{\tau }}_4{\tilde{h}}_{2,1}+{\tilde{\tau }}_3{\tilde{h}}_{2,2})=-\frac{347}{17496000},$

$\begin{array}{ccc}\hfill {\mathcal{H}}_3(x;\left\{\tilde{\tau }\right\})& \equiv & -{\left(\frac{{\tilde{h}}_{3,0}}{{\tilde{h}}_{2,0}}\right)}^2{\mathcal{H}}_1(x;\left\{\tilde{\tau }\right\})+\frac{{\tilde{h}}_{3,0}}{{\tilde{h}}_{2,0}}\left(x-\frac{{\tilde{h}}_{2,1}}{{\tilde{h}}_{2,0}}+\frac{{\tilde{h}}_{3,1}}{{\tilde{h}}_{3,0}}\right){\mathcal{H}}_2(x;\left\{\tilde{\tau }\right\})\hfill \\ & \equiv & \frac{1609}{34992000}{x}^3-\frac{347}{17496000}{x}^2-\frac{7181}{34992000}x+\frac{17}{1944000},\hfill \\ \hfill {\mathcal{H}}_4(x;\left\{\tilde{\tau }\right\})& \equiv & \frac{1379}{104976000}{x}^4+\frac{127}{17496000}{x}^3-\frac{4771}{52488000}{x}^2-\frac{13}{81000}x-\frac{379}{11664000},\hfill \\ \hfill {\mathcal{H}}_5(x;\left\{\tilde{\tau }\right\})& \equiv & \frac{403}{34992000}{x}^5+\frac{37}{1166400}{x}^4-\frac{707}{6998400}{x}^3-\frac{13}{43200}{x}^2-\frac{13}{72000}x-\frac{1}{14400},\hfill \\ \hfill {\mathcal{H}}_6(x;\left\{\tilde{\tau }\right\})& \equiv & -\frac{61}{1166400}{x}^6+\frac{1}{14400}{x}^5+\frac{181}{233280}{x}^4-\frac{17}{25920}{x}^3-\frac{841}{291600}{x}^2+\frac{1}{3600}x-\frac{1}{1620}.\hfill \end{array}$

Now we are able to compose the set of rational interpolants:

$\begin{array}{ccc}\hfill r_{0,6}\left(x\right)& =& \frac{H_7\left(\left\{\tau \right\}\right)}{{\mathcal{H}}_6(x;\left\{\tau \right\})}\hfill \\ & \equiv & -\frac{8640}{231x^6-119x^5-3067x^4+1163x^3+9532x^2-828x-8640},\hfill \\ \hfill r_{1,5}\left(x\right)& =& \frac{h_{6,0}{\mathcal{H}}_1(x;\left\{\tilde{\tau }\right\})}{{\tilde{h}}_{1,0}{\mathcal{H}}_5(x;\left\{\tau \right\})}\equiv -\frac{270(33x+17)}{109x^5-430x^4-1225x^3+4210x^2-774x-4590},\hfill \\ \hfill r_{2,4}\left(x\right)& =& \frac{h_{5,0}{\mathcal{H}}_2(x;\left\{\tilde{\tau }\right\})}{{\tilde{h}}_{2,0}{\mathcal{H}}_4(x;\left\{\tau \right\})}\equiv \frac{15(763x^2+3010x+999)}{1609x^4+694x^3-11776x^2+8804x+14985},\cdots ,\hfill \end{array}$

and finally the polynomial interpolant

$r_{6,0}\left(x\right)=\frac{h_{1,0}{\mathcal{H}}_6(x;\left\{\tilde{\tau }\right\})}{{\tilde{h}}_{6,0}}\equiv \frac{61}{720}{x}^6-\frac{9}{80}{x}^5-\frac{181}{144}{x}^4+\frac{17}{16}{x}^3+\frac{841}{180}{x}^2-\frac{9}{20}x+1.$

Formulation of Theorem 15 is due to the authors of [11]. Jacobi did not bother himself in [9] with the questions of existence or uniqueness of a solution to the interpolation problem. He just only suggested that the denominator of the (potential candidate) rational interpolant can be represented in the form of the Hankel polynomial ${\mathcal{H}}_m(x;\left\{\tau \right\})$. On its computation, the rational interpolation problem is reduced to the polynomial interpolation one for the numerator of the fraction. Jacobi did not care himself either on computational optimization of the algorithm like the one outlined in the solution to Example 4 exploiting their own Formula (24).

We next take a closer look at this JJ-identity. The linear relation of the type (4) connecting three consecutive Hankel polynomials generated by an arbitrary sequence reminds the similar looking recurrence identity for the classical orthogonal polynomials [30,31]:

(72)$P_k\left(x\right)\equiv (A_{k}x+B_k)P_{k-1}\left(x\right)+C_k{P}_{k-2}\left(x\right),k\in \{2,3,\dots \}$

for some real constants $A_k,B_k,C_k$. The proof of this identity is traditionally carried out using the orthogonality of the system ${\left\{P_j\left(x\right)\right\}}_{j=0}^k$. For instance, in case of monic polynomials $P_j\left(x\right)$, the following formula for evaluation of $B_k$ is suggested

(73)$B_k=-⟨x{P}_{k-1}\left(x\right),P_{k-1}\left(x\right)⟩/⟨P_{k-1}\left(x\right),P_{k-1}\left(x\right)⟩$

where the inner product for a pair of real polynomials $P\left(x\right),Q\left(x\right)$ is defined via

$⟨P\left(x\right),Q\left(x\right)⟩={\int }_a^{b}P\left(t\right)Q\left(t\right)w\left(t\right)dt$

with a positive and integrable on the interval $a<x<b$ function $w\left(x\right)$.

The relationship of orthogonal polynomials to the Hankel ones is evident from the form of the Gram matrix for sequence

(74)$1,x,x^2,\dots ,x^n.$

This matrix is of Hankel type: for

$P\left(x\right)=P_0+P_{1}x+\cdots +P_n{x}^n,Q\left(x\right)=Q_0+Q_{1}x+\cdots +Q_n{x}^n,$

one has

(75)$⟨P\left(x\right),Q\left(x\right)⟩=(P_0,P_1,\dots ,P_n)\mathbf{M}{(Q_0,Q_1,\dots ,Q_n)}^{\top }$

where

$\mathbf{M}={\left[{\mu }_{j+k-2}\right]}_{j,k=1}^{n+1};{\{{\mu }_j=⟨x^j,1⟩\}}_{j=0}^{2n}.$

The Gram–Schmidt orthogonalization process applied to the sequence (74) can be interpreted as that of successive computation of Hankel polynomials ${\left\{{\mathcal{H}}_k(x;\left\{\mu \right\})\right\}}_{k=1}^n$.

Though orthogonal polynomials are Hankel polynomials, the converse is not necessarily true. An immediate counterexample demonstrates this: for the generating sequence $\{1,-1,1,-1,161,999\}$, the first three Hankel polynomials are as follows

${\mathcal{H}}_1\left(x\right)\equiv x+1,{\mathcal{H}}_2\left(x\right)\equiv 0,{\mathcal{H}}_3\left(x\right)\equiv -25600(x+1),$

and ${\mathcal{H}}_1\left(x\right),{\mathcal{H}}_3\left(x\right)$ cannot be orthogonal. The necessary condition for the orthogonality of the sequence of Hankel polynomials is the positive definiteness of their Hankel matrix. In this context, the result of Theorem 6 cannot be reduced to the case of orthogonal polynomials, it indeed needs the orthogonality-free proof.

The authors of [11] take a step further. Considering the Hankel matrix $\mathbf{M}={\left[{\tau }_{j+k-2}\right]}_{j,k=1}^N$ generated by (65), they replace the requirement of its positive definiteness by that of non-vanishing of all its leading principal minors $H_1\left(\left\{\tau \right\}\right),\dots ,H_N\left(\left\{\tau \right\}\right)$. Formula (75) is kept for the definition of the inner product and it becomes equivalent to

(76)$⟨P\left(x\right),Q\left(x\right)⟩=\sum _{j=1}^N\frac{y_j}{W^{\prime }\left(x_j\right)}P\left(x_j\right)Q\left(x_j\right).$

In this way, some of the traditional properties of the classical inner product are lost. However, the relationship of the type connecting the (quasi)orthogonal Hankel polynomials ${\mathcal{H}}_1(x;\left\{\tau \right\}),{\mathcal{H}}_2(x;\left\{\tau \right\}),\dots$ remains valid with the new definition of the inner product utilized in the formulas like (73) for the involved constants computation. For finding ${\mathcal{H}}_k(x;\left\{\tau \right\})$ (provided that ${\mathcal{H}}_{k-2}(x;\left\{\tau \right\})$ and ${\mathcal{H}}_{k-1}(x;\left\{\tau \right\})$ are known), either the $2k+3$ inner products computations in ${\mathbb{R}}^{k+1}$ via (75) are needed or, equivalently, the same number of polynomial evaluations for the inner product in the form (76).

Returning to Formula (27) for computation of the constants in the JJ-identity (4), one can notice that the cost of this is equivalent to only 2 inner products in ${\mathbb{R}}^{k+1}$. We do not have any idea how to clear up the perplexity why the more general and computationally effective result is completely forgotten; we have not found any reference to it later than 1910.

We next intend to determine what assumption from those posed in Theorem 15 is responsible for the uniqueness of the solution to the rational interpolation problem, i.e., the one that prevents the cases like that in Example 3.

6. Resultant Interpolation

Let us clarify the relationship of the Hankel determinants $H_{m+1}\left(\left\{\tau \right\}\right)$ and $H_n\left(\left\{\tilde{\tau }\right\}\right)$, which appeared in Theorem 15, to the problem of existence of a nontrivial common factor for the numerator and the denominator of the rational interpolant. This is closely related to the notion mentioned in Section 2, namely the resultant $\mathcal{R}\left(p\right(x),q(x\left)\right)$ of the polynomials (6).

Given the polynomials (6) compute their values

(77)$p\left(x_j\right),q\left(x_j\right)\mathrm{for}j\in \{1,\dots ,N=m+n+1\}$

and assume that none of them is zero. Generate the Hankel matrices

$H\left(\tau \right)={\left[{\tau }_{i+j-2}\right]}_{i,j=1}^{m+1},H\left(\tilde{\tau }\right)={\left[{\tilde{\tau }}_{i+j-2}\right]}_{i,j=1}^{n+1}.$

by the sequences

(78)${\left\{{\tau }_k=\sum _{j=1}^N\frac{p\left(x_j\right)}{q\left(x_j\right)}\frac{x_j^k}{W^{\prime }\left(x_j\right)}\right\}}_{k=0}^{2m},{\left\{{\tilde{\tau }}_k=\sum _{j=1}^N\frac{q\left(x_j\right)}{p\left(x_j\right)}\frac{x_j^k}{W^{\prime }\left(x_j\right)}\right\}}_{k=0}^{2n}$

with $W\left(x\right)$ defined by (44).

Theorem 16.

The following equalities are valid

(79) $\begin{array}{ccc}\hfill H_n\left(\left\{\tilde{\tau }\right\}\right)& =& \frac{{(-1)}^{mn+n(n+1)/2}{p}_0}{{\prod }_{j=1}^{N}p\left(x_j\right)}\mathcal{R}(p\left(x\right),q\left(x\right)),\hfill \end{array}$

(80) $\begin{array}{ccc}\hfill H_{n+1}\left(\left\{\tilde{\tau }\right\}\right)& =& \frac{{(-1)}^{mn+n(n+1)/2}{q}_0}{{\prod }_{j=1}^{N}p\left(x_j\right)}\mathcal{R}(p\left(x\right),q\left(x\right)),\hfill \end{array}$

$\begin{array}{ccc}\hfill H_m\left(\left\{\tau \right\}\right)& =& \frac{{(-1)}^{m(m+1)/2}{q}_0}{{\prod }_{j=1}^{N}q\left(x_j\right)}\mathcal{R}(p\left(x\right),q\left(x\right)),\hfill \\ \hfill H_{m+1}\left(\left\{\tau \right\}\right)& =& \frac{{(-1)}^{m(m+1)/2}{p}_0}{{\prod }_{j=1}^{N}q\left(x_j\right)}\mathcal{R}(p\left(x\right),q\left(x\right)).\hfill \end{array}$

Proof. 

Will be illuminated for a particular case $n=3,m=5$. Consider first the case where $p\left(x\right)$ possesses only simple zeros; denote them by ${\lambda }_1,{\lambda }_2,{\lambda }_3$. Construct a new sequence:

${\eta }_k=\sum _{\ell =1}^3\frac{q\left({\lambda }_{\ell }\right){\lambda }_{\ell }^k}{p^{\prime }\left({\lambda }_{\ell }\right)W\left({\lambda }_{\ell }\right)}\mathrm{for}k\in \{0,1,\dots \}.$

From (49) it follows that ${\tilde{\tau }}_k=-{\eta }_k$ for $k\in \{0,\dots ,5\}$ and therefore

$H_3\left(\left\{\tilde{\tau }\right\}\right)=\left|\begin{array}{ccc}{\tilde{\tau }}_0& {\tilde{\tau }}_1& {\tilde{\tau }}_2\\ {\tilde{\tau }}_1& {\tilde{\tau }}_2& {\tilde{\tau }}_3\\ {\tilde{\tau }}_2& {\tilde{\tau }}_3& {\tilde{\tau }}_4\end{array}\right|=-\left|\begin{array}{ccc}{\eta }_0& {\eta }_1& {\eta }_2\\ {\eta }_1& {\eta }_2& {\eta }_3\\ {\eta }_2& {\eta }_3& {\eta }_4\end{array}\right|$

$=-\left|\begin{array}{ccc}1& 1& 1\\ {\lambda }_1& {\lambda }_2& {\lambda }_3\\ {\lambda }_1^2& {\lambda }_2^2& {\lambda }_3^2\end{array}\right|·\left|\begin{array}{ccc}\frac{q\left({\lambda }_1\right)}{p^{\prime }\left({\lambda }_1\right)W\left({\lambda }_1\right)}& 0& 0\\ 0& \frac{q\left({\lambda }_2\right)}{p^{\prime }\left({\lambda }_2\right)W\left({\lambda }_2\right)}& 0\\ 0& 0& \frac{q\left({\lambda }_3\right)}{p^{\prime }\left({\lambda }_3\right)W\left({\lambda }_3\right)}\end{array}\right|·\left|\begin{array}{ccc}1& {\lambda }_1& {\lambda }_1^2\\ 1& {\lambda }_2& {\lambda }_2^2\\ 1& {\lambda }_3& {\lambda }_3^2\end{array}\right|$

$=-\prod _{1\le j<k\le 3}{({\lambda }_j-{\lambda }_k)}^2\frac{{\displaystyle \prod _{\ell =1}^{3}q\left({\lambda }_{\ell }\right)}}{{\displaystyle \prod _{\ell =1}^{3}W\left({\lambda }_{\ell }\right)\prod _{\ell =1}^3{p}^{\prime }\left({\lambda }_{\ell }\right)}}.$

Using the definition of the resultant in the form (7), the last expression can be represented as

$-\prod _{1\le j<k\le 3}{({\lambda }_j-{\lambda }_k)}^2\frac{{\displaystyle \prod _{\ell =1}^{3}q\left({\lambda }_{\ell }\right)}}{{\displaystyle \frac{\mathcal{R}\left(p\right(x),W(x\left)\right)}{p_0^9}{p}_0^3{(-1)}^3\prod _{1\le j<k\le 3}{({\lambda }_j-{\lambda }_k)}^2}}$

Finally use the alternative definition of the resultant (8):

(81)$\frac{{\displaystyle p_0^6\prod _{\ell =1}^{3}q\left({\lambda }_{\ell }\right)}}{{\displaystyle \mathcal{R}\left(p\right(x),W(x\left)\right)}}=\frac{p_0\mathcal{R}(p\left(x\right),q\left(x\right))}{{\displaystyle {(-1)}^{9\times 3}\prod _{j=1}^{9}p\left(x_j\right)}}.$

Thus, the equality (79) is true.

We just proved this equality under the additional assumption that polynomial $p\left(x\right)$ has all its zeros distinct. To extend this result to the general case, the traditional trick consists in application of Theorem 4. By Corollary 1, the condition of distinction (simplicity) of zeros of polynomial with symbolic (indeterminate) coefficients can be expressed as an algebraic inequality with respect to these coefficients. In accordance with the referred theorem, the algebraic identity which is valid under an extra assumption in the form of algebraic inequality, should be valid for all values of indeterminates.

To prove the equality (80), multiply the determinant

$H_4\left(\left\{\tilde{\tau }\right\}\right)=\left|\begin{array}{cccc}{\tilde{\tau }}_0& {\tilde{\tau }}_1& {\tilde{\tau }}_2& {\tilde{\tau }}_3\\ {\tilde{\tau }}_1& {\tilde{\tau }}_2& {\tilde{\tau }}_3& {\tilde{\tau }}_4\\ {\tilde{\tau }}_2& {\tilde{\tau }}_3& {\tilde{\tau }}_4& {\tilde{\tau }}_5\\ {\tilde{\tau }}_3& {\tilde{\tau }}_4& {\tilde{\tau }}_5& {\tilde{\tau }}_6\end{array}\right|$

from the right-hand side by the determinant

$\left|\begin{array}{cccc}1& 0& 0& p_3/p_0\\ 0& 1& 0& p_2/p_0\\ 0& 0& 1& p_1/p_0\\ 0& 0& 0& 1\end{array}\right|,$

evidently equal to 1. This results in the determinant which differs from the initial one only in the last column:

$\frac{1}{p_0}\left|\begin{array}{cccc}{\tilde{\tau }}_0& {\tilde{\tau }}_1& {\tilde{\tau }}_2& p_0{\tilde{\tau }}_3+p_1{\tilde{\tau }}_2+p_2{\tilde{\tau }}_1+p_3{\tilde{\tau }}_0\\ {\tilde{\tau }}_1& {\tilde{\tau }}_2& {\tilde{\tau }}_3& p_0{\tilde{\tau }}_4+p_1{\tilde{\tau }}_3+p_2{\tilde{\tau }}_2+p_3{\tilde{\tau }}_1\\ {\tilde{\tau }}_2& {\tilde{\tau }}_3& {\tilde{\tau }}_4& p_0{\tilde{\tau }}_5+p_1{\tilde{\tau }}_4+p_2{\tilde{\tau }}_3+p_3{\tilde{\tau }}_2\\ {\tilde{\tau }}_3& {\tilde{\tau }}_4& {\tilde{\tau }}_5& p_0{\tilde{\tau }}_6+p_1{\tilde{\tau }}_5+p_2{\tilde{\tau }}_4+p_3{\tilde{\tau }}_3\end{array}\right|.$

Compute the entries of this column:

$\begin{array}{cc}\hfill p_0{\tilde{\tau }}_3+p_1{\tilde{\tau }}_2+p_2{\tilde{\tau }}_1+p_3{\tau }_0& =\sum _{j=1}^9\frac{q\left(x_j\right)}{p\left(x_j\right)}\frac{(p_0{x}_j^3+p_1{x}_j^2+p_2{x}_j+p_3)}{W^{\prime }\left(x_j\right)}\hfill \\ \hfill & =\sum _{j=1}^9\frac{q\left(x_j\right)}{p\left(x_j\right)}\frac{p\left(x_j\right)}{W^{\prime }\left(x_j\right)}=\sum _{j=1}^9\frac{q\left(x_j\right)}{W^{\prime }\left(x_j\right)}\stackrel{\left(47\right)}{=}0,\hfill \end{array}$

and similarly

$\begin{array}{cc}\hfill p_0{\tilde{\tau }}_4+p_1{\tilde{\tau }}_3+p_2{\tilde{\tau }}_2+p_3{\tilde{\tau }}_1& =0,\\ \hfill p_0{\tilde{\tau }}_5+p_1{\tilde{\tau }}_4+p_2{\tilde{\tau }}_3+p_3{\tilde{\tau }}_2& =0,\end{array}$

whereas

$p_0{\tilde{\tau }}_6+p_1{\tilde{\tau }}_5+p_2{\tilde{\tau }}_4+p_3{\tilde{\tau }}_3=\sum _{j=1}^9\frac{q\left(x_j\right)x_j^3}{W^{\prime }\left(x_j\right)}=q_0.$

Therefore,

$H_4\left(\left\{\tilde{\tau }\right\}\right)=\frac{1}{p_0}\left|\begin{array}{cccc}{\tilde{\tau }}_0& {\tilde{\tau }}_1& {\tilde{\tau }}_2& 0\\ {\tilde{\tau }}_1& {\tilde{\tau }}_2& {\tilde{\tau }}_3& 0\\ {\tilde{\tau }}_2& {\tilde{\tau }}_3& {\tilde{\tau }}_4& 0\\ {\tilde{\tau }}_3& {\tilde{\tau }}_4& {\tilde{\tau }}_5& q_0\end{array}\right|=\frac{q_0}{p_0}{H}_3\left(\left\{\tilde{\tau }\right\}\right)$

and with already deduced expression (81) for $H_3\left(\left\{\tilde{\tau }\right\}\right)$, one gets the validity of (80).  □

Remark 3.

The number of the interpolation values (77) exceeds twice that of the coefficients of both polynomials $p\left(x\right)$ and $q\left(x\right)$, i.e., the set of interpolation values is redundant for the resultant evaluation. However, from the statement of Theorem 16, one can notice that only the set ${\{p\left(x_j\right)/q\left(x_j\right)\}}_{j=1}^{n+m+1}$ of ratios (or their reciprocals) is involved in the resultant’s computation. This set is not redundant. (As for some extra multiples in the right-hand sides of Formulas (79) and (80), one may ignore the case of their vanishing as a practically improbable if the polynomials are randomly chosen over infinite fields).

7. Hankel Polynomials Computation for Extended Data Set

Let the data set (1) be extended by adding an extra point:

(82)$\begin{array}{ccccccc}x& & x_1& x_2& \dots & x_N& x_{N+1}\\ y& & y_1& y_2& \dots & y_N& y_{N+1}\end{array}\mathrm{with}{x}_{N+1}\notin \{x_1,\dots ,x_N\}.$

Is it possible to organize the computation of the rational interpolant for this set based on already computed interpolant for the initial set (1)? In other words: we are looking for an analogue of the Newton interpolation polynomial construction procedure.

In the framework of the approach discussed in Section 5, the stated problem can be reformulated as that of effective computation of Hankel polynomials generated by the sequence

(83)${\widehat{\tau }}_k=\sum _{j=1}^{N+1}\frac{x_j^k{y}_j}{{\widehat{W}}^{\prime }\left(x_j\right)};\mathrm{here}\widehat{W}\left(x\right)=\prod _{j=1}^{N+1}(x-x_j)\stackrel{\left(44\right)}{\equiv }W\left(x\right)(x-x_{N+1}).$

and a similar one obtained on replacing $y_j$ by $1/y_j$, under the assumption that the corresponding Hankel polynomials for smaller data set (1) are already at our disposal.

Theorem 17.

Hankel polynomial (we slightly simplify here the notation: $H_k=H_k\left(\left\{\tau \right\}\right)$, ${\widehat{H}}_k=H_k\left(\left\{{\widehat{\tau }}_k\right\}\right)$)

(84) ${\widehat{\mathcal{H}}}_k(x;\left\{\widehat{\tau }\right\})=\left|\begin{array}{cccc}{\widehat{\tau }}_0\hfill & {\widehat{\tau }}_1\hfill & \dots \hfill & {\widehat{\tau }}_k\hfill \\ {\widehat{\tau }}_1\hfill & {\widehat{\tau }}_2\hfill & \dots \hfill & {\widehat{\tau }}_{k+1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ {\widehat{\tau }}_{k-1}\hfill & {\widehat{\tau }}_k\hfill & \dots \hfill & {\widehat{\tau }}_{2k-1}\hfill \\ 1\hfill & x\hfill & \dots \hfill & x^k\hfill \end{array}\right|\equiv {\widehat{H}}_k{x}^k+\dots .$

generated by the sequence (83) is linked to polynomials ${\mathcal{H}}_k(x;\left\{\tau \right\})$ and ${\mathcal{H}}_{k-1}(x;\left\{\tau \right\})$ generated by the sequence (65) via the following identity

(85) $H_k{\widehat{\mathcal{H}}}_k(x;\left\{\widehat{\tau }\right\})\equiv {\widehat{H}}_k{\mathcal{H}}_k(x;\left\{\tau \right\})+{\widehat{H}}_{k+1}{\mathcal{H}}_{k-1}(x;\left\{\tau \right\}).$

Proof. 

Is based on the following equality connecting the sums (65) and (83):

(86)${\tau }_k={\widehat{\tau }}_{k+1}-x_{N+1}{\widehat{\tau }}_k$

which can be deduced in exactly the same manner as the equality (57) in the proof of Theorem 12. With its aid, represent ${\mathcal{H}}_k(x;\left\{\tau \right\})$ first as the determinant of the order $k+1$:

${\mathcal{H}}_k(x;\left\{\tau \right\})=\left|\begin{array}{cccc}{\widehat{\tau }}_1-x_{N+1}{\widehat{\tau }}_0\hfill & {\widehat{\tau }}_2-x_{N+1}{\widehat{\tau }}_1\hfill & \dots \hfill & {\widehat{\tau }}_{k+1}-x_{N+1}{\widehat{\tau }}_k\hfill \\ {\widehat{\tau }}_2-x_{N+1}{\widehat{\tau }}_1\hfill & {\widehat{\tau }}_3-x_{N+1}{\widehat{\tau }}_2\hfill & \dots \hfill & {\widehat{\tau }}_{k+2}-x_{N+1}{\widehat{\tau }}_{k+1}\hfill \\ \cdots \hfill & \cdots \hfill & & \cdots \hfill \\ {\widehat{\tau }}_k-x_{N+1}{\widehat{\tau }}_{k-1}\hfill & {\widehat{\tau }}_{k+1}-x_{N+1}{\widehat{\tau }}_k\hfill & \dots \hfill & {\widehat{\tau }}_{2k}-x_{N+1}{\widehat{\tau }}_{2k-1}\hfill \\ 1\hfill & x\hfill & \dots \hfill & x^k\hfill \end{array}\right|$

and then as the determinant of the order $k+2$:

(87)$={(-1)}^{k+1}{\left|\begin{array}{ccccc}{\widehat{\tau }}_0\hfill & {\widehat{\tau }}_1\hfill & \dots \hfill & {\widehat{\tau }}_k\hfill & \mid 1\hfill \\ {\widehat{\tau }}_1\hfill & {\widehat{\tau }}_2\hfill & \dots \hfill & {\widehat{\tau }}_{k+1}\hfill & \mid x_{N+1}\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill & ⋮\hfill \\ {\widehat{\tau }}_k\hfill & {\widehat{\tau }}_{k+1}\hfill & \dots \hfill & {\widehat{\tau }}_{2k}\hfill & \mid x_{N+1}^k\hfill \\ 1\hfill & x\hfill & \dots \hfill & x^k\hfill & \mid 0\hfill \end{array}\right|}_{(k+2)\times (k+2)}.$

Application of Sylvester’s identity (5) to the latter results in (85).  □

Corollary 1.

The following equality is valid:

(88) ${\widehat{H}}_{k+1}{\mathcal{H}}_{k-1}\left(x_{N+1}\right)={(-1)}^k{H}_k^2-{\mathcal{H}}_k\left(x_{N+1}\right){\widehat{H}}_k.$

Proof. 

Substitution $x=x_{N+1}$ into (85) results in

${\widehat{\mathcal{H}}}_k\left(x_{N+1}\right)H_k\equiv {\widehat{H}}_k{\mathcal{H}}_k\left(x_{N+1}\right)+{\widehat{H}}_{k+1}{\mathcal{H}}_{k-1}\left(x_{N+1}\right)$

Comparing the coefficients of $x^k$ in ${\mathcal{H}}_k\left(x\right)$ and in (87) yields that

$H_k={(-1)}^k{\widehat{\mathcal{H}}}_k\left(x_{N+1}\right).$

which completes the proof of (88).  □

Formula (88) gives rise to the recursive procedure for Hankel determinants ${\left\{{\widehat{H}}_{\ell }\right\}}_{\ell =1}^{k+1}$ computation, with the final result as follows:

Theorem 18.

Let ${\{{\mathcal{H}}_j\left(x_{N+1}\right)\ne 0\}}_{j=1}^k$. One has

(89) ${\widehat{H}}_{k+1}={(-1)}^k{\mathcal{H}}_k\left(x_{N+1}\right)\left(\sum _{j=1}^k\frac{H_j^2}{{\mathcal{H}}_j\left(x_{N+1}\right){\mathcal{H}}_{j-1}\left(x_{N+1}\right)}+{\widehat{\tau }}_0\right),with{\mathcal{H}}_0\left(x\right)\equiv 1.$

Proof. 

Is carried out by induction. For $k=1$, the formula evidently follows from (88). Assume that (89) holds for some k. Let us prove its validity for $k+1$. From (88) one arrives at

${\widehat{H}}_{k+2}=\frac{1}{{\mathcal{H}}_k\left(x_{N+1}\right)}\left({(-1)}^{k+1}{H}_{k+1}^2-{\mathcal{H}}_{k+1}\left(x_{N+1}\right){\widehat{H}}_{k+1}\right)$

and assuming that ${\mathcal{H}}_{k+1}\left(x_{N+1}\right)\ne 0$,

$={(-1)}^{k+1}{\mathcal{H}}_{k+1}\left(x_{N+1}\right)\left(\frac{H_{k+1}^2}{{\mathcal{H}}_{k+1}\left(x_{N+1}\right){\mathcal{H}}_k\left(x_{N+1}\right)}+\frac{{(-1)}^k}{{\mathcal{H}}_k\left(x_{N+1}\right)}{\tilde{H}}_{k+1}\right).$

Application of induction hypothesis (89) completes the proof.  □

Combination of the preceding results finalized in computational formula for

(90)$$\begin{gathered} {\widehat{\mathcal{H}}}_k\left(x\right) \\ =\frac{1}{{\mathcal{H}}_{k-1}\left(x_{N+1}\right)}\left({(-1)}^k{H}_k{\mathcal{H}}_{k-1}\left(x\right)-\frac{{\widehat{H}}_k}{H_k}\left\{{\mathcal{H}}_k\left(x_{N+1}\right){\mathcal{H}}_{k-1}\left(x\right)-{\mathcal{H}}_{k-1}\left(x_{N+1}\right){\mathcal{H}}_k\left(x\right)\right\}\right) \end{gathered}$$

under the natural assumption for the non-vanishing of the denominators. It turns out that computation of Hankel polynomials for the extended data set (82) can be performed directly via their counterparts for the initial data set (1) without evaluation of all the symmetric functions (83) (except for ${\widehat{\tau }}_0$).

Example 5.

Add the point $x=4,y=7$ to the data set of Example 4 and find the rational interpolant ${\widehat{r}}_{3,4}\left(x\right)=p\left(x\right)/q\left(x\right)$ with $degp\left(x\right)=3,degq\left(x\right)=4$ satisfying the extended set.

Solution. To find the denominator by Formula (90):

${\widehat{\mathcal{H}}}_4\left(x\right)=\frac{1}{{\mathcal{H}}_3\left(4\right)}\left[H_4{\mathcal{H}}_3\left(x\right)-\frac{{\widehat{H}}_4}{H_4}\left\{{\mathcal{H}}_4\left(4\right){\mathcal{H}}_3\left(x\right)-{\mathcal{H}}_3\left(4\right){\mathcal{H}}_4\left(x\right)\right\}\right]$

we use the already found in solution of Example 4 expressions for Hankel polynomials ${\mathcal{H}}_3\left(x\right),{\mathcal{H}}_4\left(x\right)$. The values of the leading coefficients of these polynomials, as well as those of ${\mathcal{H}}_1\left(x\right),{\mathcal{H}}_2\left(x\right)$, complemented with the additional computation of ${\widehat{\tau }}_0=-1/112$ permits one to find

${\widehat{H}}_4\stackrel{\left(89\right)}{=}-{\mathcal{H}}_3\left(4\right)\left(\frac{H_1^2}{{\mathcal{H}}_1\left(4\right)}+\frac{H_2^2}{{\mathcal{H}}_1\left(4\right){\mathcal{H}}_2\left(4\right)}+\frac{H_3^2}{{\mathcal{H}}_2\left(4\right){\mathcal{H}}_3\left(4\right)}+{\widehat{\tau }}_0\right)=\frac{17117}{27216000}.$

Thus,

${\widehat{\mathcal{H}}}_4\left(x\right)=\frac{17117}{27216000}{x}^4-\frac{249211}{54432000}{x}^3-\frac{52343}{27216000}{x}^2+\frac{1165699}{54432000}x+\frac{53}{20160}.$

Computation of the numerator of the wanted fraction is performed in a similar way using the already computed in Example 4 expressions for ${\tilde{\mathcal{H}}}_2\left(x\right)$ and ${\tilde{\mathcal{H}}}_3\left(x\right)$. Finally,

${\widehat{r}}_{3,4}\left(x\right)=\frac{-924015x^3+642060x^2+5084535x+143100}{34234x^4-249211x^3-104686x^2+1165699x+143100}.$

We do not discuss here the exceptional cases where some of polynomials ${\mathcal{H}}_k\left(x\right)$ vanish at $x_{N+1}$. As for a relationship of the obtained result to the Newton polynomial interpolant construction procedure (since the former should contain the latter as a particular case) we restrict ourselves to the reference to Theorem 13.

8. Comparison with the Barycentric Representation

We intend here to establish the relationship of Jacobi’s approach to the one based on the so-called barycentric representation.

Theorem 19.

(Ref. [6]). If a solution to the rational interpolation problem (61)–(62) with $m\le n$ exists, then $\mathbf{u}=[u_1,\dots ,u_N]$ is a vector of weights in one of its barycentric representations

$r\left(x\right)=\left(\sum _{j=1}^N\frac{u_j}{x-x_j}{y}_j\right)/\sum _{j=1}^N\frac{u_j}{x-x_j}$

if and only if $\mathbf{u}$ belongs to the kernel of the $(N-1)\times N$ stack matrix

$\mathbf{A}=\left(\begin{array}{c}\mathbf{V}\\ \mathbf{H}\end{array}\right)$

where

$\mathbf{V}=\left(\begin{array}{cccc}1& 1& \dots & 1\\ x_1& x_2& \dots & x_N\\ x_1^2& x_2^2& \dots & x_N^2\\ ⋮& ⋮& & ⋮\\ x_1^{n-1}& x_2^{n-1}& \dots & x_N^{n-1}\end{array}\right),\mathbf{H}=\left(\begin{array}{cccc}{y}_1& y_2& \dots & y_N\\ x_1{y}_1& x_2{y}_2& \dots & x_N{y}_N\\ ⋮& ⋮& & ⋮\\ x_1^{m-1}{y}_1& x_2^{m-1}{y}_2& \dots & x_N^{m-1}{y}_N\end{array}\right).$

First, we claim that $m+1$ vectors

${\Xi }_j=\left[x_1^j/W^{\prime }\left(x_1\right),\dots ,x_N^j/W^{\prime }\left(x_N\right)\right]\mathrm{for}j\in \{0,1,\dots ,m\}$

belong to the kernel of the Vandermonde matrix $\mathbf{V}$. Indeed this evidently follows from the equalities (48). Next we look for a linear combination $v_0{\Xi }_0+v_1{\Xi }_1+\cdots +v_m{\Xi }_m$, which belongs to the kernel of the matrix $\mathbf{H}$:

$\left(\begin{array}{cccc}{y}_1& y_2& \dots & y_N\\ x_1{y}_1& x_2{y}_2& \dots & x_N{y}_N\\ ⋮& ⋮& & ⋮\\ x_1^{m-1}{y}_1& x_2^{m-1}{y}_2& \dots & x_N^{m-1}{y}_N\end{array}\right)$

$\times \left(\begin{array}{cccc}1/W^{\prime }\left(x_1\right)& x_1/W^{\prime }\left(x_1\right)& \dots & x_1^m/W^{\prime }\left(x_1\right)\\ 1/W^{\prime }\left(x_2\right)& x_2/W^{\prime }\left(x_2\right)& \dots & x_2^m/W^{\prime }\left(x_2\right)\\ ⋮& ⋮& & ⋮\\ 1/W^{\prime }\left(x_N\right)& x_N/W^{\prime }\left(x_N\right)& \dots & x_N^m/W^{\prime }\left(x_N\right)\end{array}\right)\left(\begin{array}{c}{v}_0\\ v_1\\ ⋮\\ v_m\end{array}\right)={\mathbb{O}}_{m\times 1},$

or, equivalently, in the notation of symmetric functions (65) of the data set,

${\left(\begin{array}{cccc}{\tau }_0\hfill & {\tau }_1\hfill & \dots \hfill & {\tau }_m\hfill \\ {\tau }_1\hfill & {\tau }_2\hfill & \dots \hfill & {\tau }_{m+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ {\tau }_{m-1}\hfill & {\tau }_m\hfill & \dots \hfill & {\tau }_{2m-1}\hfill \end{array}\right)}_{m\times (m+1)}\left(\begin{array}{c}{v}_0\\ v_1\\ ⋮\\ v_m\end{array}\right)={\mathbb{O}}_{m\times 1}.$

Provided that the rank of the matrix standing in the left-hand side of this equation equals m, the dimension of its kernel equals 1 with the basis vector consisting of the sequence of algebraic complements to the entries of the last row of the matrix

${\left(\begin{array}{cccc}{\tau }_0\hfill & {\tau }_1\hfill & \dots \hfill & {\tau }_m\hfill \\ {\tau }_1\hfill & {\tau }_2\hfill & \dots \hfill & {\tau }_{m+1}\hfill \\ ⋮\hfill & ⋮\hfill & & ⋮\hfill \\ {\tau }_{m-1}\hfill & {\tau }_m\hfill & \dots \hfill & {\tau }_{2m-1}\hfill \\ *\hfill & *\hfill & \dots \hfill & *\hfill \end{array}\right)}_{(m+1)\times (m+1)}.$

Therefore, this vector must be proportional (collinear) to the vector of the coefficients of the denominator of the rational interpolant represented by (70).

Consequently, the Hankel polynomial approach to the rational interpolation problem can be interpreted as a method for resolution of the system of linear equations for determining the vector of weights in barycentric representation of the interpolant from Theorem 19.

9. Redundant Data Set with Systematic Errors

We now turn to a problem of data storage systems, namely the one aimed at establishing the locations of erroneous blocks in their sequence containing together with data (information) blocks and also some specially designed checksums. This problem relates to the theory of the Error Correcting Codes. To explain the idea of one of the practical procedures from this theory, we first present an artificial example.

Example 6.

Consider the data set

$\begin{array}{ccccccccc}x& & -2& -1& 0& 1& 2& 3& 4\\ y& & 30& -7& 8& 9& 11& 35& 60\end{array}$

It is known that initially it has been generated by the polynomial of degree two with integer coefficients, but later on, a suspicion arises that up to two y-values are corrupted. Find the erroneous values and restore the initial polynomial.

Solution. Compute successively Hankel polynomials from Theorem 15 and try to factorize them over $\mathbb{Q}$. We succeed at

${\mathcal{H}}_4(x;\left\{\tilde{\tau }\right\})\equiv -\frac{1}{1451520000}(x+1)(x-2)(4x^2-3x+8)$

and at

${\mathcal{H}}_2(x;\left\{\tau \right\})\equiv -\frac{77}{320}(x+1)(x-2).$

It can be easily verified that the values of polynomial $f\left(x\right)=4x^2-3x+8$ coincide with the y-values from the data set at all the given x-values except for $x=-1$ and $x=2$. Thus, the integer zeros of ${\mathcal{H}}_2(x;\left\{\tau \right\})$ give the location of the erroneous values while the square factor of ${\mathcal{H}}_4(x;\left\{\tilde{\tau }\right\})$ coincides with the “true” polynomial generating the remaining correct y-values.  □

We now concretize the problem statement. Let the data set (1) of rational numbers be redundant for the polynomial computation, i.e., the number of true values in it exceeds $k+1$ where $k=degf\left(x\right)$ is known a priori. The number of corrupted y-values is unknown, but one may expect that it is much less than those of true. We need to find x-values corresponding to the erroneous y-values. The first step in dealing with this problem should be a test for the error presence. The next result is a trivial consequence of Theorem 11.

Theorem 20.

For the polynomial interpolant generated by (1) to be of the degree $k<N-1$, it is necessary and sufficient that

(91) ${\tau }_0=0,\dots ,{\tau }_{N-k-2}=0,{\tau }_{N-k-1}\ne 0.$

Thus, the uncorrupted data set (1) passes the test ${\{{\tau }_j=0\}}_{j=0}^{N-k-2}$, and then the Formula (53) permits one to compute the polynomial. One may expect that the corrupted data set would not pass the test: generically, the degree of the interpolant for the random data set equals $N-1>k$. If the number E of erroneous values is given to one by revelation, the following result can be utilized for their location.

Theorem 21.

(Ref. [32]). Let $E\in \{1,2,\dots ,\lfloor N/2\rfloor -1\}$ and $e_1,\dots ,e_E$ be distinct numbers from $\{1,2,\dots ,N\}$. Let polynomial $f\left(x\right)$ be of a degree $k<N-2E$. Let the data set satisfy the conditions

• (a) 

  $y_j=f\left(x_j\right)$ for $j\in \{1,\dots ,N\}\backslash \{e_1,\dots ,e_E\}$,

• (b) 

  ${\widehat{y}}_{e_s}:=f\left(x_{e_s}\right)\ne y_{e_s}$ for $s\in \{1,\dots ,E\}$.

Then,

(92) ${\mathcal{H}}_E(x;\left\{\tau \right\})\equiv \frac{{\displaystyle \prod _{s=1}^E(y_{e_s}-{\widehat{y}}_{e_s})\prod _{1\le s<t\le E}{(x_{e_t}-x_{e_s})}^2}}{{\displaystyle \prod _{s=1}^E{W}^{\prime }\left(x_{e_s}\right)}}\prod _{s=1}^E(x-x_{e_s}).$

This result confirms one empirical conclusion from solution to Example 6. The other one can also be justified, i.e., under the conditions of Theorem 21, the following identity is valid, up to a sign [33]:

${\mathcal{H}}_{k+E}(x;\left\{\tilde{\tau }\right\})\equiv \frac{a_0^{N-k-2E-1}}{{\prod }_{j=1}^N{y}_j}{\mathcal{H}}_E(x;\left\{\tau \right\})f\left(x\right),$

here, $a_0$ stands for the leading coefficient of $f\left(x\right)$.

Unfortunately, in the practice of the error correcting codes, the exact number of erroneous values, as well as their location, is seldom known a priori. Therefore, in order to find the polynomial from Theorem 21, one should verify the reducibility of the successive polynomials ${\mathcal{H}}_1(x;\left\{\tau \right\}),{\mathcal{H}}_2(x;\left\{\tau \right\}),\dots$. This reason immediately refers to the results of Section 3 for the effective procedure of computation of the next polynomial in the succession when the previous ones are already known.

We now recall the idea of the Berlekamp-Welch algorithm [22] for the error correction codes. Assume the integers $y_1,\dots ,y_k$ are treated as information blocks to be protected against failures when put into computer data storage. For this aim, let us organize the checksums computation in the following way. First, chose arbitrary distinct integers $x_1,\dots ,x_k$, say ${\{x_j=j\}}_{j=1}^k$. Next, compute the polynomial interpolant $f\left(x\right)$ of the degree $\le k-1$ for the data set ${\left\{(x_j,y_j)\right\}}_{j=1}^k$. Finally, evaluate $f\left(x\right)$ at the points ${\{x_j=j\}}_{j=k+1}^N$. The obtained data set ${\left\{(x_j,f\left(x_j\right))\right\}}_{j=1}^N$ is redundant for the $f\left(x\right)$ computation, but now, in the view of preceding considerations, it might sustain some amount of data corruptions.

The real life application of the algorithm is performed in finite fields, say $\mathbf{GF}\left(2^D\right)$ for sufficiently large positive integer D. Usually, ${\{x_j={\mathfrak{a}}^{j-1}\}}_{j=1}^N$ where $\mathfrak{a}$ stands for a primitive element of the field. The polynomial (92) is addressed as the error locator polynomial. In [22], the algorithm of its construction is treated as the consequence of a linear system (64) solving. The presence of errors in the data set causes the appearance of a nontrivial common factor for the numerator and the denominator, i.e., the conditions of Theorem 15 are not fulfilled. Theorem 21 demonstrates that it is possible to make a version of the Berlekamp–Welch algorithm which is free from the rational interpolation stage.

Compared with infinite fields, several computational difficulties appear when implementing the algorithms in $\mathbf{GF}\left(2^D\right)$. We just list some of them: finding the inversion and discrete logarithm of an element, evaluation and zero finding for a polynomial over $\mathbf{GF}\left(2^D\right)$, etc. The vanishing of $H_{k-1}$ in the JJ-identity (15) is to be treated as an exception in an infinite field, but it has a non-zero probability in $\mathbf{GF}\left(2^D\right)$. A complete treatment of these issues is to be done in further publications on the subject.

10. Conclusions

We have addressed an approach for the solution of the polynomial and rational interpolation problem originated in the paper by Jacobi. It consists first in representation of the interpolant by virtue of appropriate Hankel polynomials and further recursive computation of the sequence of these polynomials via the Jacobi–Joachimsthal identity and its development. This yields one an opportunity to compute effectively the whole family of rational interpolants for the given data set.

Problems for further investigation relating the theoretical backgrounds of the Jacobi’s approach are mentioned in Section 3. Extension of the approach to the (mentioned in Introduction) rational Hermite interpolation problem [34] looks promising.

We trace some other topics of potential interest related to the approach:

• The problem of choice of the interpolant with the prescribed restrictions on its entries. For the case of polynomials over $\mathbb{R}$, this might be the demand of absence of either real zeros for the denominator or zeros lying in the right half-plane of the complex plane (in other words, the denominator must be a stable polynomial). Is it possible to resolve these types of problems in the framework of the Hankel polynomial sequence construction? The positive answer can be guessed from the title of the original Joachimsthal’s paper [13]. Indeed, it is devoted to the zero separation problem for a univariate polynomial over $\mathbb{R}$.

• Multivariate interpolation. The treatment of this problem meets more obstacles than the univariate counterpart. For instance, one can select such a data sets ${\left\{(x_j,y_j,z_j)\right\}}_{j=1}^9\subset {\mathbb{R}}^3$ with distinct pairs $(x_j,y_j)$ for which the polynomial interpolant

  $f(x,y)\in \mathbb{R}[x,y],degf=3,{\{(x_j,y_j)=z_j\}}_{j=1}^9$

  does not exist. On the other hand, there exist counterparts for the results from Section 4 for multivariate symmetric functions. For instance, an analogue of Theorem 10 is known as (what a surprise!) the Jacobi equality [35]. Therefore, the extendability of the everywhere Jacobi’s approach to the multivariate rational interpolation problem looks intriguing.

Author Contributions

Conceptualization and methodology, A.U.; investigation, I.B. and E.K.; literature search, A.U. and E.K.; algorithm implementation, I.B.; writing—original draft preparation, A.U. and I.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The authors wish to thank the anonymous reviewers for their very constructive comments.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.