when converting the FUNCTION ran (from page 1142 of NR in F90) to double precision, we find that approx 1 in every 1 billion numbers is 1.0 (exactly), which should not occur according to the code comments (exclusive of endpoints).

the code was run under win2000 (compaq fortran 6.6b) and linux (Fujitsu Fortran Compiler Driver Version 1.0 (Jun 11 1999 00:49:31)) with identical results (eg, using a seed of -460053). the behaviour does not occur under single precision.

we cannot see any limitation of the algorithm to single precision, so the behaviour seems odd. this causes problems when transforming to other deviates using the inverse cdf. We currently skip the number if it is 1.0 exactly. What causes this behaviour and is there a better fix?

example code below
------------------

module kinds
implicit none
integer,parameter::k4b=selected_int_kind(9), rk=kind(1.d0)
endmodule
!***********************************
program sample_subroutine
use kinds
implicit none
! locals
integer(k4b),parameter::n=1000000000, scale=1000000000
integer(k4b) :: i,idum,ii
real(rk)::randnum_nr,ran_nr
! Start procedure here
open(unit=10,file='bang.txt')
idum=-460053 ! cause BANG!
write(*,*) '***real_kind=',rk
write(10,*) '***real_kind=',rk
write(*,*) 'causes BANG:',idum
write(*,*) 'enter value for idum/seed:'
read(*,*) idum
write(*,*) 'idum=',idum
write(10,*) 'idum=',idum
ii=0; i=0
do
if(i>=scale)then
ii=ii+1
i=0
write(*,*)ii
write(10,*)ii
endif
i=i+1
randnum_nr=ran_nr(idum)
if(randnum_nr<=0._rk.or.randnum_nr>=1._rk)then
write(*,*)"BANG_NR",ii,i,randnum_nr
write(10,*)"BANG_NR",ii,i,randnum_nr
endif
enddo
! End main procedure here
endprogram sample_subroutine
!***********************************
FUNCTION ran_nr(idum)
use kinds
IMPLICIT NONE
INTEGER(K4B), INTENT(INOUT) :: idum
REAL(rk) :: ran_nr
INTEGER(K4B), PARAMETER :: IA=16807,IM=2147483647,IQ=127773,IR=2836
REAL(rk), SAVE :: am
INTEGER(K4B), SAVE :: ix=-1,iy=-1,k
if (idum <= 0 .or. iy < 0) then
am=nearest(1._rk,-1._rk)/IM
iy=ior(ieor(888889999,abs(idum)),1)
ix=ieor(777755555,abs(idum))
idum=abs(idum)+1
end if
ix=ieor(ix,ishft(ix,13))
ix=ieor(ix,ishft(ix,-17))
ix=ieor(ix,ishft(ix,5))
k=iy/IQ
iy=IA*(iy-k*IQ)-IR*k
if (iy < 0) iy=iy+IM
ran_nr=am*ior(iand(IM,ieor(ix,iy)),1)
END FUNCTION ran_nr
!***********************************

You can fix the code by changing the line

am=nearest(1._rk,-1._rk)/IM

to

am=nearest(1._rk/IM,-1._rk)

This doesn't affect the single precision code, and will take care of the problem in the double precision version.