4 C#include "gen/pilot.h"
5 SUBROUTINE RANLUX2(RVEC,LENV,Input_seed)
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
59 Integer Input_seed,JSDFLT_set
61 DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
62 DATA I24,J24,CARRY/24,10,0./
63 CCC Set starting seed value:
64 If(Input_seed.gt.0) then
65 JSDFLT_set = Input_seed
70 C Luxury Level 0 1 2 *3* 4
71 DATA NDSKIP/0, 24, 73, 199, 365 /
72 Corresponds to p=24 48 97 223 389
73 C time factor 1 2 3 6 10 on slow workstation
74 C 1 1.5 2 3 5 on fast mainframe
76 C NOTYET is .TRUE. if no initialization has been performed yet.
77 C Default Initialization by Multiplicative Congruential
82 WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
84 NSKIP = NDSKIP(LUXLEV)
89 WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ',
95 JSEED = 40014*(JSEED-K*53668) -K*12211
96 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
97 ISEEDS(I) = MOD(JSEED,ITWO24)
99 TWOM12 = TWOM24 * 4096.
101 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
108 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
111 C The Generator proper: "Subtract-with-borrow",
112 C as proposed by Marsaglia and Zaman,
113 C Florida State University, March, 1989
116 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
117 IF (UNI .LT. 0.) THEN
127 C small numbers (with less than 12 "significant" bits) are "padded".
128 IF (UNI .LT. TWOM12) THEN
129 RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
130 C and zero is forbidden in case someone takes a logarithm
131 IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24
133 C Skipping to luxury. As proposed by Martin Luscher.
135 IF (IN24 .EQ. 24) THEN
137 KOUNT = KOUNT + NSKIP
139 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
140 IF (UNI .LT. 0.) THEN
153 IF (KOUNT .GE. IGIGA) THEN
155 KOUNT = KOUNT - IGIGA
159 C Entry to input and float integer seeds from previous run
165 195 TWOM24 = TWOM24 * 0.5
167 TWOM12 = TWOM24 * 4096.
168 WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
169 WRITE(6,'(5X,5I12)') ISDEXT
171 SEEDS(I) = REAL(ISDEXT(I))*TWOM24
174 IF (ISDEXT(25) .LT. 0) CARRY = TWOM24
175 ISD = IABS(ISDEXT(25))
183 IF (LUXLEV .LE. MAXLEV) THEN
184 NSKIP = NDSKIP(LUXLEV)
185 WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
187 ELSE IF (LUXLEV .GE. 24) THEN
189 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
191 NSKIP = NDSKIP(MAXLEV)
192 WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
198 C Entry to ouput seeds as integers
201 ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
203 ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
204 IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25)
207 C Entry to output the "convenient" restart point
208 ENTRY RLUXAT(LOUT,INOUT,K1,K2)
215 C Entry to initialize from one or three integers
216 ENTRY RLUXGO(LUX,INS,K1,K2)
219 ELSE IF (LUX .LE. MAXLEV) THEN
221 ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
223 WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
226 DO 310 ILX= 0, MAXLEV
227 IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX
230 IF (LUXLEV .LE. MAXLEV) THEN
231 NSKIP = NDSKIP(LUXLEV)
232 WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
233 + LUXLEV,' P=', NSKIP+24
236 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
239 IF (INS .LT. 0) WRITE (6,'(A)')
240 + ' Illegal initialization by RLUXGO, negative input seed'
243 WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
247 WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
253 TWOM24 = TWOM24 * 0.5
255 JSEED = 40014*(JSEED-K*53668) -K*12211
256 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
257 ISEEDS(I) = MOD(JSEED,ITWO24)
259 TWOM12 = TWOM24 * 4096.
261 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
268 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
269 C If restarting at a break point, skip K1 + IGIGA*K2
270 C Note that this is the number of numbers delivered to
271 C the user PLUS the number skipped (if luxury .GT. 0).
274 IF (K1+K2 .NE. 0) THEN
275 DO 500 IOUTER= 1, K2+1
277 IF (IOUTER .EQ. K2+1) INNER = K1
279 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
280 IF (UNI .LT. 0.) THEN
291 C Get the right value of IN24 by direct calculation
292 IN24 = MOD(KOUNT, NSKIP+24)
293 IF (MKOUNT .GT. 0) THEN
294 IZIP = MOD(IGIGA, NSKIP+24)
295 IZIP2 = MKOUNT*IZIP + IN24
296 IN24 = MOD(IZIP2, NSKIP+24)
298 C Now IN24 had better be between zero and 23 inclusive
299 IF (IN24 .GT. 23) THEN
300 WRITE (6,'(A/A,3I11,A,I5)')
301 + ' Error in RESTARTING with RLUXGO:',' The values', INS,
302 + K1, K2, ' cannot occur at luxury level', LUXLEV