! Simulation of aeroplane problem
!
! N seats on an aeroplane and N people to fill them, each with a unique
! ticket number. Passengers proceed one at a time to their seats. The
! first person ignores their ticket number and chooses a seat at random.
! All subsequent people sit at their seat if it is available, otherwise
! they choose a seat at random. What is the probability that the last person
! can sit in their allocated seat.
!
! This simulation program is written in Fortran. To build it you'll need a
! fortran compiler. If you haven't got one get g95 from http://www.g95.org
! - it's free and is available for most OSs. If you want high performance then
! you'll find the Intel compiler (free for personal use on Linux) much better -
! just replace g95 with ifort below.
!
! To build the program type
! g95 -o aeroplane -O3 aeroplane.f90
!
! To run the program with 100 seats type at a command prompt
! aeroplane
!
! If you want to try a different number (e.g. 10) type
! aeroplane 10
!
! # seats % availability % std.dev
! 2 50.035 0.4830
! 3 49.931 0.4910
! 4 50.022 0.4979
! 5 50.026 0.5639
! 6 50.023 0.5229
! 7 50.030 0.5292
! 8 49.962 0.4597
! 9 50.021 0.5095
! 10 50.081 0.4409
!
! The % availability of all seats up to the given number is printed out along
! with the standard-deviation. If you just want the result for a particular
! number of seats use a negative number
! aeroplane -10
!
! # seats % availability % std.dev
! 10 50.020 0.4789
!
! Written by Simon Geard, 9/5/2007
! simon@whiteowl.co.uk
!
! Please feel free to adapt/use/change this program as you like.
!
!==============================================================================
module simulate
public :: hasSeat, tidy
private
integer :: nseats = 0
logical, allocatable :: seatOccupied(:)
contains
subroutine tidy
if (allocated(seatOccupied)) deallocate(seatOccupied)
end subroutine tidy
logical function hasSeat(ns)
! Run the simulation once and return .true. if the final seat
! is available
implicit none
integer, intent(in) :: ns
integer :: c
integer :: i
nseats = ns
if (.not. allocated(seatOccupied)) then
allocate(seatOccupied(ns))
end if
seatOccupied = spread(.false.,1,nseats)
! Random choice of 1st passenger
c = selectSeat()
if (c == 1) then
hasSeat = .true. ! Everyone else will be able to sit in their allocated seat
return
end if
if (c == nseats) then
hasSeat = .false. ! Last seat now occupied
return
end if
seatOccupied(c) = .true.
! Other passengers
do i=2,nseats-1
if (seatOccupied(i)) then
! Seat occupied so choose another one at random
c = selectSeat()
if (c == nseats) then
hasSeat = .false. ! Last seat now occupied
return
end if
seatOccupied(c) = .true.
else
! Seat available
seatOccupied(i) = .true.
end if
end do
hasSeat = .not. seatOccupied(nseats) ! Seat available if not occupied
end function hasSeat
integer function selectSeat()
! Choose a random unoccupied seat
implicit none
real :: h
do
call random_number(h)
selectSeat = 1+nint((nseats-1)*h)
if (.not. seatOccupied(selectSeat)) exit
end do
end function selectSeat
end module simulate
program aeroplane
use simulate
implicit none
integer, parameter :: N = 10000 ! Number of trials in a simulation
integer, parameter :: L = 100 ! Number simulations
integer :: c ! Count of the number of times in a trial the seat
! was available
integer :: ns ! Number of seats in the trial
integer :: nargs ! Number of arguments on the command-line
character(len=10) :: arg ! Holder for the first argument
real*8 :: mean ! Mean of results
real*8 :: stdev ! Standard-deviation of results
real*8 :: results(L)! Results vector, 1 result for each simulation
integer :: i, j, k ! Dummy loop counters
! Get the number of seats from the command line (default to 100)
nargs = command_argument_count()
if (nargs >= 1) then
call get_command_argument(1,arg)
read(arg,*) ns
else
ns = 100
end if
write(*,'(a)') '# seats % availability % std.dev'
if (ns < 0) then
k = -ns
call runSimulation
else
do k=2,ns
call runSimulation
end do
end if
contains
subroutine runSimulation
! Do the L simulations
do j=1,L
c = 0
! Do the N trials
do i=1,N
if (hasSeat(k)) then
c = c+1
end if
end do
results(j) = dble(c)/N
end do
mean = sum(results)/L
stdev = sqrt(dot_product(results,results)/L - mean**2)
write(*,'(3X,i3,8X,f7.3,8X,f6.4)') k,mean*100,stdev*100
call tidy
end subroutine runSimulation
end program aeroplane