subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5, 6 and 9, which
!    are commonly reserved for console I/O).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    18 September 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT, the free unit number.
!
  implicit none

  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
function i4_uniform ( a, b, seed )

!*****************************************************************************80
!
!! I4_UNIFORM returns a scaled pseudorandom I4.
!
!  Discussion:
!
!    An I4 is an integer ( kind = 4 ) value.
!
!    The pseudorandom number will be scaled to be uniformly distributed
!    between A and B.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    12 November 2006
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Springer Verlag, pages 201-202, 1983.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley Interscience, page 95, 1998.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, pages 362-376, 1986.
!
!    Peter Lewis, Allen Goodman, James Miller
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, pages 136-143, 1969.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) A, B, the limits of the interval.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, integer ( kind = 4 ) I4_UNIFORM, a number between A and B.
!
  implicit none

  integer ( kind = 4 ) a
  integer ( kind = 4 ) b
  integer ( kind = 4 ) i4_uniform
  integer ( kind = 4 ) k
  real    ( kind = 4 ) r
  integer ( kind = 4 ) seed
  integer ( kind = 4 ) value

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'I4_UNIFORM - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop
  end if

  k = seed / 127773

  seed = 16807 * ( seed - k * 127773 ) - k * 2836

  if ( seed < 0 ) then
    seed = seed + 2147483647
  end if

  r = real ( seed, kind = 4 ) * 4.656612875E-10
!
!  Scale R to lie between A-0.5 and B+0.5.
!
  r = ( 1.0E+00 - r ) * ( real ( min ( a, b ), kind = 4 ) - 0.5E+00 ) & 
    +             r   * ( real ( max ( a, b ), kind = 4 ) + 0.5E+00 )
!
!  Use rounding to convert R to an integer between A and B.
!
  value = nint ( r, kind = 4 )

  value = max ( value, min ( a, b ) )
  value = min ( value, max ( a, b ) )

  i4_uniform = value

  return
end
function i4_bit_hi1 ( n )

!*****************************************************************************80
!
!! I4_BIT_HI1 returns the position of the high 1 bit base 2 in an integer.
!
!  Discussion:
!
!    This routine uses the default integer precision, which is
!    presumed to correspond to a KIND of 4.
!
!  Example:
!
!       N    Binary    Hi 1
!    ----    --------  ----
!       0           0     0
!       1           1     1
!       2          10     2
!       3          11     2 
!       4         100     3
!       5         101     3
!       6         110     3
!       7         111     3
!       8        1000     4
!       9        1001     4
!      10        1010     4
!      11        1011     4
!      12        1100     4
!      13        1101     4
!      14        1110     4
!      15        1111     4
!      16       10000     5
!      17       10001     5
!    1023  1111111111    10
!    1024 10000000000    11
!    1025 10000000001    11
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    28 May 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the integer to be measured.
!    N should be nonnegative.  If N is nonpositive, I4_BIT_HI1
!    will always be 0.
!
!    Output, integer I4_BIT_HI1, the number of bits base 2.
!
  implicit none

  integer bit
  integer i4_bit_hi1
  integer i
  integer n

  i = n
  bit = 0

  do

    if ( i <= 0 ) then
      exit
    end if

    bit = bit + 1
    i = i / 2

  end do

  i4_bit_hi1 = bit

  return
end
function i4_bit_lo0 ( n )

!*****************************************************************************80
!
!! I4_BIT_LO0 returns the position of the low 0 bit base 2 in an integer.
!
!  Discussion:
!
!    This routine uses the default integer precision, which is
!    presumed to correspond to a KIND of 4.
!
!  Example:
!
!       N    Binary    Lo 0
!    ----    --------  ----
!       0           0     1
!       1           1     2
!       2          10     1
!       3          11     3 
!       4         100     1
!       5         101     2
!       6         110     1
!       7         111     4
!       8        1000     1
!       9        1001     2
!      10        1010     1
!      11        1011     3
!      12        1100     1
!      13        1101     2
!      14        1110     1
!      15        1111     5
!      16       10000     1
!      17       10001     2
!    1023  1111111111     1
!    1024 10000000000     1
!    1025 10000000001     1
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    28 May 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the integer to be measured.
!    N should be nonnegative.
!
!    Output, integer I4_BIT_LO0, the position of the low 1 bit.
!
  implicit none

  integer bit
  integer i
  integer i2
  integer i4_bit_lo0
  integer n

  bit = 0
  i = n

  do

    bit = bit + 1
    i2 = i / 2

    if ( i == 2 * i2 ) then
      exit
    end if

    i = i2

  end do

  i4_bit_lo0 = bit

  return
end
subroutine i4_sobol ( dim_num, seed, quasi )

!*****************************************************************************80
!
!! I4_SOBOL generates a new quasirandom Sobol vector with each call.
!
!  Discussion:
!
!    The routine adapts the ideas of Antonov and Saleev.
!
!    This routine uses the default integer precision, which is
!    presumed to correspond to a KIND of 4.
!
!    Thanks to Francis Dalaudier for pointing out that the range of allowed
!    values of DIM_NUM should start at 1, not 2!  17 February 2009.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    17 February 2009
!
!  Author:
!
!    Original FORTRAN77 version by Bennett Fox.
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    Antonov, Saleev,
!    USSR Computational Mathematics and Mathematical Physics,
!    Volume 19, 1980, pages 252 - 256.
!
!    Paul Bratley, Bennett Fox,
!    Algorithm 659:
!    Implementing Sobol's Quasirandom Sequence Generator,
!    ACM Transactions on Mathematical Software,
!    Volume 14, Number 1, pages 88-100, 1988.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom 
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, pages 362-376, 1986.
!
!    Ilya Sobol,
!    USSR Computational Mathematics and Mathematical Physics,
!    Volume 16, pages 236-242, 1977.
!
!    Ilya Sobol, Levitan, 
!    The Production of Points Uniformly Distributed in a Multidimensional 
!    Cube (in Russian),
!    Preprint IPM Akad. Nauk SSSR, 
!    Number 40, Moscow 1976.
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!    DIM_NUM must satisfy 1 <= DIM_NUM <= 40.
!
!    Input/output, integer SEED, the "seed" for the sequence.
!    This is essentially the index in the sequence of the quasirandom
!    value to be generated.  On output, SEED has been set to the
!    appropriate next value, usually simply SEED+1.
!    If SEED is less than 0 on input, it is treated as though it were 0.
!    An input value of 0 requests the first (0-th) element of the sequence.
!
!    Output, real ( kind = 4 ) QUASI(DIM_NUM), the next quasirandom vector.
!
  implicit none

  integer dim_num
  integer, parameter :: dim_max = 40
  integer, parameter :: log_max = 30

  integer atmost
  integer, save :: dim_num_save = 0
  integer i
  integer i4_bit_hi1
  integer i4_bit_lo0
  integer inc
  logical includ(8)
  logical, save :: initialized = .false.
  integer j
  integer j2
  integer k
  integer l
  integer, save, dimension(dim_max) :: lastq
  integer m
  integer, save :: maxcol
  integer newv
  integer, save, dimension(dim_max) :: poly
  real ( kind = 4 ) quasi(dim_num)
  real ( kind = 4 ), save :: recipd
  integer seed
  integer, save :: seed_save = - 1
  integer seed_temp
  integer, save, dimension(dim_max,log_max) :: v

  if ( .not. initialized .or. dim_num /= dim_num_save ) then

    initialized = .true.
!
!  Initialize (part of) V.
!
    v(1:40,1) = (/ &
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1 /)

    v(3:40,2) = (/ &
            1, 3, 1, 3, 1, 3, 3, 1, &
      3, 1, 3, 1, 3, 1, 1, 3, 1, 3, &
      1, 3, 1, 3, 3, 1, 3, 1, 3, 1, &
      3, 1, 1, 3, 1, 3, 1, 3, 1, 3 /)

    v(4:40,3) = (/ &
               7, 5, 1, 3, 3, 7, 5, & 
      5, 7, 7, 1, 3, 3, 7, 5, 1, 1, &
      5, 3, 3, 1, 7, 5, 1, 3, 3, 7, &
      5, 1, 1, 5, 7, 7, 5, 1, 3, 3 /)

    v(6:40,4) = (/ &
                     1, 7, 9,13,11, &
      1, 3, 7, 9, 5,13,13,11, 3,15, &
      5, 3,15, 7, 9,13, 9, 1,11, 7, &
      5,15, 1,15,11, 5, 3, 1, 7, 9 /)
  
    v(8:40,5) = (/ &
                           9, 3,27, &
     15,29,21,23,19,11,25, 7,13,17, &
      1,25,29, 3,31,11, 5,23,27,19, &
     21, 5, 1,17,13, 7,15, 9,31, 9 /)

    v(14:40,6) = (/ &
              37,33, 7, 5,11,39,63, &
     27,17,15,23,29, 3,21,13,31,25, &
      9,49,33,19,29,11,19,27,15,25 /)

    v(20:40,7) = (/ &
                                         13, &
     33,115, 41, 79, 17, 29,119, 75, 73,105, &
      7, 59, 65, 21,  3,113, 61, 89, 45,107 /)

    v(38:40,8) = (/ &
                                  7, 23, 39 /)
!
!  Set POLY.
!
    poly(1:40)= (/ &
        1,   3,   7,  11,  13,  19,  25,  37,  59,  47, &
       61,  55,  41,  67,  97,  91, 109, 103, 115, 131, &
      193, 137, 145, 143, 241, 157, 185, 167, 229, 171, &
      213, 191, 253, 203, 211, 239, 247, 285, 369, 299 /)
!
!  Check parameters.
!
    if ( dim_num < 1 .or. dim_max < dim_num ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'I4_SOBOL - Fatal error!'
      write ( *, '(a)' ) '  The spatial dimension DIM_NUM should satisfy:'
      write ( *, '(a,i8)' ) '    1 <= DIM_NUM <= ', dim_max
      write ( *, '(a,i8)' ) '  But this input value is DIM_NUM = ', dim_num
      stop
    end if

    dim_num_save = dim_num
!
!  The original line
!    ATMOST = 2**LOG_MAX - 1 
!  will not execute properly in the double integer arithmetic.  
!  It produces the result ATMOST = -1.
!
!  The following lines do produce the correct result.
!
!  They could also be replaced by the FORTRAN90 specific
!
!    atmost = huge ( atmost )
!
    atmost = 1
    inc = 1
    do i = 1, log_max
      inc = 2 * inc
      atmost = atmost + inc
    end do
!
!  Find the number of bits in ATMOST.
!
    maxcol = i4_bit_hi1 ( atmost )
!
!  Initialize row 1 of V.
!
    v(1,1:maxcol) = 1
!
!  Initialize the remaining rows of V.
!
    do i = 2, dim_num
!
!  The bit pattern of the integer POLY(I) gives the form
!  of polynomial I.
!
!  Find the degree of polynomial I from binary encoding.
!
      j = poly(i)
      m = 0

      do

        j = j / 2

        if ( j <= 0 ) then
          exit
        end if

        m = m + 1

      end do
!
!  We expand this bit pattern to separate components
!  of the logical array INCLUD.
!
      j = poly(i)
      do k = m, 1, -1
        j2 = j / 2
        includ(k) = ( j /= ( 2 * j2 ) )
        j = j2
      end do
!
!  Calculate the remaining elements of row I as explained
!  in Bratley and Fox, section 2.
!
      do j = m + 1, maxcol

        newv = v(i,j-m)
        l = 1

        do k = 1, m

          l = 2 * l

          if ( includ(k) ) then
            newv = ieor ( newv, l * v(i,j-k) )
          end if

        end do

        v(i,j) = newv

      end do

    end do
!
!  Multiply columns of V by appropriate power of 2.
!
    l = 1
    do j = maxcol - 1, 1, -1
      l = 2 * l
      v(1:dim_num,j) = v(1:dim_num,j) * l
    end do
!
!  RECIPD is 1/(common denominator of the elements in V) = 1 / ( 2 * L ).
!
    recipd = real ( l )
    recipd = 0.5E+00 / recipd

  end if

  if ( seed < 0 ) then
    seed = 0
  end if

  if ( seed == 0 ) then

    l = 1
    lastq(1:dim_num) = 0

  else if ( seed == seed_save + 1 ) then
!
!  Find the position of the right-hand zero in SEED.
!
    l = i4_bit_lo0 ( seed )

  else if ( seed <= seed_save ) then

    seed_save = 0
    l = 1
    lastq(1:dim_num) = 0

    do seed_temp = seed_save, seed-1
      l = i4_bit_lo0 ( seed_temp )
      lastq(1:dim_num) = ieor ( lastq(1:dim_num), v(1:dim_num,l) )
    end do

    l = i4_bit_lo0 ( seed )

  else if ( seed_save+1 < seed ) then

    do seed_temp = seed_save+1, seed-1
      l = i4_bit_lo0 ( seed_temp )
      lastq(1:dim_num) = ieor ( lastq(1:dim_num), v(1:dim_num,l) )
    end do

    l = i4_bit_lo0 ( seed )

  end if
!
!  Check that the user is not calling too many times!
!
  if ( maxcol < l ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'I4_SOBOL - Fatal error!'
    write ( *, '(a)' ) '  Too many calls!'
    write ( *, '(a,i12)' ) '  MAXCOL = ', maxcol
    write ( *, '(a,i12)' ) '  L =      ', l
    stop
  end if
!
!  Calculate the new components of QUASI.
!
  quasi(1:dim_num) = real ( lastq(1:dim_num), kind = 4 ) * recipd
  lastq(1:dim_num) = ieor ( lastq(1:dim_num), v(1:dim_num,l) )

  seed_save = seed
  seed = seed + 1

  return
end
subroutine i4_sobol_generate ( m, n, skip, r )

!*****************************************************************************80
!
!! I4_SOBOL_GENERATE generates a Sobol dataset.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    17 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points to generate.
!
!    Input, integer SKIP, the number of initial points to skip.
!
!    Output, real R(M,N), the points.
!
  implicit none

  integer m
  integer n

  integer j
  real, dimension ( m, n ) :: r
  integer seed
  integer skip

  do j = 1, n
    seed = skip + j - 1
    call i4_sobol ( m, seed, r(1:m,j) )
  end do

  return
end
subroutine i4_sobol_write ( m, n, skip, r, file_out_name )

!*****************************************************************************80
!
!! I4_SOBOL_WRITE writes a Sobol dataset to a file.
!
!  Discussion:
!
!    The initial lines of the file are comments, which begin with a
!    '#' character.
!
!    Thereafter, each line of the file contains the M-dimensional
!    components of the SKIP+I-1 entry of the Sobol sequence.
!
!    For the Sobol sequence, the value of SKIP is the same
!    as the value of SEED used to generate the first point.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    17 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of (successive) points.
!
!    Input, integer  SKIP, the number of skipped points.
!
!    Input, real R(M,N), the points.
!
!    Input, character ( len = * ) FILE_OUT_NAME, the name of
!    the output file.
!
  implicit none

  integer m
  integer n

  character ( len = * ) file_out_name
  integer file_out_unit
  integer ios
  integer j
  real r(m,n)
  integer skip
  character ( len = 40 ) string

  call get_unit ( file_out_unit )

  open ( unit = file_out_unit, file = file_out_name, status = 'replace', &
    iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'I4_SOBOL_WRITE - Fatal error!'
    write ( *, '(a)' ) '  Could not open the output file.'
    stop
  end if

  call timestring ( string )

  write ( file_out_unit, '(a)'       ) '#  ' // trim ( file_out_name )
  write ( file_out_unit, '(a)'       ) '#  created by I4_SOBOL_WRITE.F90.'
  write ( file_out_unit, '(a)'       ) '#'
  write ( file_out_unit, '(a)'       ) '#  File generated on ' &
    // trim ( string )
  write ( file_out_unit, '(a)'       ) '#'
  write ( file_out_unit, '(a,i8)'    ) '#  Spatial dimension M = ', m
  write ( file_out_unit, '(a,i8)'    ) '#  Number of points N = ', n
  write ( file_out_unit, '(a,g14.6)' ) '#  Epsilon (unit roundoff) = ', &
    epsilon ( r(1,1) )
  write ( file_out_unit, '(a,i8)'    ) '#  Initial values skipped = ', skip
  write ( file_out_unit, '(a)'       ) '#'

  write ( string, '(a,i3,a)' ) '(', m, '(2x,f10.6))'
  do j = 1, n
    write ( file_out_unit, string ) r(1:m,j)
  end do

  close ( unit = file_out_unit )

  return
end
function i4_xor ( i, j )

!*****************************************************************************80
!
!! I4_XOR calculates the exclusive OR of two integers.
!
!  Discussion:
!
!    This function is NOT needed in FORTRAN90, which supplies the
!    intrinsic IEOR function for this purpose.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    16 February 2005
!
!  Author:
!
!   John Burkardt
!
!  Parameters:
!
!    Input, integer I, J, two values whose exclusive OR is needed.
!
!    Output, integer I4_XOR, the exclusive OR of I and J.
!
  implicit none

  integer i
  integer i1
  integer i2
  integer i4_xor
  integer j
  integer j1
  integer j2
  integer k
  integer l

  i1 = i
  j1 = j
  k = 0
  l = 1

  do while ( i1 /= 0 .or. j1 /= 0 )

    i2 = i1 / 2
    j2 = j1 / 2

    if ( &
      ( ( i1 == 2 * i2 ) .and. ( j1 /= 2 * j2 ) ) .or. &
      ( ( i1 /= 2 * i2 ) .and. ( j1 == 2 * j2 ) ) ) then
      k = k + l
    end if

    i1 = i2
    j1 = j2
    l = 2 * l

  end do

  i4_xor = k

  return
end
function i8_bit_hi1 ( n )

!*****************************************************************************80
!
!! I8_BIT_HI1 returns the position of the high 1 bit base 2 in an integer.
!
!  Discussion:
!
!    This routine uses the integer precision corresponding to a KIND of 8.
!
!  Example:
!
!       N    Binary    Hi 1
!    ----    --------  ----
!       0           0     0
!       1           1     1
!       2          10     2
!       3          11     2 
!       4         100     3
!       5         101     3
!       6         110     3
!       7         111     3
!       8        1000     4
!       9        1001     4
!      10        1010     4
!      11        1011     4
!      12        1100     4
!      13        1101     4
!      14        1110     4
!      15        1111     4
!      16       10000     5
!      17       10001     5
!    1023  1111111111    10
!    1024 10000000000    11
!    1025 10000000001    11
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    28 May 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 8 ) N, the integer to be measured.
!    N should be nonnegative.  If N is nonpositive, I8_BIT_HI1
!    will always be 0.
!
!    Output, integer ( kind = 8 ) I8_BIT_HI1, the number of bits base 2.
!
  implicit none

  integer ( kind = 8 ) :: bit
  integer ( kind = 8 ) :: i8_bit_hi1
  integer ( kind = 8 ) :: i
  integer ( kind = 8 ) :: n

  i = n
  bit = 0

  do

    if ( i <= 0 ) then
      exit
    end if

    bit = bit + 1
    i = i / 2

  end do

  i8_bit_hi1 = bit

  return
end
function i8_bit_lo0 ( n )

!*****************************************************************************80
!
!! I8_BIT_LO0 returns the position of the low 0 bit base 2 in an integer.
!
!  Discussion:
!
!    This routine uses the integer precision corresponding to a KIND of 8.
!
!  Example:
!
!       N    Binary    Lo 0
!    ----    --------  ----
!       0           0     1
!       1           1     2
!       2          10     1
!       3          11     3 
!       4         100     1
!       5         101     2
!       6         110     1
!       7         111     4
!       8        1000     1
!       9        1001     2
!      10        1010     1
!      11        1011     3
!      12        1100     1
!      13        1101     2
!      14        1110     1
!      15        1111     5
!      16       10000     1
!      17       10001     2
!    1023  1111111111     1
!    1024 10000000000     1
!    1025 10000000001     1
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    28 May 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 8 ) N, the integer to be measured.
!    N should be nonnegative.
!
!    Output, integer ( kind = 8 ) I8_BIT_LO0, the position of the low 1 bit.
!
  implicit none

  integer ( kind = 8 ) :: bit
  integer ( kind = 8 ) :: i
  integer ( kind = 8 ) :: i2
  integer ( kind = 8 ) :: i8_bit_lo0
  integer ( kind = 8 ) :: n

  bit = 0
  i = n

  do

    bit = bit + 1
    i2 = i / 2

    if ( i == 2 * i2 ) then
      exit
    end if

    i = i2

  end do

  i8_bit_lo0 = bit

  return
end
subroutine i8_sobol ( dim_num, seed, quasi )

!*****************************************************************************80
!
!! I8_SOBOL generates a new quasirandom Sobol vector with each call.
!
!  Discussion:
!
!    The routine adapts the ideas of Antonov and Saleev.
!
!    This routine uses the integer and real precisions corresponding 
!    to a KIND of 8.  Many variables don't need to be set to this
!    higher precision, but I'm a little uncertain of the algorithm
!    and of FORTRAN90's double word integer arithmetic, so I set
!    EVERYTHING to KIND = 8.
!
!    Thanks to Francis Dalaudier for pointing out that the range of allowed
!    values of DIM_NUM should start at 1, not 2!  17 February 2009.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    17 February 2009
!
!  Author:
!
!    Original FORTRAN77 version by Bennett Fox.
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    Antonov, Saleev,
!    USSR Computational Mathematics and Mathematical Physics,
!    Volume 19, 1980, pages 252 - 256.
!
!    Paul Bratley, Bennett Fox,
!    Algorithm 659:
!    Implementing Sobol's Quasirandom Sequence Generator,
!    ACM Transactions on Mathematical Software,
!    Volume 14, Number 1, pages 88-100, 1988.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom 
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, pages 362-376, 1986.
!
!    Ilya Sobol,
!    USSR Computational Mathematics and Mathematical Physics,
!    Volume 16, pages 236-242, 1977.
!
!    Ilya Sobol, Levitan, 
!    The Production of Points Uniformly Distributed in a Multidimensional 
!    Cube (in Russian),
!    Preprint IPM Akad. Nauk SSSR, 
!    Number 40, Moscow 1976.
!
!  Parameters:
!
!    Input, integer ( kind = 8 ) DIM_NUM, the number of spatial dimensions.
!    DIM_NUM must satisfy 1 <= DIM_NUM <= 40.
!
!    Input/output, integer ( kind = 8 ) SEED, the "seed" for the sequence.
!    This is essentially the index in the sequence of the quasirandom
!    value to be generated.  On output, SEED has been set to the
!    appropriate next value, usually simply SEED+1.
!    If SEED is less than 0 on input, it is treated as though it were 0.
!    An input value of 0 requests the first (0-th) element of the sequence.
!
!    Output, real ( kind = 8 ) QUASI(DIM_NUM), the next quasirandom vector.
!
  implicit none

  integer ( kind = 8 ) :: dim_num
  integer ( kind = 8 ), parameter :: dim_max = 40
  integer ( kind = 8 ), parameter :: log_max = 62

  integer ( kind = 8 ) :: atmost
  integer ( kind = 8 ), save :: dim_num_save = 0
  integer ( kind = 8 ) :: i
  integer ( kind = 8 ) :: i8_bit_hi1
  integer ( kind = 8 ) :: i8_bit_lo0
  integer ( kind = 8 ) :: inc
  logical includ(8)
  logical, save :: initialized = .false.
  integer ( kind = 8 ) :: j
  integer ( kind = 8 ) :: j2
  integer ( kind = 8 ) :: k
  integer ( kind = 8 ) :: l
  integer ( kind = 8 ), save, dimension(dim_max) :: lastq
  integer ( kind = 8 ) :: m
  integer ( kind = 8 ), save :: maxcol
  integer ( kind = 8 ) :: newv
  integer ( kind = 8 ), save, dimension(dim_max) :: poly
  real ( kind = 8 ), dimension ( dim_num ) :: quasi
  real ( kind = 8 ), save :: recipd
  integer ( kind = 8 ) ::  seed
  integer ( kind = 8 ), save :: seed_save = - 1
  integer ( kind = 8 ) ::  seed_temp
  integer ( kind = 8 ), save, dimension(dim_max,log_max) :: v

  if ( .not. initialized .or. dim_num /= dim_num_save ) then

    initialized = .true.
!
!  Initialize (part of) V.
!
    v(1:40,1) = (/ &
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1 /)

    v(3:40,2) = (/ &
            1, 3, 1, 3, 1, 3, 3, 1, &
      3, 1, 3, 1, 3, 1, 1, 3, 1, 3, &
      1, 3, 1, 3, 3, 1, 3, 1, 3, 1, &
      3, 1, 1, 3, 1, 3, 1, 3, 1, 3 /)

    v(4:40,3) = (/ &
               7, 5, 1, 3, 3, 7, 5, & 
      5, 7, 7, 1, 3, 3, 7, 5, 1, 1, &
      5, 3, 3, 1, 7, 5, 1, 3, 3, 7, &
      5, 1, 1, 5, 7, 7, 5, 1, 3, 3 /)

    v(6:40,4) = (/ &
                     1, 7, 9,13,11, &
      1, 3, 7, 9, 5,13,13,11, 3,15, &
      5, 3,15, 7, 9,13, 9, 1,11, 7, &
      5,15, 1,15,11, 5, 3, 1, 7, 9 /)
  
    v(8:40,5) = (/ &
                           9, 3,27, &
     15,29,21,23,19,11,25, 7,13,17, &
      1,25,29, 3,31,11, 5,23,27,19, &
     21, 5, 1,17,13, 7,15, 9,31, 9 /)

    v(14:40,6) = (/ &
              37,33, 7, 5,11,39,63, &
     27,17,15,23,29, 3,21,13,31,25, &
      9,49,33,19,29,11,19,27,15,25 /)

    v(20:40,7) = (/ &
                                         13, &
     33,115, 41, 79, 17, 29,119, 75, 73,105, &
      7, 59, 65, 21,  3,113, 61, 89, 45,107 /)

    v(38:40,8) = (/ &
                                  7, 23, 39 /)
!
!  Set POLY.
!
    poly(1:40)= (/ &
        1,   3,   7,  11,  13,  19,  25,  37,  59,  47, &
       61,  55,  41,  67,  97,  91, 109, 103, 115, 131, &
      193, 137, 145, 143, 241, 157, 185, 167, 229, 171, &
      213, 191, 253, 203, 211, 239, 247, 285, 369, 299 /)
!
!  Check parameters.
!
    if ( dim_num < 1 .or. dim_max < dim_num ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'I8_SOBOL - Fatal error!'
      write ( *, '(a)' ) '  The spatial dimension DIM_NUM should satisfy:'
      write ( *, '(a,i8)' ) '    1 <= DIM_NUM <= ', dim_max
      write ( *, '(a,i8)' ) '  But this input value is DIM_NUM = ', dim_num
      stop
    end if

    dim_num_save = dim_num
!
!  The original line 
!    ATMOST = 2**LOG_MAX - 1 
!  will not execute properly
!  in the double integer arithmetic.  It produces the result ATMOST = -1.
!
!  The following lines do produce the correct result.
!
!  They could also be replaced by the FORTRAN90 specific
!
!    atmost = huge ( atmost )
!
    atmost = 1
    inc = 1
    do i = 1, log_max
      inc = 2 * inc
      atmost = atmost + inc
    end do
!
!  Find the number of bits in ATMOST.
!
    maxcol = i8_bit_hi1 ( atmost )
!
!  Initialize row 1 of V.
!
    v(1,1:maxcol) = 1
!
!  Initialize the remaining rows of V.
!
    do i = 2, dim_num
!
!  The bit pattern of the integer POLY(I) gives the form
!  of polynomial I.
!
!  Find the degree of polynomial I from binary encoding.
!
      j = poly(i)
      m = 0

      do

        j = j / 2

        if ( j <= 0 ) then
          exit
        end if

        m = m + 1

      end do
!
!  We expand this bit pattern to separate components
!  of the logical array INCLUD.
!
      j = poly(i)
      do k = m, 1, -1
        j2 = j / 2
        includ(k) = ( j /= ( 2 * j2 ) )
        j = j2
      end do
!
!  Calculate the remaining elements of row I as explained
!  in Bratley and Fox, section 2.
!
      do j = m + 1, maxcol

        newv = v(i,j-m)
        l = 1

        do k = 1, m

          l = 2 * l

          if ( includ(k) ) then
            newv = ieor ( newv, l * v(i,j-k) )
          end if

        end do

        v(i,j) = newv

      end do

    end do
!
!  Multiply columns of V by appropriate power of 2.
!
    l = 1
    do j = maxcol - 1, 1, -1
      l = 2 * l
      v(1:dim_num,j) = v(1:dim_num,j) * l
    end do
!
!  RECIPD is 1/(common denominator of the elements in V) = 1 / ( 2 * L ).
!
    recipd = real ( l, kind = 8 )
    recipd = 0.5D+00 / recipd

  end if

  if ( seed < 0 ) then
    seed = 0
  end if

  if ( seed == 0 ) then

    l = 1
    lastq(1:dim_num) = 0

  else if ( seed == seed_save + 1 ) then
!
!  Find the position of the right-hand zero in SEED.
!
    l = i8_bit_lo0 ( seed )

  else if ( seed <= seed_save ) then

    seed_save = 0
    l = 1
    lastq(1:dim_num) = 0

    do seed_temp = seed_save, seed-1
      l = i8_bit_lo0 ( seed_temp )
      lastq(1:dim_num) = ieor ( lastq(1:dim_num), v(1:dim_num,l) )
    end do

    l = i8_bit_lo0 ( seed )

  else if ( seed_save+1 < seed ) then

    do seed_temp = seed_save+1, seed-1
      l = i8_bit_lo0 ( seed_temp )
      lastq(1:dim_num) = ieor ( lastq(1:dim_num), v(1:dim_num,l) )
    end do

    l = i8_bit_lo0 ( seed )

  end if
!
!  Check that the user is not calling too many times!
!
  if ( maxcol < l ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'I8_SOBOL - Fatal error!'
    write ( *, '(a)' ) '  Too many calls!'
    write ( *, '(a,i12)' ) '  MAXCOL = ', maxcol
    write ( *, '(a,i12)' ) '  L =      ', l
    stop
  end if
!
!  Calculate the new components of QUASI.
!
  quasi(1:dim_num) = real ( lastq(1:dim_num), kind = 8 ) * recipd
  lastq(1:dim_num) = ieor ( lastq(1:dim_num), v(1:dim_num,l) )

  seed_save = seed
  seed = seed + 1

  return
end
subroutine i8_sobol_generate ( m, n, skip, r )

!*****************************************************************************80
!
!! I8_SOBOL_GENERATE generates a Sobol dataset.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    03 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N, the number of points to generate.
!
!    Input, integer ( kind = 8 ) SKIP, the number of initial points to skip.
!
!    Output, real ( kind = 8 ) R(M,N), the points.
!
  implicit none

  integer ( kind = 8 ) m
  integer ( kind = 8 ) n

  integer ( kind = 8 ) j
  real ( kind = 8 ), dimension ( m, n ) :: r
  integer ( kind = 8 ) seed
  integer ( kind = 8 ) skip

  do j = 1, n
    seed = skip + j - 1
    call i8_sobol ( m, seed, r(1:m,j) )
  end do

  return
end
subroutine i8_sobol_write ( m, n, skip, r, file_out_name )

!*****************************************************************************80
!
!! I8_SOBOL_WRITE writes a Sobol dataset to a file.
!
!  Discussion:
!
!    The initial lines of the file are comments, which begin with a
!    '#' character.
!
!    Thereafter, each line of the file contains the M-dimensional
!    components of the SKIP+I-1 entry of the Sobol sequence.
!
!    For the Sobol sequence, the value of SKIP is the same
!    as the value of SEED used to generate the first point.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    04 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 8 ) M, the spatial dimension.
!
!    Input, integer ( kind = 8 ) N, the number of (successive) points.
!
!    Input, integer ( kind = 8 ) SKIP, the number of skipped points.
!
!    Input, real ( kind = 8 ) R(M,N), the points.
!
!    Input, character ( len = * ) FILE_OUT_NAME, the name of
!    the output file.
!
  implicit none

  integer ( kind = 8 ) m
  integer ( kind = 8 ) n

  character ( len = * ) file_out_name
  integer file_out_unit
  integer ios
  integer ( kind = 8 ) j
  real ( kind = 8 ) r(m,n)
  integer ( kind = 8 ) skip
  character ( len = 40 ) string

  call get_unit ( file_out_unit )

  open ( unit = file_out_unit, file = file_out_name, status = 'replace', &
    iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'I8_SOBOL_WRITE - Fatal error!'
    write ( *, '(a)' ) '  Could not open the output file.'
    stop
  end if

  call timestring ( string )

  write ( file_out_unit, '(a)'       ) '#  ' // trim ( file_out_name )
  write ( file_out_unit, '(a)'       ) '#  created by I8_SOBOL_WRITE.F90.'
  write ( file_out_unit, '(a)'       ) '#'
  write ( file_out_unit, '(a)'       ) '#  File generated on ' &
    // trim ( string )
  write ( file_out_unit, '(a)'       ) '#'
  write ( file_out_unit, '(a,i8)'    ) '#  Spatial dimension M = ', m
  write ( file_out_unit, '(a,i8)'    ) '#  Number of points N = ', n
  write ( file_out_unit, '(a,g14.6)' ) '#  Epsilon (unit roundoff) = ', &
    epsilon ( r(1,1) )
  write ( file_out_unit, '(a,i8)'    ) '#  Initial values skipped = ', skip
  write ( file_out_unit, '(a)'       ) '#'

  write ( string, '(a,i3,a)' ) '(', m, '(2x,f10.6))'
  do j = 1, n
    write ( file_out_unit, string ) r(1:m,j)
  end do

  close ( unit = file_out_unit )

  return
end
function i8_xor ( i, j )

!*****************************************************************************80
!
!! I8_XOR calculates the exclusive OR of two integers.
!
!  Discussion:
!
!    This function is NOT needed in FORTRAN90, which supplies the
!    intrinsic IEOR function for this purpose.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    16 February 2005
!
!  Author:
!
!   John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 8 ) I, J, two values whose exclusive OR is needed.
!
!    Output, integer ( kind = 8 ) I8_XOR, the exclusive OR of I and J.
!
  implicit none

  integer ( kind = 8 ) i
  integer ( kind = 8 ) i1
  integer ( kind = 8 ) i2
  integer ( kind = 8 ) i8_xor
  integer ( kind = 8 ) j
  integer ( kind = 8 ) j1
  integer ( kind = 8 ) j2
  integer ( kind = 8 ) k
  integer ( kind = 8 ) l

  i1 = i
  j1 = j
  k = 0
  l = 1

  do while ( i1 /= 0 .or. j1 /= 0 )

    i2 = i1 / 2
    j2 = j1 / 2

    if ( &
      ( ( i1 == 2 * i2 ) .and. ( j1 /= 2 * j2 ) ) .or. &
      ( ( i1 /= 2 * i2 ) .and. ( j1 == 2 * j2 ) ) ) then
      k = k + l
    end if

    i1 = i2
    j1 = j2
    l = 2 * l

  end do

  i8_xor = k

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
subroutine timestring ( string )

!*****************************************************************************80
!
!! TIMESTRING writes the current YMDHMS date into a string.
!
!  Example:
!
!    STRING = 'May 31 2001   9:45:54.872 AM'
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    15 March 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) STRING, contains the date information.
!    A character length of 40 should always be sufficient.
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = * ) string
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( string, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end