function csevl ( x, cs, n )

!*****************************************************************************80
!
!! CSEVL evaluates an N term Chebyshev series.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    15 April 2003
!
!  Author:
!
!    Original FORTRAN77 version by Roger Broucke.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Roger Broucke,
!    Algorithm 446:
!    Ten Subroutines for the Manipulation of Chebyshev Series,
!    Communications of the ACM,
!    Volume 16, 1973, pages 254-256.
!
!    Leslie Fox, Ian Parker,
!    Chebyshev Polynomials in Numerical Analysis,
!    Oxford Press, 1968,
!    LC: QA297.F65.
!
!    David Kahaner, Cleve Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1989,
!    ISBN: 0-13-627258-4,
!    LC: TA345.K34.
!
!  Parameters:
!
!    Input, real ( kind = rk ) X, the argument at which the series is 
!    to be evaluated.
!    X must satisfy -1.0 <= X <= 1.0.
!
!    Input, real ( kind = rk ) CS(N), the array of N terms of a Chebyshev series.
!    In evaluating CS, only half the first coefficient is summed.
!
!    Input, integer N, the number of terms in array CS.
!    N must be at least 1, and no more than 1000.
!
!    Output, real ( kind = rk ) CSEVL, the value of the Chebyshev series.
!
  implicit none

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

  integer n

  real ( kind = rk ) b0
  real ( kind = rk ) b1
  real ( kind = rk ) b2
  real ( kind = rk ) cs(n)
  real ( kind = rk ) csevl
  integer i
  real ( kind = rk ) x

  if ( n < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CSEVL - Fatal error!'
    write ( *, '(a)' ) '  Number of terms N is less than 1.'
    stop
  end if

  if ( 1000 < n ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CSEVL - Fatal error!'
    write ( *, '(a)' ) '  The number of terms is more than 1000.'
    stop
  end if

  if ( x < -1.0D+00 .or. 1.0D+00 < x ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CSEVL - Fatal error!'
    write ( *, '(a)' ) '  The input argument X is outside the interval [-1,1].'
    stop
  end if

  b1 = 0.0D+00
  b0 = 0.0D+00

  do i = n, 1, -1
    b2 = b1
    b1 = b0
    b0 = 2.0D+00 * x * b1 - b2 + cs(i)
  end do

  csevl = 0.5D+00 * ( b0 - b2 )

  return
end
function inits ( os, nos, eta )

!*****************************************************************************80
!
!! INITS estimates the order of an orthogonal series for a given accuracy.
!
!  Discussion:
!
!    Because this is a function, it is not possible to print out
!    warning error messages.  Therefore, if an error condition is
!    detected, a bogus value of INITS is returned.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    15 April 2003
!
!  Author:
!
!    Original FORTRAN77 version by Roger Broucke.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Roger Broucke,
!    Algorithm 446:
!    Ten Subroutines for the Manipulation of Chebyshev Series,
!    Communications of the ACM,
!    Volume 16, 1973, pages 254-256.
!
!    David Kahaner, Cleve Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1989,
!    ISBN: 0-13-627258-4,
!    LC: TA345.K34.
!
!  Parameters:
!
!    Input, real ( kind = rk ) OS(NOS), the coefficients in the series.
!
!    Input, integer NOS, the number of coefficients.  NOS must be
!    at least 1, or an error condition arises.
!
!    Input, real ( kind = rk ) ETA, the requested accuracy of the series.
!    Ordinarily, ETA will be chosen to be one-tenth machine precision.
!
!    Output, integer INITS, the order of the series guaranteeing 
!    the given accuracy.  However, on error, INITS will be returned
!    as a negative number.
!
  implicit none

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

  integer nos

  real ( kind = rk ) err
  real ( kind = rk ) eta
  integer i
  integer ii
  integer inits
  real ( kind = rk ) os(nos)
!
!  Fatal error.  Number of coefficients less than 1.
!  But an error message cannot be printed out because this is a function.
!
  if ( nos < 1 ) then
    inits = - huge ( i )
    return
  end if

  err = 0.0D+00

  i = 0

  do ii = 1, nos

    i = nos + 1 - ii
    err = err + abs ( os(i) )

    if ( eta < err ) then
      i = nos + 1 - ii
      exit
    end if

  end do
!
!  Eta may be too small.  Accuracy cannot be guaranteed.
!  But an error message cannot be printed out because this is a function.
!
  if ( i == 0 ) then
    i = - nos
  end if

  inits = i

  return
end
subroutine p00_c ( prob, m, c )

!*****************************************************************************80
!
!! P00_C computes a random C parameter vector for any problem.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer PROB, the problem number.
!
!    Input, integer M, the spatial dimension.
!
!    Output, real ( kind = rk ) C(M), the parameter vector.
!
  implicit none

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

  integer m

  real ( kind = rk ), dimension ( 6 ) :: b = (/ &
    1.5D+00, 0.0D+00, 1.85D+00, 7.03D+00, 20.4D+00, 4.3D+00 /)
  real ( kind = rk ) c(m)
  integer prob

  b(2) = real ( m, kind = rk )

  call r8vec_uniform_01 ( m, c )
  c(1:m) = b(prob) * c(1:m) / sum ( c )

  return
end
subroutine p00_d ( prob, m, id, c, w, n, x, d )

!*****************************************************************************80
!
!! P00_D returns a derivative component of any function.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer PROB, the index of the function.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer ID, the spatial coordinate to differentiate.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of evaluation points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) D(N), the ID-th derivative component.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) d(n)
  integer id
  integer prob
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)  

  if ( id < 0 .or. m < id ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_D - Fatal error!'
    write ( *, '(a,i4)' ) '  Illegal spatial coordinate ID = ', id
    stop
  end if

  if ( prob == 1 ) then
    call p01_d ( m, id, c, w, n, x, d )
  else if ( prob == 2 ) then
    call p02_d ( m, id, c, w, n, x, d )
  else if ( prob == 3 ) then
    call p03_d ( m, id, c, w, n, x, d )
  else if ( prob == 4 ) then
    call p04_d ( m, id, c, w, n, x, d )
  else if ( prob == 5 ) then
    call p05_d ( m, id, c, w, n, x, d )
  else if ( prob == 6 ) then
    call p06_d ( m, id, c, w, n, x, d )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_D - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal function index PROB = ', prob
    stop
  end if

  return
end
subroutine p00_f ( prob, m, c, w, n, x, f )

!*****************************************************************************80
!
!! P00_F returns the value of any function.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer PROB, the index of the function.
!
!    Input, integer M, the spatial dimension.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of evaluation points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) F(N), the function values.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) f(n)
  integer prob
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)  

  if ( prob == 1 ) then
    call p01_f ( m, c, w, n, x, f )
  else if ( prob == 2 ) then
    call p02_f ( m, c, w, n, x, f )
  else if ( prob == 3 ) then
    call p03_f ( m, c, w, n, x, f )
  else if ( prob == 4 ) then
    call p04_f ( m, c, w, n, x, f )
  else if ( prob == 5 ) then
    call p05_f ( m, c, w, n, x, f )
  else if ( prob == 6 ) then
    call p06_f ( m, c, w, n, x, f )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_F - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal function index PROB = ', prob
    stop
  end if

  return
end
subroutine p00_prob_num ( prob_num )

!*****************************************************************************80
!
!! P00_PROB_NUM returns the number of test functions available.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!   Output, integer PROB_NUM, the number of test functions.
!
  implicit none

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

  integer prob_num

  prob_num = 6

  return
end
subroutine p00_q ( prob, m, c, w, q )

!*****************************************************************************80
!
!! P00_Q returns the integral of any function.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer PROB, the index of the function.
!
!    Input, integer M, the spatial dimension.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Output, real ( kind = rk ) Q, the integral.
!
  implicit none

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

  integer m

  real ( kind = rk ) c(m)
  integer prob
  real ( kind = rk ) q
  real ( kind = rk ) w(m)

  if ( prob == 1 ) then
    call p01_q ( m, c, w, q )
  else if ( prob == 2 ) then
    call p02_q ( m, c, w, q )
  else if ( prob == 3 ) then
    call p03_q ( m, c, w, q )
  else if ( prob == 4 ) then
    call p04_q ( m, c, w, q )
  else if ( prob == 5 ) then
    call p05_q ( m, c, w, q )
  else if ( prob == 6 ) then
    call p06_q ( m, c, w, q )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_Q - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal function index PROB = ', prob
    stop
  end if

  return
end
subroutine p00_title ( prob, title )

!*****************************************************************************80
!
!! P00_TITLE returns the title for any function.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer PROB, the index of the function.
!
!    Output, character ( len = * ) TITLE, the function title.
!
  implicit none

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

  integer prob
  character ( len = * ) title

  if ( prob == 1 ) then
    call p01_title ( title )
  else if ( prob == 2 ) then
    call p02_title ( title )
  else if ( prob == 3 ) then
    call p03_title ( title )
  else if ( prob == 4 ) then
    call p04_title ( title )
  else if ( prob == 5 ) then
    call p05_title ( title )
  else if ( prob == 6 ) then
    call p06_title ( title )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_TITLE - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal function index PROB = ', prob
    stop
  end if

  return
end
subroutine p00_w ( prob, m, w )

!*****************************************************************************80
!
!! P00_W computes a random W parameter vector for any problem.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer PROB, the problem number.
!
!    Input, integer M, the spatial dimension.
!
!    Output, real ( kind = rk ) W(M), the parameter vector.
!
  implicit none

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

  integer m

  integer prob
  real ( kind = rk ) w(m)

  call r8vec_uniform_01 ( m, w )

  return
end
subroutine p01_d ( m, id, c, w, n, x, d )

!*****************************************************************************80
!
!! P01_D evaluates any derivative component for problem p01.
!
!  Discussion:
!
!    f(x) = cos ( 2 * pi * w(1) + sum ( c(1:m) * x(1:m) ) )
!
!    Default values are:
!
!    c(1:m) = 1/m
!    w(1) = 0.3
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, integer ID, the spatial coordinate to differentiate.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) D(N), the ID-th derivative component.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) d(n)
  integer i
  integer id
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  d(1:n) = 2.0D+00 * pi * w(1)
  do i = 1, m
    d(1:n) = d(1:n) + c(i) * x(i,1:n)
  end do
  d(1:n) = - c(id) * sin ( d(1:n) )

  return
end
subroutine p01_f ( m, c, w, n, x, f )

!*****************************************************************************80
!
!! P01_F evaluates the function for problem p01.
!
!  Discussion:
!
!    f(x) = cos ( 2 * pi * w(1) + sum ( c(1:m) * x(1:m) ) )
!
!    Default values are:
!
!    c(1:m) = 1/m
!    w(1) = 0.3
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) F(N), the function values.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) f(n)
  integer i
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  f(1:n) = 2.0D+00 * pi * w(1)
  do i = 1, m
    f(1:n) = f(1:n) + c(i) * x(i,1:n)
  end do
  f(1:n) = cos ( f(1:n) )

  return
end
subroutine p01_q ( m, c, w, q )

!*****************************************************************************80
!
!! P01_Q evaluates the integral for problem p01.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    29 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Output, real ( kind = rk ) Q, the integral.
!
  implicit none

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

  integer m

  real ( kind = rk ) c(m)
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  real ( kind = rk ) q
  real ( kind = rk ) w(m)

  q = 2.0D+00 ** m * &
    cos ( 2.0D+00 * pi * w(1) + 0.5D+00 * sum ( c(1:m) ) ) * &
    product ( sin ( 0.5D+00 * c(1:m) ) / c(1:m) )

  return
end
subroutine p01_title ( title )

!*****************************************************************************80
!
!! P01_TITLE returns the name of problem p01.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) TITLE, the title of the problem.
!
  implicit none

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

  character ( len = * ) title

  title = 'Oscillatory'

  return
end
subroutine p02_d ( m, id, c, w, n, x, d )

!*****************************************************************************80
!
!! P02_D evaluates an derivative component for problem p02.
!
!  Discussion:
!
!    f(x) = 1 / product ( c(1:m)^(-2) + ( x(1:m) - w(1:m) )^2 )
!
!    Default values are:
!
!    c(1:m) = 1
!    w(1:m) = 0.5
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, integer ID, the spatial coordinate to differentiate.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) D(N), the ID-th derivative component.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) d(n)
  integer i
  integer id
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  d(1:n) = 1.0D+00
  do i = 1, m
    d(1:n) = d(1:n) * ( c(i)**(-2) + ( x(i,1:n) - w(i) )**2 )
  end do
  d(1:n) = - 2.0D+00 / d(1:n) * ( x(id,1:n) - w(id) ) / &
    ( c(id) ** ( -2 ) + ( x(id,1:n) - w(id) ) ** 2 )

  return
end
subroutine p02_f ( m, c, w, n, x, f )

!*****************************************************************************80
!
!! P02_F evaluates the function for problem p02.
!
!  Discussion:
!
!    f(x) = 1 / product ( c(1:m)^(-2) + ( x(1:m) - w(1:m) )^2 )
!
!    Default values are:
!
!    c(1:m) = 1
!    w(1:m) = 0.5
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) F(N), the function values.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) f(n)
  integer i
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  f(1:n) = 1.0D+00
  do i = 1, m
    f(1:n) = f(1:n) * ( c(i)**(-2) + ( x(i,1:n) - w(i) )**2 )
  end do
  f(1:n) = 1.0D+00 / f(1:n)

  return
end
subroutine p02_q ( m, c, w, q )

!*****************************************************************************80
!
!! P02_Q evaluates the integral for problem p02.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    29 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Output, real ( kind = rk ) Q, the integral.
!
  implicit none

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

  integer m

  real ( kind = rk ) c(m)
  real ( kind = rk ) q
  real ( kind = rk ) w(m)

  q = product ( &
               (   atan ( ( 1.0D+00 - w(1:m) ) * c(1:m) ) &
                 + atan (             w(1:m)   * c(1:m) ) &
               ) * c(1:m) &
             )
  return
end
subroutine p02_title ( title )

!*****************************************************************************80
!
!! P02_TITLE returns the title of problem p02.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) TITLE, the title of the problem.
!
  implicit none

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

  character ( len = * ) title

  title = 'Product Peak'

  return
end
subroutine p03_d ( m, id, c, w, n, x, d )

!*****************************************************************************80
!
!! P03_D evaluates any derivative component for problem p03.
!
!  Discussion:
!
!    f(x) = 1 / ( 1 + sum ( c(1:m) * x(1:m) ) ) ^ ( m + 1 )
!
!    Default values are:
!
!    c(1:m) = 1/m
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, integer ID, the spatial coordinate to differentiate.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) D(N), the ID-th derivative component.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) d(n)
  integer i
  integer id
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  d(1:n) = 1.0D+00
  do i = 1, m
    d(1:n) = d(1:n) + c(i) * x(i,1:n)
  end do
  d(1:n) = - c(id) * real ( m + 1, kind = rk ) / d(1:n) ** ( m + 2 )

  return
end
subroutine p03_f ( m, c, w, n, x, f )

!*****************************************************************************80
!
!! P03_F evaluates the function for problem p03.
!
!  Discussion:
!
!    f(x) = 1 / ( 1 + sum ( c(1:m) * x(1:m) ) ) ^ ( m + 1 )
!
!    Default values are:
!
!    c(1:m) = 1/m
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) F(N), the function values.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) f(n)
  integer i
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  f(1:n) = 1.0D+00
  do i = 1, m
    f(1:n) = f(1:n) + c(i) * x(i,1:n)
  end do
  f(1:n) = 1.0D+00 / f(1:n) ** ( m + 1 )

  return
end
subroutine p03_q ( m, c, w, q )

!*****************************************************************************80
!
!! P03_Q evaluates the integral for problem p03.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    29 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Output, real ( kind = rk ) Q, the integral.
!
  implicit none

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

  integer m

  real ( kind = rk ) c(m)
  integer ivec(m)
  real ( kind = rk ) q
  integer rank
  real ( kind = rk ) r8_factorial
  real ( kind = rk ) r8_mop
  real ( kind = rk ) r8vec_i4vec_dot_product
  integer s
  real ( kind = rk ) w(m)
!
!  Here, we need to generate all possible DIM_NUM tuples with
!  values of 0 or 1.
! 
  q = 0.0D+00
  rank = 0

  do

    call tuple_next ( 0, 1, m, rank, ivec )

    if ( rank == 0 ) then
      exit
    end if

    s = sum ( ivec(1:m) )

    q = q + r8_mop ( s ) / ( 1.0D+00 + r8vec_i4vec_dot_product ( m, c, ivec ) )

  end do

  q = q / ( r8_factorial ( m ) * product ( c(1:m) ) )

  return
end
subroutine p03_title ( title )

!*****************************************************************************80
!
!! P03_TITLE returns the title of problem p03.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) TITLE, the title of the problem.
!
  implicit none

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

  character ( len = * ) title

  title = 'Corner Peak'

  return
end
subroutine p04_d ( m, id, c, w, n, x, d )

!*****************************************************************************80
!
!! P04_D evaluates any derivative component for problem p04.
!
!  Discussion:
!
!    f(x) = exp ( - sum ( c(1:m)^2 * ( x(1:m) - w(1:m) )^2 )
!
!    Default values are:
!
!    c(1:m) = 1 / m
!    w(1:m) = 0.5
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, integer ID, the spatial coordinate to differentiate.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) D(N), the ID-th derivative component.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) d(n)
  integer i
  integer id
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  d(1:n) = 0.0D+00
  do i = 1, m
    d(1:n) = d(1:n) + ( c(i) * ( x(i,1:n) - w(i) ) )**2
  end do
  d(1:n) = exp ( - d(1:n) ) * c(id)**2 * ( - 2.0D+00 ) * ( x(id,1:n) - w(id) )

  return
end
subroutine p04_f ( m, c, w, n, x, f )

!*****************************************************************************80
!
!! P04_F evaluates the function for problem p04.
!
!  Discussion:
!
!    f(x) = exp ( - sum ( c(1:m)^2 * ( x(1:m) - w(1:m) )^2 )
!
!    Default values are:
!
!    c(1:m) = 1 / m
!    w(1:m) = 0.5
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) F(N), the function values.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) f(n)
  integer i
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  f(1:n) = 0.0D+00
  do i = 1, m
    f(1:n) = f(1:n) + ( c(i) * ( x(i,1:n) - w(i) ) )**2
  end do
  f(1:n) = exp ( - f(1:n) )

  return
end
subroutine p04_q ( m, c, w, q )

!*****************************************************************************80
!
!! P04_Q evaluates the integral for problem p04.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    29 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Output, real ( kind = rk ) Q, the integral.
!
  implicit none

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

  integer m

  real ( kind = rk ) c(m)
  integer i
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  real ( kind = rk ) q
  real ( kind = rk ) r8_error
  real ( kind = rk ) w(m)

  q = 1.0D+00

  do i = 1, m

    q = q * sqrt ( pi ) &
      * ( r8_error ( c(i) * ( 1.0D+00 - w(i) ) ) &
        + r8_error ( c(i) *             w(i) ) ) &
      / ( 2.0D+00 * c(i) )

  end do

  return
end
subroutine p04_title ( title )

!*****************************************************************************80
!
!! P04_TITLE returns the title of problem p04.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) TITLE, the title of the problem.
!
  implicit none

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

  character ( len = * ) title

  title = 'Gaussian'

  return
end
subroutine p05_d ( m, id, c, w, n, x, d )

!*****************************************************************************80
!
!! P05_D evaluates any derivative component for problem p05.
!
!  Discussion:
!
!    f(x) = exp ( - sum ( c(1:m) * abs ( x(1:m) - w(1:m) ) ) )
!
!    Default values are:
!
!    c(1:m) = 2.0
!    w(1:m) = 0.5
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, integer ID, the spatial coordinate to differentiate.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) D(N), the ID-th derivative component.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) d(n)
  integer i
  integer id
  integer j
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  d(1:n) = 0.0D+00
  do i = 1, m
    d(1:n) = d(1:n) + c(i) * abs ( x(i,1:n) - w(i) )
  end do
  d(1:n) = exp ( - d(1:n) )

  do j = 1, n
    if ( x(id,j) - w(id) <= 0.0D+00 ) then
      d(j) = d(j) * c(id)
    else
      d(j) = - d(j) * c(id)
    end if
  end do

  return
end
subroutine p05_f ( m, c, w, n, x, f )

!*****************************************************************************80
!
!! P05_F evaluates the function for problem p05.
!
!  Discussion:
!
!    f(x) = exp ( - sum ( c(1:m) * abs ( x(1:m) - w(1:m) ) ) )
!
!    Default values are:
!
!    c(1:m) = 2.0
!    w(1:m) = 0.5
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) F(N), the function values.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) f(n)
  integer i
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  f(1:n) = 0.0D+00
  do i = 1, m
    f(1:n) = f(1:n) + c(i) * abs ( x(i,1:n) - w(i) )
  end do
  f(1:n) = exp ( - f(1:n) )

  return
end
subroutine p05_q ( m, c, w, q )

!*****************************************************************************80
!
!! P05_Q evaluates the integral for problem p05.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    29 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Output, real ( kind = rk ) Q, the integral.
!
  implicit none

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

  integer m

  real ( kind = rk ) c(m)
  integer i
  real ( kind = rk ) q
  real ( kind = rk ) w(m)

  q = 1.0D+00

  do i = 1, m
!
!  W < 0 < 1
!
!  | X - W | = X - W from 0 to 1.
!
    if ( w(i) < 0.0D+00 ) then

      q = q * &
      ( exp ( - c(i) * (         - w(i) ) ) &
      - exp ( - c(i) * ( 1.0D+00 - w(i) ) ) ) / c(i)
!
!  0 < W < 1
!
!  | X - W | = W - X from 0 to Z, 
!            = X - W from      Z to 1.
!
    else if ( w(i) < 1.0D+00 ) then

      q = q * ( 2.0D+00 &
          - exp ( - c(i) * (           w(i) ) ) &
          - exp ( - c(i) * ( 1.0D+00 - w(i) ) ) ) / c(i)
!
!  0 < 1 < W
!
!  | X - W | = W - X from 0 to 1.
!
    else

      q = q * &
      ( exp ( - c(i) * ( w(i) - 1.0D+00 ) ) &
      - exp ( - c(i) * ( w(i)           ) ) ) / c(i)

    end if

  end do

  return
end
subroutine p05_title ( title )

!*****************************************************************************80
!
!! P05_TITLE returns the title of problem p05.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) TITLE, the title of the problem.
!
  implicit none

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

  character ( len = * ) title

  title = 'Continuous'

  return
end
subroutine p06_d ( m, id, c, w, n, x, d )

!*****************************************************************************80
!
!! P06_D evaluates any derivative component for problem p06.
!
!  Discussion:
!
!    f(x) = exp ( c(1:m) * x(1:m) ) if x(1) <= w(1) and x(2) <= w(2).
!           0                          otherwise
!
!    Default values are:
!
!    c(1:m) = 0.5^(1/m)
!    w(1:2) = 0.5^(1/m)
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, integer ID, the spatial coordinate to differentiate.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) D(N), the ID-th derivative component.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) d(n)
  integer i
  integer id
  integer j
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  if ( m == 1 ) then

    d(1:n) = c(1) * exp ( c(1) * x(1,1:n) )

    do j = 1, n
      if ( w(1) < x(1,j) ) then
        d(j) = 0.0D+00
      end if
    end do

  else

    d(1:n) = 0.0D+00
    do i = 1, m
      d(1:n) = d(1:n) + c(i) * x(i,1:n)
    end do
    d(1:n) = c(id) * exp ( d(1:n) )

    do j = 1, n
      if ( w(1) < x(1,j) .or. w(2) < x(2,j) ) then
        d(j) = 0.0D+00
      end if
    end do

  end if

  return
end
subroutine p06_f ( m, c, w, n, x, f )

!*****************************************************************************80
!
!! P06_F evaluates the function for problem p06.
!
!  Discussion:
!
!    f(x) = exp ( c(1:m) * x(1:m) ) if x(1) <= w(1) and x(2) <= w(2).
!           0                          otherwise
!
!    Default values are:
!
!    c(1:m) = 0.5^(1/m)
!    w(1:2) = 0.5^(1/m)
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Input, integer N, the number of points.
!
!    Input, real ( kind = rk ) X(M,N), the evaluation points.
!
!    Output, real ( kind = rk ) F(N), the function values.
!
  implicit none

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

  integer m
  integer n

  real ( kind = rk ) c(m)
  real ( kind = rk ) f(n)
  integer i
  integer j
  real ( kind = rk ) w(m)
  real ( kind = rk ) x(m,n)

  if ( m == 1 ) then

    f(1:n) = exp ( c(1) * x(1,1:n) )

    do j = 1, n
      if ( w(1) < x(1,j) ) then
        f(j) = 0.0D+00
      end if
    end do

  else

    f(1:n) = 0.0D+00
    do i = 1, m
      f(1:n) = f(1:n) + c(i) * x(i,1:n)
    end do
    f(1:n) = exp ( f(1:n) )

    do j = 1, n
      if ( w(1) < x(1,j) .or. w(2) < x(2,j) ) then
        f(j) = 0.0D+00
      end if
    end do

  end if

  return
end
subroutine p06_q ( m, c, w, q )

!*****************************************************************************80
!
!! P06_Q evaluates the integral for problem p06.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    29 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Alan Genz,
!    A Package for Testing Multiple Integration Subroutines,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    Reidel, 1987, pages 337-340,
!    ISBN: 9027725144,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer M, the dimension of the argument.
!
!    Input, real ( kind = rk ) C(M), W(M), the problem parameters.
!
!    Output, real ( kind = rk ) Q, the integral.
!
  implicit none

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

  integer m

  real ( kind = rk ) c(m)
  integer i
  real ( kind = rk ) q
  real ( kind = rk ) w(m)
!
!  To simplify the calculation, force W(3:M) to be at least 1.0.
!
  w(3:m) = 1.0D+00

  q = 1.0D+00

  do i = 1, m

    if ( w(i) <= 0.0D+00 ) then

      q = q * 0.0D+00

    else if ( w(i) <= 1.0D+00 ) then

      if ( c(i) == 0.0D+00 ) then
        q = q * w(i)
      else
        q = q * ( exp ( c(i) * w(i) ) - 1.0D+00 ) / c(i)
      end if

    else

      if ( c(i) /= 0.0D+00 ) then
        q = q * ( exp ( c(i) * w(i) ) - 1.0D+00 ) / c(i)
      end if

    end if

  end do

  return
end
subroutine p06_title ( title )

!*****************************************************************************80
!
!! P06_TITLE returns the title of problem p06.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    04 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) TITLE, the title of the problem.
!
  implicit none

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

  character ( len = * ) title

  title = 'Discontinuous'

  return
end
function r8_error ( x )

!*****************************************************************************80
!
!! R8_ERROR computes the error function.
!
!  Discussion:
!
!    This function was renamed "R8_ERROR" from "ERF", to avoid a conflict
!    with the name of a corresponding routine often, but not always,
!    supplied as part of the math support library.
!
!    The definition of the error function is:
!
!      erf(x) = ( 2 / sqrt ( pi ) ) * Integral ( 0 <= t <= x ) exp ( -t^2 ) dt
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    15 April 2003
!
!  Author:
!
!    Original FORTRAN77 version by Wayne Fullerton.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    David Kahaner, Cleve Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1989,
!    ISBN: 0-13-627258-4,
!    LC: TA345.K34.
!
!  Parameters:
!
!    Input, real ( kind = rk ) X, the argument of the error function.
!
!    Output, real ( kind = rk ) R8_ERROR, the value of the error function at X.
!
  implicit none

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

  real ( kind = rk ) csevl
  real ( kind = rk ) r8_error
  real ( kind = rk ) r8_errorc
  real ( kind = rk ), parameter, dimension ( 13 ) :: erfcs = (/ &
    -0.049046121234691808D+00, -0.14226120510371364D+00, &
     0.010035582187599796D+00, -0.000576876469976748D+00, &
     0.000027419931252196D+00, -0.000001104317550734D+00, &
     0.000000038488755420D+00, -0.000000001180858253D+00, &
     0.000000000032334215D+00, -0.000000000000799101D+00, &
     0.000000000000017990D+00, -0.000000000000000371D+00, &
     0.000000000000000007D+00 /)
  integer inits
  integer, save :: nterf = 0
  real ( kind = rk ), save :: sqeps = 0.0D+00
  real ( kind = rk ), parameter :: sqrtpi = 1.7724538509055160D+00
  real ( kind = rk ) x
  real ( kind = rk ), save :: xbig = 0.0D+00
  real ( kind = rk ) y
!
!  Initialize the Chebyshev series.
!
  if ( nterf == 0 ) then
    nterf = inits ( erfcs, 13, 0.1D+00 * epsilon ( erfcs ) )
    xbig = sqrt ( - log ( sqrtpi * epsilon ( xbig ) ) )
    sqeps = sqrt ( 2.0D+00 * epsilon ( sqeps ) )
  end if

  y = abs ( x )

  if ( y <= sqeps ) then
    r8_error = 2.0D+00 * x / sqrtpi
  else if ( y <= 1.0D+00 ) then
    r8_error = x &
      * ( 1.0D+00 + csevl ( 2.0D+00 * x**2 - 1.0D+00, erfcs, nterf ) )
  else if ( y <= xbig ) then
    r8_error = sign ( 1.0D+00 - r8_errorc ( y ), x )
  else
    r8_error = sign ( 1.0D+00, x )
  end if

  return
end
function r8_errorc ( x )

!*****************************************************************************80
!
!! R8_ERRORC computes the complementary error function.
!
!  Discussion:
!
!    This function was renamed "R8_ERRORC" from "ERFC", to avoid a conflict
!    with the name of a corresponding routine often, but not always,
!    supplied as part of the math support library.
!
!    The definition of the complementary error function is:
!
!      ERFC(X) = 1 - ERF(X)
!
!    where ERF(X) is the error function.
!
!  Modified:
!
!    26 August 2001
!
!  Author:
!
!    Original FORTRAN77 version by Wayne Fullerton.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    David Kahaner, Cleve Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1989,
!    ISBN: 0-13-627258-4,
!    LC: TA345.K34.
!
!  Parameters:
!
!    Input, real ( kind = rk ) X, the argument of the error function.
!
!    Output, real ( kind = rk ) R8_ERRORC, the value of the complementary 
!    error function at X.
!
  implicit none

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

  real ( kind = rk ) csevl
  real ( kind = rk ), parameter, dimension ( 13 ) :: erfcs = (/ &
   -0.049046121234691808D+00, &
   -0.14226120510371364D+00,  &
    0.010035582187599796D+00, &
   -0.000576876469976748D+00, &
    0.000027419931252196D+00, &
   -0.000001104317550734D+00, &
    0.000000038488755420D+00, &
   -0.000000001180858253D+00, &
    0.000000000032334215D+00, &
   -0.000000000000799101D+00, &
    0.000000000000017990D+00, &
   -0.000000000000000371D+00, &
    0.000000000000000007D+00 /)
  real ( kind = rk ), parameter, dimension ( 24 ) :: erfccs = (/ &
    0.0715179310202925D+00, &
   -0.026532434337606719D+00, &
    0.001711153977920853D+00, &
   -0.000163751663458512D+00, &
    0.000019871293500549D+00, &
   -0.000002843712412769D+00, &
    0.000000460616130901D+00, &
   -0.000000082277530261D+00, &
    0.000000015921418724D+00, &
   -0.000000003295071356D+00, &
    0.000000000722343973D+00, &
   -0.000000000166485584D+00, &
    0.000000000040103931D+00, &
   -0.000000000010048164D+00, &
    0.000000000002608272D+00, &
   -0.000000000000699105D+00, &
    0.000000000000192946D+00, &
   -0.000000000000054704D+00, &
    0.000000000000015901D+00, &
   -0.000000000000004729D+00, &
    0.000000000000001432D+00, &
   -0.000000000000000439D+00, &
    0.000000000000000138D+00, &
   -0.000000000000000048D+00 /)
  real ( kind = rk ), parameter, dimension ( 23 ) :: erc2cs = (/ &
   -0.069601346602309501D+00, &
   -0.041101339362620893D+00, &
    0.003914495866689626D+00, &
   -0.000490639565054897D+00, &
    0.000071574790013770D+00, &
   -0.000011530716341312D+00, &
    0.000001994670590201D+00, &
   -0.000000364266647159D+00, &
    0.000000069443726100D+00, &
   -0.000000013712209021D+00, &
    0.000000002788389661D+00, &
   -0.000000000581416472D+00, &
    0.000000000123892049D+00, &
   -0.000000000026906391D+00, &
    0.000000000005942614D+00, &
   -0.000000000001332386D+00, &
    0.000000000000302804D+00, &
   -0.000000000000069666D+00, &
    0.000000000000016208D+00, &
   -0.000000000000003809D+00, &
    0.000000000000000904D+00, &
   -0.000000000000000216D+00, &
    0.000000000000000052D+00 /)
  real ( kind = rk ) eta
  integer inits
  integer, save :: nterc2 = 0
  integer, save :: nterf = 0
  integer, save :: nterfc = 0
  real ( kind = rk ) r8_errorc
  real ( kind = rk ), save :: sqeps = 0.0D+00
  real ( kind = rk ), parameter :: sqrtpi = 1.7724538509055160D+00
  real ( kind = rk ) x
  real ( kind = rk ), save :: xmax = 0.0D+00
  real ( kind = rk ), save :: xsml = 0.0D+00
  real ( kind = rk ) y

  if ( nterf == 0 ) then

    eta = 0.1D+00 * epsilon ( eta )
    nterf = inits ( erfcs, 13, eta )
    nterfc = inits ( erfccs, 24, eta )
    nterc2 = inits ( erc2cs, 23, eta )

    xsml = -sqrt ( - log ( sqrtpi * epsilon ( xsml ) ) )
    xmax = sqrt ( - log ( sqrtpi * tiny ( xmax ) ) )
    xmax = xmax - 0.5D+00 * log ( xmax ) / xmax - 0.01D+00
    sqeps = sqrt ( 2.0D+00 * epsilon ( sqeps ) )

  end if

  if ( x <= xsml ) then
    r8_errorc = 2.0D+00
    return
  end if
!
!  X so big that ERFC will underflow.
!
  if ( xmax < x ) then
    r8_errorc = 0.0D+00
    return
  end if

  y = abs ( x )
!
!  erfc(x) = 1.0D+00 - erf(x) for -1 <= x <= 1.
!
  if ( y <= 1.0D+00 ) then

    if ( y < sqeps ) then
      r8_errorc = 1.0D+00 - 2.0D+00 * x / sqrtpi
    else if ( sqeps <= y ) then
      r8_errorc = 1.0D+00 - x * ( 1.0D+00 + &
        csevl ( 2.0D+00 * x * x - 1.0D+00, erfcs, nterf ) )
    end if

    return

  end if
!
!  For 1 < |x| <= xmax, erfc(x) = 1.0D+00 - erf(x)
!
  y = y * y

  if ( y <= 4.0D+00 ) then
    r8_errorc = exp ( -y ) / abs ( x ) * ( 0.5D+00 &
      + csevl ( ( 8.0D+00 / y - 5.0D+00 ) / 3.0D+00, erc2cs, nterc2 ) )
  else
    r8_errorc = exp ( -y ) / abs ( x ) * ( 0.5D+00 &
      + csevl ( 8.0D+00 / y - 1.0D+00, erfccs, nterfc ) )
  end if

  if ( x < 0.0D+00 ) then
    r8_errorc = 2.0D+00 - r8_errorc
  end if

  return
end
subroutine tuple_next ( m1, m2, n, rank, x )

!*****************************************************************************80
!
!! TUPLE_NEXT computes the next element of a tuple space.
!
!  Discussion:
!
!    The elements are N vectors.  Each entry is constrained to lie
!    between M1 and M2.  The elements are produced one at a time.
!    The first element is
!      (M1,M1,...,M1),
!    the second element is
!      (M1,M1,...,M1+1),
!    and the last element is
!      (M2,M2,...,M2)
!    Intermediate elements are produced in lexicographic order.
!
!    For example, if our scalar input was
!
!      N = 2,
!      M1 = 1,
!      M2 = 3
!
!    then the sequence of input and output vectors would be:
!
!      INPUT        OUTPUT
!      -------      -------
!      Rank  X      Rank   X
!      ----  ---    -----  ---
!      0     * *    1      1 1
!      1     1 1    2      1 2
!      2     1 2    3      1 3
!      3     1 3    4      2 1
!      4     2 1    5      2 2
!      5     2 2    6      2 3
!      6     2 3    7      3 1
!      7     3 1    8      3 2
!      8     3 2    9      3 3
!      9     3 3    0      0 0
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    18 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M1, M2, the minimum and maximum entries.
!
!    Input, integer N, the number of components.
!
!    Input/output, integer RANK, counts the elements.
!    On first call, set RANK to 0.  Thereafter, the output value of RANK
!    will indicate the order of the element returned.  When there are no
!    more elements, RANK will be returned as 0.
!
!    Input/output, integer X(N), on input the previous tuple.
!    On output, the next tuple.
!
  implicit none

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

  integer n

  integer i
  integer m1
  integer m2
  integer rank
  integer x(n)

  if ( m2 < m1 ) then
    rank = 0
    return
  end if

  if ( rank <= 0 ) then

    x(1:n) = m1
    rank = 1

  else

    rank = rank + 1
    i = n

    do

      if ( x(i) < m2 ) then
        x(i) = x(i) + 1
        exit
      end if

      x(i) = m1

      if ( i == 1 ) then
        rank = 0
        x(1:n) = m1
        exit
      end if

      i = i - 1

    end do

  end if

  return
end