Re: random ints on symmetric interval
  Home FAQ Contact Sign in
comp.lang.fortran only
 
Advanced search
POPULAR GROUPS

more...

 Up
Re: random ints on symmetric interval         

Group: comp.lang.fortran · Group Profile
Author: Ron Ford
Date: Sep 18, 2008 16:02

On Wed, 17 Sep 2008 06:38:36 -0700 (PDT), e p chandler posted:
> On Sep 17, 2:52 am, Ron Ford wrote:
>
>> Why do I believe that
>> b = min(max(1,ceiling(harvest*sides)),sides)
>> is better than
>>   b = -i_max + int( harvest * i_tot )
>
> That's a matter of "taste". :-). For each of us our choice illustrates
> our own thought process.

I think your taste is more legible than mine. Unsurprisingly, Glen is
right about the ultimate calculation here:

! calculates percent of UB for b-a on a symmetric interval
implicit none
integer, parameter :: i13 = selected_int_kind(13)
integer(i13), parameter:: i_max = 32767
integer(i13), parameter:: i_tot = i_max*2 +1
integer(i13), parameter:: sides = i_tot
integer(i13), parameter:: trials = 1000000
integer(i13) k(-i_max:i_max)
integer(i13):: a, b, i, g, inbounds, outbounds

real:: harvest, ratio, percent

! seed random num generator

CALL init_random_seed

! prime the pump
call random_number(harvest)
b = 3 + nint(10*harvest)
do i=1,b
call random_number(harvest)
print *, i, harvest
end do

!main control
inbounds = 0
outbounds = 0

do i = 1, trials

call random_number(harvest)
b = -i_max + int( harvest * i_tot )
call random_number(harvest)
a = -i_max + int( harvest * i_tot )
g= b-a
! count outcomes
if (abs (g).le.i_max) then
inbounds = inbounds + 1
else
outbounds =outbounds + 1
end if

end do

ratio = real(inbounds)/ real(outbounds)
percent = real(outbounds)/real(trials)
print *, "inbounds = ", inbounds
print *, "outbounds = ", outbounds
print *, "ratio is ", ratio
print *, "percent =" , 100.0 *percent, " %%"

contains
SUBROUTINE init_random_seed()
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed

CALL RANDOM_SEED(size = n)
print *, "n=", n
ALLOCATE(seed(n))

CALL SYSTEM_CLOCK(COUNT=clock)
print *, "clock=", clock

seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)
print *, "seed= ", seed

DEALLOCATE(seed)
END SUBROUTINE
end program
! gfortran -o pop -Wall freeformat59.f95

n= 1
clock= 30134791
seed= 30134791
1 0.197780
2 0.764005
3 0.471637
4 0.507431
5 0.156869
6 0.739084
7 0.341144
inbounds = 749970
outbounds = 250030
ratio is 2.99952
percent = 25.0030 %%

Press RETURN to close window . . .

I wonder how big the interval can get before the numbers start to get
"grainy" with b = -i_max + int( harvest * i_tot )?

How would a person stich together 2 pseudorandoms?

--
A Sunday school is a prison in which children do penance for the evil
conscience of their parents.
H. L. Mencken
no comments
diggit! del.icio.us! reddit!