function w = weights_ls ( d, a, b, n, x )

%*****************************************************************************80
%
%% weights_ls() computes weights for a least squares quadrature rule.
%
%  Discussion:
%
%    We suppose that we wish to estimate integrals of functions f(x)
%    over the interval [a,b], given values of the function at N points X,
%    and that we wish to do this in such a way that our estimate is
%    based on the least squares polynomial approximant p(x) of degree D.
%
%    The resulting quadrature rule will only be exact if N = D + 1
%    and f(x) is a polynomial of degree D or less.
%
%    This function is based in part on ideas of Mac Hyman, of Tulane
%    University.
%
%  Method:
%
%    The approximant p(x) to f(x) has coefficients C(1:D + 1).
%    The values of p(x) at the abscissas are found by V * C = P.
%    The values of f(x) at the abscissas are in the vector F.
%    We want to minimize |P-F| = |V * C - F|.
%    We set up the D + 1 x D + 1 system V' * V * C = V' * F.
%
%    Now, for quadrature, we have Q * C = W * F,
%    so formally Q * inv ( V' * V ) * V' = W.
%
%    Hence, to solve for the weights W, we solve two Vandermonde systems:
%
%    1) Solve   Q' = V' * S  for S, an N vector.
%    2) Solve   S  = V  * R  for R, a D + 1 vector.
%    3) Compute W  = R' * V' for W, an N vector.
%
%    so that, by construction (and writing inv() for matrices which
%    actually don't have inverses):
%
%      W  = R' * V'
%      W' = V  * R
%      W' = V  * inv ( V )      * S
%      W' = V  * inv ( V )      * inv ( V' ) * Q'
%      W  = Q  * inv ( V )      * inv ( V' ) * V'
%      W  = Q  * inv ( V' * V )              * V'
%
%    Thus, knowing the values Q and the abscissas X, we can construct
%    V and solve for W, and the value of W depends (as it should) only
%    on N and X, but not on f(x).
%
%    Having W, we can estimate the integral of any f(x)
%    if we have its values F(1:n) at the points X(1:n), by the formula:
%
%      I(f) \approx Q(f) = sum ( 1 <= i <= n ) W(i) * F(i)
%
%    A simple test is to use N points associated with a known quadrature
%    rule, and specify D = N-1.  Assuming the rule is exact for polynomials
%    of degree D, this procedure should determine the same weights as 
%    those used for the given quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    11 March 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer D, the desired degree of the polynomial approximant.
%    0 <= D.
%
%    Input, real A, B, the left and right endpoints of the interval.
%
%    Input, integer N, the number of data points or abscissas.
%
%    Input, real X(N), the data points.
%
%    Output, real W(N), the quadrature weights.
%

%
%  Integrate the basis functions x^0, x^1, ..., x^d over [a,b].
%
  q = zeros ( 1, d + 1 );
  bi = b;
  ai = a;
  for i = 1 : d + 1
    q(i) = ( bi - ai ) / i;
    ai = ai * a;
    bi = bi * b;
  end
%
%  Form the Vandermonde matrix V.
%
  v = zeros ( n, d + 1 );

  v(1:n,1) = 1.0;
  for j = 2 : d + 1
    v(1:n,j) = v(1:n,j-1) .* x(1:n);
  end
%
%  Solve Q' = V' * S for vector S of length N.
%
  s = v' \ q';
%
%  Solve S = V * R for vector R, of length D + 1.
%
  r = v \ s;
%
%  Compute W = R' * V';
%
  w = r' * v';
%
%  Make W a column vector.
%
  w = w';

  return
end

