k

cides, after this particular choice of gi; j , with

1 @Fj

g= det : i; j = 1; : : : ; d ;

deg f · · · deg f @i

1 d

which is independent of k now; in other words, there exist polynomials gk ; k = 1; : : : ; N , such that

ˆ

deg gk ¡ deg g and det Gk = g + gk . Moreover, g is a homogeneous polynomial of degree at most

ˆ ˆ

deg f + · · · + deg f ’ d. Now, let p be any polynomial, then

1 d

N N N

det Gj (·) p(zj ) p(zj )

Kp = p(zj ) =g + g:

ˆ (8)

det Gj (zj ) j

det Gj (zj ) det Gj (zj )

j=1 j=1 j=1

Combining (8) with (6) then yields the existence of polynomials q1 ; : : : ; qd such that

N N d

p(zj ) p(zj )

p=g + g+

ˆ qj fj

det Gj (zj ) j

det Gj (zj )

j=1 j=1 j=1

and comparing homogeneous terms of degree deg g Kronecker realized that either, for any p such

that deg p ¡ deg g,

N

p(zj )

=0 (9)

det Gj (zj )

j=1

or there exist homogeneous polynomials h1 ; : : : ; hd such that

d

g= hj det F : (10)

j

j=1

The latter case, Eq. (10), says (in algebraic terminology) that there is a syzygy among the leading

terms of the polynomials F ; j = 1; : : : ; d, and is equivalent to the fact that N ¡ deg f · · · deg f ,

j 1 d

while (9) describes and even characterizes the complete intersection case that N = deg f · · · deg f .

1 d

In his paper, Kronecker also mentions that the condition (10) has been overlooked in [21]. Jacobi

dealt there with the common zeros of two bivariate polynomials and derived explicit representations

for the functional

N

f(zj )

[z1 ; : : : ; zN ]f:= ; (11)

det Gj (zj )

j=1

M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35 27

which behaves very much like a divided di erence, since it is a combination of point evaluations

d

which, provided that (9) hold true, annihilates deg g’1 .

In addition, Kronecker refers to a paper [6] which he says treats the case of symmetric functions,

probably elementary symmetric polynomials. Unfortunately, this paper is unavailable to us so far.

3. Bivariate tables, the natural approach

Only very few research papers on multivariate polynomial interpolation were published during the

ÿrst part of this century. In the classical book Interpolation [45], where one section (Section 19) is

devoted to this topic, the author only refers to two related papers, recent at that time (1927), namely

[27,28]. The latter one [28], turned out to be inaccessible to us, unfortunately, but it is not di cult

to guess that it might have pursued a tensor product approach, because this is the unique point of

view of [45] (see also [31]).

The formulas given in [27] are Newton formulas for tensor product interpolation in two variables,

and the author, Narumi, claims (correctly) that they can be extended to “many variables”. Since it is

a tensor product approach, the interpolation points are of the form (xi ; yj ), 06i6m; 06j6n, with

xi ; yj arbitrarily distributed on the axes OX and OY , respectively. Bivariate divided di erences for

these sets of points are obtained in [27], by recurrence, separately for each variable. With the usual

notations, the interpolation formula from [27] reads as

j’1

m n i’1

p(x; y) = f[x0 ; : : : ; xi ; y0 ; : : : ; yj ] (x ’ xh ) (y ’ xk ); (12)

i=0 j=0 h=0 k=0

where empty products have the value 1. Remainder formulas based on the mean value theorem are

also derived recursively from the corresponding univariate error formulas in [27]. For f su ciently

smooth there exist values ; ; Á; Á such that

@m+1 f( ; y) m (x ’ xh ) @n+1 f(x; Á) n (y ’ yk )

h=0 k=0

R(x; y) = +

@x (m + 1)! @y (n + 1)!

m+1 n+1

@m+n+2 f( ; Á ) m (x ’ xh ) n (y ’ yk )

h=0 k=0

’ : (13)

@x m+1 @y n+1 (m + 1)! (n + 1)!

The special case of equidistant points on both axes is particularly considered in [27], and since

the most popular formulas at that time were based on ÿnite di erences with equally spaced argu-

ments, Narumi shows how to extend Gauss, Bessel and Stirling univariate interpolation formulas for

equidistant points to the bivariate case by tensor product. He also applies the formulas he obtained

to approximate the values of bivariate functions, but he also mentions that some of his formulas had

been already used in [49].

In [45], the Newton formula (12) is obtained in the same way, with the corresponding remainder

formula (13). Moreover, Ste ensen considers a more general case, namely when for each i; 06i6m,

the interpolation points are of the form y0 ; : : : ; yni , with 06ni 6n. Now with a similar argument the

interpolating polynomial becomes

j’1

ni

m i’1

p(x; y) = f[x0 ; : : : ; xi ; y0 ; : : : ; yj ] (x ’ xh ) (y ’ xk ) (14)

i=0 j=0 h=0 k=0

28 M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35

with a slightly more complicated remainder formula. The most interesting particular cases occur

when ni = n, which is the Cartesian product considered above, and when ni = m ’ i. This triangular

case (triangular not because of the geometrical distribution of the interpolation points, but of the

indices (i; j)), gives rise to the interpolating polynomial

j’1

m m’i i’1

p(x; y) = f[x0 ; : : : ; xi ; y0 ; : : : ; yj ] (x ’ xh ) (y ’ xk ); (15)

i=0 j=0 h=0 k=0

that is

j’1

i’1

p(x; y) = f[x0 ; : : : ; xi ; y0 ; : : : ; yj ] (x ’ xh ) (y ’ xk ): (16)

06i+j6m h=0 k=0

Ste ensen refers for this formula to Biermann™s lecture notes [4] from 1905, and actually it seems

that Biermann has been the ÿrst who considered polynomial interpolation on the triangular grid in

a paper [3] from 1903 (cf. [44]) in the context of cubature.

Since the triangular case corresponds to looking at the “lower triangle” of the tensor product

situation only, this case can be resolved by tensor product methods. In particular, the respective

error formula can be written as

m+1 i’1 m’i

@m+1 f( i ; Ái ) h=0 (x ’ xh ) k=0 (y’ yk )

R(x; y) = : (17)

@xi @ym+1’i i! (m ’ i + 1)!

i=0

In the case of Cartesian product Ste ensen also provides the Lagrange formula for (12), which

can be obviously obtained by tensor product of univariate formulas.

Remainder formulas based on intermediate points ( i ; Ái ) can be written in many di erent forms.

For them we refer to Stancu™s paper [44] which also contains a brief historical introduction where the

author refers, among others, to [3,15,27,40,41]. Multivariate remainder formulas with Peano (spline)

kernel representation, however, have not been derived until very recently in [42] and, in particular,

in [43] which treats the triangular situation.

4. Salzer™s papers: from bivariate tables to general sets

In 1944, Salzer [33] considered the interpolation problem at points of the form (x1 + s1 h1 ; : : : ; x n +

sn hn ) where

(i) (x1 ; : : : ; x n ) is a given point in Rn ,

(ii) h1 ; : : : ; hn are given real numbers,

(iii) s1 ; : : : ; sn are nonnegative integers summing up to m.

This is the multivariate extension of the triangular case (16) for equally spaced arguments, where

ÿnite di erences can be used. Often, di erent names are used for the classical Newton interpolation

formula in the case of equally spaced arguments using forward di erences: Newton“Gregory, Harriot“

Briggs, also known by Mercator and Leibnitz, etc. See [18] for a nice discussion of this issue. In

[33], Salzer takes the natural multivariate extension of this formula considering the polynomial

q(t1 ; : : : ; tn ) := p(x1 + t1 h1 ; : : : ; x n + tn hn ) of total degree not greater than m in the variables t1 ; : : : ; tn ,

which interpolates a function f(x1 +t1 h1 ; : : : ; x n +tn hn ) at the points corresponding to ti =si , i=1; : : : ; n,

M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35 29

where the si are all nonnegative integers such that 06s1 + · · · + sn 6m. The formula, which is called

in [33] a multiple Gregory“Newton formula, is rewritten there in terms of the values of the function

f at the interpolation points, i.e., in the form

t1 tn m ’ t1 ’ · · · ’ tn

q(t1 ; : : : ; tn ) = ··· f(x1 + s1 h1 ; : : : ; x n + sn hn ): (18)

s1 sn m ’ s1 ’ · · · ’ sn

s1 +···+sn 6m

Note that (18) is the Lagrange formula for this interpolation problem. Indeed, each function

t1 tn m ’ t1 ’ · · · ’ tn

··· (19)

s1 sn m ’ s1 ’ · · · ’ sn

is a polynomial in t1 ; : : : ; tn of total degree m which vanishes at all points (t1 ; : : : ; tn ) with ti non-

negative integers 06t1 + · · · + tn 6m, except at the point (s1 ; : : : ; sn ), where it takes the value 1. In

particular, for n = 1 we get the well-known univariate Lagrange polynomials

t’i

t m’t

˜s (t) = =

s m’s s’i

06i6m;

i=s

for s = 0; : : : ; m.

Salzer used these results in [34] to compute tables for the polynomials (18) and, some years later

in [35], he studied in a similar form how to get the Lagrange formula for the more general case of

formula (16), even starting with this formula. He obtained the multivariate Lagrange polynomials

by a rather complicated expression involving the univariate ones.

It should be noted that several books related to computations and numerical methods published

around this time include parts on multivariate interpolation to some extent, surprisingly, more than

most of the recent textbooks in Numerical Analysis. We have already mentioned Ste ensen™s book

[45], but we should also mention Whittaker and Robinson [51, pp. 371“374], Mikeladze [26, Chapter

XVII] and especially Kunz [23, pp. 248“274], but also Isaacson and Keller [20, pp. 294 “299] and

Berezin and Zhidkov [2, pp. 156 “194], although in any of them not really much more than in [45]

is told.

In [36,37], Salzer introduced a concept of bivariate divided di erences abandoning the idea of

iteration for each variable x and y taken separately. Apparently, this was the ÿrst time (in spite

of the similarity with (11)), that bivariate divided di erences were explicitly deÿned for irregularly

distributed sets of points. Divided di erences with repeated arguments are also considered in [37] by

coalescence of the ones with di erent arguments. Since [36] was just a ÿrst attempt of [37], we only

explain the latter one. Salzer considers the set of monomials {xi yj }, with i; j nonnegative integers,

ordered in a graded lexical term order, that is,

(i; j) ¡ (h; k) ” i + j ¡ h + k or i + j = h + k; i ¿ h: (20)

Hence, the monomials are listed as

{1; x; y; x2 ; xy; y2 ; x3 ; : : :}: (21)

For any set of n + 1 points (xi ; yi ), Salzer deÿnes the associated divided di erence

n

[01 : : : n]f:= Ak f(xk ; yk ); (22)

k=0

30 M. Gasca, T. Sauer / Journal of Computational and Applied Mathematics 122 (2000) 23“35

choosing the coe cients Ak in such a form that (22) vanishes when f is any of the ÿrst n monomials

of list (21) and takes the value 1 when f is the (n + 1)st monomial of that list. In other words, the

coe cients Ak are the solution of the linear system

n

ij

xi yj any of the ÿrst n monomials of (21);

Ak xk yk = 0;

k=0

n

ij

xi yj the (n + 1)th monomial of (21):

Ak xk yk = 1; (23)

k=0

These generalized divided di erences share some of the properties of the univariate ones but not

all. Moreover, they have some limitations, for example, they exist only if the determinant of the

coe cients in (23) is di erent from zero, and one has no control of that property in advance. On

the other hand, observe that for example the simple divided di erence with two arguments (x0 ; y0 )

and (x; y), which is

f(x; y) ’ f(x0 ; y0 )

;

x ’ x0

gives, when applied to f(x; y) = xy, the rational function

xy ’ x0 y0

x ’ x0

and not a polynomial of lower degree. In fact, Salzer™s divided di erences did not have great success.

Several other deÿnitions of multivariate divided di erences had appeared since then, trying to keep

as many as possible of the good properties of univariate divided di erences, cf. [16].

5. Reduction of a problem to other simpler ones

Around the 1950s an important change of paradigm happened in multivariate polynomial inter-

polation, as several people began to investigate more general distributions of points, and not only

(special) subsets of Cartesian products. So, when studying cubature formulae [32], Radon observed

the following in 1948: if a bivariate interpolation problem with respect to a set T ‚ R2 of ( k+2 )

2

2

interpolation points is unisolvent in k , and U is a set of k + 2 points on an arbitrary straight

line ˜ ‚ R2 such that ˜ © T = …, then the interpolation problem with respect to T ∪ U is unisolvent

2

in k+1 . Radon made use of this observation to build up point sets which give rise to unisolvent

interpolation problems for m recursively by degree. Clearly, these interpolation points immediately

yield interpolatory cubature formulae.

The well-known BÃ zout theorem, cf. [50], states that two planar algebraic curves of degree m and

e