ccccccccccccccccccccccccc     Program 2.2     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE PFIT(N,M,X,F,A,U)
C
C Subroutine generating orthonormal polynomials U(M,N) up to
C (M-1)th order and coefficients A(M), for the least squares
C approximation of the function F(N) at X(N). Other variables
C used: G(K) for g_k, H(K) for h_k, S(K) for <u_k|u_k>.
C
      PARAMETER (NMAX = 101)
      PARAMETER (MMAX = 101)
      DIMENSION G(MMAX),H(MMAX),S(MMAX),X(N),F(N),U(M,N),A(M)
C
      IF(N.GT.NMAX) STOP 'Too many points'
      IF(M.GT.MMAX) STOP 'Order too high'
C
C Set up the zeroth order polynomial u_0
C
      DO      100  I = 1, N
        U(1,I) = 1.
  100 CONTINUE
      DO    200   I = 1, N
        TMP  = U(1,I)*U(1,I)
        S(1) = S(1) + TMP
        G(1) = G(1) + X(I)*TMP
        A(1) = A(1) + U(1,I)*F(I)
  200 CONTINUE
      G(1) = G(1)/S(1)
      H(1) = 0.0
      A(1) = A(1)/S(1)
C
C Set up the first order polynomial u_1
C
      DO      300  I = 1, N
        U(2,I) = X(I)*U(1,I) - G(1)*U(1,I)
        S(2)   = S(2) + U(2,I)**2
        G(2)   = G(2) + X(I)*U(2,I)**2
        H(2)   = H(2) + X(I)*U(2,I)*U(1,I)
        A(2)   = A(2) + U(2,I)*F(I)
  300 CONTINUE
      G(2) = G(2)/S(2)
      H(2) = H(2)/S(1)
      A(2) = A(2)/S(2)
C
      IF(M.LT.3) RETURN
C
C Higher order polynomials u_k from the recursive relation
C
      DO      700  I = 2, M-1
        DO    500  J = 1, N
          U(I+1,J) = X(J)*U(I,J) - G(I)*U(I,J) - H(I)*U(I-1,J)
          S(I+1)   = S(I+1) + U(I+1,J)**2
          G(I+1)   = G(I+1) + X(J)*U(I+1,J)**2
          H(I+1)   = H(I+1) + X(J)*U(I+1,J)*U(I,J)
          A(I+1)   = A(I+1) + U(I+1,J)*F(J)
  500   CONTINUE
        G(I+1) = G(I+1)/S(I+1)
        H(I+1) = H(I+1)/S(I)
        A(I+1) = A(I+1)/S(I+1)
  700 CONTINUE
      RETURN
      END