4.5.3 Viskovatov's Algorithm for the General Case

In order to generalize Viskovatov's algorithm to degenerate cases, one must make an inspired choice for the representation of the particular \L/M] Pade aDDroximant as a continued fraction. We consider

which should be compared with the general C-fraction (5.12) and Magnus P-fraction (5.16). In (5.17), each p,(z) is a polynomial with ^(0)^=0 and each a; is a positive integer: each p^z) and a, are fully defined in (5.21) and (5.22). As usual, we suppose that we are given the series expansion

and that we are requested to find the [L/M] Pade approximant of (5.18) with L>M. If L<M, we should remove the factor zk from (5.18), where k is the highest power such that zk factors (5.18), and start with its reciprocal series. We choose to consider an infinite series in (5.18) for convenience of representation; in fact only L + M+ 1 terms of (5.18) are required, and questions of convergence do not arise. Following Werner [1979, 1980], we define the generalized Viskovatov reduction:

We may then define the polynomial Pj(z) of degree strictly less than a by

Initialization. Let

and

Iteration. For j=0,1,2,..., proceed recursively. Define a] to be the least integer for which

Our aim is to find the [/-/w -] Pade approximant of the left-hand side of (5.22). If at >lj + m , then let n=j. In this case, (5.17) is the Pade approxi­mant required, as is proved in Theorem 4.5.1. If

the [L/M] Pade approximant to the given series (5.18) does not exist, as is proved in Theorem 4.5.1.

Next, we consider the reciprocal series defined by

so that, formally at least,

In order to find the [lj/m ] Pade approximant of the left-hand side of (5.25), we seek the [my//y—ay] Pade approximant of the left-hand side of (5.24V Hence we define

and

The requirement that

is guaranteed by (5.21), and also

in view of (5.23V Hence we need to find the

of the series

defined by (5.24) and (5.26), and the process

becomes recursive.

Having defined the iterative process for the construction of (5.17), we must verify that it fulfills its specification. From (5.26), we note that

and therefore

where deg{/?(z)} is the degree of p(z). We see that the fraction (5.17) has at most L + M divisors, which occurs if a,— 1 for / = 0,1,2,..., n, and hence the process terminates.

Theorem 4.5.1. The generalized Viskovatov reduction defined in (5.19)— (5.26) is reliable: if the [L/M] Pade approximant exists, the process finds if, if not, (5.23) is satisfied and indicates this fact.

Proof Suppose that the [L/M] approximant to (5.20) exists, yet the algorithm fails. Let the approximant have the Maclaurin expansion

for

and

derived by the reduction is identical when applied to

itself. Initially, with

we find

so that

Since the reduction process (5.25) for calculating the [L/M] approximant is based on the series-expansion coefficients only, the representation (5.17)

This division of

yields 0n(z) as the quotient anc

as the remainder. Unless on(z) itself is the [L/ Wm] Padé approximant of

/(z), (5.21) ensures that

and so we may use the Maclaurin

expansion

We sunnose that

in (5.31) is reduced sequentially by the

process indicated by (5.19)-(5.26) as far as possible. Let

fory=0,1k. Hence we reduce the Padé approximant sequentially by

where we define

and

Suppose that the process fails at stage k, namely because (5.23) is satisfied with j=k. From (5.33),

where lk >ak >lk +mk. Because the left-hand side of (5.37) has degree at most lk, then either ak *Zlk or else T(k\z) = 0. This is incompatible with a failure at stage k. Consequently the reduction process cannot fail if the [L/M] Padé approximant exists.

The proof of the converse result, that R(0)(z) constructed by (5.17)-(5.26 > is in fact the [L/M] Padé approximant, uses the definitions of these equations. From (5.22) and (5.28), we see that

where

we define

and for

We define

We see directly from (5.40) that

Let us suppose that

It then follows from (5.25) and (5.41) that (5.42) is true with k^k~\. Using (5.39), we have proved that

by induction. This completes the proof of the theorem.

variable

is also a factor ol

but it is obvious that each common factor of

Hence we see

that we may view the reduction process and Theorem 4.5.1 as a proof, using

In the course of the proof we can see the connection between (5.31). (5.331 and the Euclidean algorithm. The connection is direct using the

the Euclidean algorithm, that P(0\z) and (?<0|(z ) have no common factors. Of course, the proof in Theorem 1.4.3 is more direct.

There is also a connection between (5.33) and the recurrence relations (2.14). From (5.33), we find that

for

We may suppose that the relative normalization

is fixed by defining

From (5.35), (5.36), and (5.44), it follows that

and

which are remarkably similar to (5.2.14) when viewed as a recurrence for j=n—\,n — 2,...,l. In view of (5.34), (5.41), and (5.45), the recurrence defines the numerator and denominator of

However, the recurrence

(5.48)

initialized by (5.45), is probably the more convenient

form.

It is now clear how to define a modified Viskovatov algorithm, based on (5.17)-(5.26) and incorporating a generalization of (3.1) for efficiency. In fact, we may as well consider the problem of finding an [L/M] Padé approximant, with L>M, for

In (5.49), the coefficients gj°\ dj0) for / = 0.1,2           L + M are given, and

^ 7^0. This problem is solved using bigradients in Section 1.6. We outline a generalized Viskovatov algorithm, emphasizing only the alterations to (5.17)—(5.26) [Bultheel, 1980b],

Initialization. Define l0—L, m0=M.

Recurrence. The process is iterative for7 = 0,1,2,.... We construct the polynomial pj{z) of degree at most lJ — mJ such that

(5.50)

with a ■ > lj —mr In fact, we need to construct the coefficients of

(5.511

coefficients of the series

(5.521

The construction of the coefficients of (5.51) and (5.52) is performed iteratively, according to the algebra implied by the following analysis. Let

For k = 0,1,...,/—m define

and

HencePj(z) is specified by (5.51), and for the coefficients of (5.52) we have

(5.53)

We define a; as the least integer such that (5.53) holds with h\tn 7^0, and hence the coefficients of (5.52) are fully specified.

If

Otherwise, let

or if (5.53) vanishes exactly, the iteration is terminated

the Padé approximant required is nonexistent.

and

for

In this case,

the /th stage of the reduction

process (5.54) may be iterated. The algorithm is now fully specified.

The generalized Viskovatov algorithm in (5.49)-(5.54) is designed to yield an [L/M] Padé approximant reliably. It follows a path on the Padé table which avoids nonexistent entries in blocks in the table. We conclude that the path that the convergents of the fraction (5.17) follows in the Padé table depends on the values of the coefficients in (5.49) as well as the values of L and M. As a computational algorithm, it can be recommended wherever exact arithmetic is available and the values of a0, a,, ...,«„_, are unambigu­ously defined. The techniques of the generalized Viskovatov algorithm can be adapted to rational interpolation using a generalized Thiele fraction (see Part II, Section 1.1, and Werner [ 1979a, b]) except that the algorithm must incorporate a facility for reordering the interpolation points if necessary [Thacher and Tukey, 1960; Graves-Morris and Hopkins, 1978; Graves- Morris, 1980a], Similar techniques based on the Euclidean algorithm allow Kronecker's algorithm to be generalized to avoid inessential degeneracies (see Section 2.4; Graves-Morris [1980a]).