Dummy implementation of Alarms
[u/mrichter/AliRoot.git] / MEVSIM / ranlux.F
CommitLineData
c28b8f8d 1C*
2C* $Id$
3C*
c28b8f8d 4C#include "gen/pilot.h"
5 SUBROUTINE RANLUX(RVEC,LENV)
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 LOGICAL NOTYET
60 DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
61 DATA I24,J24,CARRY/24,10,0./
62C default
63C Luxury Level 0 1 2 *3* 4
64 DATA NDSKIP/0, 24, 73, 199, 365 /
65Corresponds to p=24 48 97 223 389
66C time factor 1 2 3 6 10 on slow workstation
67C 1 1.5 2 3 5 on fast mainframe
68C
69C NOTYET is .TRUE. if no initialization has been performed yet.
70C Default Initialization by Multiplicative Congruential
71 IF (NOTYET) THEN
72 NOTYET = .FALSE.
73 JSEED = JSDFLT
74 INSEED = JSEED
75 WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
76 LUXLEV = LXDFLT
77 NSKIP = NDSKIP(LUXLEV)
78 LP = NSKIP + 24
79 IN24 = 0
80 KOUNT = 0
81 MKOUNT = 0
82 WRITE(6,'(A,I2,A,I4)') ' RANLUX DEFAULT LUXURY LEVEL = ',
83 + LUXLEV,' p =',LP
84 TWOM24 = 1.
85 DO 25 I= 1, 24
86 TWOM24 = TWOM24 * 0.5
87 K = JSEED/53668
88 JSEED = 40014*(JSEED-K*53668) -K*12211
89 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
90 ISEEDS(I) = MOD(JSEED,ITWO24)
91 25 CONTINUE
92 TWOM12 = TWOM24 * 4096.
93 DO 50 I= 1,24
94 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
95 NEXT(I) = I-1
96 50 CONTINUE
97 NEXT(1) = 24
98 I24 = 24
99 J24 = 10
100 CARRY = 0.
101 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
102 ENDIF
103C
104C The Generator proper: "Subtract-with-borrow",
105C as proposed by Marsaglia and Zaman,
106C Florida State University, March, 1989
107C
108 DO 100 IVEC= 1, LENV
109 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
110 IF (UNI .LT. 0.) THEN
111 UNI = UNI + 1.0
112 CARRY = TWOM24
113 ELSE
114 CARRY = 0.
115 ENDIF
116 SEEDS(I24) = UNI
117 I24 = NEXT(I24)
118 J24 = NEXT(J24)
119 RVEC(IVEC) = UNI
120C small numbers (with less than 12 "significant" bits) are "padded".
121 IF (UNI .LT. TWOM12) THEN
122 RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
123C and zero is forbidden in case someone takes a logarithm
124 IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24
125 ENDIF
126C Skipping to luxury. As proposed by Martin Luscher.
127 IN24 = IN24 + 1
128 IF (IN24 .EQ. 24) THEN
129 IN24 = 0
130 KOUNT = KOUNT + NSKIP
131 DO 90 ISK= 1, NSKIP
132 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
133 IF (UNI .LT. 0.) THEN
134 UNI = UNI + 1.0
135 CARRY = TWOM24
136 ELSE
137 CARRY = 0.
138 ENDIF
139 SEEDS(I24) = UNI
140 I24 = NEXT(I24)
141 J24 = NEXT(J24)
142 90 CONTINUE
143 ENDIF
144 100 CONTINUE
145 KOUNT = KOUNT + LENV
146 IF (KOUNT .GE. IGIGA) THEN
147 MKOUNT = MKOUNT + 1
148 KOUNT = KOUNT - IGIGA
149 ENDIF
150 RETURN
151C
152C Entry to input and float integer seeds from previous run
153 ENTRY RLUXIN(ISDEXT)
154 NOTYET = .FALSE.
155 TWOM24 = 1.
156 DO 195 I= 1, 24
157 NEXT(I) = I-1
158 195 TWOM24 = TWOM24 * 0.5
159 NEXT(1) = 24
160 TWOM12 = TWOM24 * 4096.
161 WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
162 WRITE(6,'(5X,5I12)') ISDEXT
163 DO 200 I= 1, 24
164 SEEDS(I) = REAL(ISDEXT(I))*TWOM24
165 200 CONTINUE
166 CARRY = 0.
167 IF (ISDEXT(25) .LT. 0) CARRY = TWOM24
168 ISD = IABS(ISDEXT(25))
169 I24 = MOD(ISD,100)
170 ISD = ISD/100
171 J24 = MOD(ISD,100)
172 ISD = ISD/100
173 IN24 = MOD(ISD,100)
174 ISD = ISD/100
175 LUXLEV = ISD
176 IF (LUXLEV .LE. MAXLEV) THEN
177 NSKIP = NDSKIP(LUXLEV)
178 WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
179 + LUXLEV
180 ELSE IF (LUXLEV .GE. 24) THEN
181 NSKIP = LUXLEV - 24
182 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
183 ELSE
184 NSKIP = NDSKIP(MAXLEV)
185 WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
186 LUXLEV = MAXLEV
187 ENDIF
188 INSEED = -1
189 RETURN
190C
191C Entry to ouput seeds as integers
192 ENTRY RLUXUT(ISDEXT)
193 DO 300 I= 1, 24
194 ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
195 300 CONTINUE
196 ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
197 IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25)
198 RETURN
199C
200C Entry to output the "convenient" restart point
201 ENTRY RLUXAT(LOUT,INOUT,K1,K2)
202 LOUT = LUXLEV
203 INOUT = INSEED
204 K1 = KOUNT
205 K2 = MKOUNT
206 RETURN
207C
208C Entry to initialize from one or three integers
209 ENTRY RLUXGO(LUX,INS,K1,K2)
210 IF (LUX .LT. 0) THEN
211 LUXLEV = LXDFLT
212 ELSE IF (LUX .LE. MAXLEV) THEN
213 LUXLEV = LUX
214 ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
215 LUXLEV = MAXLEV
216 WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
217 ELSE
218 LUXLEV = LUX
219 DO 310 ILX= 0, MAXLEV
220 IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX
221 310 CONTINUE
222 ENDIF
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
227 ELSE
228 NSKIP = LUXLEV - 24
229 WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
230 ENDIF
231 IN24 = 0
232 IF (INS .LT. 0) WRITE (6,'(A)')
233 + ' Illegal initialization by RLUXGO, negative input seed'
234 IF (INS .GT. 0) THEN
235 JSEED = INS
236 WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
237 + JSEED, K1,K2
238 ELSE
239 JSEED = JSDFLT
240 WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
241 ENDIF
242 INSEED = JSEED
243 NOTYET = .FALSE.
244 TWOM24 = 1.
245 DO 325 I= 1, 24
246 TWOM24 = TWOM24 * 0.5
247 K = JSEED/53668
248 JSEED = 40014*(JSEED-K*53668) -K*12211
249 IF (JSEED .LT. 0) JSEED = JSEED+ICONS
250 ISEEDS(I) = MOD(JSEED,ITWO24)
251 325 CONTINUE
252 TWOM12 = TWOM24 * 4096.
253 DO 350 I= 1,24
254 SEEDS(I) = REAL(ISEEDS(I))*TWOM24
255 NEXT(I) = I-1
256 350 CONTINUE
257 NEXT(1) = 24
258 I24 = 24
259 J24 = 10
260 CARRY = 0.
261 IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
262C If restarting at a break point, skip K1 + IGIGA*K2
263C Note that this is the number of numbers delivered to
264C the user PLUS the number skipped (if luxury .GT. 0).
265 KOUNT = K1
266 MKOUNT = K2
267 IF (K1+K2 .NE. 0) THEN
268 DO 500 IOUTER= 1, K2+1
269 INNER = IGIGA
270 IF (IOUTER .EQ. K2+1) INNER = K1
271 DO 450 ISK= 1, INNER
272 UNI = SEEDS(J24) - SEEDS(I24) - CARRY
273 IF (UNI .LT. 0.) THEN
274 UNI = UNI + 1.0
275 CARRY = TWOM24
276 ELSE
277 CARRY = 0.
278 ENDIF
279 SEEDS(I24) = UNI
280 I24 = NEXT(I24)
281 J24 = NEXT(J24)
282 450 CONTINUE
283 500 CONTINUE
284C 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)
290 ENDIF
291C 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
296 IN24 = 0
297 ENDIF
298 ENDIF
299 RETURN
300 END