function px = least_val ( nterms, b, c, d, x )

%*****************************************************************************80
%
%% LEAST_VAL evaluates a least squares polynomial defined by LEAST_SET.
%
%  Discussion:
%
%    The least squares polynomial is assumed to be defined as a sum
%
%      P(X) = SUM ( I = 1 to NTERMS ) D(I) * P(I-1,X)
%
%    where the orthogonal basis polynomials P(I,X) satisfy the following
%    three term recurrence:
%
%      P(-1,X) = 0
%      P(0,X) = 1
%      P(I,X) = ( X - B(I-1) ) * P(I-1,X) - C(I-1) * P(I-2,X)
%
%    Therefore, the least squares polynomial can be evaluated as follows:
%
%    If NTERMS is 1, then the value of P(X) is D(1) * P(0,X) = D(1).
%
%    Otherwise, P(X) is defined as the sum of NTERMS > 1 terms.  We can
%    reduce the number of terms by 1, because the polynomial P(NTERMS,X)
%    can be rewritten as a sum of polynomials;  Therefore, P(NTERMS,X)
%    can be eliminated from the sum, and its coefficient merged in with
%    those of other polynomials.  Repeat this process for P(NTERMS-1,X)
%    and so on until a single term remains.
%    P(NTERMS,X) of P(NTERMS-1,X)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Samuel Conte, Carl deBoor,
%    Elementary Numerical Analysis,
%    Second Edition,
%    McGraw Hill, 1972,
%    ISBN: 07-012446-4,
%    LC: QA297.C65.
%
%  Parameters:
%
%    Input, integer NTERMS, the number of terms in the least squares
%    polynomial.  NTERMS must be at least 1.  The input value of NTERMS
%    may be reduced from the value given to LEAST_SET.  This will
%    evaluate the least squares polynomial of the lower degree specified.
%
%    Input, real B(NTERMS), C(NTERMS), D(NTERMS), the information
%    computed by LEAST_SET.
%
%    Input, real X, the point at which the least squares polynomial
%    is to be evaluated.
%
%    Output, real PX, the value of the least squares 
%    polynomial at X.
%
  px = d(nterms);
  prev = 0.0;

  for i = nterms-1 : -1 : 1

    prev2 = prev;
    prev = px;

    if ( i == nterms-1 )
      px = d(i) + ( x - b(i) ) * prev;
    else
      px = d(i) + ( x - b(i) ) * prev - c(i+1) * prev2;
    end

  end

  return
end