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