APPENDIX

A FORTRAN FUNCTION

Specification. PADE(C,IC,A,IA,B,IB,L,M,X,W,Wl,W2,IK,rW) is a FORTRAN IV FUNCTION for calculating Pade approximants.

Method. The method is based on the solution of the linear equations (1.1.6, 7). A successful FUNCTION call causes the numerator and denominator coefficients of the Pade approximant to be calculated and returns the value of the approximant at some prespecified point X as the value of PADE. If the FUNCTION call is not successful, the approximant has failed a simple empirical degeneracy test which is built into the Gauss-Jordan elimination. The method is designed for data coefficients with full single precision accuracy. Higher precision data coefficients are often warranted: in this case, the code must be modified correspondingly.

Parameters.

C a REAL ARRAY of dimension IC. On entry, C contains the given series coefficients and is unchanged on exit.

A a REAL ARRAY of dimension IA. On exit, A(1),A(2),... ,A(L + 1) contain the numerator parameters a0,au...,aL respectively.

B a REAL ARRAY of dimension IB. On exit, B(1),B(2),...,B(M+ 1) contain the numerator parameters b0,bl,...,bM respectively.

L an INTEGER specifying the degree of the numerator.

M an INTEGER specifying the degree of the denominator.

W a two-dimensional DOUBLE PRECISION REAL ARRAY of dimensions (IW,IW) which must be provided for workspace.

W1,W2 two DOUBLE PRECISION REAL ARRAYs of dimension IW which must be provided for workspace.

IK an INTEGER ARRAY of dimension IW which must be provided for work­space.

IC,IA,IB,IW INTEGERS specifying various array dimensions, to be set on entry and un­changed on exit.

X X is the REAL variable of the given power series which must be set with some value on entry. X is unchanged on exit.

Error indicators. If the integers IC, IA, IB, and IW are inconsistently chosen (we require that L»0, M>0, IW>M, IB>M, IA>L and IC>L + M), or the approximant is degenerate, an error message is printed, and control is returned to the calling program. If the approximant

ENCYCLOPEDIA OF MATHEMATICS and Its Applications, Gian-Carlo Rota (ed.). Vol. 13: George A. Baker, Jr., and Peter R. Graves-Morris, Pade Approximants: Basic Theory, Parti       ISBN 0-201-13512-4

Copyright © 1981 by Addison-Wesley Publishing Company, Inc., Advanced Book Program. All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior permission of the publisher.

is evaluated at a pole, an error message is printed. The arrays A, B are set correctly, and the artificial value PADE( )=0 is returned.

5. Example. The following program calculates the [4/4] Pade approximant of exp(x) to estimate e. (> = 2.71828182....) Channel 7 is used for input and channel 2 is used for output.

PROGRAM BOOKPROG

DIMENSION CC(10),AA(5),BB(5),WW(10,10),WX(10),WY(10),IK(10) DOUBLE PRECISION WW,WX,WY CC(1) = 1.0 DO 1 1 = 1,8

CC(I + 1) = CC(I)/I

FORMAT(1 X.33HTHE GIVEN SERIES COEFFICIENTS ARE/9F13.8) WRITE(2,2) (CC(I),I = 1,9)

FORMAT (2I3.F13.5) READ (7,3) L,M,X LP1 = L + 1 MP1=M + 1

FORMAT(1X,//11H FORM THE [,14,1 H/,I4,13H] PADE ATX =,F13.8//) WRITE (2,4) L,M,X

Y = PADE(CC,10,AA,5,BB,5,L,M,X,WW,WX,WY,IK,10)

FORMAT(1 X.30HTHE NUMERATOR COEFFICIENTS ARE//9F13.8) WRITE(2,5) (AA(I),I = 1,LP1)

FORMAT(1 X/33H THE DENOMINATOR COEFFICIENTS ARE//9F13.8) WRITE(2,6) (BB(I),I = 1,MP1)

FORMAT(1 X/42H THE VALUE OF THE PADE APPROXIMANT AT X IS.F13.8) WRITE(2,7) Y

STOP END

FUNCTION PADE (C,IC,A,IA,B,IB,L,M,X,W,W1 ,W2,IK,IW)

C

C THE REAL FUNCTION PADE RETURNS THE VALUE OF THE [L/M] PADE C APPROXIMANT EVALUATED AT X. THE PADE NUMERATOR COEFFICIENTS C ARE STORED IN THE FIRST L + 1 LOCATIONS OF THE ARRAY A. C THE DENOMINATOR COEFFICIENTS ARE STORED IN THE FIRST M + 1 C LOCATIONS OF THE ARRAY B.

C DOUBLE PRECISION ARRAYS W(IW,IW), W1 (IW), W2(IW) AND AN C INTEGER ARRAY IK(IW) ARE REQUIRED FOR WORKSPACE. THE STRICT C INEQUALITIES IW>0 AND IW .GE. M, IB>M, IA>L, IC>L + M ARE C NECESSARY TO FULFIL THE STORAGE REQUIREMENTS AT ENTRY. C

DIMENSION C(IC),B(IB),A(IA),W(IW,IW),W1(IW),W2(IW),IK(IW) DOUBLE PRECISION W,W1,W2 DOUBLE PRECISION ONE,OUGHT,DET,T DATA ONE,OUGHT/0.1 D 01,0.0D 00/ C NEXT STATEMENT APPLIES TO COMPUTERS WITH 12 FIGURE PRECISION EPS=0.1E-10 LP1 = L + 1 B(1) = 1.0 IF(L) 40,1,1

IF(M) 40,6,2

IF(IW-M) 40,3,3

3

4

5

6

7

8

9

10

11

12

13

C

C

14

15

16

17

18

19

20

21

22

23

IF(IB-M) 40,40,4 IF(IA-L) 40,40,5 IF(IC-L-M) 40,40,8 DO 7 1 = 1,LP1 A(I) = C(I) D = 1.0 GO TO 35 DO 13 1 = 1,M DO 12 J = 1,M IF (L + l-J) 11,10,10 11 =L+I-J + 1 W(I,J) = C(I1) GO TO 12 W(l,J)=OUGHT CONTINUE 11 =L + I + 1 W1(I)=-C(I1) CONTINUE

SOLVE DENOMINATOR EQUATIONS USING GAUSS JORDAN ELIMINATION

WITH FULL PIVOTING AND DOUBLE PRECISION.

DET = ONE

DO 14 J = 1,M

IK(J) = 0

DO 30 1 = 1,M

T=OUGHT

DO 19 J = 1,M

IF(IK(J)-1) 15,19,15

DO 18 K = 1,M

IF(IK(K)-1) 16,18,30

I F(D ABS(T)-D ABS(W( J, K))) 17,17,18

IROW = J

ICOL = K

T = W(J,K)

CONTINUE

CONTINUE

IK(ICOL)=IK(ICOL) + 1

IF(IROW-ICOL) 20,22,20

DET=-DET

DO 21 N = 1,M

T = W(IROW,N)

W(IROW,N) = W(ICOL,N)

W(ICOL,N) = T

T=W1(IROW)

W1 (I ROW) = W1 (ICOL)

W1(ICOL)=T

W2(l) = W(ICOL,ICOL)

IM1 =1-1

IF(I.EQ.1) GO TO 25

T=DEXP(DLOG(DABS(DET))/DBLE(FLOAT(IM1))) IF(DABS(W2(l))-T*DBLE(FLOAT(l))*EPS) 23,25,25 DET = OUGHT

FORMAT(1X,5HTHE [,I4,1H/,I4,

*45H] PADE APPROXIMANT APPARENTLY IS REDUCIBLE OR/ *49H DOES NOT EXIST. IF THIS IS SO, TRY USING A LOWER/ *50H ORDER APPROXIMANT. OTHERWISE TRY HIGHER PRECISION/ *41H THROUGHOUT OR THE NAG LIBRARY ALGORITHM.). WRITE(2,24) L,M RETURN

DET = DET*W2(I) W(ICOL,ICOL)=ONE DO 26 N = 1,M

W(ICOL,N)=W(ICOL,N)/W2(l) W1(ICOL)=W1(ICOL)/W2(l) DO 29 LI = 1,M IF(LI-ICOL) 27,29,27

T = W(LI,ICOL) W(LI,ICOL) = OUGHT DO 28 N = 1,M

W(LI,N)=W(LI,N)-W(ICOL,N)*T W1(LI)=W1(LI)-W1(ICOL)*T

CONTINUE

CONTINUE

С EVALUATE DENOMINATOR DO 31 1 = 1,M

B(I + 1)=W1(I) D = B(M + 1) DO 32 1 = 1,M I1=M + 1-I

D = X*D + B(I1)

С EVALUATE NUMERATOR DO 34 I = 1 ,LP1 K = M + 1 IF(K.GT.I) K = l

= 0.0

DO 33 J = 1,K I1=l-J + 1

Y = Y + B(J)*C(I1) A(I) = Y

CONTINUE

= ABS(D)/M IF(Y.LT.EPS) GO TO 42

IF(L) 40,36,37

PADE=C(1)/D RETURN

Y = A(L+1) DO 38 1 = 1,L I1=L-I + 1

Y = Y*X + A(I1) PADE=Y/D RETURN

FORMAT(4H IC = ,I4,4H IA = ,I4,4H IB = ,I4,3H L = ,I4,3H M = ,I4,4H IW = , *I4/ 43H AND THERE IS AN ERROR AMONG THESE INTEGERS)

WRITE(2,39) IC,IA,IB,L,M,IW PADE = 0.0

RETURN

FORMAT(1 X,E16.8,40H APPARENTLY IS A POLE OF THE APPROXIMANT)

WRITE(2,41) X PADE = 0.0 RETURN END

The output from the program BOOKPROG is:

THE GIVEN SERIES COEFFICIENTS ARE

1.00000000 1.00000000 0.50000000 0.16666667 0.04166667 FORM THE [ 4/ 4] PADE AT X = 1.00000000 THE NUMERATOR COEFFICIENTS ARE

1.00000000 0.50000000 0.10714286 0.01190476 0.00059524 THE DENOMINATOR COEFFICIENTS ARE

1.00000000 -0.50000000 0.10714286 -0.01190476 0.00059524 THE VALUE OF THE PADE APPROXIMANT AT X IS 2.71828172