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 RANLUX2(RVEC,LENV,Input_seed)
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
68 Integer Input_seed,JSDFLT_set
70 DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
71 DATA I24,J24,CARRY/24,10,0./
72 CCC Set starting seed value:
73 If(Input_seed.gt.0) then
74 JSDFLT_set = Input_seed
79 C Luxury Level 0 1 2 *3* 4
80 DATA NDSKIP/0, 24, 73, 199, 365 /
81 Corresponds to p=24 48 97 223 389
82 C time factor 1 2 3 6 10 on slow workstation
83 C 1 1.5 2 3 5 on fast mainframe
85 C NOTYET is .TRUE. if no initialization has been performed yet.
86 C Default Initialization by Multiplicative Congruential
91 WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
93 NSKIP = NDSKIP(LUXLEV)
98 WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ',
102 TWOM24 = TWOM24 * 0.5
104 JSEED = 40014*(JSEED-K*53668) -K*12211
105 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
106 ISEEDS(I) = MOD(JSEED,ITWO24)
108 TWOM12 = TWOM24 * 4096.
110 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
117 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
120 C The Generator proper: "Subtract-with-borrow",
121 C as proposed by Marsaglia and Zaman,
122 C Florida State University, March, 1989
125 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
126 IF (UNI .LT. 0.) THEN
136 C small numbers (with less than 12 "significant" bits) are "padded".
137 IF (UNI .LT. TWOM12) THEN
138 RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
139 C and zero is forbidden in case someone takes a logarithm
140 IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24
142 C Skipping to luxury. As proposed by Martin Luscher.
144 IF (IN24 .EQ. 24) THEN
146 KOUNT = KOUNT + NSKIP
148 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
149 IF (UNI .LT. 0.) THEN
162 IF (KOUNT .GE. IGIGA) THEN
164 KOUNT = KOUNT - IGIGA
168 C Entry to input and float integer seeds from previous run
174 195 TWOM24 = TWOM24 * 0.5
176 TWOM12 = TWOM24 * 4096.
177 WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
178 WRITE(6,'(5X,5I12)') ISDEXT
180 SEEDS(I) = REAL(ISDEXT(I))*TWOM24
183 IF (ISDEXT(25) .LT. 0) CARRY = TWOM24
184 ISD = IABS(ISDEXT(25))
192 IF (LUXLEV .LE. MAXLEV) THEN
193 NSKIP = NDSKIP(LUXLEV)
194 WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
196 ELSE IF (LUXLEV .GE. 24) THEN
198 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
200 NSKIP = NDSKIP(MAXLEV)
201 WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
207 C Entry to ouput seeds as integers
210 ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
212 ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
213 IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25)
216 C Entry to output the "convenient" restart point
217 ENTRY RLUXAT(LOUT,INOUT,K1,K2)
224 C Entry to initialize from one or three integers
225 ENTRY RLUXGO(LUX,INS,K1,K2)
228 ELSE IF (LUX .LE. MAXLEV) THEN
230 ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
232 WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
235 DO 310 ILX= 0, MAXLEV
236 IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX
239 IF (LUXLEV .LE. MAXLEV) THEN
240 NSKIP = NDSKIP(LUXLEV)
241 WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
242 + LUXLEV,' P=', NSKIP+24
245 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
248 IF (INS .LT. 0) WRITE (6,'(A)')
249 + ' Illegal initialization by RLUXGO, negative input seed'
252 WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
256 WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
262 TWOM24 = TWOM24 * 0.5
264 JSEED = 40014*(JSEED-K*53668) -K*12211
265 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
266 ISEEDS(I) = MOD(JSEED,ITWO24)
268 TWOM12 = TWOM24 * 4096.
270 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
277 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
278 C If restarting at a break point, skip K1 + IGIGA*K2
279 C Note that this is the number of numbers delivered to
280 C the user PLUS the number skipped (if luxury .GT. 0).
283 IF (K1+K2 .NE. 0) THEN
284 DO 500 IOUTER= 1, K2+1
286 IF (IOUTER .EQ. K2+1) INNER = K1
288 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
289 IF (UNI .LT. 0.) THEN
300 C Get the right value of IN24 by direct calculation
301 IN24 = MOD(KOUNT, NSKIP+24)
302 IF (MKOUNT .GT. 0) THEN
303 IZIP = MOD(IGIGA, NSKIP+24)
304 IZIP2 = MKOUNT*IZIP + IN24
305 IN24 = MOD(IZIP2, NSKIP+24)
307 C Now IN24 had better be between zero and 23 inclusive
308 IF (IN24 .GT. 23) THEN
309 WRITE (6,'(A/A,3I11,A,I5)')
310 + ' Error in RESTARTING with RLUXGO:',' The values', INS,
311 + K1, K2, ' cannot occur at luxury level', LUXLEV