home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 18 REXX / 18-REXX.zip / mathapps.zip / CurveFit2.cmd < prev    next >
OS/2 REXX Batch file  |  1998-08-28  |  24KB  |  606 lines
  1. /* Typing "CurveFit2" without parameters will provide a help screen.         */
  2. /*                                                                           */
  3. /* Read a file of X,Y values and obtain a least squares estimate of          */
  4. /* nth degree polynomial fit to the data.  I have not extensively tested the */
  5. /* program so I don't know how high a degree it will take or how long it will*/
  6. /* take to compute higher degree fits.  I've left the max matrix size in the */
  7. /* Gauss-Jordan routine at 50.  This gives a 49th degree fits.               */
  8. /* It is left as an exercise for the student to evaluate the regression      */
  9. /* at different values of x, compute residuals, or determine the             */
  10. /* appropriate degree to use in the fit.                                     */
  11. /*                                                                           */
  12. /* The format of the printed solution lists coefficients in descending       */
  13. /* powers of x, i.e.  d*x^3 + c*x^2 + b*x + a                                */
  14. /* The user may reverse this by appropriate moving of the comments in the    */
  15. /* section of code starting "/* Print the solution matrix */"                */
  16. /*                                                                           */
  17. /* If the format of the output table is not satisfactory, changing it is left*/
  18. /* as an exercise for the student.                                           */
  19. /*                                                                           */
  20. /* If this routine is called as an external routine a normal termination     */
  21. /* will return a value of 1.  A 0 is an error return.                        */
  22. /*                                                                           */
  23. /* I wish to thank Mike Montenegro for his criticisms and assistance!        */
  24. /*                                                                           */
  25. /* The program makes use of the REXXLIB routines.  It is left to the user to */
  26. /* obtain and install this shareware package.  It is readily availble at     */
  27. /* Hobbes and other ftp sites, as well as from Quercus.                      */
  28. /* The package is referenced in more than 4 places in this code, for example:*/
  29. /*    "if rxfuncquery('rexxlibregister') ...."                               */
  30. /*    "if dosisfile(in)<>1 then do ..."                                      */
  31. /*    "SumX.i=SumX.i+pow(x.n,i) ..."                                         */
  32. /*    "SumYX.i=SumYX.i+(y.n*pow(x.n,i)) ..."                                 */
  33. /* Such lines would need to be modified or replaced in REXXLIB is not used.  */
  34. /*                                                                           */
  35. /* This program may be copied and modified as desired provided appropriate   */
  36. /* credit is given.  Hopefully others will see fit to improve on it.         */
  37. /*                                                                           */
  38. /* Doug Rickman March 31, 1998 Marshall Space Flight Center, NASA.           */
  39. /*                                                                           */
  40. /* Programmer's notes:                                                       */
  41. /* I have done only one thing that is at all tricky, the rest is very ho-hum.*/
  42. /* I wanted the linear algebra routines to be generic subroutines, so they   */
  43. /* could be called and used in other programs.  To do this they must be able */
  44. /* to accept arrays of any name.  Since a subroutine in REXX can not return  */
  45. /* an array or even a compound variable the array must be created by the     */
  46. /* calling routine and the subroutine must work with it.  The Gauss-Jordan   */
  47. /* routine does its matrix operations in place!  In other words it writes    */
  48. /* the solution into the original, source matrices.  The two facts caused me */
  49. /* to set up the subroutines as follows -                                    */
  50. /*       The calling routine sets a variable called "elist" (for "edit list")*/
  51. /*    prior to calling the subroutine.  elist will contain the names of all  */
  52. /*    the matrices and other parameters needed by the subroutine.   The      */
  53. /*    subroutine is then called using only elist as a parameter.             */
  54. /*       In the subroutine elist is opened and the external variable names   */
  55. /*    are set to internal variable names.  The routine then uses the         */
  56. /*    "interpret" instruction to get the name of and work with the external  */
  57. /*    variable name.                                                         */
  58. /*                                                                           */
  59. /* This is decidedly a pain.  I could have also accomplished the objective   */
  60. /* using queues but I thought this might be quicker.  There is a decided     */
  61. /* performance penalty, but at this time it is not a real issue.             */
  62. /*                                                                           */
  63.  
  64. signal on Halt
  65. signal on NotReady
  66.  
  67. Numeric Digits 9
  68.  
  69. if rxfuncquery('rexxlibregister') then do         /* this will start rexxlib */
  70.     call rxfuncadd 'rexxlibregister', 'rexxlib', 'rexxlibregister'  
  71.     call rexxlibregister
  72.     end
  73.  
  74. if rxfuncquery('sysloadfuncs') then do           /* this will start rexxutil */
  75.     CALL RxFuncAdd 'SysLoadFuncs', 'RexxUtil', 'SysLoadFuncs' 
  76.     CALL SysLoadFuncs
  77.     end
  78.  
  79. arg in degree print
  80. in    =strip(in)
  81. degree=strip(degree)
  82. print =strip(print)
  83.  
  84. if in='' | in='?' | in='-?' | in='/?' then call Help
  85. if dosisfile(in)<>1 then do
  86.    say 'The input file: ' in' is not a valid file.'
  87.     exit
  88.     end /* do */
  89.  
  90. /* --------------------------------------------------------------------------*/
  91. /* --- begin MAIN                                               -------------*/
  92. /* Read the data                                                             */
  93. /* Data are assumed to be x y pairs, one per line, with or without delimiting*/
  94. /* commas.                                                                   */
  95. rc=stream(in,'c','open read')
  96. n=0
  97. do while lines(in)\=0
  98.    n=n+1
  99.    data=linein(in)
  100.    data=translate(data,' ',',')
  101.    parse var data x.n y.n
  102.    end
  103. x.0=n
  104. rc=lineout(in)
  105.  
  106. /* What degree polynomial? */
  107. if degree='' then do
  108.    say 'What degree polynomial should be fit to the data?'
  109.    pull degree
  110.    end
  111.  
  112. /* Compute the individual sums */
  113. call Sums
  114.  
  115. /* Build the two matrices */
  116. do i=1 to degree+1 /*row */
  117.    do j=1 to degree+1 /* column */
  118.       m=(j-1)+(i-1)
  119.       A.i.j=SumX.m
  120.       end j
  121.    end i
  122. NRowsA=i-1
  123. NColsA=NRowsA
  124.  
  125. do i=1 to degree+1
  126.    ii=i-1
  127.    b.i.1=SumYX.ii
  128.    end i
  129. NColsB=1
  130. NRowsB=NRowsA
  131.  
  132. /* Backup the original A array. */
  133. do i=1 to NRowsA
  134.    do j=1 to NRowsA
  135.       AA.i.j=A.i.j
  136.       end j
  137.    end i
  138.  
  139. /* now call the Gauss-Jordan routine */
  140. elist= 'A. NRowsA NColsA B. NRowsB NColsB'
  141. call GaussJordan elist
  142.  
  143. /* Critical quality check!!!!!                                               */
  144. /* Multiply original A. by A inverse to check for identity.                  */
  145. elist= 'A. NRowsA NRowsB AA. NRowsA NColsA C. NRowsC NColsC'
  146. call MatrixMultiply    elist
  147.  
  148. elist= 'C. NRowsC NColsC'
  149. rc=IndentityCheck()
  150. if rc=1 then do
  151.    /* say 'The result is close to the identity matrix.  All is well.'  */
  152.    end
  153. else do
  154.    say 
  155.    say 'Checking A x A(inverse):  Oh Oh!!'
  156.    say 
  157.    say 'The inverse of A. may not be numerically precise enough.  You should '
  158.    say 'examine the values in A. x A.inverse, in theory this should be the '
  159.    say 'identity matrix.  You may need to either increase the precision of '
  160.    say 'Numeric Digits (set at the beginning of the program) or increase the '
  161.    say 'amount of Numeric Fuzz set in the subroutine IndentityCheck:.  Of '
  162.    say 'course you might also consider putting in better data.  :-)   You can'
  163.    say 'also change the value of the variable "QualityChecks" internally in '
  164.    say 'the code to have the software dump additional information.'
  165.    say
  166.    end
  167.  
  168. /* Optional quality checks.  Set "QualityChecks=YES" to use.                 */
  169. QualityChecks='NO'
  170. if QualityChecks=YES then call QualityChecks
  171.  
  172.  
  173. /* Print the solution matrix */ /* There are two formats available */
  174. say
  175. say 'Solution: '
  176.  
  177. /* This format will print the solution with ascending exponents */
  178. /*
  179. Solution=b.1.1
  180. do i=2 to Degree+1
  181.    Solution=Solution' + 'b.i.1'*x^'i-1
  182.    end i
  183. say 'y = 'Solution
  184. */
  185.  
  186. /* This format will print the solution with descending exponents */
  187. solution=''
  188. do i= Degree+1 to 2 by -1
  189.    Solution=Solution b.i.1'*x^'i-1 '+'
  190.    end i
  191. Solution=solution  b.1.1
  192. say 'y = 'Solution
  193.  
  194. /* Determine the standard error */
  195. call StandardError
  196.  
  197. /* To print the table or not to print! That is the question. */
  198. select
  199.    when print='Y' then call PrintTable
  200.    when print='N' then nop
  201.    otherwise do
  202.       say
  203.       say 'Do you wish to see a table of x, y, y estimate, delta, delta^2?'
  204.       say 'Enter a "y" for yes, any other response will exit the program.'
  205.       key=translate(sysgetkey())
  206.       say
  207.       if key='Y' then call PrintTable
  208.       end
  209.    end /* select */
  210.  
  211. return 1
  212.  
  213.  
  214. /* --- end MAIN                                                 -------------*/
  215. /* --------------------------------------------------------------------------*/
  216.  
  217. /* --------------------------------------------------------------------------*/
  218. /* --- begin subroutine - Help:                                 -------------*/
  219. Help:
  220. rc= charout(,'1b'x||'[31;7m'||'CurveFit2:'||'1b'x||'[0m'||'0d0a'x)
  221. say 'Fit an nth degree polynomial to x,y data by the method of least squares.'
  222.  
  223. rc= charout(,'1b'x||'[33;1m'||'usage:'||'1b'x||'[0m')
  224. say ' CurveFit2 in degree print'
  225. say ''
  226.  
  227. rc= charout(,'1b'x||'[33;1m'||'where:'||'1b'x||'[0m')
  228. say ' in     = the file holding x,y pairs, one pair per line, comma or blank '
  229. say '                delimited.'
  230. say '       degree = degree of polynomial to fit to data.'
  231. say '       print  = if "y" a table of x, y, y estimate, delta, delta^2 is printed.'
  232. say '              = if "n" then the table will not be printed.'
  233. say ''
  234.  
  235. rc= charout(,'1b'x||'[33;1m'||'Exam: '||'1b'x||'[0m')
  236. say ' convolution Conv.data Conv.kernel Conv.out '
  237. say ' CurveFit2.cmd Curve.data 2 y'
  238.  
  239. rc= charout(,'1b'x||'[33;1m'||'notes:'||'1b'x||'[0m')
  240. say ' If the parameters "degree" and "print" are not given on the command line'
  241. say ' the user will be prompted for the polynomial degree to fit to the data'
  242. say ' and then will be asked if the table is to be printed.'
  243. say ' WARNINGS: There are relatively few checks on user input.  For example, if one'
  244. say ' attempts a 2nd degree polynomial with only 2 points, or if the input '
  245. say ' file has a blank line CurveFit will probably go Boom.  There is a check of'
  246. say ' the inverse matrix computed (an essential step).  If quality is not very'
  247. say ' high the user is notified and suggestions made to address the problem.'
  248. say ''
  249.  
  250. say 'Doug Rickman  March 21, 1998'
  251. exit
  252. return
  253.  
  254. /* --- end  subroutine - Help:                                  -------------*/
  255. /* --------------------------------------------------------------------------*/
  256.  
  257. /* --------------------------------------------------------------------------*/
  258. /* --- begin subroutine - Sums:                                 -------------*/
  259. /* Computes the sum of X, Y, X^2, Y^2, .....                                 */
  260. Sums:
  261. do i=1 to degree*2
  262.    SumX.i=0
  263.    do n=1 to x.0
  264.       SumX.i=SumX.i+pow(x.n,i) 
  265.       end n
  266.    end i
  267. SumX.0=n-1
  268.  
  269. SumY=0
  270. do n=1 to x.0
  271.    SumY=SumY+y.n
  272.    end n
  273. SumYX.0=SumY
  274.  
  275. do i=1 to degree
  276.    SumYX.i=0
  277.    do n=1 to x.0
  278.       SumYX.i=SumYX.i+(y.n*pow(x.n,i))
  279.       end n
  280.    end i
  281. return
  282. /* --- end subroutine   - Sums:                                 -------------*/
  283. /* --------------------------------------------------------------------------*/
  284.  
  285. /* --------------------------------------------------------------------------*/
  286. /* --- begin subroutine - GaussJordan:                          -------------*/
  287. /* Solve a square array using the Gauss-Jordan algorithm as outlined in      */
  288. /* numerical recipes.                                                        */
  289. /* Implemented by Doug Rickman March 13, 1998                                */
  290. /* elist holds:                                                              */
  291. /*      name of the first array,                                             */
  292. /*      the variable holding the number of rows in the first array,          */
  293. /*      the variable holding the number columns in the first array.          */
  294. /*      name of the second array,                                            */
  295. /*      the variable holding the number of rows in the second array,         */
  296. /*      the variable holding the number columns in the second array.         */
  297. /*                                                                           */
  298. /*      A. NRowsA NColsA B. NRowsB NColsB                                    */
  299. GaussJordan:
  300. procedure expose (elist)
  301.  
  302. parse var elist VArrayName1 VNRows1 VNCols1 VArrayName2 VNRows2 VNCols2 
  303.  
  304. NRows1=value(VNRows1)
  305. NCols1=value(VNCols1)
  306. ArrayName1=strip(VArrayName1,'t','.')
  307.  
  308. NRows2=value(VNRows2)
  309. NCols2=value(VNCols2)
  310. ArrayName2=strip(VArrayName2,'t','.')
  311.  
  312. N=NRows1     /* number of elements and rows */
  313. M=NCols2     /* right hand vactors is an array N by M */
  314.  
  315. NMax=50
  316. do j=1 to N
  317.    IPIV.j=0
  318.    end j
  319.  
  320. do i=1 to N
  321.    BIG=0
  322.    do j=1 to N
  323.       if IPIV.j <> 1 then do k=1 to N
  324.          interpret 'Temp='Arrayname1'.j.k'
  325.          if IPIV.k=0 & abs(Temp)>=BIG then do
  326.             interpret 'BIG=abs('Arrayname1'.j.k)'
  327.             irow=j
  328.             icol=k
  329.             end 
  330.          else if IPIV.k > 1 then do
  331.             say 'Singular matrix! Stop 1'
  332.             exit
  333.             end
  334.          end k
  335.       end j
  336.  
  337.    IPIV.icol=IPIV.icol+1
  338.  
  339.    if irow<>icol then do
  340.       do /*14*/ L=1 to N
  341.          interpret 'DUM='Arrayname1'.irow.L'
  342.          interpret Arrayname1'.irow.L='Arrayname1'.icol.L'
  343.          interpret Arrayname1'.icol.L=DUM'
  344.          end L
  345.       do L=1 to M
  346.          interpret 'DUM='Arrayname2'.irow.L'
  347.          interpret Arrayname2'.irow.L='Arrayname2'.icol.L'
  348.          interpret Arrayname2'.icol.L=DUM'
  349.          end L
  350.       end /* if irow<>icol then do ... */
  351.  
  352.    INDXR.i=irow
  353.    INDXC.i=icol
  354.    interpret 'Temp='Arrayname1'.icol.icol'
  355.    if Temp=0 then do
  356.       say 'Singular matrix! Stop 2.'
  357.       exit
  358.       end
  359.    
  360.    interpret 'PIVINV=1/'Arrayname1'.icol.icol'
  361.    interpret Arrayname1'.icol.icol=1'
  362.    do L=1 to N
  363.       interpret Arrayname1'.icol.L='Arrayname1'.icol.L*PIVINV'
  364.       end L
  365.    do L=1 to M
  366.       interpret Arrayname2'.icol.L='Arrayname2'.icol.L*PIVINV'
  367.       end L
  368.    do LL=1 to N
  369.       if LL \= icol then do
  370.          interpret 'DUM='Arrayname1'.LL.icol'
  371.          interpret Arrayname1'.LL.icol=0'
  372.          do L=1 to N
  373.             interpret Arrayname1'.LL.L='Arrayname1'.LL.L-'Arrayname1'.icol.L*DUM'
  374.             end L
  375.          do L=1 to M
  376.             interpret Arrayname2'.LL.L='Arrayname2'.LL.L-'Arrayname2'.icol.L*DUM'
  377.             end L
  378.          end /* if LL \= icol then do ... */      
  379.       end LL
  380.  
  381.    end i
  382.  
  383. /* Unscramble */
  384. do L=N to 1 by -1
  385.    if INDXR.L \= INDXC.L then do K=1 to N
  386.       INDXRL=INDXR.L
  387.       INDXCL=INDXC.L
  388.       interpret 'DUM='Arrayname1'.K.INDXRL'
  389.       interpret Arrayname1'.K.INDXRL='Arrayname1'.K.INDXCL'
  390.       interpret Arrayname1'.K.INDXCL=DUM'
  391.       end K
  392.    end L
  393.  
  394. return
  395. /* --- end subroutine   - GaussJordan:                          -------------*/
  396. /* --------------------------------------------------------------------------*/
  397.  
  398. /* --------------------------------------------------------------------------*/
  399. /* --- begin subroutine - MatrixMultiply:                       -------------*/
  400. /* Returns a 1 if successful, result is in C.                                */
  401. /* Returns a 2 if the matrices are not the correct size.                     */
  402. /* elist holds:                                                              */
  403. /*      name of the first array,                                             */
  404. /*      the variable holding the number of rows in the first array,          */
  405. /*      the variable holding the number columns in the first array.          */
  406. /*      name of the second array,                                            */
  407. /*      the variable holding the number of rows in the second array,         */
  408. /*      the variable holding the number columns in the second array.         */
  409. /*      name of the result array,                                            */
  410. /*      the variable holding the number of rows in the result array,         */
  411. /*      the variable holding the number columns in the result array.         */
  412. /*                                                                           */
  413. /*      A. NRowsA NColsA B. NRowsB NColsB C. NColsC NRowsC                   */
  414.  
  415. MatrixMultiply:
  416. procedure expose (elist)
  417.  
  418. parse var elist VArrayName1 VNRows1 VNCols1 ,
  419.                 VArrayName2 VNRows2 VNCols2 ,
  420.                 VArrayNameR VNRowsR VNColsR
  421.  
  422. NRows1=value(VNRows1)
  423. NCols1=value(VNCols1)
  424. ArrayName1=strip(VArrayName1,'t','.')
  425.  
  426. NRows2=value(VNRows2)
  427. NCols2=value(VNCols2)
  428. ArrayName2=strip(VArrayName2,'t','.')
  429.  
  430. ArrayNameR=strip(VArrayNameR,'t','.')
  431.  
  432. if NCols1=NRows2 then do
  433.    NRowsR=Nrows1
  434.    NColsR=NCols2
  435.    do i=1 to NRowsR
  436.       do j=1 to NColsR
  437.          interpret ArrayNameR'.i.j=0'
  438.          do k=1 to NCols1
  439.             interpret ArrayNameR'.i.j='ArraynameR'.i.j+('ArrayName1'.i.k*'ArrayName2'.k.j)'
  440.             end k
  441.          end j
  442.       end i
  443.    end /* if ... */
  444.  
  445.    interpret VNRowsR'=NRowsR'
  446.    interpret VNColsR'=NColsR'
  447.  
  448.    return 1
  449.    end
  450. else return 2
  451.  
  452. /* --- end subroutine   - MatrixMultiply:                       -------------*/
  453. /* --------------------------------------------------------------------------*/
  454.  
  455. /* --------------------------------------------------------------------------*/
  456. /* --- begin subroutine - IndentityCheck:                       -------------*/
  457. /* Check to see if matrix is an identity matrix (diagonal=1, other=0).       */
  458. /* 1 is returned for and identity matrix, a 0 otherwise.                     */
  459. /* Precision can be increased by increasing NUMERIC DIGITS for the Gauss-    */
  460. /* Jordan subroutine.  Also the tolerance for error can be increased by      */
  461. /* changing the NUMERIC FUZZ value in this routine.                          */
  462. /* elist holds:                                                              */
  463. /*          name of the array,                                               */
  464. /*          the variable holding the number of rows,                         */
  465. /*          the variable holding the number columns.                         */
  466. IndentityCheck: procedure expose (elist)
  467.  
  468. parse var elist ArrayName1 VNRows1 VNCols1
  469.  
  470. NRows=value(VNRows1)
  471. NCols=value(VNCols1)
  472. Arrayname=strip(Arrayname1,'t','.')
  473.  
  474. n=digits()
  475. numeric fuzz n-3
  476. do i=1 to NRows
  477.    do j=1 to NCols
  478.       interpret 'test=1+'Arrayname'.i.j'
  479.       if i=j & test=2 then iterate
  480.       if test=1 then iterate
  481.       else do
  482.          numeric fuzz 0
  483.          return 0
  484.          end 
  485.       end j
  486.    end i
  487. numeric fuzz 0
  488. return 1
  489. /* --- end subroutine   - IndentityCheck:                       -------------*/
  490. /* --------------------------------------------------------------------------*/
  491.  
  492. /* --------------------------------------------------------------------------*/
  493. /* --- begin subroutine - ShowMatrix:                           -------------*/
  494. /* elist holds:                                                              */
  495. /*          name of the array,                                               */
  496. /*          the variable holding the number of rows,                         */
  497. /*          the variable holding the number columns.                         */
  498.  
  499. ShowMatrix: procedure expose (elist)
  500.  
  501. parse var elist ArrayName1 VNRows1 VNCols1
  502.  
  503. NRows=value(VNRows1)
  504. NCols=value(VNCols1)
  505. Arrayname=strip(Arrayname1,'t','.')
  506. /* say 'Array: 'ArrayName */
  507. do i=1 to NRows
  508.    row=''
  509.    do j=1 to NCols
  510.       interpret 'row=row 'ArrayName'.i.j'
  511.       end j
  512.    say row
  513.    end i
  514. return
  515.  
  516. /* --- end subroutine   - ShowMatrix:                           -------------*/
  517. /* --------------------------------------------------------------------------*/
  518.  
  519. /* --------------------------------------------------------------------------*/
  520. /* --- begin subroutine -                                       -------------*/
  521. QualityChecks:
  522.  
  523. say
  524. /* Print the inverse of A. */
  525. say 'Inverse of A is '
  526. elist= 'A. NRowsA NColsA'
  527. call ShowMatrix elist
  528.  
  529. say
  530. /* Multiply A inverse times the solution vector.  The result should equal the*/
  531. /* original B. matrix.                                                       */
  532. elist= 'AA. NRowsA NRowsB B. NRowsB NColsB C. NRowsC NColsC'
  533. say 'Inverse of A times the solution B: (should equal the original matrix B.)'
  534. call MatrixMultiply    elist
  535. elist= 'C. NRowsC NColsC'
  536. call ShowMatrix elist
  537.  
  538. say
  539. /* Multiply original by inverse to check for identity.                       */
  540. elist= 'A. NRowsA NRowsB AA. NRowsA NColsA C. NRowsC NColsC'
  541. say 'Original A times the inverse of A: (should equal identity matrix)'
  542. call MatrixMultiply    elist
  543. elist= 'C. NRowsC NColsC'
  544. call ShowMatrix elist
  545.  
  546. return
  547.  
  548. /* --- end subroutine   -                                       -------------*/
  549. /* --------------------------------------------------------------------------*/
  550.  
  551. /* --------------------------------------------------------------------------*/
  552. /* --- begin subroutine - StandardError:                        -------------*/
  553. StandardError:
  554.  
  555. SE=0
  556. do n=1 to x.0
  557.    yHat=b.1.1
  558.    do i=2 to Degree+1
  559.       exponent=i-1
  560.       yHat=yHat + b.i.1*pow(x.n,exponent)
  561.       end i
  562.    SE=SE+(y.n-yHat)*(y.n-yHat)
  563.    end n
  564. SE=sqrt(SE/x.0)
  565.  
  566. say
  567. say 'Standard Error of Estimate = 'SE
  568.  
  569. return
  570. /* --- end  subroutine  - StandardError:                        -------------*/
  571. /* --------------------------------------------------------------------------*/
  572.  
  573. /* --------------------------------------------------------------------------*/
  574. /* --- begin subroutine - PrintTable:                           -------------*/
  575. PrintTable:
  576. say
  577. say '         x          y y estimate    delta    delta^2'
  578. say '__________ __________ __________ __________ __________'
  579.  
  580. do n=1 to x.0
  581.    yHat=b.1.1
  582.    do i=2 to Degree+1
  583.       exponent=i-1
  584.       yHat=yHat + b.i.1*pow(x.n,exponent)
  585.       end i
  586.    say right(x.n,10) right(y.n,10) right(yHat,10) right(format(y.n-yHat,,3,3,3),10) right(format(pow(y.n-yHat,2),,3,3,3),10)
  587.    end n
  588.  
  589. return
  590. /* --- end  subroutine  - PrintTable:                           -------------*/
  591. /* --------------------------------------------------------------------------*/
  592.  
  593. /* --------------------------------------------------------------------------*/
  594. /* --- begin subroutine - Halt:                                 -------------*/
  595. Halt:
  596. say 'This is a graceful exit from a Cntl-C'
  597. return 0
  598. /* --- end  subroutine - Halt:                                  -------------*/
  599. /* --------------------------------------------------------------------------*/
  600. /* --- begin subroutine - NotReady:                             -------------*/
  601. NotReady:
  602. say 'It would seem that you are pointing at non-existant data.  Oops.  Bye!'
  603. return 0
  604. /* --- end  subroutine - NotReady:                              -------------*/
  605. /* --------------------------------------------------------------------------*/
  606.