function c = least_squares_approximant_coef ( nd, xd, yd, m )

%*****************************************************************************80
%
%% least_squares_approximant_coef() determines the least squares polynomial.
%
%  Discussion:
%
%    The least squares polynomial approximant of order M to ND data values
%    (XD(I),YD(I)) is the unique polynomial P(X) of degree at most M-1 which
%    minimizes the square root of the sum of squares of the errors
%
%      err = sqrt ( sum ( 1 <= I <= ND ) ( P(XD(I)) - YD(I) ) ^ 2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer ND, the number of data points.
%    ND must be at least 1.
%
%    real XD(ND,1), the data points.
%
%    real YD(ND,1), the data values.
%
%    integer M, the order of the desired polynomial.
%
%  Output:
%
%    real C(M,1), the polynomial coefficients.
%    C(1,1) is the constant term and C(M,1) is the coefficient of X^(M-1).
%
  if ( nd < m ) 

    A = least_squares_approximant_matrix ( nd, xd, nd );

    c1 = A \ yd;

    [ u, s, v ] = svd ( A );

    for i = 1 : nd
      if ( s(i,i) < sqrt ( eps ) )
        s(i,i) = 0.0;
      else
        s(i,i) = 1.0 / s(i,i);
      end
    end
    c2 = v * s' * u' * yd;

    c(1:nd,1) = c1;
    c(nd+1:m,1) = 0.0;

  else

    A = least_squares_approximant_matrix ( nd, xd, m );

    c1 = A \ yd;

    [ u, s, v ] = svd ( A );

    for i = 1 : min ( m, nd )
      if ( s(i,i) < sqrt ( eps ) )
        s(i,i) = 0.0;
      else
        s(i,i) = 1.0 / s(i,i);
      end
    end
    c2 = v * s' * u' * yd;

    c = c1;

  end

  return
end
