> 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.\n`):
> fi:
> if LM = 0 and LN = 0 then
> OK := FALSE:
> print(`Not both m and n can be zero\n`):
> 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:\\name.ext`):
> print(`for example:   A:\\DATA.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:
> 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,`\n\n`):
> # 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: \n`):
> 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,`\n`):
> od:
> fprintf(OUP,`\n`):
> # 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: \n`):
> 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,`\n`):
> od:
> fprintf(OUP,`\n`):
> 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:\n`):
> print(`1. Output to screen\n`):
> print(`2. Output to text file\n`):
> print(`Enter 1 or 2\n`):
> FLAG := scanf(`%d`)[1]:print(`Your input is `):print(FLAG):
> if FLAG = 2 then
> print(`Input the file name in the form - drive:\\name.ext\n`):
> print(`for example:   A:\\OUTPUT.DTA\n`):
> NAME := scanf(`%s`)[1]:print(`Output file is `):print(NAME):
> OUP := fopen(NAME,WRITE,TEXT):
> else
> OUP := default:
> fi:
> fprintf(OUP, `PADE' RATIONAL APPROXIMATION\n\n`):
> fprintf(OUP, `Denominator Coefficients Q[0], ..., Q[M] \n`):
> for I1 from 0 to LM do
> fprintf(OUP, ` %11.8f`, Q[I1]):
> od:
> fprintf(OUP, `\n`):
> fprintf(OUP, `Numerator Coefficients P[0], ..., P[N]\n`):
> for I1 from 0 to LN do
> fprintf(OUP, ` %11.8f`, P[I1]):
> od:
> fprintf(OUP, `\n`):
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi:
> fi:
> if OK = FALSE then
> print(`System has no unique solution\n`):
> fi:
> fi: