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






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








The requirement that


is guaranteed by (5.21), and also

Pade approximant

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










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






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






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




 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.




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







We may suppose that the relative normalization


is fixed by defining





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






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




initialized by (5.45), is probably the more convenient


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


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



coefficients of the series



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






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


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.





Otherwise, let

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

the Padé approximant required is nonexistent.









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]).