# Computation of the c-Table Related to the Padé Approximation

Radoslaw Jedynak, Jacek Gilewicz

<figure><img src="image_1.png" /></figure>

## To cite this version:

Radoslaw Jedynak, Jacek Gilewicz. Computation of the c-Table Related to the Padé Approximation. Journal of Applied Mathematics, 2013, 2013, pp.185648. 10.1155/2013/185648. hal-01278031

HAL Id: hal-01278031

https://hal.science/hal-01278031

Submitted on 2 May 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Research Article

Computation of the c-Table Related to the Padé Approximation

Radosław Jedynak<fnref n="1" /> and Jacek Gilewicz<fnref n="2" />,<fnref n="3" />

Correspondence should be addressed to Jacek Gilewicz; gilewicz@cpt.univ-mrs.fr

Received 27 June 2013; Revised 13 September 2013; Accepted 24 September 2013

Academic Editor: Nicola Mastronardi

Copyright © 2013 R. Jedynak and J. Gilewicz. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

The aim of this paper is to give a complete and practical method for numerical application of Padé approximation with the help of the c-table analysis. We present an exhaustive list of useful formulas to compute a c-table related to a formal power series C(z) = Σ∞n=0 c_n z^n. Some of these formulas are not widely known, because they were presented in publications of limited circulation. Some others were never published, as three symmetric Paszkowski-like formulas to overcome the blocks in a c-table or an extension of local error formula for Padé approximants in the blocks. All formulas are given in two versions: in terms of Toeplitz determinants (c-table) and in the version of Hankel determinants (c-table). We compare the theory with numerical observations by reproducing different computational aspects of software producing the c-tables with the presence of blocks and their evolution following the evolution of computer environment.

1. Introduction

The Padé approximation method is largely used to solve
many problems of numerical analysis such as the convergence
acceleration, analytic continuation of complex functions,
moment problems, and numerical integration [1–3] and, in
general, to approximate functions of the complex variable
represented by a truncated power series and the detection
of their zeros and poles. This method is also commonly
applied to solve numerous problems in physical modeling [4].
We can find some historical background and fundamental
concepts of this method in [5] and the recently published
book by Trefethen [6]. In [7] the authors present the methods
used to construct matrix Padé approximants if the coefficient
matrices of the input matrix polynomial are triangular. The
extended Euclidean algorithm is applied to solve this prob-
lem. The application of Padé approximants in computational
problems starts frequently by the computation of the auxiliary
table, called c-table [8]. The entries of c-table are the Toeplitz
determinants of matrices of linear systems defining the
denominators of Padé approximants. The square blocks of
zeros in the c-table indicate the existence of corresponding
blocks in the table of Padé approximants. The so-called

valleys in the c-table are the lines of minimal absolute values
of entries, which indicate the lines of best Padé approximants
(BPA) in the Padé table [9]. Because the computation of the
c-table is simpler, it is recommended to begin by this com-
putation to obtain global preliminary information about the
interesting Padé approximants before their own calculations.

Let us recall some definitions and results related to Padé approximation.

Let

$$
C(z) = \sum_{n=0}^{\infty} c_n z^n, \quad c_0 \neq 0 \qquad (1)
$$

be a formal power series of a function *f*. It can be a Maclaurin
series of an analytic function at *z* = 0 or an asymptotic series,
such as a Stieltjes series having a zero radius of convergence
[10, 11].

Definition 1. Let

$$
\begin{gather*}
\left\{ \frac{m}{n} \right\} (z) = \frac{P_m(z)}{Q_n(z)} = \frac{p_0 + p_1 z + \cdots + p_m z^m}{q_0 + q_1 z + \cdots + q_n z^n}, \quad (2) \\
\deg P_m \le m, \quad \deg Q_n \le n
\end{gather*}
$$

<footnote n="1">¹ Uniwersytet Technologiczno-Humanistyczny im. K. Pulaskiego, ul. Malczewskiego 20a, 26-600 Radom, Poland</footnote>

<footnote n="2">² Aix Marseille Université, CPT CNRS, UMR 7332, Luminy, Case 907, 13288 Marseille Cedex 9, France</footnote>

<footnote n="3">³ Université de Toulon, CNRS, CPT, UMR 7332, 83957 La Garde, France</footnote>
<table><caption>TABLE 1</caption><thead><tr><th>m</th><th colspan="3">n</th></tr><tr><th></th><th>0</th><th>1</th><th>2</th><th>...</th></tr></thead><tbody><tr><th>0</th><td>{0/0}</td><td>{0/1}</td><td>{0/2}</td><td>...</td></tr><tr><th>1</th><td>{1/0}</td><td>{1/1}</td><td>...</td><td>...</td></tr><tr><th>2</th><td>{2/0}</td><td>...</td><td>...</td><td>...</td></tr><tr><th>...</th><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

denote in general a rational function. The Padé table (p-table) of the formal power series (1) is a doubly infinite array of irreducible rational functions {$m/n$} called reduced Padé forms, determined in such a manner that the Maclaurin expansion of {$m/n$} agrees with $C(z)$ as far as possible.

The first column in Table 1 contains the truncated series (1). Because a reduced Padé form can be normalized by dividing all coefficients by one coefficient $\neq 0$, then {$m/n$} contains only $m+n+1$ unknowns which can be determined by $m+n+1$ first coefficients of $C(z)$ or, equivalently, by $m+n+1$ first derivatives of $f$ at $z=0$:

$$ \left\{ \frac{m}{n} \right\} (0) = f(0), \qquad (3) $$

$$ \left\{ \frac{m}{n} \right\}' (0) = f'(0), \dots, \quad \left\{ \frac{m}{n} \right\}^{(m+n)} (0) = f^{(m+n)} (0). \qquad (4) $$

This is usually expressed by the condition

$$ f(z) - \frac{P_m(z)}{Q_n(z)} = O(z^{m+n+1}) \qquad (5) $$

which means that the power series expansion of $P_m(z)/Q_n(z)$ matches the power series expansion of $f$ at $z=0$ at least up to the power $m+n$.

**Definition 2** (see Gilewicz [11, page 171]). If the reduced Padé form {$m/n$} satisfies (5), then it is a Padé approximant and is denoted by $[m/n]: \{m/n\} = [m/n]$.

**Definition 3.** The power series and associated p-table are said to be normal if every element of the table is different from any other element. Each Padé approximant $[m/n]$ reproduces the series $C(z)$ exactly up to the term in $z^{m+n}$.

Multiplying (5) by $Q_n$ we obtain

$$ Q_n(z) f(z) - P_m(z) = O(z^{m+n+1}) \qquad (6) $$

which leads to the separated linear system for the coefficients $q_j$ and $p_j$ ($c_i \equiv 0$ if $i < 0$):

$$ -p_k + \sum_{j=0}^{n} c_{k-j} q_j = 0 \quad k = 0, 1, \dots, m, \qquad (7) $$

$$ \sum_{j=0}^{n} c_{k-j} q_j = 0 \quad k = m+1, \dots, m+n. \qquad (8) $$

The system (8) of $n$ equations for ($n+1$) unknowns $q_j$ always a nontrivial solution (Frobenius) which determines the reduced Padé form. Let us note that (6) is equivalent to (5) only if $Q_n$ is reversible, that is, if $Q_n(0) \neq 0$. Following this remark Baker presented an equivalent definition of the Padé approximant as follows.

**Definition 4** (see Baker [10, page 5]). The Padé approximant is a solution of linear systems (8) and (7) with the condition $Q_n(0) = q_0 = 1$ ($c_i \equiv 0$ if $i < 0$):

$$ -p_k + \sum_{j=1}^{n} c_{k-j} q_j = -c_k \quad k = 0, 1, \dots, m, \qquad (9) $$

$$ \sum_{j=1}^{n} c_{k-j} q_j = -c_k \quad k = m+1, \dots, m+n. \qquad (10) $$

The power series converges only inside the convergence circle; the Padé approximants reproduce zeros and poles of meromorphic functions and give a good approximation of these functions also outside the circle of convergence. Moreover, in particular cases where the series represents a rational function, the Padé approximant method finds this function automatically. In other cases the BPA algorithms allow finding excellent approximations. Padé approximant is the best local rational approximation of $f$ in the vicinity of the point of expansion of $f$ in the power series. Of course, we are not limited by the point $z=0$ and if we are interested in other regions in $\mathbb{C}$ we can compute the Padé approximants at arbitrary point $z_0$ of analyticity of $f$ knowing the coefficients of the series $\sum a_j(z-z_0)^j$.

A determinant of the matrix of the system (10) is the Toeplitz determinant:

$$ m \ge 0, n \ge 1: \quad C_n^m = \begin{vmatrix} c_m & c_{m-1} & \cdots & c_{m-n+1} \\ c_{m+1} & \cdot & \cdots & \cdot \\ \cdot & \cdot & \cdots & \cdot \\ \cdot & \cdot & \cdots & \cdot \\ c_{m+n-1} & \cdot & \cdots & c_m \end{vmatrix}, \qquad (11) $$

where we put $c_k \equiv 0$ if $k < 0$. Defining $C_0^m := 1$ we can build an infinite table of $C_n^m$'s, called c-table.

The second column of Table 2 contains the coefficients of the power series $C_1^m = c_m$. The first row contains the powers of $c_0$, $C_n^0 = (c_0)^n$ and the second row can be computed recursively by expanding $C_n^1$:

$$ C_n^1 = \sum_{j=1}^{n} (-c_0)^{j-1} c_j C_{n-j}^1. \qquad (12) $$

The c-table was first introduced by Gragg [8]. Baker [10] used an alternative definition permuting rows in (11), calling the
<table><caption>TABLE 2</caption><thead><tr><th>m</th><th colspan="5">n</th><th>m</th><th colspan="4">n</th></tr><tr><th></th><th>0</th><th>1</th><th>2</th><th>3</th><th>...</th><th></th><th>0</th><th>1</th><th>2</th><th>...</th></tr></thead><tbody><tr><th>0</th><td>1</td><td>c<sub>0</sub></td><td>(c<sub>0</sub>)<sup>2</sup></td><td>(c<sub>0</sub>)<sup>3</sup></td><td>...</td><td></td><td>0</td><td>C<sub>0</sub><sup>0</sup></td><td>C<sub>1</sub><sup>0</sup></td><td>C<sub>2</sub><sup>0</sup></td></tr><tr><th>1</th><td>1</td><td>c<sub>1</sub></td><td>...</td><td>...</td><td>...</td><td>≡</td><td>1</td><td>C<sub>0</sub><sup>1</sup></td><td>C<sub>1</sub><sup>1</sup></td><td>C<sub>2</sub><sup>1</sup></td></tr><tr><th>2</th><td>1</td><td>c<sub>2</sub></td><td>...</td><td>...</td><td>...</td><td></td><td>2</td><td>C<sub>0</sub><sup>2</sup></td><td>C<sub>1</sub><sup>2</sup></td><td>C<sub>2</sub><sup>2</sup></td></tr><tr><th>...</th><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td></td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

resulting Hankel determinants $C(m/n)$ and the corresponding table a c-table. An obvious relation between $C(m/n)$ and $C_n^m$ is

$$C_n^m = (-1)^{n(n-1)/2} C\left(\frac{m}{n}\right),$$

$$C\left(\frac{m}{n}\right) = \begin{vmatrix} c_{m-n+1} & \cdots & c_m \\ c_{m-n+2} & \cdots & c_{m+1} \\ \vdots & \ddots & \vdots \\ c_m & \cdots & c_{m+n-1} \end{vmatrix} \quad (13)$$

Interest in the c-table is due to its direct relation to the Padé table [8, 11] and to the particular ease of its computation. The c-table contains information about the block structure of the Padé table. Each square block of zeros in the c-table defines a corresponding block in the Padé table. Moreover, a dominant term of the error of Padé approximant located outside a block is given by a ratio of two Toeplitz determinants:

$$f(z) - \frac{P_m(z)}{Q_n(z)} = (-1)^n \frac{C_{n+1}^{m+1}}{C_m^n} z^{m+n+1} + \cdots, \quad (14)$$

A theory of valleys in the c-table [9] leading to a numerical algorithm of choice of the best Padé approximant [11] is based on property (14).

In general a p-table presents a square block structure. We denote by $(m, n; k)$ a $k \times k$ square block of all equal entries with a west-north corner $\{m/n\}$. In each block

$$\forall i, j < k: \left\{ m + \frac{i}{n} + j \right\}_{i+j \le k-1} = \left[ \frac{m}{n} \right] \quad (15)$$

are Padé approximants and satisfy

$$f(z) - \left[ \frac{m}{n} \right](z) = O\left(z^{m+n+k}\right), \quad (16)$$

but the reduced forms located under the antidiagonal of a block are not Padé approximants:

$$\forall i, j < k: \left\{ m + \frac{i}{n} + j \right\}_{i+j \ge k} \neq \left[ \frac{m}{n} \right] \quad (17)$$

because $(m+i)+(n+j)+1 > m+n+k$ (see (16)). If the p-table presents a $(m, n; k)$ block, the c-table presents a $(k-1) \times (k-1)$ block of zeros with the west-north corner on the position $(m+1), (n+1)$. Now the following property becomes evident.

**Property 1.** The p-table is normal if all $C_n^m \neq 0$.

Let us complete the error formula (14) by the corresponding formula for $[m/n]$ belonging to the block $(m, n; k)$:

$$f(z) - \left[ \frac{m}{n} \right](z) = rz^{m+n+k} + \cdots, \quad (18)$$

where

$$\begin{align*}
r &= \frac{C_{n+1}^{m+k}}{C_n^{m+k-1}} (-1)^n \\
&= \frac{C_{n+k}^{m+1}}{C_{n+k-1}^m} (-1)^{n+k-1} \tag{19} \\
&= \frac{C(m+k/n+1)}{C(m+k-1/n)} \\
&= \frac{C(m+1/n+k)}{C(m/n+k-1)}. \notag
\end{align*}$$

This formula given in [10, 12] is given here in terms of Toeplitz and Hankel determinants. An extension of representations of $r$ using other determinants around the zeros block in the c-table will be presented in the next section.

The existence of an infinite block $(m, n; \infty)$ means that the series $C(z)$ represents a rational function $[m/n]$. If a p-table presents a block $(m, n; k)$, then the $[m/n]$ Padé approximant computed with $m+n+1$ coefficients $c_0, c_1, \dots, c_{m+n}$ located on the antidiagonal $m+n=\text{const}$ in the p-table is clearly the best Padé approximant (BPA) on this antidiagonal, because it is the only approximant which reproduces exactly the series $C(z)$ up to the power $z^{m+n+k-1}$. A number of algorithms of choice of BPA are presented in [11]. The detailed modern theory of Padé approximants and the exhaustive list of algorithms of their calculation are presented in [10, 11]. The entries of p-table can be calculated directly by solving the linear system (10) ((9) is restricted only to the substitutions) or by more efficient recursive algorithms adapted to the Toeplitz systems. One of these is based on the Wynn formula

$$W \begin{array}{ccc} & N & \\ C & E & \\ S & & \end{array} \frac{1}{W-C} + \frac{1}{E-C} = \frac{1}{N-C} + \frac{1}{S-C} \quad (20)$$

which permits computing the $E$ elements from four others (ascending algorithm). This method is too laborious and so it is simpler to discover the block structure of the p-table by analyzing the block structure of the c-table. In numerical practice we are always interested in BPA and a more rapid way to determine the position of BPA passes through the analysis of the c-table.
## 2. Computation of the c-Table in the Normal Case: Numerical Recommendations

The particular form of Toeplitz (resp. Hankel) determinants permits computing them recursively using the Sylvester formula

$$ W \begin{array}{ccc} & N & \\ C & E & NS + EW = C^2 \quad (\text{for } C_n^m), \\ & S & \end{array} \qquad (21) $$

$$ NS - EW = C^2 \quad (\text{for } C\left(\frac{m}{n}\right)) $$

avoiding the laborious direct computation by (11) (resp. (13)). Starting from the two first columns one can compute the “East” $C_n^m$ elements by the *ascending algorithm* [13]:

$$ E = \frac{C}{W}C - \frac{S}{W}N. \qquad (22) $$

Alternatively, starting from the first two rows one can compute the “South” elements by the *descending algorithm*:

$$ S = \frac{C}{N}C - \frac{E}{N}W. \qquad (23) $$

The above order of arithmetic operations diminishes a risk of overflows or underflows. We have studied in [13] the complexity of both the above algorithms counting multiplications and divisions (including that of (12) needed for the descending algorithm) necessary to obtain the $C_n^m$ from the normalized sequence ($c_j$), that is, with $c_0 = 1$. The cost of the normalization of each coefficient $c_1, c_2, \dots$ in the second column of the c-table is equal to 1. Then the costs in the second row are 0, 0, 1, 2, ..., $n-1$, ... (12). Then the cost $p_E$ of $E$ computed by (22) is $p_E = p_C + p_W + p_S + p_N + 4$, where $p_C, p_W, p_S$, and $p_N$ are the previous costs of $C, W, S$, and $N$, respectively. In fact it is the cost of the calculation of all $n^2$ determinants disposed in the triangle pointed on the last $C_n^m$ if it is situated under or on the diagonal. Analogously, the global cost of the calculation of all $m^2$ determinants with the last $C_n^m$ by the descending algorithm is $p_S = p_C + p_N + p_E + p_W + 4$. The following separation line optimizes the cost of each algorithm.

In the following example we give the number of operations needed by the descending algorithm followed by the number of operations needed by the ascending one:

$$ C_{11}^8: \quad 105\,266, 159\,995 \qquad C_{11}^9: \quad 349\,092, 160\,070. \qquad (24) $$

The precision of $E$ computed by (22) is, crudely estimating, proportional to the global cost $p_E = 2p_C + 2p_W + p_S + p_N + 5$, where $p_C, p_W, p_S$, and $p_N$ are the previous global costs. We count twice $P_C$ and $P_W$ because each presence of an element in the arithmetic expression introduces its error. On the other hand if we wish to optimize the precision, we must take into account all arithmetic operations, including additions and subtractions. The errors of differences play a dominant role, but they are not considered here. Consequently the separation line will be modified. The precision observed

<table><caption>TABLE 3</caption><thead><tr><th>Row</th><th>AA</th><th>DA</th><th>Exact values</th></tr></thead><tbody><tr><td>8</td><td><u>.1345 - 54</u></td><td>.1345 - 54</td><td>.1344 - 54</td></tr><tr><td>9</td><td><u>.3457 - 64</u></td><td>.3445 - 64</td><td>.3465 - 64</td></tr><tr><td>10</td><td><u>-.1041 - 73</u></td><td>-.1073 - 73</td><td>-.1035 - 73</td></tr><tr><td>11</td><td><u>-.6476 - 83</u></td><td>-.5203 - 83</td><td>-.6054 - 83</td></tr><tr><td>12</td><td><u>-.4552 - 92</u></td><td>.1140 - 91</td><td>.1443 - 91</td></tr><tr><td>13</td><td><u>.1653 - 97</u></td><td>.9137 - 98</td><td>.1387 - 98</td></tr><tr><td>14</td><td><u>-.9270 - 103</u></td><td>-.7552 - 103</td><td>.9683 - 105</td></tr><tr><td>15</td><td><u>-.2095 - 108</u></td><td>-.6141 - 108</td><td>.2715 - 110</td></tr><tr><td>16</td><td><u>-.2044 - 111</u></td><td>-.1911 - 111</td><td>.2241 - 115</td></tr></tbody></table>

numerically using multiprecision calculus (see details in [14]) showed that the separation line mounts a few places giving a clear advantage to the ascended algorithm.

*Example 5.* Fragment of the 13th column of the c-table of the series of the Stieltjes function $-(1/x) \log(1-x)$ is computed by ascending algorithm (AA) and descending algorithm (DA) in double precision (see Table 3).

The underlined elements (Example 5) computed by the ascending algorithm situated in the region of accuracy of the descending algorithm are a little better than those computed by descending algorithm. The exact values were determined by multiprecision calculus with 32 digits.

## 3. Overflow, Underflow, and Detection of Blocks: Numerical Recommendations

No problems arise in a numerical computation if all coefficients $c_j$ are integers [8, 12]. If a block exists, the integer arithmetic produces zeros exactly. Unfortunately in a real case the computational problem is generally much more difficult. Numerical experiments show that $C_n^m$'s decrease or increase quite rapidly. The risk of the overflow or underflow appears approximately for $m + n \ge 3p$, where $p$ is the number of decimal digits of the floating point representation ([11, page 344]). Fortran subroutines given in [13, 15] use masks allowing the determination of exponents of real constants composing of the arithmetic expression (22) or (23) before their evaluation. After that the program checks if a global exponent is inside an authorized range or not and, if so, if it authorizes the evaluation or not. Knowing the computer used it is possible to overcome overflows without stopping an execution of the program, but the use of this technique embarrasses the software portability.

The underflow problem is related to another important question: is the computed value of $C_n^m$ equal to zero or not? The Vignes permutation-perturbation method to detect the so-called “zero informatique” in French [16] gives a satisfactory statistical estimation, allowing one to decide if a determinant vanishes or not, but in our practical computational problems its cost is too great. The crude estimation is the following:

$$ \text{if } |C_n^m| < \text{base}^{-(m+n)}, \quad \text{then } C_n^m \approx 0, \qquad (25) $$
<table><caption>TABLE 4</caption><thead><tr><th colspan="5" style="text-align: center;">(a)</th></tr><tr><th style="text-align: center;">m</th><th style="text-align: center;">4</th><th style="text-align: center;">5</th><th style="text-align: center;">6</th><th style="text-align: center;">7</th></tr></thead><tbody><tr><th style="text-align: center;">3</th><td style="text-align: center;">1. - 8</td><td style="text-align: center;">1. - 10</td><td style="text-align: center;">1. - 11</td><td style="text-align: center;">1. - 12</td></tr><tr><th style="text-align: center;">4</th><td style="text-align: center;">1. - 10</td><td style="text-align: center;">1. - 11</td><td style="text-align: center;">1. - 12</td><td style="text-align: center;">1. - 13</td></tr><tr><th style="text-align: center;">5</th><td style="text-align: center;">1. - 11</td><td style="text-align: center;">1. - 12</td><td style="text-align: center;">1. - 13</td><td style="text-align: center;">1. - 14</td></tr><tr><th style="text-align: center;">6</th><td style="text-align: center;">1. - 12</td><td style="text-align: center;">1. - 13</td><td style="text-align: center;">1. - 14</td><td style="text-align: center;">1. - 16</td></tr><tr><th style="text-align: center;">7</th><td style="text-align: center;">1. - 13</td><td style="text-align: center;">1. - 14</td><td style="text-align: center;">1. - 16</td><td style="text-align: center;">1. - 17</td></tr></tbody><thead><tr><th colspan="5" style="text-align: center;">(b)</th></tr><tr><th style="text-align: center;">m</th><th style="text-align: center;">4</th><th style="text-align: center;">5</th><th style="text-align: center;">6</th><th style="text-align: center;">7</th></tr></thead><tbody><tr><th style="text-align: center;">3</th><td style="text-align: center;">6.28 + 03</td><td style="text-align: center;">-5.60 + 04</td><td style="text-align: center;">4.99 + 05</td><td style="text-align: center;">-4.45 + 06</td></tr><tr><th style="text-align: center;">4</th><td style="text-align: center;">6.61 - 01</td><td style="text-align: center;">2.19 + 00</td><td style="text-align: center;">7.18 + 00</td><td style="text-align: center;">2.34 + 01</td></tr><tr><th style="text-align: center;">5</th><td style="text-align: center;">6.61 - 07</td><td style="text-align: center;">-6.61 - 07</td><td style="text-align: center;">6.61 - 07</td><td style="text-align: center;">-6.61 - 07</td></tr><tr><th style="text-align: center;">6</th><td style="text-align: center;">9.06 - 10</td><td style="text-align: center;">6.55 - 26</td><td style="text-align: center;">-2.26 - 26</td><td style="text-align: center;">7.48 - 27</td></tr><tr><th style="text-align: center;">7</th><td style="text-align: center;">1.24 - 12</td><td style="text-align: center;">-3.11 - 29</td><td style="text-align: center;">3.48 - 47</td><td style="text-align: center;">9.31 - 47</td></tr></tbody></table>

where base denotes the base of representation of the real constant in the computer (essentially 2, 10, or 16). The second method, purely empirical, is due to Guziński [17]. Guziński follows the monotonicity of the ascending or descending sequences: ($|C_{n-j}^{m+j}|$, $j = n, n-1, \dots, 1, 0$ or $j = -m, -m+1, \dots, 0$, respectively. Using the 32-bit arithmetic (single precision) he estimates that the true value of determinant (let us call it “new”) is equal to zero if the previous monotonicity changes significantly, that is, if the new determinant with respect to the previous determinant satisfies the following inequality:

$$ \text{new} < \text{previous} * 5 * 10^{-12}. \qquad (26) $$

We recommend this efficient method which allows the detection of the zeros blocks in c-tables very quickly. Of course, it must be readapted to the used real constant representation. These two criteria, (25) and (26), are well-illustrated by Example 6, where the second table represents a fragment of c-table for a series produced by a rational function:

$$ f(x) = \left[ \frac{5}{4} \right] = \frac{1}{(1 - (x/3))^2} - \frac{9x^3}{(1 - (x/9))^2}. \qquad (27) $$

The numerical constants are denoted shortly by decimal exponents: for instance $0.9 \times 10^{-5}$ is denoted by $.9 - 5$. The first table represents the approximative limits (25) for the base $= 16 : 16^{-(m+n)} \approx 10^{-(6/5)(m+n)}$. The underlined elements correspond to the beginning of the infinite block of zeros starting at element $C_5^6$. As shown here both criteria work: the underlined elements on the right (c-table) are less than the corresponding limits on the left and also Guziński criterion is satisfied for these elements.

Example 6 (zeros block in the c-table). If the first zero is detected by the ascending algorithm, it corresponds to the west-north corner of the block of zeros. Then, we stop at the

moment of the ascending computation of this antidiagonal and we return to the next antidiagonal to verify the value of the determinant situated under this first zero. Suppose that we find $k-1$ zeros starting from $C_{n+1}^{m+1}$ which indicate the existence of the block $(m, n; k)$. The next section is devoted to the strategy of following the calculus.

The “normal” monotonicity, by opposition to the jumps (26), can be tested by reference to the normal c-table corresponding, for instance, to the Stieltjes function (see Table 4)

$$ f(x) = -\frac{1}{x} \log(1-x) = 1 + \frac{1}{2}x + \frac{1}{3}x^2 + \dots \qquad (28) $$

Example 7 presents a fragment of the c-table for a series produced by the mentioned function.

Example 7 (fragment of the c-table of the Stieltjes function $f - (1/x) \log(1-x)$). The minima on the antidiagonals are underlined. They correspond to the positions of the best approximants on the corresponding antidiagonals in the Padé table. Joining these minima we obtain the valley structure in the c-table mentioned in Section 6 (see Table 5).

Using the 64-bit arithmetic (double precision) the risk of overflows or underflows is rather marginal [14]. Both standards of mentioned representations of numerical values are specified by the IEEE Standard for Floating-Point Arithmetic (IEEE 754).

## 4. Nonnormal Case: Froissart-Gilewicz and Paszkowski Identities to Overcome Blocks

The proof of geometrical progression of elements surrounding the zeros blocks in the c-table and an analytic proof of the identity allowing the computation of the east or south elements of these blocks when the Sylvester identity fails are presented in [11]. The purely algebraic proof of an equivalent identity is given by Paszkowski in [12]. Three new Paszkowski-like identities and an extension of formula (19) are proved and are given in this section. Figure 2 illustrates the c-table corresponding to the block $(m, n; k)$ with $k = 8$ and two “shells” adjacent to the $7 \times 7$ zeros block.

Elements surrounding the blocks of zeros form geometrical sequences (attention: the sense of west and north arrows and the superscripts $+$, $-$ on the south and east sides are inverted with respect to [11, 12] and now correspond to the sense of calculation of the elements of the c-table by presented algorithms) from the corners as indicated by arrows and with respective ratios $\hat{S}$, $\hat{W}$, $\hat{N}$, and $\hat{E}$. Then, it suffices to know four entries $C_n^m$, $C_{n+1}^m$, $C_n^{m+1}$ (these three are already computed before the detection of the block), and $C_{n+1}^{m+k}$ (first nonzero element after $k-1$ zeros on the $n+1$ column which identifies the block) to calculate all elements surrounding the block of zeros in the following way:

$$ \hat{S} = \frac{C_{n+1}^{m+k}}{C_n^{m+k}}, \qquad C_{n+i}^{m+k} = C_n^{m+k} * \hat{S}^i, $$

$$ \hat{W} = \frac{C_n^{m+1}}{C_n^m}, \qquad C_n^{m+i} = C_n^m * \hat{W}^i, $$
<figure><img src="image_3.png" /><figcaption>FIGURE 1: Separation line optimizing the cost of each algorithm.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2: Location of elements in (m, n; 8) block for Gilewicz-Froissart and Paszkowski identities.</figcaption></figure>

row which allows following the use of respective, ascending, or descending algorithm. In fact, the Sylvester algorithms fail only in the rare situations illustrated in Figure 5. Let A be an intersection of an east shadow region produced by a block (m, n; 4) and a south shadow region produced by the block (M, N; 5). Let $E_4$ be an intersection of A with the second east 3-element column on the right side of the block (m, n; 4) and let $S_5$ be an intersection of A with the second south 4-element row in the bottom of the block (M, N; 5). The combined Sylvester algorithm fails only if both $E_4$ and $S_5$ (in general: $E_K$ and $S_K$) are nonempty.

In this particular case it is necessary to compute $E_4 = C_{n+5}^{m+3}$ by the Froissart-Gilewicz or Paszkowski identity, to be able to follow the recurrence computation of the c-table by the combined Sylvester algorithms. If not, we can compute the south black elements $S_5$ at the bottom of the blocks of zeros and follow the computation by Sylvester algorithms. The last strategy is recommended if $card(S_K) < card(E_k)$, that is, in the opposite situation, as seen in Figure 5.

However this “theoretical” strategy based on an extensive application of the combined Sylvester algorithm can be only used near the diagonal of a c-table. Indeed, the complexity diagram (Figure 1, Section 2) shows that outside this region the cost of the computation of elements needed is much lower

<figure><img src="image_5.png" /><figcaption>FIGURE 3: Location of determinants giving different forms of error $r$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4: East and south shadow regions in the c-table.</figcaption></figure>

if they are computed by the Froissart-Gilewicz or Paszkowski formulas.

The first identity allows going around a block (m, n; k) in the c-table (i.e., for $C_n^m$). It was proven by Froissart and Gilewicz [11]. Using the notation of Figure 2 it can be presented in the following form of a determinant:

$$k \ge 3: \quad (-1)^k \left| \begin{array}{cc} \frac{N'}{N} & \frac{1}{\tilde{W}} \\ \frac{S'}{S} & \hat{E} \end{array} \right| = (-1)^i \left| \begin{array}{cc} \frac{W'}{W} & \tilde{N} \\ \frac{E'}{E} & \frac{1}{\tilde{S}} \end{array} \right|, \quad (35)$$ $$k = 2: \quad \frac{N'}{N^2} + \frac{S'}{S^2} = \frac{W'}{W^2} + \frac{E'}{E^2}.$$

After computing elements $E'$ or/and $S'$ one can continue the computing using the Sylvester formulas. This identity written for $C(m/n)$ (now $N'$, $N$, ... denote Hankel determinants) is

$$k \ge 3: \quad (-1)^k \left| \begin{array}{cc} \frac{N'}{N} & \frac{1}{\tilde{W}} \\ \frac{S'}{S} & \tilde{E} \end{array} \right| = (-1)^i \left| \begin{array}{cc} -\frac{W'}{W} & \tilde{N} \\ \frac{E'}{E} & (-1)^k \frac{1}{\tilde{S}} \end{array} \right|, \quad (36)$$ $$k = 2: \quad \frac{N'}{N^2} + \frac{S'}{S^2} = -\frac{W'}{W^2} - \frac{E'}{E^2}.$$
<figure><img src="image_7.png" /><figcaption>FIGURE 5: Intersection of two shadow regions in the c-table.</figcaption></figure>

Paszkowski identity for $C(m/n)$, equivalent to the previous, gives also $E'$ or/and $S'$:

$$E' = (-1)^{ki} S^+ \left[ \frac{N'}{W^-} + \frac{N^+}{W^+} \left( (-1)^k \frac{S'}{S^+} + (-1)^i \frac{W'}{W^-} \right) \right]. \quad (37)$$

Paszkowski presented this identity in the elegant form of a determinant:

$$\begin{vmatrix} N' & N^+ & E^+ \\ W^- & S^+ \\ W' & W^+ & S' \end{vmatrix} = 0. \qquad (38)$$

Following (32) we have $W^- * E^+ = N^- * S^+$ and $N^+ * S^- = E^- * W^+$ which give three new Paszkowski-like identities:

$$\begin{align} \begin{vmatrix} N' & N^+ & E^+ \\ N^- & E^+ \\ W' & W^+ & S' \end{vmatrix} &= \begin{vmatrix} N' & E^- & E^+ \\ W^- & S^- & S' \\ W' & S^- & S' \end{vmatrix} \tag{39} \\ &= \begin{vmatrix} N' & N^+ & E^+ \\ N^- & E^+ \\ W' & W^+ & S' \end{vmatrix} = 0. \notag \end{align}$$

These identities for $C_n^m$ become

$$\begin{align} \begin{vmatrix} N' & N^+(-1)^k & E^+(-1)^{kn} \\ W^- & S^+ \\ -W' & W^+ & S'(-1)^{k+n} \end{vmatrix} &= 0, \tag{40} \\ \begin{vmatrix} N'(-1)^{k+1} & N^+ & E^+(-1)^n \\ N^- & E^+ \\ W'(-1)^n & W^+ & S'(-1)^{k+1} \end{vmatrix} &= \begin{vmatrix} N'(-1)^{kn} & E^- & E^+ \\ W^- & S^+ \\ W'(-1)^n & S^- & S'(-1)^{k+1} \end{vmatrix} \notag \\ &= \begin{vmatrix} N'(-1)^n & N^+ & E^+(-1)^{k+1} \\ N^- & E^+ \\ W'(-1)^{kn} & W^+ & S' \end{vmatrix} = 0. \notag \end{align}$$

Achuthan and Ponnuswamy have published [18] a detailed analysis of a table of approximants of McCabe M-fractions in the nonnormal case. To this so-called M-table corresponds a table of Toeplitz (or Hankel) determinants, like the c-table (or c-table). These tables have a square block structure analogous to that presented in this work. Consequently, all our formulas can be applied to compute the elements around the blocks in the table of Hankel determinants related to the M-table.

## 5. Computation of Padé Approximants around the Blocks: Numerical Recommendations

Our goal consists of indicating the location of the blocks in the p-table and then giving some information about the location of the BPA by means of the analysis of the structure of the c-table. The next step in the work of approximation consists of computing the relevant Padé approximant. Many methods of computation of Padé approximants use the relation between the convergents of continued fractions and Padé approximants [19]. However if we wish to compute an individual Padé approximant $[m/n]$ solving the linear system (10) followed by the substitutions (9), we can use Cholesky method, the cost of which is

$$\frac{1}{6}n(n^2 + 2) + mn. \qquad (41)$$

The cost of Cholesky method or other methods adapted to this particular symmetric matrix (see [11, page 370] or [20]) represents one-half of the cost of the solution by the Gauss method. For instance the cost of the computation of [6/4] by Cholesky is 36 but by Longman algorithm is 465. Using this last algorithm we compute all 45 Padé approximants located in the triangle ended by the antidiagonal $m+n=10$. If we need to know all p-tables it is recommended to use any efficient recursive algorithm (see [10, 11, 20]) computing the so-called L-triangle in the p-table knowing $L+1$ first coefficients of the power series of $f$. To do this it is also useful to know the generalization of the Wynn cross identity (20) obtained by Florent Cordellier (see [11, page 369]) for the Padé approximants surrounding a block $(m, n; k)$ in the p-table:

$$i = 1, 2, \dots, k: \quad (S_i - C)^{-1} + (N_i - C)^{-1} \\ = (W_i - C)^{-1}(E_i - C)^{-1}, \qquad (42)$$

where $C = [m/n]$, $S_i = [m + k/n + k - i]$, $N_i = [m - 1/n + i - 1]$, $W_i = [m + i - 1/n - 1]$, and $E_i = [m + k - i/n + k]$ occupied symmetric positions around the block as indicated in Figure 6.

## 6. Best Padé Approximant (BPA)

The Padé approximant to a function $f$ is the best local approximation of $f$ with respect to the norm of uniform convergence. However in practice we are commonly interested in the nonlocal quality of this approximation as, for instance, an analytic continuation of $f$. Knowing $L+1$ first coefficients $c_j$ of power expansion (1) of $f$ we can compute all Padé approximants $[m/n]$ with $m+n \le L$, that is, the so-called L-triangle
<figure><img src="image_8.png" /><figcaption>FIGURE 6: Location of Padé approximants around the block in the Cordellier identity.</figcaption></figure>

in the *p*-table. The entries on the last antidiagonal $m + n = L$ contain complete initial information. The natural question is: where is the BPA on each antidiagonal? The answer to this question is given in [11] where four empirical algorithms of choice of BPA, based on the numerical experiments and justified for few classes of functions, are proposed:

(i) method $\rho$: analysis of the behavior of the sequence $\{|c_n/c_{n+1}|\}$ which also allows the elimination of badly computed coefficients;

(ii) method of valleys: the minima of $|C_n^m|$ on the antidiagonals indicate the positions of BPA and/or the positions of the blocks (see Example 7);

(iii) method of coefficients of PA: estimation of coefficients which can be neglected allowing the detection of the blocks;

(iv) method of Gram matrix: analysis of the eigenvalues of $G = AA^T$, where $A$ is the matrix of the system (10); this shows if $|C_n^m| = (\prod_{i=1}^n \lambda_i)^{1/2}$ is equal to zero or not.

In numerous scientific papers where the Padé approximation method is applied, the authors automatically use the diagonal Padé approximants $[n/n]$ without any justification. This frequently leads to erroneous conclusions as those of Van Dyke (reported in [11, page 414]) who, after selecting for his problem the [2/2] PA, claimed that PA method is not satisfactory in his case. Curiously, the [2/2] PA is the worst choice among all possible elements in the *p*-table and the BPA gives a result ten times better than the better approximation proposed by Van Dyke!

We repeat that our goal is to help the authors wishing to apply correctly the Padé approximation method detecting the BPA quickly.

## 7. Numerical Observations and Computational Suggestions

In numerical practice the coefficients $c_j$ used for the computation of Padé approximants are first computed numerically

or determined experimentally. Then they are affected by some noise and some errors, in general increasing with the index. The $\rho$-method of choice of best Padé approximant (BPA) [11] allows the elimination of the last badly computed coefficients. This empirical method and the recommendation do not compute the approximants. Many papers (see [11, 21–24]) were devoted to analyzing the effect of noise on the Padé approximants. The so-called *Froissart doublets* (zero-pole doublets) appearing in the Padé approximants are created by the noise. Unfortunately one has not obtained at the present moment an efficient method of filtering the noise from the coefficients: the simple elimination of Froissart doublets spoils the quality of the resulting fraction obtained by the cleaning of the Padé approximant.

## 8. Conclusion

This text contains a complete list of useful identities and empirical numerical rules. Certain of these are new or were never published in current journals. We hope that the presented strategy of analysis of the *c*-table before the calculation of the BPA or of the *p*-table will allow the improvement of the efficient use of Padé approximation method.

## Acknowledgments

The second author wishes to thank Professor Luc Wuytack who has encouraged him many years ago to write the unpublished note [25] which is in the origin of this review. The authors would like to thank the referees for their very valuable corrections.

## References

[1] D. Bahuguna, A. Ujlayan, and D. N. Pandey, "A comparative study of numerical methods for solving an integro-differential equation," *Computers & Mathematics with Applications*, vol. 57, no. 9, pp. 1485-1493, 2009.

[2] R. Jedynak and J. Gilewicz, "Approximation of the integrals of the gaussian distribution of asperity heights in the Greenwood-Tripp contact model of two rough surfaces revisited," *Journal of Applied Mathematics*, vol. 2013, Article ID 459280, 7 pages, 2013.

[3] R. Jedynak and J. Gilewicz, "Approximation of smooth functions by weighted means of N-point Padé approximants," *Ukrainian Mathematical Journal*, vol. 65, no. 10, pp. 1410-1419, 2013.

[4] D. Roy, "Global approximation for some functions," *Computer Physics Communications*, vol. 180, no. 8, pp. 1315-1337, 2009.

[5] A. I. Aptekarev, V. I. Buslaev, A. Martines-Finkel'shtein, and S. P. Suetin, "Padé approximants, continued fractions, and orthogonal polynomials," *Russian Mathematical Surveys*, vol. 66, no. 6(402), pp. 37-122, 2011.

[6] L. N. Trefethen, *Approximation Theory and Approximation Practice*, SIAM, Philadelphia, Pa, USA, 2013.

[7] M. Kaliyappan, S. Ponnusamy, and S. Sundar, "Recursive formulation of the matrix Padé approximation in packed storage," *Computers & Mathematics with Applications*, vol. 59, no. 4, pp. 1532-1540, 2010.

[8] W. B. Gragg, "The Padé table and its relation to certain algorithms of numerical analysis," *SIAM Review*, vol. 14, pp. 1-16, 1972.
[9] J. Gilewicz and A. Magnus, “Valleys in c-table,” in *Padé Approximation and Its Applications (Proc. Conf., Univ. Antwerp, Antwerp, 1979)*, L. Wuytack, Ed., vol. 765 of *Lecture Notes in Mathematics*, pp. 135–149, Springer, Berlin, Germany, 1979.

[10] G. A. Baker, Jr., *Essentials of Padé Approximants*, Academic Press, New York, NY, USA, 1975.

[11] J. Gilewicz, *Approximants de Padé*, vol. 667 of *Lecture Notes in Mathematics*, Springer, Berlin, Germany, 1978.

[12] S. Paszkowski, “Evaluation of a C-table,” *Journal of Computational and Applied Mathematics*, vol. 44, no. 2, pp. 219–233, 1992.

[13] J. Gilewicz and E. Leopold, “Subroutine CTABLE for best Padé approximant detection,” preprint CPT-81/P.1324, 1981.

[14] J. Gilewicz and M. Pindor, “Computation of c-table and detection of valleys,” *Ukrainian Journal of Mathematics*, vol. 62, pp. 762–772, 2010.

[15] M. H. Lallemand, *Logiciel pour la détéction du Meilleur Approximant de Padé* [Ph.D. thesis], DEA, CPT-CNRS, Marseille, France, 1985.

[16] J. Vignes, “Review on stochastic approach to round-off error analysis and its applications,” *Mathematics and Computers in Simulation*, vol. 30, no. 6, pp. 481–491, 1988.

[17] W. Guziński, PADELIB: Program Library for Padé Approximation, Inst. of Nuclear Research, Warsaw, Poland, 1978, (in Polish).

[18] P. Achuthan and S. Ponnuswamy, “On general two-point continued fraction expansions and Padé tables,” *Journal of Approximation Theory*, vol. 64, no. 3, pp. 291–314, 1991.

[19] J. Gilewicz and R. Jedynak, “Compatibility of continued fraction convergents with Padé approximants,” in *Approximation and Computation. In Honor of Gradimir V. Milovanovic*, W. Gautschi, G. Mastroianni, and Th. M. Rassias, Eds., vol. 42 of Springer *Optimization and Its Application*, pp. 135–144, Springer, New York, NY, USA, 2011.

[20] C. Brezinski, “Computation of Padé approximants and continued fractions,” *Journal of Computational and Applied Mathematics*, vol. 2, no. 2, pp. 113–123, 1976.

[21] B. Beckermann, D. Bessis, L. Perotti, and D. Vrinceanu, “Computing high precision matrix Padé approximants,” *Numerical Algorithms*, vol. 61, no. 2, pp. 189–208, 2012.

[22] D. Bessis and L. Perotti, “Universal analytic properties of noise: introducing the J-matrix formalism,” *Journal of Physics A*, vol. 42, no. 36, Article ID 365202, 15 pages, 2009.

[23] J. Gilewicz and M. Pindor, “Padé approximants and noise: a case of geometric series,” *Journal of Computational and Applied Mathematics*, vol. 87, no. 2, pp. 199–214, 1997.

[24] J. Gilewicz and M. Pindor, “Padé approximants and noise: rational functions,” *Journal of Computational and Applied Mathematics*, vol. 105, no. 1-2, pp. 285–297, 1999.

[25] J. Gilewicz, “Computation of the c-table related to the Padé approximation,” preprint CPT-91/P.2901, 1991.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

Hindawi

Submit your manuscripts at
http://www.hindawi.com

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>