function [ w, x ] = gm_unit_rule_set ( rule, m, n )

%*****************************************************************************80
%
%% gm_unit_rule_set() sets a Grundmann-Moeller rule for a unit simplex.
%
%  Discussion:
%
%    This is a revised version of the calculation which seeks to compute
%    the value of the weight in a cautious way that avoids intermediate
%    overflow.  Thanks to John Peterson for pointing out the problem on
%    26 June 2008.
%
%    This rule returns weights and abscissas of a Grundmann-Moeller
%    quadrature rule for the M-dimensional unit simplex.
%
%    The dimension N can be determined by calling GM_RULE_SIZE.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 March 2017
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Axel Grundmann, Michael Moeller,
%    Invariant Integration Formulas for the N-Simplex
%    by Combinatorial Methods,
%    SIAM Journal on Numerical Analysis,
%    Volume 15, Number 2, April 1978, pages 282-290.
%
%  Input:
%
%    integer RULE, the index of the rule.
%    0 <= RULE.
%
%    integer M, the spatial dimension.
%    1 <= M.
%
%    integer N, the number of points in the rule.
%
%  Output:
%
%    real W(N), the weights.
%
%    real X(M,N), the abscissas.
%
  w = zeros ( n, 1 );
  x = zeros ( m, n );

  s = rule;
  d = 2 * s + 1;
  k = 0;
  one_pm = 1;

  for i = 0 : s

    weight =  one_pm;

    for j = 1 : max ( m, max ( d, d + m - i ) )

      if ( j <= m )
        weight = weight *  ( j );
      end
      if ( j <= d )
        weight = weight * ( d + m - 2 * i );
      end
      if ( j <= 2 * s )
        weight = weight / 2.0;
      end
      if ( j <= i )
        weight = weight / j;
      end
      if ( j <= d + m - i )
        weight = weight / j;
      end

    end

    one_pm = - one_pm;

    beta_sum = s - i;
    more = 0;
    beta = [];
    h = 0;
    t = 0;

    while ( true )

      [ beta, more, h, t ] = comp_next ( beta_sum, m + 1, ...
        beta, more, h, t );

      k = k + 1;

      w(k) = weight;

      x(1:m,k) = ( 2 * beta(2:m+1)' + 1 ) / ( d + m - 2 * i );

      if ( ~ more )
        break
      end

    end

  end
%
%  Normalize.
%
  volume1 = simplex_unit_volume ( m );
  w(1:n) = w(1:n) * volume1;

  return
end
