> restart:
> # CHEBYSHEV RATIONAL APPROXIMATION ALGORITHM 8.2
> #
> # To obtain the rational approximation
> #
> # rT(x) = (p0*T0 + p1*T1 +...+ pn*Tn) / (q0*T0 + q1*T1 +...+ qm*Tm)
> #
> # for a given function f(x):
> #
> # INPUT  nonnegative integers m and n.
> #
> # OUTPUT  coefficients q0, q1, ... , qm, p0, p1, ... , pn.
> #
> # The coefficients of the Chebyshev expansion a0, a1, ..., aN could
> # be calculated instead of input as is assumed in this program.
> print(`This is Chebyshev Approximation.\n\n`):
> OK := FALSE:
> while OK = FALSE do
> print(`Input m and n separated by a space\n`):
> 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:
> OK := FALSE:
> while OK = FALSE do
> print(`The Chebyshev coefficients a[0], a[1], ... , a[N+m]\n`):
> print(`are to be input.\n`):
> print(`Choice of input method:\n`):
> print(`1. Input entry by entry from keyboard\n`):
> print(`2. Input data from a text file\n`):
> print(`Choose 1 or 2 please\n`):
> FLAG := scanf(`%d`)[1]:print(`Input is `):print(FLAG):
> if FLAG = 1 or FLAG = 2 then
> OK := TRUE:
> fi:
> od:
> if FLAG = 1 then
> print(`Input in order a[0] to a[N+m]\n`):
> for I2 from 0 to BN+LM do
> print(`Input A[I] for I = `):printS(I2):
> AA[I2] := scanf(`%f`)[1]:print(`a[I] = `):print(AA[I2]):
> od:
> fi:
> if FLAG = 2 then
> print(`The text file may contain as many entries\n`):
> print(`per line as desired each separated by blank.\n`):
> print(`Has such a text file been created?\n`):
> print(`Enter 1 for yes or 2 for no`):
> ANS := scanf(`%d`)[1]:PRINT(`Input is `):print(ANS):
> if ANS = 1 then
> print(`Input the file name in the form - `):
> print(`drive:\\name.ext\n`):
> print(`for example:   A:\\DATA.DTA\n`):
> NAME := scanf(`%s`)[1]:
> INP := fopen(NAME,READ,TEXT):
> for I2 from 0 to BN+LM do
> AA[I2] := fscanf(INP, `%f`)[1]:
> od:fclose(INP):
> else
> print(`The program will end so the input file can `):
> print(`be created.\n`):
> OK := FALSE:
> fi:
> fi:
> if OK = TRUE then
> OUP := default:
> fprintf(OUP,`Coefficients a[I] for I = 0, 1, 2, ... are`):
> for I2 from 0 to BN do
> fprintf(OUP, ` %11.8f`, AA[I2]):
> od:
> fprintf(OUP,`\n\n`):
> # Step 1
> N := BN:
> M := N+1:
> # Step 2 - performed during input
> for I2 from 0 to N do
> NROW[I2] := I2:
> od:
> # Initialize row pointer
> NN := N-1:
> # Step 3
> Q[0] := 1.0:
> # Step 4
> # Set up a linear system with A used in place of B
> for I2 from 0 to N do
> # Step 5
> for J from 0 to I2 do
> if J <= LN then
> A[I2,J] := 0:
> fi:
> od:
> # Step 6
> if I2 <= LN then
> A[I2,I2] := 1.0:
> fi:
> # Step 7
> for J from I2+1 to LN do
> A[I2,J] := 0:
> od:
> # Step 8
> for J from LN+1 to N do
> if I2 <> 0 then
> PP := I2-J+LN:
> if PP < 0 then
> PP := -PP:
> fi:
> A[I2,J] := -(AA[I2+J-LN]+AA[PP])/2.0:
> else
> A[I2,J] := -AA[J-LN]/2.0:
> fi:
> od:
> A[I2,N+1] := AA[I2]:
> od:
> # Step 9
> A[0,N+1] := A[0,N+1]/2.0:
> fprintf(OUP,`The Linear System is: \n`):
> for I1 from 1 to N+1 do
> for J1 from 1 to N+2 do
> fprintf(OUP, ` %11.8f`, A[I1-1,J1-1]):
> od:
> fprintf(OUP,`\n`):
> od:
> fprintf(OUP,`\n`):
> # Steps 10 - 21 solve the linear system using partial pivoting
> I2 := LN+1:
> # Step 10
> while OK = TRUE and I2 <= NN do
> # Step 11
> IMAX := NROW[I2]:
> AMAX := abs(A[IMAX,I2]):
> IMAX := I2:
> JJ := I2+1:
> for IP from JJ to N do
> JP := NROW[IP]:
> if abs(A[JP,I2]) > AMAX then
> AMAX := abs(A[JP,I2]):
> IMAX := IP:
> fi:
> od:
> # Step 12
> if AMAX <= 1.0e-20 then
> OK := false:
> else
> # Step 13
> # Simulate row interchange
> if NROW[I2] <> NROW[IMAX] then
> NCOPY := NROW[I2]:
> NROW[I2] := NROW[IMAX]:
> NROW[IMAX] := NCOPY:
> fi:
> I1 := NROW[I2]:
> # Step 14
> # Perform elimination
> for J from JJ to M-1 do
> J1 := NROW[J]:
> # Step 15
> XM := A[J1,I2]/A[I1,I2]:
> # Step 16
> for K from JJ to M do
> A[J1,K] := A[J1,K]-XM*A[I1,K]:
> od:
> # Step 17
> A[J1,I2] := 0:
> od:
> fi:
> I2 := I2+1:
> od:
> if OK = TRUE then
> # Step 18
> N1 := NROW[N]:
> if abs(A[N1,N]) <= 1.0e-20 then
> OK := false:
> # System has no unique solution
> else
> # Step 19
> fprintf(OUP,`The Reduced Linear System is: \n`):
> for I1 from 1 to N+1 do
> for J1 from 1 to N+2 do
> fprintf(OUP, ` %11.8f`, A[I1-1,J1-1]):
> od:
> fprintf(OUP,`\n`):
> od:
> fprintf(OUP,`\n`):
> # Start backward substitution
> if LM > 0 then
> Q[LM] := A[N1,M]/A[N1,N]:
> A[N1,M] := Q[LM]:
> fi:
> PP := 1:
> # Step 20
> for K from LN+1 to NN do
> I2 := NN-K+LN+1:
> JJ := I2+1:
> N2 := NROW[I2]:
> SUM := A[N2,N+1]:
> for KK from JJ to N do
> LL := NROW[KK]:
> SUM := SUM - A[N2,KK] * A[LL,M]:
> od:
> A[N2,M] := SUM / A[N2,I2]:
> Q[LM-PP] := A[N2,M]:
> PP := PP+1:
> od:
> # Step 21
> for K from 0 to LN do
> I2 := LN-K:
> N2 := NROW[I2]:
> SUM := A[N2,N+1]:
> for KK from LN+1 to N do
> LL := NROW[KK]:
> SUM := SUM-A[N2,KK]*A[LL,M]:
> od:
> A[N2,M] := SUM:
> P[LN-K] := A[N2,M]:
> od:
> # Step 22
> # Procedure completed
> 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(`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, `CHEBYSHEV RATIONAL APPROXIMATION\n\n`):
> fprintf(OUP, `Denominator Coefficients Q[0], ..., Q[M] \n`):
> for I2 from 0 to LM do
> fprintf(OUP, ` %11.8f`, Q[I2]):
> od:
> fprintf(OUP, `\n`):
> fprintf(OUP, `Numerator Coefficients P[0], ..., P[N]\n`):
> for I2 from 0 to LN do
> fprintf(OUP, ` %11.8f`, P[I2]):
> 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: