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.\134n\134n`):

OK := FALSE:

while OK = FALSE do

print(`Input m and n separated by a space\134n`):

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:

OK := FALSE:

while OK = FALSE do

print(`The Chebyshev coefficients a[0], a[1], ... , a[N+m]\134n`):

print(`are to be input.\134n`):

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

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

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

print(`Choose 1 or 2 please\134n`):

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]\134n`):

for I2 from 0 to BN+LM do

print(`Input A[I] for I = `):print(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\134n`):

print(`per line as desired each separated by blank.\134n`):

print(`Has such a text file been created?\134n`):

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:\134\134name.ext\134n`):

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

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.\134n`):

OK := FALSE:

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,`\134n\134n`):

# 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: \134n`):

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,`\134n`):

od:

fprintf(OUP,`\134n`):

# 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: \134n`):

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,`\134n`):

od:

fprintf(OUP,`\134n`):

# 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:\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(`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, `CHEBYSHEV RATIONAL APPROXIMATION\134n\134n`):

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

for I2 from 0 to LM do

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

od:

fprintf(OUP, `\134n`):

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

for I2 from 0 to LN do

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

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:

SUNUaGlzfmlzfkNoZWJ5c2hldn5BcHByb3hpbWF0aW9uLnwrfCtHNiI=

SURJbnB1dH5tfmFuZH5ufnNlcGFyYXRlZH5ieX5hfnNwYWNlfCtHNiI=

SSVtfj1+RzYi

IiIj

SSVufj1+RzYi

IiIk

SVVUaGV+Q2hlYnlzaGV2fmNvZWZmaWNpZW50c35hWzBdLH5hWzFdLH4uLi5+LH5hW04rbV18K0c2Ig==

STJhcmV+dG9+YmV+aW5wdXQufCtHNiI=

STlDaG9pY2V+b2Z+aW5wdXR+bWV0aG9kOnwrRzYi

SUcxLn5JbnB1dH5lbnRyeX5ieX5lbnRyeX5mcm9tfmtleWJvYXJkfCtHNiI=

SUAyLn5JbnB1dH5kYXRhfmZyb21+YX50ZXh0fmZpbGV8K0c2Ig==

STZDaG9vc2V+MX5vcn4yfnBsZWFzZXwrRzYi

SSpJbnB1dH5pc35HNiI=

IiIi

ST9JbnB1dH5pbn5vcmRlcn5hWzBdfnRvfmFbTittXXwrRzYi