subroutine monomial_value ( m, n, e, x, v )

!*****************************************************************************80
!
!! monomial_value() evaluates a monomial.
!
!  Discussion:
!
!    This routine evaluates a monomial of the form
!
!      product ( 1 <= i <= m ) x(i)^e(i)
!
!    where the exponents are nonnegative integers.  Note that
!    if the combination 0^0 is encountered, it should be treated
!    as 1.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    20 April 2014
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points at which the
!    monomial is to be evaluated.
!
!    Input, integer E(M), the exponents.
!
!    Input, real ( kind = rk ) X(M,N), the point coordinates.
!
!    Output, real ( kind = rk ) V(N), the value of the monomial.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m
  integer n

  integer e(m)
  integer i
  real ( kind = rk ) v(n)
  real ( kind = rk ) x(m,n)

  v(1:n) = 1.0D+00

  do i = 1, m
    if ( 0 /= e(i) ) then
      v(1:n) = v(1:n) * x(i,1:n) ** e(i)
    end if
  end do

  return
end
subroutine r8ge_det ( n, a, pivot, det )

!*****************************************************************************80
!
!! r8ge_det() computes the determinant of a matrix factored by R8GE_FA.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    29 March 2003
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Jack Dongarra, James Bunch, Cleve Moler, Pete Stewart,
!    LINPACK User's Guide,
!    SIAM, 1979
!
!  Parameters:
!
!    Input, integer N, the order of the matrix.
!    N must be positive.
!
!    Input, real ( kind = rk ) A(N,N), the LU factors computed by R8GE_FA.
!
!    Input, integer PIVOT(N), as computed by R8GE_FA.
!
!    Output, real ( kind = rk ) DET, the determinant of the matrix.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) a(n,n)
  real ( kind = rk ) det
  integer i
  integer pivot(n)

  det = 1.0D+00

  do i = 1, n
    det = det * a(i,i)
    if ( pivot(i) /= i ) then
      det = - det
    end if
  end do

  return
end
subroutine r8ge_fa ( n, a, pivot, info )

!*****************************************************************************80
!
!! r8ge_fa() factors a general matrix.
!
!  Discussion:
!
!    R8GE_FA is a simplified version of the LINPACK routine DGEFA.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Jack Dongarra, James Bunch, Cleve Moler, Pete Stewart,
!    LINPACK User's Guide,
!    SIAM, 1979
!
!  Parameters:
!
!    Input, integer N, the order of the matrix.
!    N must be positive.
!
!    Input/output, real ( kind = rk ) A(N,N), the matrix to be factored.
!    On output, A contains an upper triangular matrix and the multipliers
!    which were used to obtain it.  The factorization can be written
!    A = L * U, where L is a product of permutation and unit lower
!    triangular matrices and U is upper triangular.
!
!    Output, integer PIVOT(N), a vector of pivot indices.
!
!    Output, integer INFO, singularity flag.
!    0, no singularity detected.
!    nonzero, the factorization failed on the INFO-th step.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) a(n,n)
  integer i
  integer info
  integer pivot(n)
  integer j
  integer k
  integer l
  real ( kind = rk ) t

  info = 0

  do k = 1, n - 1
!
!  Find L, the index of the pivot row.
!
    l = k
    do i = k+1, n
      if ( abs ( a(l,k) ) < abs ( a(i,k) ) ) then
        l = i
      end if
    end do

    pivot(k) = l
!
!  If the pivot index is zero, the algorithm has failed.
!
    if ( a(l,k) == 0.0D+00 ) then
      info = k
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'R8GE_FA - Warning!'
      write ( *, '(a,i8)' ) '  Zero pivot on step ', info
      return
    end if
!
!  Interchange rows L and K if necessary.
!
    if ( l /= k ) then
      t      = a(l,k)
      a(l,k) = a(k,k)
      a(k,k)= t
    end if
!
!  Normalize the values that lie below the pivot entry A(K,K).
!
    a(k+1:n,k) = - a(k+1:n,k) / a(k,k)
!
!  Row elimination with column indexing.
!
    do j = k+1, n

      if ( l /= k ) then
        t      = a(l,j)
        a(l,j) = a(k,j)
        a(k,j) = t
      end if

      a(k+1:n,j) = a(k+1:n,j) + a(k+1:n,k) * a(k,j)

    end do

  end do

  pivot(n) = n

  if ( a(n,n) == 0.0D+00 ) then
    info = n
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'R8GE_FA - Warning!'
    write ( *, '(a,i8)' ) '  Zero pivot on step ', info
  end if

  return
end
subroutine simplex_general_sample ( m, n, t, x )

!*****************************************************************************80
!
!! simplex_general_sample() samples a general simplex in M dimensions.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    02 March 2017
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Reuven Rubinstein,
!    Monte Carlo Optimization, Simulation, and Sensitivity 
!    of Queueing Networks,
!    Krieger, 1992,
!    ISBN: 0894647644,
!    LC: QA298.R79.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) T(M,M+1), the simplex vertices.
!
!    Output, real ( kind = rk ) X(M,N), the points.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m
  integer n

  real ( kind = rk ) t(m,m+1)
  real ( kind = rk ) x(m,n)
  real ( kind = rk ) x1(m,n)

  call simplex_unit_sample ( m, n, x1 )

  call simplex_unit_to_general ( m, n, t, x1, x )

  return
end
subroutine simplex_general_volume ( m, t, volume )

!*****************************************************************************80
!
!! simplex_general_volume() computes the volume of a simplex in N dimensions.
!
!  Discussion:
!
!    The formula is: 
!
!      volume = 1/M! * det ( B )
!
!    where B is the M by M matrix obtained by subtracting one
!    vector from all the others.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    29 March 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the dimension of the space.
!
!    Input, real ( kind = rk ) T(M,M+1), the vertices.
!
!    Output, real ( kind = rk ) VOLUME, the volume of the simplex.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m

  real ( kind = rk ) t(m,m+1)
  real ( kind = rk ) b(m,m)
  real ( kind = rk ) det
  integer i
  integer info
  integer j
  integer pivot(m)
  real ( kind = rk ) volume

  b(1:m,1:m) = t(1:m,1:m)
  do j = 1, m
    b(1:m,j) = b(1:m,j) - t(1:m,m+1)
  end do

  call r8ge_fa ( m, b, pivot, info )

  if ( info /= 0 ) then

    volume = 0.0D+00

  else

    call r8ge_det ( m, b, pivot, det )

    volume = abs ( det )
    do i = 1, m
      volume = volume / real ( i, kind = rk )
    end do

  end if

  return
end
subroutine simplex_unit_monomial_integral ( m, e, integral )

!*****************************************************************************80
!
!! simplex_unit_monomial_integral(): integral in unit simplex in M dimensions.
!
!  Discussion:
!
!    The monomial is F(X) = product ( 1 <= I <= M ) X(I)^E(I)
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    13 January 2014
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer E(3), the exponents.  
!    Each exponent must be nonnegative.
!
!    Output, real ( kind = rk ) INTEGRAL, the integral.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m

  integer e(m)
  integer i
  real ( kind = rk ) integral
  integer j
  integer k

  if ( any ( e(1:m) < 0 ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SIMPLEX_UNIT_MONOMIAL_INTEGRAL - Fatal error!'
    write ( *, '(a)' ) '  All exponents must be nonnegative.'
    stop 1
  end if

  k = 0
  integral = 1.0D+00

  do i = 1, m

    do j = 1, e(i)
      k = k + 1
      integral = integral * real ( j, kind = rk ) / real ( k, kind = rk )
    end do

  end do

  do i = 1, m
    k = k + 1
    integral = integral / real ( k, kind = rk )
  end do

  return
end
subroutine simplex_unit_sample ( m, n, x )

!*****************************************************************************80
!
!! simplex_unit_sample() samples the unit simplex in M dimensions.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    13 January 2014
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Reuven Rubinstein,
!    Monte Carlo Optimization, Simulation, and Sensitivity 
!    of Queueing Networks,
!    Krieger, 1992,
!    ISBN: 0894647644,
!    LC: QA298.R79.
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points.
!
!    Output, real ( kind = rk ) X(M,N), the points.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m
  integer n

  real ( kind = rk ) e(m+1)
  real ( kind = rk ) e_sum
  integer j
  real ( kind = rk ) x(m,n)

  do j = 1, n

    call random_number ( harvest = e(1:m+1) )

    e(1:m+1) = - log ( e(1:m+1) )

    e_sum = sum ( e(1:m+1) )

    x(1:m,j) = e(1:m) / e_sum

  end do

  return
end
subroutine simplex_unit_to_general ( m, n, t, ref, phy )

!*****************************************************************************80
!
!! simplex_unit_to_general() maps the unit simplex to a general simplex.
!
!  Discussion:
!
!    Given that the unit simplex has been mapped to a general simplex
!    with vertices T, compute the images in T, under the same linear
!    mapping, of points whose coordinates in the unit simplex are REF.
!
!    The vertices of the unit simplex are listed as suggested in the
!    following:
!
!      (0,0,0,...,0)
!      (1,0,0,...,0)
!      (0,1,0,...,0)
!      (0,0,1,...,0)
!      (...........)
!      (0,0,0,...,1)
!
!    Thanks to Andrei ("spiritualworlds") for pointing out a mistake in the
!    previous implementation of this routine, 02 March 2008.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    02 March 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points to transform.
!
!    Input, real ( kind = rk ) T(M,M+1), the vertices of the
!    general simplex.  
!
!    Input, real ( kind = rk ) REF(M,N), points in the 
!    reference triangle.
!
!    Output, real ( kind = rk ) PHY(M,N), corresponding points 
!    in the physical triangle.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m
  integer n

  integer dim
  real ( kind = rk ) phy(m,n)
  real ( kind = rk ) ref(m,n)
  real ( kind = rk ) t(m,m+1)
  integer vertex
!
!  The image of each point is initially the image of the origin.
!
!  Insofar as the pre-image differs from the origin in a given vertex
!  direction, add that proportion of the difference between the images
!  of the origin and the vertex.
!
  do dim = 1, m

    phy(dim,1:n) = t(dim,1)

    do vertex = 2, m + 1

      phy(dim,1:n) = phy(dim,1:n) &
        + ( t(dim,vertex) - t(dim,1) ) * ref(vertex-1,1:n)

    end do

  end do

  return
end
function simplex_unit_volume ( m )

!*****************************************************************************80
!
!! simplex_unit_volume() computes the volume of the unit simplex in M dimensions.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    13 January 2014
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Output, real ( kind = rk ) SIMPLEX_UNIT_VOLUME, the volume.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer i
  integer m
  real ( kind = rk ) simplex_unit_volume
  real ( kind = rk ) value

  value = 1.0D+00
  do i = 1, m
    value = value / real ( i, kind = rk )
  end do

  simplex_unit_volume = value

  return
end