EXPLANATION FILE OF PROGRAM LSQPLY
                  ==================================


  Least-Squares Fitting With Orthogonal Polynomials
  -------------------------------------------------
   
    The major problem associated with applying the one-dimensional least-squares 
  polynomial algorithm given before (see file Lstsqr.txt) is round-off error.
  This error is due to the large number of mathematical operations associated
  with implementing equation (1.4.4):
   
                              t  -1 t
                        D = (X X)  X Y 

    By using Forsythe polynomials instead of simple powers of X, however, you
  can still apply the same general matrix equation while drastically reducing
  the number of calculations. (For a discussion of this technique, see Applied
  Numerical Methods, Ref. 10. It is interesting to note that the very signifi-
  cant round-off error implications of this alternative method are fairly
  obvious, but appear to have escaped forceful attention in the common litera-
  ture). 

    We will denote the Forsythe polynomial of degree j as Pj(x). The associated
  approximation equation which is equivalent to equation (1.4.1) is 

            y(X) = D + D P (x) + D P (X) + ... + D P (X)     (1.6.1) 
                    0   1 1       2 2             n n

  We will define the matrix P as


                   | 1  P1(x1)  P2(x1) ... Pn(x1) |
                   | 1                            |
               P = | 1                            |
                   | 1                            |
                   | 1                            |
                   | 1  P1(xm)  P2(xm) ... Pn(xm) |


  This has the sa me appearance as the X matrix shown in section 1.4, and we can 
  thereby write the approximation equation as 

                            Y = PD                             (1.6.2)
                             
  Again, following the analysis of that section, we will define the error func-
  tion for the fit to be 
                                      t
                      Error = (Y - PD) (Y - PD)                (1.6.3)
                       
  The corresponding solution obtained for the coefficient vector, 0, is 

                                t  -1 t
                          D = (P P)  P Y                       (1.6.4)

  The savings in computation (and thus in round-off error) is associated with
  the form of P^tP. If all the off-diagonal terms of that matrix could somehow
  be forced to be zero, the total number of calculations would be greatly re-
  duced. (P^tP)^-1 would then not require a full matrix inversion that is very
  prone to round-off error. ln fact, if that could be accomplished, the coef-
  ficients would simply be

                                m
                               Sum P (xi) yi
                               i=1  j
                          D  = --------------                  (1.6.5)  
                           j    m        2
                               Sum P (xi)
                               i=1  j

  The ability to force this condition rests in the choice of the Pj(x) poly-
  nomials. Forsythe found just a set. The recursion relation for calculating
  them is 

                     P (x) = (x - g ) P   (x) - d  P   (x)     (1.6.6)
                      j            j   j-1       j  j-2

               with  P  (x) = 0  and  P (x) = 1     
                      -1               1

                            m               2
                           Sum xi |P   (xi)|
                           i=1      j-1
                      g  = -----------------                   (1.6.7)
                       j    m            2
                           sum |P   (xi)|
                           i=1   j-1

                            m             
                           Sum xi P   (xi) P   (xi)
                           i=1     j-1      j-2
                      d  = ------------------------            (1.6.8)
                       j        m            2
                               sum |P   (xi)|
                               i=1   j-2

  Besides having the nice round-off error property, this formulation also has
  the advantage that the order of the fit can be increased by one degree without
  extensive recalculation. The calculation for the additional degree can simply
  be added on. 

    Program LSQPLY performs these calculations. The user inputs the degree of
  fit (N), the number of data points (M), and the data pairs [X(I),Y(I)]. Also,
  an optional error reduction factor (E) can be included. If E = 0, the subrou-
  tine will provide a fit to the degree N. If E > 0, an increasing sequence of
  fits will be tried until the standard deviation either is reduced by less
  than that factor or increases, or until the Nth-degree fit is reached. ln the
  case of termination before the Nth degree, the actual degree of fit is retur-
  ned in L. 

    It is apparent from the examples given that this algorithm is fairly resis-
  tant to round-off error.
   
    The precautions involved in using this subroutine are simple. The number of
  distinct data points must exceed the degree of fit (M >= N + 1). E must be
  nonnegative. This method is highly recommended for polynomial regression.

  From [BIBLI 01].
-----------------------------------------------------
End of file lsqply.txt