On Tue, 16 Sep 2008 20:45:39 -0700 (PDT), e p chandler posted:
>> integer(i13), parameter:: i_max = 63
>> integer(i13), parameter:: i_tot = i_max*2 +1
>
>> integer(i13), parameter:: sides = i_tot
>> integer(i13), parameter:: trials = 10
>>
>> integer(i13), dimension(trials)::C
>> integer(i13), dimension(trials)::F
>>
>> integer(i13):: b, i, &
>> Â Â Â Â Â Â diff2, c3, c4
>> real:: harvest, diff, t3, t4
>>
>> ! seed random num generator
> [snip]
>>
>> do i = 1, trials
>>
>> Â Â call random_number(harvest)
>> Â Â b = min(max(1,ceiling(harvest*sides)),sides)
>>
>> C(i) = b
>>
>> end do
>> print *, 'C ='
>> print *, C
>>
>> print *, 'F ='
>>
>> do i = 1, trials
>>
>> F(i) = C(i) - i_max
>> end do
>> print *, F
>
> [snip]
>
>> does this prog give
>> equiprobable integers on the interval of [-i_max , i_max] ?
>
> No. It looks like you have an off by one bug. Consider this test
> program:
>
> C:\Users\epc\temp>type rf.f90
> implicit none
>
> integer,parameter :: i_max = 63
> integer,parameter :: i_tot = 2 * i_max + 1
> integer,parameter :: trials = 10 * i_tot
>
> integer k(-i_max:i_max)
> integer b, i
>
> real harvest
>
> k = 0
>
> do i=1,trials
> call random_number(harvest)
> b = -i_max + int( harvest * i_tot )
> k(b) = k(b) + 1
> end do
>
> print *,k
>
> end
>
> C:\Users\epc\temp>g95 rf.f90
>
> C:\Users\epc\temp>a
> 11 9 9 8 12 8 10 10 11 10 13 13 10 6 4 9 13 12 6 10 6 12 9 11 10 8 16
> 8 6 13 10
> 7 5 7 10 9 6 11 15 9 16 12 7 10 11 9 9 11 6 18 8 8 11 12 11 12 9 15 9
> 14 10 6 1
> 0 10 7 12 12 15 11 10 9 10 8 11 14 8 11 12 11 11 9 9 7 10 6 8 10 16 14
> 7 9 2 12
> 8 11 11 11 14 5 12 8 4 12 10 11 10 12 11 10 7 9 11 14 8 11 9 7 12 13
> 12 7 8 10 1
> 0 15 11 9
>
> int( harvest * i_tot ) gives an integer from 0 to i_tot-1
> OR from 0 to 2*i_max
>
> subtracting i_max gives an integer from -i_max to i_max
>
> This is one traditional way of picking a random integer over a range.
>
> 1 + int ( 6 * harvest), for example, is the classical way of rolling a
> die.
>
> - e
I *do* have off by one bug, but I think we're both wrong at this point, as
your classical way to roll a die has uneven outcomes.
The big step forward for this endeavor though was
integer(i13) k(-i_max:i_max)
, which allows me to count outcomes.
n= 1
clock= 269500907
seed= 269500907
1 0.224650
2 0.820390
3 0.150279
4 0.285857
C =
F =
k is 0 347 353
354 398
361 384 376
379 402
355 379 381
343 388
371 403 340
341 383
390 355 370
387 411
366 319
system clock was 0.00000
cpu time was 0
Press RETURN to close window . . .
k(-i_max) is clearly not getting any business. I thought I had the nuts on
how to roll a die in fortran, but I'm a litle perplexed here.
implicit none
integer, parameter :: i13 = selected_int_kind(13)
integer(i13), parameter:: i_max = 13
integer(i13), parameter:: i_tot = i_max*2 +1
! integer(i13), parameter:: sides = 8
integer(i13), parameter:: sides = i_tot
integer(i13), parameter:: trials = 10000
integer(i13), dimension(trials)::C
integer(i13), dimension(trials)::F
integer(i13) k(-i_max:i_max)
integer(i13):: b, i, &
diff2, c3, c4
real:: harvest, diff, t3, t4
! 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
do i = 1, trials
call random_number(harvest)
b = min(max(1,ceiling(harvest*sides)),sides)
C(i) = b
end do
print *, 'C ='
! print *, C
print *, 'F ='
do i = 1, trials
F(i) = C(i) - i_max
end do
! print *, F
! count outcomes
k = 0
do i = 1, trials
k(F(i)) = k(F(i)) + 1
end do
print *, " k is ", k
??
--
Opera in English is, in the main, just about as sensible as baseball in
Italian.
H. L. Mencken
