C********************************************************************C C The material in this file is provided to you "as-is" and without C C warranty of any kind, express, implied or otherwise, including C C without limitation, any warranty of fitness for a particular C C purpose. C C This file may be distributed freely for non-commercial purposes, C C provided this notice, as well as the copyright notice, included C C as part of the file. C C Distribution of this file for any other purpose is unlawful. C C Copyright (c) 1999 by Tamar Schlick & Xiaoliang Qian C C********************************************************************C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Combined Linear Congruential Generator (CLFG) C C Adapted from Pierre L'Ecuyer & Terry H Andres' C version code. C C References C C [1] P. L'Ecuyer and T. H. Andres, C C ``A Random Number Generator Based on the Combination C C of Four LCGs'', Mathematics and Computers in Simulation, C C 44 (1997), 99--107. C C [2] http://www.iro.umontreal.ca/~lecuyer/ C C C C For further information, please contact C C Tamar Schlick C C schlick@nyu.edu C C Converted to FORTRAN by C C Xiaoliang Qian 10/7/99 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C This pseudorandom number generator combines 4 linear C C congruential generators (LCGs) to get the long period C C of about 2^121 for the resulting sequence. This sequence has C C passed many statistical tests of `randomness'. C C Essentially, the four LCGs defined as C C C C X{j,n} = a{j}X{j,n-1} mod m{j} (1) C C j=1,2,3,4 C C n=1,2,3, ... C C C C with m{1}= 2147483647, m{2}=2147483543, m{3}=2147483423, C C m{4}=2147483323 and a{1}=45991, a{2}=207707, a{3}=138556, C C a{4}=49689. C C C C The construct C C Z{n} = (Sum [(-1)^{j+1}*X{j,n}/m{j}]) mod 1 (2) C C C C for n=1,2,3, ... is then a uniformly distributed random sequence C C in (0,1). It can be proved that the LCG corresponding to the C C combined generator has modulus, multiplier, and period length of C C 21267641435849934371830464348413044909, C C 5494569482908719143153333426731027229, C C (2^{31}-2)(2^{31}-106)(2^{31}-226)(2^{31}-326) ~ 2^{121}, C C respectively. C C C C The default initial seed is the vector {11111111, 22222222, C C 33333333, 44444444} and can be changed by calling SetIniSD C C after calling Init. C C C C This RNG can be used under parallel conditions to give C C independent random number sequence when each processor C C calls with different stream number g (e.g., GenVal(g)). C C C C To use these RNG routines, the user should proceed as follows: C C C C 1. Call routine Init() to initialize all Maxgen (100) streams C C using four default initial seeds. C C 2. [Optional] Call SetiniSD(sd) with desired seed array C C (4 values) to override the default values specified in Init(). C C 3. Call function GenVal(k), where k is an integer from 1 to 100 C C specifying the stream number. For parallel codes, k can be set C C to a processor id number. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccC C A test program for the CLCG below. C C f77 -o clcgf clcgf.f C C C Program main integer ng,nsteps,sd(4) C ng number of independent random number streams C nsteps number of random numbers in one sequence C sd(4) the 4 initial seeds for the first stream parameter (ng = 3, nsteps = 1000) double precision rnd(ng),GenVal data sd /12345,67890,54321,9876/ integer i,j,k C Set up the Maxgen number of streams with default settings call Init() C Optional: reset the initial seeds of each stream call SetiniSD(sd) C Generate nsteps of random numbers from each stream do 100 i = 1,nsteps do 50 k = 1,ng rnd(k) = GenVal (k) 50 continue write(*,'(f20.16,1x,f20.16,1x,f20.16)') (rnd(j),j = 1,ng) 100 continue end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine Init ( ) C------------------------------------------------------------------------C C Initialize the RNG with seed values in vector sd of dimension 4. C C Each such initial vector generates one stream of random variates C C combining the 4 LCG sequences with resulting sequence length C C 2^{121}. C C The aw and avw arrays of dimension 4 have default values C C aw{j} = a{j}^{2^31} mod m{j} and avw{j} = a{j}^{2^41} mod m{j}, C C for j=1, ..., 4 corresponding to the 4 LCGs. C integer sd(4) integer v,w,j,i parameter (v = 31, w = 41) data sd/11111111,22222222,33333333,44444444/ integer IniSD,LstSD,NewSD,Maxgen parameter (Maxgen = 100, IniSD = 1, LstSD = 2, NewSD = 3) integer a(4),m(4),aw(4),avw(4) integer Ig(4,Maxgen),Lg(4,Maxgen),Cg(4,Maxgen) data a/45991,207707,138556,49689/ data m/2147483647,2147483543,2147483423,2147483323/ common /RANDOM/Ig,Lg,Cg,aw,avw SAVE /RANDOM/ do 100 j = 1,4 aw(j) = a(j) do 85 i = 1,w aw(j) = MulMod (aw(j),aw(j),m(j)) 85 continue avw(j) = aw(j) do 90 i = 1,v avw(j) = MulMod (avw(j),avw(j),m(j)) 90 continue 100 continue call SetiniSD (sd) end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine SetiniSD (s) C------------------------------------------------------------------------C C Set initial seed values for all 100 (Maxgen) streams using the C C initial seeds from the first stream. C integer g,s(4),j integer IniSD,LstSD,NewSD,Maxgen parameter (Maxgen = 100, IniSD = 1, LstSD = 2, NewSD =3) integer m(4),aw(4),avw(4) integer Ig(4,Maxgen),Lg(4,Maxgen),Cg(4,Maxgen) data m/2147483647,2147483543,2147483423,2147483323/ common /RANDOM/Ig,Lg,Cg,aw ,avw SAVE /RANDOM/ do 70 j = 1,4 Ig(j,1) = s(j) 70 continue call IniGen (1,IniSD) do 80 g = 2,Maxgen do 75 j = 1,4 Ig(j,g) = MulMod (avw(j),Ig(j,g-1),m(j)) 75 continue call IniGen (g,IniSD) 80 continue end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC double precision function GenVal (g) C------------------------------------------------------------------------C C Return a double precision uniformly distributed random number in C C (0,1) from the gth stream and reset the current seed Cg accordingly C C (i.e., using one of the 100 initial seed vectors generated in the C C SetiniSD routine). C integer g,k,s,j integer IniSD,LstSD,NewSD,Maxgen double precision u(4),one,zero parameter (Maxgen = 100, IniSD = 1, LstSD = 2, NewSD = 3, & one = 1.0d0, zero = 0.0d0) integer a(4),m(4),aw(4),avw(4) integer Ig(4,Maxgen),Lg(4,Maxgen),Cg(4,Maxgen) data a/45991,207707,138556,49689/ data m/2147483647,2147483543,2147483423,2147483323/ common /RANDOM/Ig,Lg,Cg,aw,avw SAVE /RANDOM/ integer dv(4),mv(4) data dv/46693,10339,15499,43218/ data mv/25884,870,3979,24121/ data u/4.65661287524579692d-10,-4.65661310075985993d-10, & 4.65661336096842131d-10,-4.65661357780891134d-10/ if (g .gt. Maxgen) write(*,*) 'ERROR:GenVal With g>Maxgen!' GenVal = zero do 110 j = 1,4 s = Cg(j,g) k = s / dv(j) s = a(j) * (s - k *dv(j)) - k * mv(j) if (s .lt. 0) s = s + m(j) Cg(j,g) = s GenVal = GenVal + u(j) * s if (GenVal .lt. zero) GenVal = GenVal + one if (GenVal .ge. one) GenVal = GenVal - one 110 continue return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine SetSeed (g,s) C------------------------------------------------------------------------C C This optional routine uses the input seed value s for stream g C C instead of the default settings (routine SetiniSD). C integer g,s(4),j integer IniSD,LstSD,NewSD,Maxgen parameter (Maxgen = 100, IniSD = 1, LstSD = 2, NewSD = 3) integer Ig(4,Maxgen) common /RANDOM/Ig SAVE /RANDOM/ if (g .gt. Maxgen) write(*,*) 'ERROR:SetSeed With g>Maxgen!' do 50 j = 1,4 Ig(j,g) = s(j) 50 continue call IniGen (g,IniSD) end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine IniGen (g,type) C------------------------------------------------------------------------C C This optional routine resets the gth stream so that the initial seed C C is either the original initial seed (if type = 1) or the last seed C C (if type = 3). C integer g,type,j integer IniSD,LstSD,NewSD,Maxgen parameter (Maxgen = 100, IniSD = 1, LstSD = 2, NewSD = 3) integer m(4),aw(4),avw(4) integer Ig(4,Maxgen),Lg(4,Maxgen),Cg(4,Maxgen) data m/2147483647,2147483543,2147483423,2147483323/ common /RANDOM/Ig,Lg,Cg,aw ,avw SAVE /RANDOM/ if (g .gt. Maxgen) write(*,*) 'ERROR:IniGen With g>Maxgen!' do 60 j = 1,4 if (type .eq. IniSD) then Lg(j,g) = Ig(j,g) else if (type .eq. NewSD) Lg(j,g) = MulMod (aw(j),Lg(j,g),m(j)) endif Cg(j,g) = Lg(j,g) 60 continue end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC integer function MulMod (s,t,M) C--------------------------------------------------------------------C C Return s*t mod M. All numbers out of range are truncated C C before the mod operation to avoid overflow. See Park & Miller, C C Comm. ACM 31, 1192 (1988) for this multiplication procedure. C integer s,t,M,H parameter (H = 32768) integer S0,S1,q,qh,rh,k if (s .lt. 0) s = s + M if (t .lt. 0) t = t + M if (s .lt. H) then S0 = s MulMod = 0 else S1 = s / H S0 = s - H * S1 qh = M / H rh = M - H * qh if (S1 .ge. H) then S1 = S1 - H k = t / qh MulMod = H * (t - k * qh) - k * rh 10 if (MulMod .lt. 0) then MulMod = MulMod + M goto 10 endif else MulMod = 0 endif if (S1 .ne. 0) then q =M / S1 k = t / q MulMod = MulMod - k * (M - S1 * q) if (MulMod .gt. 0) MulMod = MulMod - M MulMod = MulMod + S1 * (t - k * q) 20 if (MulMod .lt. 0) then MulMod = MulMod + M goto 20 endif endif k = MulMod / qh MulMod = H * (MulMod - k * qh) - k * rh 30 if (MulMod .lt. 0) then MulMod = MulMod + M goto 30 endif endif if (S0 .ne. 0) then Q =M / S0 k = t / q MulMod = MulMod - k * (M - S0 * q) if (MulMod .gt. 0) MulMod = MulMod - M MulMod = MulMod + S0 * (t - k * q) 40 if (MulMod .lt. 0) then MulMod = MulMod + M goto 40 endif endif return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC