4 C#include "gen/pilot.h"
5 SUBROUTINE RANLUX(RVEC,LENV)
6 C Subtract-and-borrow random number generator proposed by
7 C Marsaglia and Zaman, implemented by F. James with the name
8 C RCARRY in 1991, and later improved by Martin Luescher
9 C in 1993 to produce "Luxury Pseudorandom Numbers".
10 C Fortran 77 coded by F. James, 1993
13 C ------ ------ The available luxury levels are:
15 C level 0 (p=24): equivalent to the original RCARRY of Marsaglia
16 C and Zaman, very long period, but fails many tests.
17 C level 1 (p=48): considerable improvement in quality over level 0,
18 C now passes the gap test, but still fails spectral test.
19 C level 2 (p=97): passes all known tests, but theoretically still
21 C level 3 (p=223): DEFAULT VALUE. Any theoretically possible
22 C correlations have very small chance of being observed.
23 C level 4 (p=389): highest possible luxury, all 24 bits chaotic.
25 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
26 C!!! Calling sequences for RANLUX: ++
27 C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++
28 C!!! 32-bit random floating point numbers between ++
29 C!!! zero (not included) and one (also not incl.). ++
30 C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++
31 C!!! one 32-bit integer INT and sets Luxury Level LUX ++
32 C!!! which is integer between zero and MAXLEV, or if ++
33 C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++
34 C!!! should be set to zero unless restarting at a break++
35 C!!! point given by output of RLUXAT (see RLUXAT). ++
36 C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
37 C!!! which can be used to restart the RANLUX generator ++
38 C!!! at the current point by calling RLUXGO. K1 and K2++
39 C!!! specify how many numbers were generated since the ++
40 C!!! initialization with LUX and INT. The restarting ++
41 C!!! skips over K1+K2*E9 numbers, so it can be long.++
42 C!!! A more efficient but less convenient way of restarting is by: ++
43 C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++
44 C!!! ISVEC of 25 32-bit integers (see RLUXUT) ++
45 C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++
46 C!!! 32-bit integer seeds, to be used for restarting ++
47 C!!! ISVEC must be dimensioned 25 in the calling program ++
48 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
50 DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
51 PARAMETER (MAXLEV=4, LXDFLT=3)
52 DIMENSION NDSKIP(0:MAXLEV)
54 PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
55 PARAMETER (ITWO24=2**24, ICONS=2147483563)
56 SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
57 SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
60 DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
61 DATA I24,J24,CARRY/24,10,0./
63 C Luxury Level 0 1 2 *3* 4
64 DATA NDSKIP/0, 24, 73, 199, 365 /
65 Corresponds to p=24 48 97 223 389
66 C time factor 1 2 3 6 10 on slow workstation
67 C 1 1.5 2 3 5 on fast mainframe
69 C NOTYET is .TRUE. if no initialization has been performed yet.
70 C Default Initialization by Multiplicative Congruential
75 WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
77 NSKIP = NDSKIP(LUXLEV)
82 WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ',
88 JSEED = 40014*(JSEED-K*53668) -K*12211
89 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
90 ISEEDS(I) = MOD(JSEED,ITWO24)
92 TWOM12 = TWOM24 * 4096.
94 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
101 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
104 C The Generator proper: "Subtract-with-borrow",
105 C as proposed by Marsaglia and Zaman,
106 C Florida State University, March, 1989
109 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
110 IF (UNI .LT. 0.) THEN
120 C small numbers (with less than 12 "significant" bits) are "padded".
121 IF (UNI .LT. TWOM12) THEN
122 RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
123 C and zero is forbidden in case someone takes a logarithm
124 IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24
126 C Skipping to luxury. As proposed by Martin Luscher.
128 IF (IN24 .EQ. 24) THEN
130 KOUNT = KOUNT + NSKIP
132 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
133 IF (UNI .LT. 0.) THEN
146 IF (KOUNT .GE. IGIGA) THEN
148 KOUNT = KOUNT - IGIGA
152 C Entry to input and float integer seeds from previous run
158 195 TWOM24 = TWOM24 * 0.5
160 TWOM12 = TWOM24 * 4096.
161 WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
162 WRITE(6,'(5X,5I12)') ISDEXT
164 SEEDS(I) = REAL(ISDEXT(I))*TWOM24
167 IF (ISDEXT(25) .LT. 0) CARRY = TWOM24
168 ISD = IABS(ISDEXT(25))
176 IF (LUXLEV .LE. MAXLEV) THEN
177 NSKIP = NDSKIP(LUXLEV)
178 WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
180 ELSE IF (LUXLEV .GE. 24) THEN
182 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
184 NSKIP = NDSKIP(MAXLEV)
185 WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
191 C Entry to ouput seeds as integers
194 ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
196 ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
197 IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25)
200 C Entry to output the "convenient" restart point
201 ENTRY RLUXAT(LOUT,INOUT,K1,K2)
208 C Entry to initialize from one or three integers
209 ENTRY RLUXGO(LUX,INS,K1,K2)
212 ELSE IF (LUX .LE. MAXLEV) THEN
214 ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
216 WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
219 DO 310 ILX= 0, MAXLEV
220 IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX
223 IF (LUXLEV .LE. MAXLEV) THEN
224 NSKIP = NDSKIP(LUXLEV)
225 WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
226 + LUXLEV,' P=', NSKIP+24
229 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
232 IF (INS .LT. 0) WRITE (6,'(A)')
233 + ' Illegal initialization by RLUXGO, negative input seed'
236 WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
240 WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
246 TWOM24 = TWOM24 * 0.5
248 JSEED = 40014*(JSEED-K*53668) -K*12211
249 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
250 ISEEDS(I) = MOD(JSEED,ITWO24)
252 TWOM12 = TWOM24 * 4096.
254 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
261 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
262 C If restarting at a break point, skip K1 + IGIGA*K2
263 C Note that this is the number of numbers delivered to
264 C the user PLUS the number skipped (if luxury .GT. 0).
267 IF (K1+K2 .NE. 0) THEN
268 DO 500 IOUTER= 1, K2+1
270 IF (IOUTER .EQ. K2+1) INNER = K1
272 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
273 IF (UNI .LT. 0.) THEN
284 C Get the right value of IN24 by direct calculation
285 IN24 = MOD(KOUNT, NSKIP+24)
286 IF (MKOUNT .GT. 0) THEN
287 IZIP = MOD(IGIGA, NSKIP+24)
288 IZIP2 = MKOUNT*IZIP + IN24
289 IN24 = MOD(IZIP2, NSKIP+24)
291 C Now IN24 had better be between zero and 23 inclusive
292 IF (IN24 .GT. 23) THEN
293 WRITE (6,'(A/A,3I11,A,I5)')
294 + ' Error in RESTARTING with RLUXGO:',' The values', INS,
295 + K1, K2, ' cannot occur at luxury level', LUXLEV