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