restart:

# PADE' RATIONAL APPROXIMATION ALGORITHM 8.1

# To obtain the rational approximation

# r(x) = p(x) / q(x)

# = (p0 + p1*x + ... + pn*x^n) / (q0 + q1*x + ... + qm*x^m)

# for a given function f(x):

# INPUT nonnegative integers m and n.

# OUTPUT coefficients q0, q1, ... , qm, p0, p1, ... , pn.

# The coefficients of the Maclaurin polynomial a0, a1, ... , aN could

# be calculated instead of input as is assumed in this program.

print(`This is Pade Approximation.`):

OK := FALSE:

while OK = FALSE do

print(`Input m and n separated by blanks`):

LM := scanf(`%d`)[1]:

LN := scanf(`%d`)[1]:print(`m = `):print(LM):print(`n = `):print(LN):

BN := LM+LN:

if LM >= 0 and LN >= 0 then

OK := TRUE:

else

print(`m and n must both be nonnegative.\134n`):

fi:

if LM = 0 and LN = 0 then

OK := FALSE:

print(`Not both m and n can be zero\134n`):

fi:

od:

print(`The MacLaurin coefficients a[0], a[1], ... , a[N]`):

print(`are to be input.`):

print(`Choice of input method:`):

print(`1. Input entry by entry from keyboard`):

print(`2. Input data from a text file`):

print(`Choose 1 or 2 please`):

FLAG := scanf(`%d`)[1]:print(`Your input is `):print(FLAG):

if FLAG = 1 then

print(`Input in order a[0] to a[N] `):

for I1 from 0 to BN do

print(`Input coefficient a[I] for I = `):print(I1):

AA[I1] := scanf(`%f`)[1]:print(`Your input is `):print(AA[I1]):

od:

OK := TRUE:

else

print(`As many entries as desired can be placed`):

print(`on each line of the file each separated by blank.`):

print(`Has the input file been created? - enter 1 for yes or 2 for no.`):

ANS := scanf(`%d`)[1]: print(`Your response is`): print(ANS):

if ANS = 1 then

print(`Input the file name in the form - `):

print(`drive:\134\134name.ext`):

print(`for example: A:\134\134DATA.DTA`):

NAME := scanf(`%s`)[1]:

INP := fopen(NAME,READ,TEXT):

for I1 from 0 to BN do

AA[I1] := fscanf(INP, `%f`)[1]:

od:

fclose(INP):

OK := TRUE:

else

print(`The program will end so the input file can `):

print(`be created.`):

OK := FALSE:

fi:

if OK = TRUE then

OUP := default:

fprintf(OUP,`Coefficients a[I] for I = 0, 1, 2, ... are`):

for I1 from 0 to BN do

fprintf(OUP, ` %11.8f`, AA[I1]):

od:

fprintf(OUP,`\134n\134n`):

# Step 1

N := BN:

M := N+1:

# Step 2 - performed during input

for I3 from 1 to N do

NROW[I3-1] := I3:

od:

# Initialize row pointer for linear system

NN := N-1:

# Step 3

Q[0] := 1:

P[0] := AA[0]:

# Step 4

# Set up a linear system but use A[i,j] in place of B[i,j]

for I3 from 1 to N do

# Step 5

for J from 1 to I3-1 do

if J <= LN then

A[I3-1,J-1] := 0:

fi:

od:

# Step 6

if I3 <= LN then

A[I3-1,I3-1] := 1:

fi:

# Step 7

for J from I3+1 to N do

A[I3-1,J-1] := 0:

od:

# Step 8

for J from 1 to I3 do

if J <= LM then

A[I3-1,LN+J-1] := -AA[I3-J]:

fi:

od:

# Step 9

for J from LN+I3+1 to N do

A[I3-1,J-1] := 0:

od:

# Step 10

A[I3-1,N] := AA[I3]:

od:

fprintf(OUP,`The Linear System is: \134n`):

for I1 from 1 to N do

for J1 from 1 to N+1 do

fprintf(OUP, ` %11.8f`, A[I1-1,J1-1]):

od:

fprintf(OUP,`\134n`):

od:

fprintf(OUP,`\134n`):

# Solve the linear system using partial pivoting.

I3 := LN+1:

# Step 11

while OK = TRUE and I3 <= NN do

# Step 12

IMAX := NROW[I3-1]:

AMAX := abs(A[IMAX-1,I3-1]):

IMAX := I3:

JJ := I3+1:

for IP from JJ to N do

JP := NROW[IP-1]:

if abs(A[JP-1,I3-1]) > AMAX then

AMAX := abs(A[JP-1,I3-1]):

IMAX := IP:

fi:

od:

# Step 13

if AMAX <= 1.0e-20 then

OK := false:

else

# Step 14

# Simulate row interchange

if NROW[I3-1] <> NROW[IMAX-1] then

NCOPY := NROW[I3-1]:

NROW[I3-1] := NROW[IMAX-1]:

NROW[IMAX-1] := NCOPY:

fi:

I1 := NROW[I3-1]:

# Step 15

# Perform elimination

for J from JJ to M-1 do

J1 := NROW[J-1]:

# Step 16

XM := A[J1-1,I3-1]/A[I1-1,I3-1]:

# Step 17

for K from JJ to M do

A[J1-1,K-1] := A[J1-1,K-1]-XM * A[I1-1,K-1]:

od:

# Step 18

A[J1-1,I3-1] := 0:

od:

fi:

I3 := I3+1:

od:

fprintf(OUP,`The Reduced Linear System is: \134n`):

for I4 from 1 to N do

for J1 from 1 to N+1 do

fprintf(OUP, ` %11.8f`, A[I4-1,J1-1]):

od:

fprintf(OUP,`\134n`):

od:

fprintf(OUP,`\134n`):

if OK = TRUE then

# Step 19

N1 := NROW[N-1]:

if abs(A[N1-1,N-1]) <= 1.0e-20 then

OK := FALSE:

# System has no unique solution

else

# Step 20

# Start backward substitution

if LM > 0 then

Q[LM] := A[N1-1,M-1]/A[N1-1,N-1]:

A[N1-1,M-1] := Q[LM]:

fi:

PP := 1:

# Step 21

for K from LN+1 to NN do

I3 := NN-K+LN+1:

JJ := I3+1:

N2 := NROW[I3-1]:

SUM := A[N2-1,N]:

for KK from JJ to N do

LL := NROW[KK-1]:

SUM := SUM-A[N2-1,KK-1]*A[LL-1,M-1]:

od:

# Step 22

A[N2-1,M-1] := SUM/A[N2-1,I3-1]:

Q[LM-PP] := A[N2-1,M-1]:

PP := PP+1:

od:

for K from 1 to LN do

I3 := LN-K+1:

N2 := NROW[I3-1]:

SUM := A[N2-1,N]:

for KK from LN+1 to N do

LL := NROW[KK-1]:

SUM := SUM-A[N2-1,KK-1]*A[LL-1,M-1]:

od:

# Step 23

# Procedure is completed.

A[N2-1,M-1] := SUM:

P[LN-K+1] := A[N2-1,M-1]:

od:

print(`Choice of output method:\134n`):

print(`1. Output to screen\134n`):

print(`2. Output to text file\134n`):

print(`Enter 1 or 2\134n`):

FLAG := scanf(`%d`)[1]:print(`Your input is `):print(FLAG):

if FLAG = 2 then

print(`Input the file name in the form - drive:\134\134name.ext\134n`):

print(`for example: A:\134\134OUTPUT.DTA\134n`):

NAME := scanf(`%s`)[1]:print(`Output file is `):print(NAME):

OUP := fopen(NAME,WRITE,TEXT):

else

OUP := default:

fi:

fprintf(OUP, `PADE' RATIONAL APPROXIMATION\134n\134n`):

fprintf(OUP, `Denominator Coefficients Q[0], ..., Q[M] \134n`):

for I1 from 0 to LM do

fprintf(OUP, ` %11.8f`, Q[I1]):

od:

fprintf(OUP, `\134n`):

fprintf(OUP, `Numerator Coefficients P[0], ..., P[N]\134n`):

for I1 from 0 to LN do

fprintf(OUP, ` %11.8f`, P[I1]):

od:

fprintf(OUP, `\134n`):

if OUP <> default then

fclose(OUP):

print(`Output file `,NAME,` created successfully`):

fi:

if OK = FALSE then

print(`System has no unique solution\134n`):

fi:

STxUaGlzfmlzflBhZGV+QXBwcm94aW1hdGlvbi5HNiI=

SUJJbnB1dH5tfmFuZH5ufnNlcGFyYXRlZH5ieX5ibGFua3NHNiI=

SSVtfj1+RzYi

IiIj

SSVufj1+RzYi

IiIk

SVJUaGV+TWFjTGF1cmlufmNvZWZmaWNpZW50c35hWzBdLH5hWzFdLH4uLi5+LH5hW05dRzYi

STFhcmV+dG9+YmV+aW5wdXQuRzYi

SThDaG9pY2V+b2Z+aW5wdXR+bWV0aG9kOkc2Ig==

SUYxLn5JbnB1dH5lbnRyeX5ieX5lbnRyeX5mcm9tfmtleWJvYXJkRzYi

ST8yLn5JbnB1dH5kYXRhfmZyb21+YX50ZXh0fmZpbGVHNiI=

STVDaG9vc2V+MX5vcn4yfnBsZWFzZUc2Ig==

SS9Zb3VyfmlucHV0fmlzfkc2Ig==

IiIi

ST1JbnB1dH5pbn5vcmRlcn5hWzBdfnRvfmFbTl1+RzYi