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