/
durand.f
42 lines (41 loc) · 1.15 KB
/
durand.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
FUNCTION ran1(idum)
INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
DOUBLE PRECISION ran1,AM,EPS,RNMX
PARAMETER (IA=16807,IM=2147483647,AM=1.d0/IM,IQ=127773,IR=2836,
*NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=3.d-16,RNMX=1.d0-EPS)
INTEGER j,k,iv(NTAB),iy
SAVE iv,iy
DATA iv /NTAB*0/, iy /0/
if (idum.le.0.or.iy.eq.0) then
idum=max(-idum,1)
do 11 j=NTAB+8,1,-1
k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k
if (idum.lt.0) idum=idum+IM
if (j.le.NTAB) iv(j)=idum
11 continue
iy=iv(1)
endif
k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k
if (idum.lt.0) idum=idum+IM
j=1+iy/NDIV
iy=iv(j)
iv(j)=idum
ran1=min(AM*iy,RNMX)
return
END
subroutine durand(seed, npts, x)
implicit none
integer npts, i, idum
real*8 seed, ran1, x(npts)
if (seed .lt. 0.0d0 .or. seed .gt. 2147483648.0) then
write(6,*) 'seed must be a positive integer < 2.1474*10**9'
seed = 135791113.0
end if
idum = -int(seed)
do i = 1, npts
x(i) = ran1(idum)
end do
seed = dble(idum)
end