subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! get_unit() returns a free FORTRAN unit number.
!
!  Discussion:
!
!    A "free" FORTRAN unit number is a value 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 a value 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 MIT license.
!
!  Modified:
!
!    26 October 2008
!
!  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
subroutine high_card_probability ( n, p )

!*****************************************************************************80
!
!! HIGH_CARD_PROBABILITY: winning probabilities for the high card game.
!
!  Discussion:
!
!    The high card game presents the player with a deck of cards, each
!    having an unknown value.  The player is allowed to go throught the
!    deck once, looking at the cards one at a time.  At any time, the player
!    may decide to take a particular card, winning that amount and stopping
!    the game.  If the player continues to the end, by default the last card
!    indicates the amount won.
!
!    An optimal strategy for selecting the highest card is as follows:
!    * look at, but do not select, the first k-1 cards;
!    * stop at the first card, from k to n, that is higher than the first 
!      k-1 cards.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2014
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of cards.  
!
!    Output, real ( kind = rk ) P(N).  P(K) is the probability that a strategy 
!    that skips K-1 cards will win, given that the deck has N cards.
!
  implicit none

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

  integer n

  integer i
  integer j
  real ( kind = rk ) p(n)
  real ( kind = rk ) t

  do i = 1, n
    t = 0.0D+00
    do j = i, n - 1
      t = t + 1.0D+00 / real ( j, kind = rk )
    end do
    p(i) = ( 1.0D+00 + real ( i - 1, kind = rk ) * t ) / real ( n, kind = rk )
  end do

  return
end
subroutine high_card_shuffle ( n, seed, sequence )

!*****************************************************************************80
!
!! HIGH_CARD_SHUFFLE generates a sequence of numeric "cards" for a game.
!
!  Discussion:
!
!    In this game, you know that the deck contains N cards.  You win by 
!    choosing the highest card in the deck.  You don't know what this card
!    is, and you must choose your card by saying "stop" as, one by one,
!    the cards of the deck are exposed.  
!
!    A random guesser would get the high card with probability 1/N.
!
!    An intelligent guesser can do much better.
!
!    It is the goal of this program so "shuffle" a deck of cards suitable
!    for this game.  The problem is that we know the highest card in an
!    ordinary deck.  Let's replace the cards by integers.  Then if we know
!    in advance the range of the cards (say, they must lie between 1 and 
!    1,000), it may be true that we can guess the card that is the maximum.
!
!    However, this program produces a sequence of integer card values for
!    which no information can be gained from the values.  It does this
!    by regarding the card values as binary integers between 1 and 2^N - 1.
!    We can make a perfectly information-free sequence as follows:
!
!      Card 1 sets bit N-1 to 1.
!      Card 2 sets bit N-2 to 1, bit  N-1 randomly.
!      ...
!      Card I sets bit N-I to 1, bits N-1 down to N-I+1 randomly.
!      ...
!      Card N sets bit N-N to 1, bits N-1 down to 1 randomly.
!
!    The I-th card has equal probability to fall in any of the I intervals
!    defined by the I-1 previous cards.  So, knowing the previous cards tells
!    you absolutely nothing about where the next card will fall, and each
!    card is, at the moment you see it, as good a guess for the maximum as
!    the unseen cards.
!
!    For example, the command "high_card_shuffle(7)" might yield
!
!      64    96    80     8     4    82    29
!    or
!      64    32    16    24    12    58    73
!    or
!      64    96    48     8   100    26    93
!    or
!      64    96    16    56    76   122    71
!
!    in which the highest card is #2, #7, #5, or #6 respectively.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2014
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of cards.  N probably needs to 
!    be less than 32.
!
!    Input/output, integer SEED, a seed for the random
!    number generator.
!
!    Output, integer SEQUENCE(N), a set of N integer values 
!    that can be used as the cards in the high card guessing game.
!
  implicit none

  integer n

  integer c
  integer i
  integer i4_uniform_ab
  integer j
  integer k
  integer seed
  integer sequence(n)

  if ( 32 <= n ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HIGH_CARD_SHUFFLE - Fatal error!'
    write ( *, '(a)' ) '  This program can only handle N < 32.'
    stop 1
  end if

  do i = 1, n
    c = 2 ** ( n - i )
    do j = 1, i - 1
      k = i4_uniform_ab ( 0, 1, seed )
      c = c + k * 2 ** ( n - i + j )
    end do
    sequence(i) = c
  end do

  return
end
subroutine high_card_simulation ( deck_size, skip_num, trial_num, seed, p )

!*****************************************************************************80
!
!! HIGH_CARD_SIMULATION simulates a game of choosing the highest card in a deck.
!
!  Discussion:
!
!    You are given a deck of DECK_SIZE cards.
!
!    Your goal is to select the high card.  For convenience, we can assume 
!    the cards are a permutation of the integers from 1 to DECK_SIZE, but in
!    fact the user mustn't see such values or else it's obvious which is the
!    largest card.
!
!    However, your choice is made under the following rules:  You may turn over
!    one card at a time.  When a card is turned over, you may declare that to be
!    your choice, or else turn over another card.  If you have not chosen a card
!    by the end, then your choice is the final card.
!
!    If you have no idea what to do, and simply decide in advance to pick
!    a card "at random", that is, for example, you decide to pick the 15th card
!    before having seen any cards, then your probability of winning is 
!    1/DECK_SIZE.
!
!    The question is, can you do better than that?
!
!    Your strategy is as follows: always look at the first SKIP_NUM cards 
!    without choosing them.  Then choose the very next card you encounter 
!    that is larger than the cards you skipped.
!
!    Using this program, you can easily see that skipping 5 cards is much better
!    than picking one at random, skipping 10 is even better, and so on.  
!    Of course, you can't skip too many cards, and in fact, the results seem 
!    to be best for somewhere around 30 to 35 cards skipped.  For problems 
!    like this, the optimal value is somewhere around 1 / e, where E is the 
!    base of the natural logarithm system.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2014
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DECK_SIZE, the number of cards in the deck.
!    2 <= DECK_SIZE.  Default value is 52;
!
!    Input, integer SKIP_NUM, the number of initial cards you plan 
!    to examine but will NOT select.  If SKIP_NUM is 0, you don't look at any 
!    cards first.  0 <= SKIP_NUM < DECK_SIZE.
!
!    Input, integer TRIAL_NUM, the number of times we will 
!    simulate this process.
!
!    Input/output, integer SEED, a seed for the random 
!    number generator.
!
!    Output, real ( kind = rk ) P, the estimated probability that your strategy 
!    of skipping SKIP_NUM cards and then selecting the next card that is 
!    bigger, will result in choosing the highest card.
!
  implicit none

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

  integer deck_size

  integer card
  integer cards(deck_size)
  integer choice
  integer correct
  integer, parameter :: i4_huge = 2147483647
  real ( kind = rk ) p
  integer seed
  integer skip_max
  integer skip_num
  integer trial
  integer trial_num
  integer true_max
!
!  Check values.
!
  if ( deck_size < 2 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'HIGH_CARD_SIMULATION - Fatal error!'
    write ( *, '(a)' ) '  DECK_SIZE must be at least 2.'
    write ( *, '(a,i6)' ) '  Your value was ', deck_size
    stop 1
  end if

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

  if ( deck_size <= skip_num ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'HIGH_CARD_SIMULATION - Fatal error!'
    write ( *, '(a)' ) '  SKIP_NUM must be less than DECK_SIZE.'
    write ( *, '(a,i6)' ) '  Your DECK_SIZE = ', deck_size
    write ( *, '(a,i6)' ) '  Your SKIP_NUM = ', skip_num
    stop 1
  end if

  if ( trial_num < 1 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'HIGH_CARD_SIMULATION - Fatal error!'
    write ( *, '(a)' ) '  TRIAL_NUM must be at least 1.'
    write ( *, '(a,i6)' ) '  Your TRIAL_NUM was ', trial_num 
    stop 1
  end if

  correct = 0

  do trial = 1, trial_num

    call permutation_random ( deck_size, seed, cards )

    if ( 1 <= skip_num ) then
      skip_max = maxval ( cards(1:skip_num) )
    else
      skip_max = - i4_huge
    end if

    true_max = maxval ( cards(1:deck_size) )
!
!  In case you don't encounter a card larger than SKIP_MAX,
!  we'll assume you pick the last card in the deck, even though
!  you know it's a loser.
!
    choice = cards(deck_size)
!
!  Turn over the remaining cards in the deck, but stop
!  immediately when you find one bigger than SKIP_MAX.
!
    do card = skip_num + 1, deck_size
      if ( skip_max < cards(card) ) then
        choice = cards(card)
        exit
      end if
    end do
!
!  Record successful choices.
!
    if ( choice == true_max ) then
      correct = correct + 1
    end if

  end do
!
!  Estimate the probability.
!
  p = real ( correct, kind = rk ) / real ( trial_num, kind = rk )

  return
end
function i4_uniform_ab ( a, b, seed )

!*****************************************************************************80
!
!! I4_UNIFORM_AB returns a scaled pseudorandom I4 between A and B.
!
!  Discussion:
!
!    An I4 is an integer value.
!
!    The pseudorandom number will be scaled to be uniformly distributed
!    between A and B.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    02 October 2012
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Second Edition,
!    Springer, 1987,
!    ISBN: 0387964673,
!    LC: QA76.9.C65.B73.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, December 1986, pages 362-376.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley, 1998,
!    ISBN: 0471134031,
!    LC: T57.62.H37.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, Number 2, 1969, pages 136-143.
!
!  Parameters:
!
!    Input, integer A, B, the limits of the interval.
!
!    Input/output, integer SEED, the "seed" value, which
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, integer I4_UNIFORM_AB, a number between A and B.
!
  implicit none

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

  integer a
  integer b
  integer, parameter :: i4_huge = 2147483647
  integer i4_uniform_ab
  integer k
  real ( kind = rk ) r
  integer seed
  integer value

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'I4_UNIFORM_AB - 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 + i4_huge
  end if

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

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

  i4_uniform_ab = value

  return
end
subroutine i4vec_print ( n, a, title )

!*****************************************************************************80
!
!! I4VEC_PRINT prints an I4VEC.
!
!  Discussion:
!
!    An I4VEC is a vector of I4's.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    02 May 2010
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, integer A(N), the vector to be printed.
!
!    Input, character ( len = * ) TITLE, a title.
!
  implicit none

  integer n

  integer a(n)
  integer i
  character ( len = * ) title

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( title )
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,i8,a,2x,i12)' ) i, ':', a(i)
  end do

  return
end
subroutine permutation_random ( n, seed, p )

!*****************************************************************************80
!
!! PERMUTATION_RANDOM returns a random permutation.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2014
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of objects to permute.
!
!    Input, integer SEED, a seed for the random 
!    number generator.
!
!    Output, integer P(N), a permutation of the integers 
!    from 1 to N.
!
  implicit none

  integer n

  integer i
  integer i4_uniform_ab
  integer k
  integer p(n)
  integer p1
  integer seed

  do i = 1, n
    p(i) = i
  end do

  do i = 1, n - 1

    k = i4_uniform_ab ( i, n, seed )

    p1   = p(i)
    p(i) = p(k)
    p(k) = p1

  end do

  return
end
subroutine timestamp ( )

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

  character ( len = 8 ) ampm
  integer d
  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
  integer values(8)
  integer y

  call date_and_time ( values = 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 ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end