double precision function ranu () c The function ranu returns a pseudo-random number from a uniform dis- c tribution between 0 and 1 exclusive. The subroutine ranset (see be- c low) should be called (once) before the first time this function is c called. implicit none integer seed1, seed2, iz, k double precision minv parameter (minv=1.d0/2147483563.d0) common /random/seed1, seed2 save /random/ c data seed1,seed2 /12345,67890/ k = seed1/53668 seed1 = 40014*(seed1-k*53668)-k*12211 if (seed1 .lt. 0) seed1 = seed1 + 2147483563 k = seed2/52774 seed2 = 40692*(seed2-k*52774)-k*3791 if (seed2 .lt. 0) seed2 = seed2 + 2147483399 iz = seed1 - seed2 if(iz .lt. 1) iz = iz + 2147483562 ranu = iz * minv return end c******************************************************************** subroutine ranuv (n, vec) c The subroutine ranuv generates a vector of n pseudo-random numbers, c each of which is from a uniform distribution between 0 and 1 exclu- c sive, storing this vector in the parameter vec. The subroutine c ranset (see below) should be called (once) before the first time c this subroutine is called. implicit none integer seed1, seed2, iz, n, i, k double precision minv, vec(n) parameter (minv=1.d0/2147483563.d0) common /random/seed1, seed2 save /random/ data seed1,seed2 /12345,67890/ if (n .lt. 1) return do 10 i = 1, n k = seed1/53668 seed1 = 40014*(seed1-k*53668)-k*12211 if (seed1 .lt. 0) seed1 = seed1 + 2147483563 k = seed2/52774 seed2 = 40692*(seed2-k*52774)-k*3791 if (seed2 .lt. 0) seed2 = seed2 + 2147483399 iz = seed1 - seed2 if(iz .lt. 1) iz = iz + 2147483562 vec(i) = iz * minv 10 continue return end c******************************************************************** subroutine ranget (iseed1,iseed2) c The subroutine ranget sets its (only) parameter cseed to the current c value of the seed used in ranu and ranuv to generate random numbers. integer seed1, seed2, iseed1, iseed2 common /random/ seed1, seed2 save /random/ iseed1 = seed1 iseed2 = seed2 return end c*********************************************************************** subroutine ranset (iseed1,iseed2) c The subroutine ranset sets the seed (used in ranu and ranuv to gener- c ate random numbers) to the value of its (only) parameter iseed, pro- c vided iseed is between 0 and the modulus m (see the parameter state- c ment below) exclusive. Otherwise, ranset sets the seed to 1. This c subroutine should be called (once) before the first time ranu or ranuv c is called. integer iseed1,iseed2,seed1,seed2 common /random/ seed1,seed2 save /random/ if ((iseed1 .gt. 0) .and. (iseed1 .lt. 2147483563)) then seed1 = iseed1 else seed1 = 12345 endif if ((iseed2 .gt. 0) .and. (iseed2 .lt. 2147483399)) then seed2 = iseed2 else seed2 = 67890 endif return end c*********************************************************************** subroutine rannv1 (n, vec) c The subroutine rannv1 generates a vector of n pseudo-random numbers, c each of which is from a standard normal distribution (mean zero, vari- c ance one), based on the method of Odeh and Evans: "The Percentage c Points of the Normal Distribution," Applied Statistics, 23 (1974), c pp. 96-97. For a mean mu other than zero and/or a variance vector c sigma other than one, set: vec(i) = mu + sqrt(sigma(i))*vec(i). c See also the subroutine rannv2 below. The subroutine ranset (above) c should be called once before the first time this subroutine is called. integer n double precision vec(n),temp,p0,p1,p2,p3,p4,q0,q1,q2,q3,q4 parameter (p0=-.322232431088d0, p1=-1d0, p2=-.342242088547d0, * p3=-.204231210245d-1, p4=-.453642210148d-4, * q0=.99348462606d-1, q1=.588581570495d0, q2=.531103462366d0, * q3=.10353775285d0, q4=.38560700634d-2) if (n .lt. 1) return call ranuv(n, vec) do 10 i = 1, n temp = vec(i) if (temp .gt. 0.5d0) vec(i) = 1d0 - vec(i) vec(i) = sqrt(log(1d0/vec(i)**2)) vec(i) = vec(i) + * ((((vec(i) * p4 + p3) * vec(i) + p2) * * vec(i) + p1) * vec(i) + p0) / * ((((vec(i) * q4 + q3) * vec(i) + q2) * * vec(i) + q1) * vec(i) + q0) if (temp .lt. 0.5d0) vec(i) = -vec(i) 10 continue return end c*********************************************************************** subroutine rannv2 (n, vec, mean, var) c The subroutine rannv2 generates a vector of n pseudo-random numbers, c the ith component of which is from a normal distribution with the c specified mean and variance var(i). This vector is stored in the c parameter vec. The subroutine ranset (see above) should be called c (once) before the first time this subroutine is called. This subrou- c tine offers an alternative to rannv1 above. integer n double precision vec(n), mean, var(n) call rannv1(n, vec) do 10 i = 1, n vec(i) = sqrt(var(i)) * vec(i) 10 continue if (mean .ne. 0d0) then do 20 i = 1, n vec(i) = mean + vec(i) 20 continue end if return end c*********************************************************************** subroutine rannv3 (n, vec, mean, var) c This routine assigns to each component i of the vector vec(n) a c random number from a Gaussian distribution with prescribed mean and c variance var(i). vec(i) is assigned as a weighted average of m irv's c {r1, r2, r3, ... rm} chosen from a uniform dist. between -1 and 1: c vec(i) <-- alpham * (r1 + r2 + .... + rm), where alpham is chosen c so that vec(i) (which has a Gaussian dist. by the central limit thm. c with mean 0) will have a variance satisfying < vec(i)**2 > = var(i). c From < vec(i)**2 > = alpham **2 * < (r1 + r2 + ... + rm) **2 > c = alpham **2 * m * 1/3 = var(i), c alpham = sqrt ( 3 * var(i) / m ). Since the sum of irv's with c normal distributions has a normal density with a mean and variance c that are the sums of the individual means and variances, c < sum {vec(i)**2} > = sum {var(i)}. implicit double precision(a-h,o-z) parameter(zero=0.d0,one=1.d0,two=2.d0,three=3.d0,mlge=20) integer n double precision mean,vec(n),var(n) do 10 i = 1, n sumran = zero do 5 j = 1, mlge sumran = sumran + (one - two * ranu()) 5 continue alpham = sqrt ( (var(i) * three) / (dble(mlge)) ) vec(i) = alpham * sumran 10 continue if (mean.ne.zero) then do 20 i = 1, n vec(i) = mean + vec(i) 20 continue endif return end