5 C* Revision 1.2 1997/09/22 13:45:47 mclareni
6 C* Correct error in initializing RANLUX by using RLUXIN with the output of
7 C* RLUXUT from a previous run.
9 C* Revision 1.1.1.1 1996/04/01 15:02:55 mclareni
13 C#include "gen/pilot.h"
14 SUBROUTINE RANLUX(RVEC,LENV)
15 C Subtract-and-borrow random number generator proposed by
16 C Marsaglia and Zaman, implemented by F. James with the name
17 C RCARRY in 1991, and later improved by Martin Luescher
18 C in 1993 to produce "Luxury Pseudorandom Numbers".
19 C Fortran 77 coded by F. James, 1993
22 C ------ ------ The available luxury levels are:
24 C level 0 (p=24): equivalent to the original RCARRY of Marsaglia
25 C and Zaman, very long period, but fails many tests.
26 C level 1 (p=48): considerable improvement in quality over level 0,
27 C now passes the gap test, but still fails spectral test.
28 C level 2 (p=97): passes all known tests, but theoretically still
30 C level 3 (p=223): DEFAULT VALUE. Any theoretically possible
31 C correlations have very small chance of being observed.
32 C level 4 (p=389): highest possible luxury, all 24 bits chaotic.
34 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35 C!!! Calling sequences for RANLUX: ++
36 C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++
37 C!!! 32-bit random floating point numbers between ++
38 C!!! zero (not included) and one (also not incl.). ++
39 C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++
40 C!!! one 32-bit integer INT and sets Luxury Level LUX ++
41 C!!! which is integer between zero and MAXLEV, or if ++
42 C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++
43 C!!! should be set to zero unless restarting at a break++
44 C!!! point given by output of RLUXAT (see RLUXAT). ++
45 C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
46 C!!! which can be used to restart the RANLUX generator ++
47 C!!! at the current point by calling RLUXGO. K1 and K2++
48 C!!! specify how many numbers were generated since the ++
49 C!!! initialization with LUX and INT. The restarting ++
50 C!!! skips over K1+K2*E9 numbers, so it can be long.++
51 C!!! A more efficient but less convenient way of restarting is by: ++
52 C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++
53 C!!! ISVEC of 25 32-bit integers (see RLUXUT) ++
54 C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++
55 C!!! 32-bit integer seeds, to be used for restarting ++
56 C!!! ISVEC must be dimensioned 25 in the calling program ++
57 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
59 DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
60 PARAMETER (MAXLEV=4, LXDFLT=3)
61 DIMENSION NDSKIP(0:MAXLEV)
63 PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
64 PARAMETER (ITWO24=2**24, ICONS=2147483563)
65 SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
66 SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
69 DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
70 DATA I24,J24,CARRY/24,10,0./
72 C Luxury Level 0 1 2 *3* 4
73 DATA NDSKIP/0, 24, 73, 199, 365 /
74 Corresponds to p=24 48 97 223 389
75 C time factor 1 2 3 6 10 on slow workstation
76 C 1 1.5 2 3 5 on fast mainframe
78 C NOTYET is .TRUE. if no initialization has been performed yet.
79 C Default Initialization by Multiplicative Congruential
84 WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
86 NSKIP = NDSKIP(LUXLEV)
91 WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ',
97 JSEED = 40014*(JSEED-K*53668) -K*12211
98 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
99 ISEEDS(I) = MOD(JSEED,ITWO24)
101 TWOM12 = TWOM24 * 4096.
103 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
110 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
113 C The Generator proper: "Subtract-with-borrow",
114 C as proposed by Marsaglia and Zaman,
115 C Florida State University, March, 1989
118 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
119 IF (UNI .LT. 0.) THEN
129 C small numbers (with less than 12 "significant" bits) are "padded".
130 IF (UNI .LT. TWOM12) THEN
131 RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
132 C and zero is forbidden in case someone takes a logarithm
133 IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24
135 C Skipping to luxury. As proposed by Martin Luscher.
137 IF (IN24 .EQ. 24) THEN
139 KOUNT = KOUNT + NSKIP
141 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
142 IF (UNI .LT. 0.) THEN
155 IF (KOUNT .GE. IGIGA) THEN
157 KOUNT = KOUNT - IGIGA
161 C Entry to input and float integer seeds from previous run
167 195 TWOM24 = TWOM24 * 0.5
169 TWOM12 = TWOM24 * 4096.
170 WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
171 WRITE(6,'(5X,5I12)') ISDEXT
173 SEEDS(I) = REAL(ISDEXT(I))*TWOM24
176 IF (ISDEXT(25) .LT. 0) CARRY = TWOM24
177 ISD = IABS(ISDEXT(25))
185 IF (LUXLEV .LE. MAXLEV) THEN
186 NSKIP = NDSKIP(LUXLEV)
187 WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
189 ELSE IF (LUXLEV .GE. 24) THEN
191 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
193 NSKIP = NDSKIP(MAXLEV)
194 WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
200 C Entry to ouput seeds as integers
203 ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
205 ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
206 IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25)
209 C Entry to output the "convenient" restart point
210 ENTRY RLUXAT(LOUT,INOUT,K1,K2)
217 C Entry to initialize from one or three integers
218 ENTRY RLUXGO(LUX,INS,K1,K2)
221 ELSE IF (LUX .LE. MAXLEV) THEN
223 ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
225 WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
228 DO 310 ILX= 0, MAXLEV
229 IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX
232 IF (LUXLEV .LE. MAXLEV) THEN
233 NSKIP = NDSKIP(LUXLEV)
234 WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
235 + LUXLEV,' P=', NSKIP+24
238 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
241 IF (INS .LT. 0) WRITE (6,'(A)')
242 + ' Illegal initialization by RLUXGO, negative input seed'
245 WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
249 WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
255 TWOM24 = TWOM24 * 0.5
257 JSEED = 40014*(JSEED-K*53668) -K*12211
258 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
259 ISEEDS(I) = MOD(JSEED,ITWO24)
261 TWOM12 = TWOM24 * 4096.
263 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
270 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
271 C If restarting at a break point, skip K1 + IGIGA*K2
272 C Note that this is the number of numbers delivered to
273 C the user PLUS the number skipped (if luxury .GT. 0).
276 IF (K1+K2 .NE. 0) THEN
277 DO 500 IOUTER= 1, K2+1
279 IF (IOUTER .EQ. K2+1) INNER = K1
281 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
282 IF (UNI .LT. 0.) THEN
293 C Get the right value of IN24 by direct calculation
294 IN24 = MOD(KOUNT, NSKIP+24)
295 IF (MKOUNT .GT. 0) THEN
296 IZIP = MOD(IGIGA, NSKIP+24)
297 IZIP2 = MKOUNT*IZIP + IN24
298 IN24 = MOD(IZIP2, NSKIP+24)
300 C Now IN24 had better be between zero and 23 inclusive
301 IF (IN24 .GT. 23) THEN
302 WRITE (6,'(A/A,3I11,A,I5)')
303 + ' Error in RESTARTING with RLUXGO:',' The values', INS,
304 + K1, K2, ' cannot occur at luxury level', LUXLEV