Plane efficiency framework (Guiseppe)
[u/mrichter/AliRoot.git] / MEVSIM / ranlux2.F
CommitLineData
c28b8f8d 1C*
2C* $Id$
3C*
c28b8f8d 4C#include "gen/pilot.h"
5 SUBROUTINE RANLUX2(RVEC,LENV,Input_seed)
6C Subtract-and-borrow random number generator proposed by
7C Marsaglia and Zaman, implemented by F. James with the name
8C RCARRY in 1991, and later improved by Martin Luescher
9C in 1993 to produce "Luxury Pseudorandom Numbers".
10C Fortran 77 coded by F. James, 1993
11C
12C LUXURY LEVELS.
13C ------ ------ The available luxury levels are:
14C
15C level 0 (p=24): equivalent to the original RCARRY of Marsaglia
16C and Zaman, very long period, but fails many tests.
17C level 1 (p=48): considerable improvement in quality over level 0,
18C now passes the gap test, but still fails spectral test.
19C level 2 (p=97): passes all known tests, but theoretically still
20C defective.
21C level 3 (p=223): DEFAULT VALUE. Any theoretically possible
22C correlations have very small chance of being observed.
23C level 4 (p=389): highest possible luxury, all 24 bits chaotic.
24C
25C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
26C!!! Calling sequences for RANLUX: ++
27C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++
28C!!! 32-bit random floating point numbers between ++
29C!!! zero (not included) and one (also not incl.). ++
30C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++
31C!!! one 32-bit integer INT and sets Luxury Level LUX ++
32C!!! which is integer between zero and MAXLEV, or if ++
33C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++
34C!!! should be set to zero unless restarting at a break++
35C!!! point given by output of RLUXAT (see RLUXAT). ++
36C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
37C!!! which can be used to restart the RANLUX generator ++
38C!!! at the current point by calling RLUXGO. K1 and K2++
39C!!! specify how many numbers were generated since the ++
40C!!! initialization with LUX and INT. The restarting ++
41C!!! skips over K1+K2*E9 numbers, so it can be long.++
42C!!! A more efficient but less convenient way of restarting is by: ++
43C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++
44C!!! ISVEC of 25 32-bit integers (see RLUXUT) ++
45C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++
46C!!! 32-bit integer seeds, to be used for restarting ++
47C!!! ISVEC must be dimensioned 25 in the calling program ++
48C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
49 DIMENSION RVEC(LENV)
50 DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
51 PARAMETER (MAXLEV=4, LXDFLT=3)
52 DIMENSION NDSKIP(0:MAXLEV)
53 DIMENSION NEXT(24)
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
58 INTEGER LUXLEV
59 Integer Input_seed,JSDFLT_set
60 LOGICAL NOTYET
61 DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
62 DATA I24,J24,CARRY/24,10,0./
63CCC Set starting seed value:
64 If(Input_seed.gt.0) then
65 JSDFLT_set = Input_seed
66 Else
67 JSDFLT_set = JSDFLT
68 End If
69C default
70C Luxury Level 0 1 2 *3* 4
71 DATA NDSKIP/0, 24, 73, 199, 365 /
72Corresponds to p=24 48 97 223 389
73C time factor 1 2 3 6 10 on slow workstation
74C 1 1.5 2 3 5 on fast mainframe
75C
76C NOTYET is .TRUE. if no initialization has been performed yet.
77C Default Initialization by Multiplicative Congruential
78 IF (NOTYET) THEN
79 NOTYET = .FALSE.
80 JSEED = JSDFLT_set
81 INSEED = JSEED
82 WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
83 LUXLEV = LXDFLT
84 NSKIP = NDSKIP(LUXLEV)
85 LP = NSKIP + 24
86 IN24 = 0
87 KOUNT = 0
88 MKOUNT = 0
89 WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ',
90 + LUXLEV,' p =',LP
91 TWOM24 = 1.
92 DO 25 I= 1, 24
93 TWOM24 = TWOM24 * 0.5
94 K = JSEED/53668
95 JSEED = 40014*(JSEED-K*53668) -K*12211
96 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
97 ISEEDS(I) = MOD(JSEED,ITWO24)
98 25 CONTINUE
99 TWOM12 = TWOM24 * 4096.
100 DO 50 I= 1,24
101 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
102 NEXT(I) = I-1
103 50 CONTINUE
104 NEXT(1) = 24
105 I24 = 24
106 J24 = 10
107 CARRY = 0.
108 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
109 ENDIF
110C
111C The Generator proper: "Subtract-with-borrow",
112C as proposed by Marsaglia and Zaman,
113C Florida State University, March, 1989
114C
115 DO 100 IVEC= 1, LENV
116 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
117 IF (UNI .LT. 0.) THEN
118 UNI = UNI + 1.0
119 CARRY = TWOM24
120 ELSE
121 CARRY = 0.
122 ENDIF
123 SEEDS(I24) = UNI
124 I24 = NEXT(I24)
125 J24 = NEXT(J24)
126 RVEC(IVEC) = UNI
127C small numbers (with less than 12 "significant" bits) are "padded".
128 IF (UNI .LT. TWOM12) THEN
129 RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
130C and zero is forbidden in case someone takes a logarithm
131 IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24
132 ENDIF
133C Skipping to luxury. As proposed by Martin Luscher.
134 IN24 = IN24 + 1
135 IF (IN24 .EQ. 24) THEN
136 IN24 = 0
137 KOUNT = KOUNT + NSKIP
138 DO 90 ISK= 1, NSKIP
139 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
140 IF (UNI .LT. 0.) THEN
141 UNI = UNI + 1.0
142 CARRY = TWOM24
143 ELSE
144 CARRY = 0.
145 ENDIF
146 SEEDS(I24) = UNI
147 I24 = NEXT(I24)
148 J24 = NEXT(J24)
149 90 CONTINUE
150 ENDIF
151 100 CONTINUE
152 KOUNT = KOUNT + LENV
153 IF (KOUNT .GE. IGIGA) THEN
154 MKOUNT = MKOUNT + 1
155 KOUNT = KOUNT - IGIGA
156 ENDIF
157 RETURN
158C
159C Entry to input and float integer seeds from previous run
160 ENTRY RLUXIN(ISDEXT)
161 NOTYET = .FALSE.
162 TWOM24 = 1.
163 DO 195 I= 1, 24
164 NEXT(I) = I-1
165 195 TWOM24 = TWOM24 * 0.5
166 NEXT(1) = 24
167 TWOM12 = TWOM24 * 4096.
168 WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
169 WRITE(6,'(5X,5I12)') ISDEXT
170 DO 200 I= 1, 24
171 SEEDS(I) = REAL(ISDEXT(I))*TWOM24
172 200 CONTINUE
173 CARRY = 0.
174 IF (ISDEXT(25) .LT. 0) CARRY = TWOM24
175 ISD = IABS(ISDEXT(25))
176 I24 = MOD(ISD,100)
177 ISD = ISD/100
178 J24 = MOD(ISD,100)
179 ISD = ISD/100
180 IN24 = MOD(ISD,100)
181 ISD = ISD/100
182 LUXLEV = ISD
183 IF (LUXLEV .LE. MAXLEV) THEN
184 NSKIP = NDSKIP(LUXLEV)
185 WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
186 + LUXLEV
187 ELSE IF (LUXLEV .GE. 24) THEN
188 NSKIP = LUXLEV - 24
189 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
190 ELSE
191 NSKIP = NDSKIP(MAXLEV)
192 WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
193 LUXLEV = MAXLEV
194 ENDIF
195 INSEED = -1
196 RETURN
197C
198C Entry to ouput seeds as integers
199 ENTRY RLUXUT(ISDEXT)
200 DO 300 I= 1, 24
201 ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
202 300 CONTINUE
203 ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
204 IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25)
205 RETURN
206C
207C Entry to output the "convenient" restart point
208 ENTRY RLUXAT(LOUT,INOUT,K1,K2)
209 LOUT = LUXLEV
210 INOUT = INSEED
211 K1 = KOUNT
212 K2 = MKOUNT
213 RETURN
214C
215C Entry to initialize from one or three integers
216 ENTRY RLUXGO(LUX,INS,K1,K2)
217 IF (LUX .LT. 0) THEN
218 LUXLEV = LXDFLT
219 ELSE IF (LUX .LE. MAXLEV) THEN
220 LUXLEV = LUX
221 ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
222 LUXLEV = MAXLEV
223 WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
224 ELSE
225 LUXLEV = LUX
226 DO 310 ILX= 0, MAXLEV
227 IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX
228 310 CONTINUE
229 ENDIF
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
234 ELSE
235 NSKIP = LUXLEV - 24
236 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
237 ENDIF
238 IN24 = 0
239 IF (INS .LT. 0) WRITE (6,'(A)')
240 + ' Illegal initialization by RLUXGO, negative input seed'
241 IF (INS .GT. 0) THEN
242 JSEED = INS
243 WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
244 + JSEED, K1,K2
245 ELSE
246 JSEED = JSDFLT_set
247 WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
248 ENDIF
249 INSEED = JSEED
250 NOTYET = .FALSE.
251 TWOM24 = 1.
252 DO 325 I= 1, 24
253 TWOM24 = TWOM24 * 0.5
254 K = JSEED/53668
255 JSEED = 40014*(JSEED-K*53668) -K*12211
256 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
257 ISEEDS(I) = MOD(JSEED,ITWO24)
258 325 CONTINUE
259 TWOM12 = TWOM24 * 4096.
260 DO 350 I= 1,24
261 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
262 NEXT(I) = I-1
263 350 CONTINUE
264 NEXT(1) = 24
265 I24 = 24
266 J24 = 10
267 CARRY = 0.
268 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
269C If restarting at a break point, skip K1 + IGIGA*K2
270C Note that this is the number of numbers delivered to
271C the user PLUS the number skipped (if luxury .GT. 0).
272 KOUNT = K1
273 MKOUNT = K2
274 IF (K1+K2 .NE. 0) THEN
275 DO 500 IOUTER= 1, K2+1
276 INNER = IGIGA
277 IF (IOUTER .EQ. K2+1) INNER = K1
278 DO 450 ISK= 1, INNER
279 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
280 IF (UNI .LT. 0.) THEN
281 UNI = UNI + 1.0
282 CARRY = TWOM24
283 ELSE
284 CARRY = 0.
285 ENDIF
286 SEEDS(I24) = UNI
287 I24 = NEXT(I24)
288 J24 = NEXT(J24)
289 450 CONTINUE
290 500 CONTINUE
291C 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)
297 ENDIF
298C 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
303 IN24 = 0
304 ENDIF
305 ENDIF
306 RETURN
307 END