* README for zufall random number package
* ------ --- ------ ------ ------ -------
* This package contains a portable random number generator set
* for: uniform (u in [0,1)), normal (<g> = 0, <g^2> = 1), and
* Poisson distributions. The basic module, the uniform generator,
* uses a lagged Fibonacci series generator:
*
* t = u(n-273) + u(n-607)
* u(n) = t - float(int(t))
*
* where each number generated, u(k), is floating point. Since
* the numbers are floating point, the left end boundary of the
* range contains zero. This package is nearly portable except
* for the following. (1) It is written in lower case, (2) the
* test package contains a timer (second) which is not portable,
* and (3) there are cycle times (in seconds) in data statements
* for NEC SX-3, Fujitsu VP2200, and Cray Y-MP. Select your
* favorite and comment out the others. Replacement functions
* for 'second' are included - comment out the others. Otherwise
* the package is portable and returns the same set of floating
* point numbers up to word precision on any machine. There are
* compiler directives ($cdir for Cray, *vdir for SX-3, and VOCL
* for Fujitsu VP2200) which should be otherwise ignored.
*
* To compile this beast, note that all floating point numbers
* are declared 'double precision'. On Cray X-MP, Y-MP, and C-90
* machines, use the cft77 (cf77) option -dp to run this in 64
* bit mode (not 128 bit double).
*
* External documentation, "Lagged Fibonacci Random Number Generators
* for the NEC SX-3," is to be published in the International
* Journal of High Speed Computing (1994). Otherwise, ask the
* author:
*
* W. P. Petersen
* IPS, RZ F-5
* ETHZ
* CH 8092, Zurich
* Switzerland
*
* e-mail: wpp@ips.ethz.ch.
*
* The package contains the following routines:
*
* ------------------------------------------------------
* UNIFORM generator routines:
*
* subroutine zufalli(seed)
* integer seed
* c initializes common block containing seeds. if seed=0,
* c the default value is 1802.
*
* subroutine zufall(n,u)
* integer n
* double precision u(n)
* c returns set of n uniforms u(1), ..., u(n).
*
* subroutine zufallsv(zusave)
* double precision zusave(608)
* c saves buffer and pointer in zusave, for later restarts
*
* subroutine zufallrs(zusave)
* double precision zusave(608)
* c restores seed buffer and pointer from zusave
* ------------------------------------------------------
*
* NORMAL generator routines:
*
* subroutine normalen(n,g)
* integer n
* double precision g(n)
* c returns set of n normals g(1), ..., g(n) such that
* c mean <g> = 0, and variance <g**2> = 1.
*
* subroutine normalsv(normsv)
* double precision normsv(1634)
* c saves zufall seed buffer and pointer in normsv
* c buffer/pointer for normalen restart also in normsv
*
* subroutine normalrs(normsv)
* double precision normsv(1634)
* c restores zufall seed buffer/pointer and
* c buffer/pointer for normalen restart from normsv
* ------------------------------------------------------
*
* POISSON generator routine:
*
* subroutine fische(n,mu,q)
* integer n,q(n)
* double precision mu
* c returns set of n integers q, with poisson
* c distribution, density p(q,mu) = exp(-mu) mu**q/q!
* c
* c USE zufallsv and zufallrs for stop/restart sequence
* c
program fibotest
implicit none
integer n
parameter(n=5000)
c
c this program tests three random number generators:
c------------------------------------------------------------
c IMPORTANT!! be sure to compile with -dp option on
c the Cray cft77 compiler. This disables 128-bit double
c precision, converting it to 64-bit single. Likewise,
c comment out the 'second' function. On other machines,
c uncomment lines for appropriate version of 'second'.
c------------------------------------------------------------
c zufallt = a uniform distribution of r.v.'s test
c
c normalt = a gaussian distribution of r. v.'s test
c
c fischet = poisson distribution of random numbers test
c
double precision a(n)
double precision t1,t2,second
integer p(n)
integer seed
c
c initialize the seeds for the uniform random number generator
c
seed = 0
t1 = second()
call zufalli(seed)
t2 = second()
t1 = t2 - t1
print 100
print 200
print 300,t1
c
c first do the uniform distribution test
c
print 400
call zufallt(n,a)
c
c next, the box-muller gaussian generator test
c
print 500
call normalt(n,a)
c
c finally, the poisson distribution test yields integers
c
print 600
call fischet(n,p)
c
print 700
c
* 15x,'===== BEGIN TESTS ====='/,
* 15x,24('='))
* 1x,'***** Initialization of zufall *****'/,
* 1x,36('-'))
* 1x,'***** Uniform distribution test *****'/,
* 1x,37('-')/)
* 1x,'***** Gaussian distribution test *****'/,
* 1x,37('-')/)
* 1x,'***** Poisson distribution test *****'/,
* 1x,37('-')/)
* 15x,'===== END OF TESTS ====='/,
* 15x,24('=')/)
stop
end
c
subroutine zufallt(n,a)
implicit none
integer n
c
double precision a(n)
integer ia(20)
double precision diff,t0,t1,t2,t3,second
double precision svblk(608)
double precision b(607),buff(607)
double precision CYCLE
integer i,ii,k,nits,ptr
common /klotz0/buff,ptr
c
c clock cycle for machine
c
c cycle for SX-3:
c data CYCLE/2.9E-9/
c cycle for Y-MP:
c data CYCLE/6.0E-9/
c cycle for VP2200
c data CYCLE/3.2E-9/
c cycle for SGI Indigo
c data CYCLE/1.0E-8/
c cycle for Sparc 51
data CYCLE/2.0E-8/
c
c number of iterations of test: nits
c
nits = 128
c
do 1 k=1,20
ia(k) = 0
t0 = 100.
do 2 k=1,nits
t1 = second()
t2 = second()
t1 = t2 - t1
t0 = min(t0,t1)
t1 = 100.
do 3 k=1,nits
t2 = second()
call zufall(n,a)
t3 = second()
t2 = t3 - t2
t1 = min(t2,t1)
do 4 i=1,n
ii = int(a(i)*20.)+1
ia(ii) = ia(ii) + 1
c
c last time, save klotz0 for save/resore test
c
if(k.eq.nits-1)then
call zufallsv(svblk)
endif
c
c
c test save/restore sequence
c
call zufallrs(svblk)
call zufall(607,b)
diff = 0.
do 5 i=1,min(n,607)
diff = diff + abs(b(i) - a(i))
if(diff.ne.0.) then
print *,' ERROR in start/restart: diff = ',diff
else
print *,' zufall save/restore test OK'
endif
c
t1 = (t1 - t0)/float(n)
print 100,t1
print 200,t1/CYCLE
print 300,(k,ia(k),k=1,20)
* 1x,' ------- ---------',/,
* 20(/1x,' bin(',i2,') = ',i9))
return
end
c
subroutine normalt(n,x)
implicit none
integer n
double precision x(n),y(128),boxsv(1634)
double precision x1,x2,x3,x4,x5,x6,xx2,xx4
double precision diff,t0,t1,t2,t3,second
integer bin(21),i,k,kk,nits
c
double precision CYCLE
c
c clock cycle for machine
c
c cycle for SX-3:
c data CYCLE/2.9E-9/
c cycle for Y-MP:
c data CYCLE/6.0E-9/
c cycle for VP2200
c data CYCLE/3.2E-9/
c cycle for SGI Indigo
c data CYCLE/1.0E-8/
c cycle for Sparc 51
data CYCLE/2.0E-8/
c
c number of iterations of test
c
nits = 128
c
c initialize moments
c
x1 = 0.
x2 = 0.
x3 = 0.
x4 = 0.
x5 = 0.
x6 = 0.
c
call normalen(n,x)
do 1 i = 1,21
bin(i) = 0
c
t0 = 10.
do 2 k = 1,nits
t1 = second()
t2 = second()
t1 = t2 - t1
t0 = min(t1,t0)
t1 = 100.
do 3 k = 1,nits
c
c save seeds and pointers for save/restore test
c
if(k.eq.nits) call normalsv(boxsv)
c
t2 = second()
call normalen(n,x)
t3 = second()
t2 = t3 - t2
t1 = min(t1,t2)
c
do 4 i=1,n
kk = int(2.0*(x(i)+5.25)) + 1
bin(kk) = bin(kk) + 1
do 5 i=1,n
x1 = x1 + x(i)
xx2 = x(i)*x(i)
x2 = x2 + xx2
x3 = x3 + xx2*x(i)
xx4 = xx2*xx2
x4 = x4 + xx4
x5 = x5 + xx4*x(i)
x6 = x6 + xx4*xx2
c
c restore previous seeds and pointers for save/restore test
c
if(k.eq.nits) then
call normalrs(boxsv)
call normalen(128,y)
endif
c
c
c save/restore check:
c
do 6 i=1,128
diff = diff + abs(y(i) - x(i))
if(diff.ne.0.) then
print *,' ERROR in normalsv/normalrs: diff = ',diff
else
print *,' normalen save/restore test OK'
endif
c
x1 = x1/float(n*nits)
x2 = x2/float(n*nits)
x3 = x3/float(n*nits)
x4 = x4/float(n*nits)
x5 = x5/float(n*nits)
x6 = x6/float(n*nits)
c
t1 = (t1 - t0)/float(n)
print 100,t1
print 200,t1/CYCLE
print 300,x1,x2,x3,x4,x5,x6
print 400
do 7 k=1,21
print 500,k,bin(k)
c
* 4x,'Compare to: (0.0)',18x,'(1.0)'/
* 4x,' <x> = ',e12.5,', <x**2> = ',e12.5,//
* 4x,'Compare to: (0.0)',18x,'(3.0)'/
* 4x,' <x**3> = ',e12.5,', <x**4> = ',e12.5,//
* 4x,'Compare to: (0.0)',17x,'(15.0)'/
* 4x,' <x**5> = ',e12.5,', <x**6> = ',e12.5)
* 1x,' --------- -- -------- ------------'/)
c
return
end
c
subroutine fischet(n,p)
implicit none
double precision mu,fp
double precision p1,p2,p3,p4
double precision x1,x2,x3,x4
double precision t0,t1,t2,t3,second
integer bin(20)
integer i,k,kk,n,nits
integer p(n)
c
double precision CYCLE
c
c clock cycle for machine
c
c cycle for SX-3:
c data CYCLE/2.9E-9/
c cycle for Y-MP:
c data CYCLE/6.0E-9/
c cycle for VP2200
c data CYCLE/3.2E-9/
c cycle for SGI Indigo
c data CYCLE/1.0E-8/
c cycle for Sparc 51
data CYCLE/2.0E-8/
c
mu = 2.0
nits = 128
c
do 1 k=1,20
bin(k) = 0
c
c moment comparison values
c
p1 = mu
p2 = mu + mu*mu
p3 = mu + 3.*mu*mu + mu*mu*mu
p4 = mu + 7.*mu*mu + 6.*mu*mu*mu + mu**4
c
x1 = 0.
x2 = 0.
x3 = 0.
x4 = 0.
c
t0 = 10.
do 2 k=1,nits
t1 = second()
t2 = second()
t1 = t2 - t1
t0 = min(t0,t1)
t1 = 10.
do 3 k=1,nits
c
t2 = second()
call fische(n,mu,p)
t3 = second()
t2 = t3 - t2
t1 = min(t1,t2)
c
do 4 i=1,n
kk = p(i)+1
bin(kk) = bin(kk) + 1
c
do 5 i=1,n
fp = float(p(i))
x1 = x1 + fp
x2 = x2 + fp*fp
x3 = x3 + fp*fp*fp
x4 = x4 + fp*fp*fp*fp
c
c
x1 = x1/float(n*nits)
x2 = x2/float(n*nits)
x3 = x3/float(n*nits)
x4 = x4/float(n*nits)
c
t1 = (t1 - t0)/float(n)
print 100,t1
print 200,t1/CYCLE
print 300,p1,p2,x1,x2,p3,p4,x3,x4
print 400,mu
do 6 k=1,20
print 500,k,bin(k)
c
* 3x,'Compare: (',e12.5,') (',e12.5,')'/
* 3x,' <p> = ',e12.5,', <p**2> = ',e12.5,//
* 3x,'Compare: (',e12.5,') (',e12.5,')'/
* 3x,' <p**3> = ',e12.5,', <p**4> = ',e12.5/)
* 1x,' --------- -- ------- ------------'/)
return
end
c
double precision function second()
c
c portable version of Cray function second(), comment out
c entire function for Y-MP. For SX-3, Fujitsu VP2200,
c or generic Unix, uncomment the version you want:
c
c NEC SX-3 version
c double precision xx(2)
c call clock(xx)
c
c VP2200 version
c double precision xx(2)
c call clockv(xx(2),xx(1),0,2)
c
c Generic Unix version
real xx(2)
call etime(xx)
c
second = xx(1)
c
return
end
c
c ---------------- end of test programs -------------
c
subroutine zufall(n,a)
implicit none
c
c portable lagged Fibonacci series uniform random number
c generator with "lags" -273 und -607:
c
c t = u(i-273)+buff(i-607) (floating pt.)
c u(i) = t - float(int(t))
c
c W.P. Petersen, IPS, ETH Zuerich, 19 Mar. 92
c
double precision a(*)
double precision buff(607)
double precision t
integer i,k,ptr,VL,k273,k607
integer buffsz,nn,n,left,q,qq
integer aptr,aptr0,bptr
c
common /klotz0/buff,ptr
data buffsz/607/
c
aptr = 0
nn = n
c
c
if(nn .le. 0) return
c
c factor nn = q*607 + r
c
q = (nn-1)/607
left = buffsz - ptr
c
if(q .le. 1) then
c
c only one or fewer full segments
c
if(nn .lt. left) then
do 2 i=1,nn
a(i+aptr) = buff(ptr+i)
ptr = ptr + nn
return
else
do 3 i=1,left
a(i+aptr) = buff(ptr+i)
ptr = 0
aptr = aptr + left
nn = nn - left
c buff -> buff case
VL = 273
k273 = 334
k607 = 0
do 4 k=1,3
cdir$ ivdep
*vdir nodep
*VOCL LOOP, TEMP(t), NOVREC(buff)
do 5 i=1,VL
t = buff(k273+i) + buff(k607+i)
buff(k607+i) = t - float(int(t))
k607 = k607 + VL
k273 = k273 + VL
VL = 167
if(k.eq.1) k273 = 0
c
goto 1
endif
else
c
c more than 1 full segment
c
do 6 i=1,left
a(i+aptr) = buff(ptr+i)
nn = nn - left
ptr = 0
aptr = aptr+left
c
c buff -> a(aptr0)
c
VL = 273
k273 = 334
k607 = 0
do 7 k=1,3
if(k.eq.1)then
*VOCL LOOP, TEMP(t)
do 8 i=1,VL
t = buff(k273+i) + buff(k607+i)
a(aptr+i) = t - float(int(t))
k273 = aptr
k607 = k607 + VL
aptr = aptr + VL
VL = 167
else
cdir$ ivdep
*vdir nodep
*VOCL LOOP, TEMP(t)
do 9 i=1,VL
t = a(k273+i) + buff(k607+i)
a(aptr+i) = t - float(int(t))
k607 = k607 + VL
k273 = k273 + VL
aptr = aptr + VL
endif
nn = nn - 607
c
c a(aptr-607) -> a(aptr) for last of the q-1 segments
c
aptr0 = aptr - 607
VL = 607
c
*vdir novector
do 10 qq=1,q-2
k273 = 334 + aptr0
cdir$ ivdep
*vdir nodep
*VOCL LOOP, TEMP(t), NOVREC(a)
do 11 i=1,VL
t = a(k273+i) + a(aptr0+i)
a(aptr+i) = t - float(int(t))
nn = nn - 607
aptr = aptr + VL
aptr0 = aptr0 + VL
c
c a(aptr0) -> buff, last segment before residual
c
VL = 273
k273 = 334 + aptr0
k607 = aptr0
bptr = 0
do 12 k=1,3
if(k.eq.1) then
*VOCL LOOP, TEMP(t)
do 13 i=1,VL
t = a(k273+i) + a(k607+i)
buff(bptr+i) = t - float(int(t))
k273 = 0
k607 = k607 + VL
bptr = bptr + VL
VL = 167
else
cdir$ ivdep
*vdir nodep
*VOCL LOOP, TEMP(t), NOVREC(buff)
do 14 i=1,VL
t = buff(k273+i) + a(k607+i)
buff(bptr+i) = t - float(int(t))
k607 = k607 + VL
k273 = k273 + VL
bptr = bptr + VL
endif
goto 1
endif
end
c
subroutine zufalli(seed)
implicit none
c
c generates initial seed buffer by linear congruential
c method. Taken from Marsaglia, FSU report FSU-SCRI-87-50
c variable seed should be 0 < seed <31328
c
integer seed
integer ptr
double precision s,t
double precision buff(607)
integer ij,kl,i,ii,j,jj,k,l,m
common /klotz0/buff,ptr
data ij/1802/,kl/9373/
c
if(seed.ne.0) ij = seed
c
i = mod(ij/177,177) + 2
j = mod(ij,177) + 2
k = mod(kl/169,178) + 1
l = mod(kl,169)
do 1 ii=1,607
s = 0.0
t = 0.5
do 2 jj=1,24
m = mod(mod(i*j,179)*k,179)
i = j
j = k
k = m
l = mod(53*l+1,169)
if(mod(l*m,64).ge.32) s = s+t
t = .5*t
buff(ii) = s
return
end
c
subroutine zufallsv(svblk)
implicit none
c
c saves common blocks klotz0, containing seeds and
c pointer to position in seed block. IMPORTANT: svblk must be
c dimensioned at least 608 in driver. The entire contents
c of klotz0 (pointer in buff, and buff) must be saved.
c
double precision buff(607)
integer ptr,i
double precision svblk(*)
common /klotz0/buff,ptr
c
svblk(1) = ptr
do 1 i=1,607
svblk(i+1) = buff(i)
c
return
end
subroutine zufallrs(svblk)
implicit none
c
c restores common block klotz0, containing seeds and pointer
c to position in seed block. IMPORTANT: svblk must be
c dimensioned at least 608 in driver. The entire contents
c of klotz0 must be restored.
c
double precision buff(607)
integer i,ptr
double precision svblk(*)
common /klotz0/buff,ptr
c
ptr = svblk(1)
do 1 i=1,607
buff(i) = svblk(i+1)
c
return
end
subroutine normalen(n,x)
implicit none
c
c Box-Muller method for Gaussian random numbers
c
double precision x(*)
double precision xbuff(1024)
integer i,ptr,xptr,first
integer buffsz,nn,n,left
common /klotz1/xbuff,first,xptr
data buffsz/1024/
c
nn = n
if(nn .le. 0) return
if(first.eq.0)then
call normal00
first = 1
endif
ptr = 0
c
left = buffsz - xptr
if(nn .lt. left) then
do 2 i=1,nn
x(i+ptr) = xbuff(xptr+i)
xptr = xptr + nn
return
else
do 3 i=1,left
x(i+ptr) = xbuff(xptr+i)
xptr = 0
ptr = ptr+left
nn = nn - left
call normal00
goto 1
endif
end
subroutine normal00
implicit none
double precision pi,twopi
parameter(pi=3.141592653589793)
double precision xbuff(1024),r1,r2,t1,t2
integer first,xptr,i
common /klotz1/xbuff,first,xptr
c
twopi = 2.*pi
call zufall(1024,xbuff)
*VOCL LOOP, TEMP(r1,r2,t1,t2), NOVREC(xbuff)
do 1 i=1,1024,2
r1 = twopi*xbuff(i)
t1 = cos(r1)
t2 = sin(r1)
r2 = sqrt(-2.*log(1.-xbuff(i+1)))
xbuff(i) = t1*r2
xbuff(i+1) = t2*r2
c
return
end
subroutine normalsv(svbox)
implicit none
c
c saves common block klotz0 containing buffers
c and pointers. IMPORTANT: svbox must be dimensioned at
c least 1634 in driver. The entire contents of blocks
c klotz0 (via zufallsv) and klotz1 must be saved.
c
double precision buff(607)
integer i,k,ptr
double precision xbuff(1024)
integer xptr,first
double precision svbox(*)
common /klotz0/buff,ptr
common /klotz1/xbuff,first,xptr
c
if(first.eq.0)then
print *,' ERROR in normalsv, save of unitialized block'
endif
c
c save zufall block klotz0
c
call zufallsv(svbox)
c
svbox(609) = first
svbox(610) = xptr
k = 610
do 1 i=1,1024
svbox(i+k) = xbuff(i)
c
return
end
subroutine normalrs(svbox)
implicit none
c
c restores common blocks klotz0, klotz1 containing buffers
c and pointers. IMPORTANT: svbox must be dimensioned at
c least 1634 in driver. The entire contents
c of klotz0 and klotz1 must be restored.
c
double precision buff(607)
integer ptr
double precision xbuff(1024)
integer i,k,xptr,first
double precision svbox(*)
common /klotz0/buff,ptr
common /klotz1/xbuff,first,xptr
c
c restore zufall blocks klotz0 and klotz1
c
call zufallrs(svbox)
first = svbox(609)
if(first.eq.0)then
print *,' ERROR in normalsv, restoration of unitialized block'
endif
xptr = svbox(610)
k = 610
do 1 i=1,1024
xbuff(i) = svbox(i+k)
c
return
end
subroutine fische(n,mu,p)
implicit none
integer p(*)
integer indx(1024)
integer n,i,ii,jj,k,left,nl0,nsegs,p0
double precision u(1024),q(1024)
double precision q0,pmu,mu
c
c Poisson generator for distribution function of p's:
c
c q(mu,p) = exp(-mu) mu**p/p!
c
c initialize arrays, pointers
c
if (n.le.0) return
c
pmu = exp(-mu)
p0 = 0
c
nsegs = (n-1)/1024
left = n - nsegs*1024
nsegs = nsegs + 1
nl0 = left
c
do 2 k = 1,nsegs
c
do 3 i=1,left
indx(i) = i
p(p0+i) = 0
q(i) = 1.0
c
c Begin iterative loop on segment of p's
c
c
c Get the needed uniforms
c
call zufall(left,u)
c
jj = 0
c
cdir$ ivdep
*vdir nodep
*VOCL LOOP, TEMP(ii,q0), NOVREC(indx,p,q)
do 4 i=1,left
ii = indx(i)
q0 = q(ii)*u(i)
q(ii) = q0
if( q0.gt.pmu ) then
jj = jj + 1
indx(jj) = ii
p(p0+ii) = p(p0+ii) + 1
endif
c
c any left in this segment?
c
left = jj
if(left.gt.0)then
goto 1
endif
c
p0 = p0 + nl0
nl0 = 1024
left = 1024
c
c
return
end
c
block data
implicit none
c
c globally accessable, compile-time initialized data
c
integer ptr,xptr,first
double precision buff(607),xbuff(1024)
common /klotz0/buff,ptr
common /klotz1/xbuff,first,xptr
data ptr/0/,xptr/0/,first/0/
end
.
