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

n=0 c n z . 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.


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][2][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 problem.The application of Padé approximants in computational problems starts frequently by the computation of the auxiliary table, called -table [8].The entries of -table are the Toeplitz determinants of matrices of linear systems defining the denominators of Padé approximants.The square blocks of zeros in the -table indicate the existence of corresponding blocks in the table of Padé approximants.The so-called valleys in the -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 -table is simpler, it is recommended to begin by this computation 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 be a formal power series of a function .It can be a Maclaurin series of an analytic function at  = 0 or an asymptotic series, such as a Stieltjes series having a zero radius of convergence [10,11].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 ̸ = 0, then {/} contains only ++1 unknowns which can be determined by  +  + 1 first coefficients of () or, equivalently, by  +  + 1 first derivatives of  at  = 0: This is usually expressed by the condition which means that the power series expansion of   ()/  () matches the power series expansion of  at  = 0 at least up to the power  + .
which leads to the separated linear system for the coefficients   and   (  ≡ 0 if  < 0): The system (8) of  equations for ( + 1) unknowns   always a nontrivial solution (Frobenius) which determines the reduced Padé form.Let us note that ( 6) is equivalent to (5) only if   is reversible, that is, if   (0) ̸ = 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 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.
where we put   ≡ 0 if  < 0. Defining   0 := 1 we can build an infinite table of    's, called -table.The second column of Table 2 contains the coefficients of the power series   1 =   .The first row contains the powers of  0 ,  0  = ( 0 )  and the second row can be computed recursively by expanding  1  : The -table was first introduced by Gragg [8].Baker [10] used an alternative definition permuting rows in (11), calling the Interest in the -table is due to its direct relation to the Padé table [8,11] and to the particular ease of its computation A theory of valleys in the -table [9] leading to a numerical algorithm of choice of the best Padé approximant [11] is based on property (14).
In general a -table presents a square block structure.We denote by (, ; ) a  ×  square block of all equal entries with a west-north corner {/}.In each block are Padé approximants and satisfy but the reduced forms located under the antidiagonal of a block are not Padé approximants: because (+)+(+)+1 > ++ (see (16)). where This formula given in [10,12] is given here in terms of Toeplitz and Hankel determinants.An extension of representations of  using other determinants around the zeros block in the table will be presented in the next section.The existence of an infinite block (, ; ∞) means that the series () represents a rational function [/].If a -table presents a block (, ; ), then the [/] Padé approximant computed with  +  + 1 coefficients  0 ,  1 , . . .,  + located on the antidiagonal  +  = const in the -table is clearly the best Padé approximant (BPA) on this antidiagonal, because it is the only approximant which reproduces exactly the series () up to the power  ++−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 -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 which permits computing the  elements from four others (ascending algorithm).This method is too laborious and so it is simpler to discover the block structure of the -table by analyzing the block structure of the -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 -table.

Computation of the 𝑐-Table in the Normal Case: Numerical Recommendations
The particular form of Toeplitz (resp.Hankel) determinants permits computing them recursively using the Sylvester formula avoiding the laborious direct computation by (11) (resp.( 13)).
Starting from the two first columns one can compute the "East"    elements by the ascending algorithm [13]: Alternatively, starting from the first two rows one can compute the "South" elements by the descending algorithm: 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    from the normalized sequence (  ), that is, with  0 = 1.The cost of the normalization of each coefficient  1 ,  2 , . . . in the second column of the -table is equal to 1. Then the costs in the second row are 0, 0, 1, 2, . . .,  − 1, . . .(12).Then the cost   of  computed by (22) is   =   +   +   +   + 4, where   ,   ,   , and   are the previous costs of , , , and , respectively.In fact it is the cost of the calculation of all  2 determinants disposed in the triangle pointed on the last    if it is situated under or on the diagonal.Analogously, the global cost of the calculation of all  2 determinants with the last    by the descending algorithm is   =   +   +   +   + 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: The precision of  computed by ( 22) is, crudely estimating, proportional to the global cost   = 2  + 2  +   +   + 5, where   ,   ,   , and   are the previous global costs.We count twice   and   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  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.

Overflow, Underflow, and Detection of Blocks: Numerical Recommendations
No problems arise in a numerical computation if all coefficients   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    's decrease or increase quite rapidly.The risk of the overflow or underflow appears approximately for  +  ≥ 3, where  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    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: 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: (| + − |),  = ,  − 1, . . ., 1, 0 or  = −, − + 1, . . ., 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: new < previous * 5 * 10 −12 . (26) We recommend this efficient method which allows the detection of the zeros blocks in -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 -table for a series produced by a rational function: The numerical constants are denoted shortly by decimal exponents: for instance .9× 10 −5 is denoted by .9− 5.The first table represents the approximative limits (25) for the base = 16 : 16 −(+) ≈ 10 −(6/5)(+) .The underlined elements correspond to the beginning of the infinite block of zeros starting at element  6 5 .As shown here both criteria work: the underlined elements on the right (-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 -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  − 1 zeros starting from  +1 +1 which indicate the existence of the block (, ; ).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 -table corresponding, for instance, to the Stieltjes function (see Table 4)  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).

Nonnormal Case: Froissart-Gilewicz and Paszkowski Identities to Overcome Blocks
The proof of geometrical progression of elements surrounding the zeros blocks in the -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 -table corresponding to the block (, ; ) with  = 8 and two "shells" adjacent to the 7 × 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 -table by presented algorithms) from the corners as indicated by arrows and with respective ratios Ŝ, Ŵ, N, and Ê.Then, it suffices to know four entries    ,   +1 ,  +1  (these three are already computed before the detection of the block), and  + +1 (first nonzero element after  − 1 zeros on the  + 1 column which identifies the block) to calculate all elements surrounding the block of zeros in the following way: Two additional relations were proved in [11]: where here  = (/ + ),  = ( + / +  − ),  = ( + /), and  = ( +  − / + ).Using the above relations we can easily extend the formulas for  in (18) to other ratios of Toeplitz determinants located around the zeros block.Let us simplify the notation of the left part of the formula (19) as follows: and denote the next determinants by   ,   ,  *  , and  *  as indicated in Figure 3 corresponding to the block (, ; ).
Then, we have the following equivalent formulas: where [/2] denote the remainder of division of  by 2.
The ascending Sylvester algorithm fails for the second east column because of the division by  = 0, that is, for  + ++1 with  = 1, . . .,  − 1.The descending algorithm fails for the second south row, that is, for  ++1 + with  = 1, . . .,  − 1.This creates the "east shadow" region and the "south shadow" region where the elements can not be computed by respective algorithms.These regions are shown in Figure 4.
If no block exists in the way of descending algorithm then the east shadow region can be computed by this algorithm.Analogously the ascending algorithm can be used to compute the south shadow region.Then, the two combined algorithms allow the easy computation of many -tables with blocks without the use of the special identities which will be presented in this section.Of course, these identities can be used to compute the black east column or the black south Journal of Applied Mathematics   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  be an intersection of an east shadow region produced by a block (, ; 4) and a south shadow region produced by the block (, ; 5).Let  4 be an intersection of  with the second east 3-element column on the right side of the block (, ; 4) and let  5 be an intersection of  with the second south 4-element row in the bottom of the block (, ; 5).The combined Sylvester algorithm fails only if both  4 and  5 (in general:   and   ) are nonempty.
In this particular case it is necessary to compute  4 =  +3 +5 by the Froissart-Gilewicz or Paszkowski identity, to be able to follow the recurrence computation of the -table by the combined Sylvester algorithms.If not, we can compute the south black elements  5 at the bottom of the blocks of zeros and follow the computation by Sylvester algorithms.The last strategy is recommended if card(  ) < card(  ), 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 -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   The first identity allows going around a block (, ; ) in the -table (i.e., for    ).It was proven by Froissart and Gilewicz [11].Using the notation of Figure 2   Paszkowski identity for (/), equivalent to the previous, gives also   or/and   : Paszkowski presented this identity in the elegant form of a determinant:  [19].However if we wish to compute an individual Padé approximant [/] solving the linear system (10) followed by the substitutions (9), we can use Cholesky method, the cost of which is 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  +  = 10.
If we need to know all -tables it is recommended to use any efficient recursive algorithm (see [10,11,20]) computing the so-called -triangle in the -

Best Padé Approximant (BPA)
The Padé approximant to a function  is the best local approximation of  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 .Knowing  + 1 first coefficients   of power expansion (1) of  we can compute all Padé approximants [/] with  +  ≤ , that is, the so-called -triangle in the -table.The entries on the last antidiagonal  +  =  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: (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  =   , where  is the matrix of the system (10); this shows if |   | = (Π  =1   ) 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 [/] 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 -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.

Numerical Observations and Computational Suggestions
In numerical practice the coefficients   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 -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][22][23][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.

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 -table before the calculation of the BPA or of the -table will allow the improvement of the efficient use of Padé approximation method.

Figure 1 :
Figure 1: Separation line optimizing the cost of each algorithm.

Figure 3 :
Figure 3: Location of determinants giving different forms of error .

Figure 4 :
Figure 4: East and south shadow regions in the -table.
if they are computed by the Froissart-Gilewicz or Paszkowski formulas.

4 Figure 5 :
Figure 5: Intersection of two shadow regions in the -table.

Figure 6 :
Figure 6: Location of Padé approximants around the block in the Cordellier identity.
(i) method : analysis of the behavior of the sequence {|  / +1 |} which also allows the elimination of badly computed coefficients;(ii) method of valleys: the minima of |   | on the antidiagonals indicate the positions of BPA and/or the positions of the blocks (seeExample 7); in general a rational function.The Padé table (table) of the formal power series (1) is a doubly infinite array of irreducible rational functions {/} called reduced Padé forms, determined in such a manner that the Maclaurin expansion of {/} agrees with () as far as possible.

Table 2
. The -table contains information about the block structure of the Padé table.Each square block of zeros in the -table defines a corresponding block in the Padé table.

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

Table 4 (
Example 7 presents a fragment of the -table for a series produced by the mentioned function.Example 7 (fragment of the -table of the Stieltjes function  − (1/) log(1 − )).The minima on the antidiagonals are underlined.They correspond to the positions of the best padé approximants on the corresponding antidiagonals in the Padé table.Joining these minima we obtain the valley structure in the -table mentioned in Section 6 (see Table it can be presented in the following form of a determinant:

Computation of Padé Approximants around the Blocks: Numerical Recommendations Our
+  +  −  +    +                − *  + =  − *  + and  + *  − =  − *  + which give three new Paszkowski-like identities:   (−1)   −   (−1) +1              =                (−1)   +  + (−1) +1  −  +   (−1)   +                Achuthan and Ponnuswamy have published [18] a detailed analysis of a table of approximants of McCabe -fractions in the nonnormal case.To this so-called -table corresponds a table of Toeplitz (or Hankel) determinants, like the -table (or -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 -table.goal consists of indicating the location of the blocks in the -table and then giving some information about the location of the BPA by means of the analysis of the structure of the -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 5.