65a8c847932f6168ff28be575e841d59a48ee3f5
[u/mrichter/AliRoot.git] / MEVSIM / ranlux2.F
1 C*
2 C* $Id$
3 C*
4 C* $Log$
5 C* Revision 1.2  1997/09/22 13:45:47  mclareni
6 C* Correct error in initializing RANLUX by using RLUXIN with the output of
7 C* RLUXUT from a previous run.
8 C*
9 C* Revision 1.1.1.1  1996/04/01 15:02:55  mclareni
10 C* Mathlib gen
11 C*
12 C*
13 C#include "gen/pilot.h"
14       SUBROUTINE RANLUX2(RVEC,LENV,Input_seed)
15 C         Subtract-and-borrow random number generator proposed by
16 C         Marsaglia and Zaman, implemented by F. James with the name
17 C         RCARRY in 1991, and later improved by Martin Luescher
18 C         in 1993 to produce "Luxury Pseudorandom Numbers".
19 C     Fortran 77 coded by F. James, 1993
20 C
21 C   LUXURY LEVELS.
22 C   ------ ------      The available luxury levels are:
23 C
24 C  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
25 C           and Zaman, very long period, but fails many tests.
26 C  level 1  (p=48): considerable improvement in quality over level 0,
27 C           now passes the gap test, but still fails spectral test.
28 C  level 2  (p=97): passes all known tests, but theoretically still
29 C           defective.
30 C  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
31 C           correlations have very small chance of being observed.
32 C  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
33 C
34 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35 C!!!  Calling sequences for RANLUX:                                  ++
36 C!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
37 C!!!                   32-bit random floating point numbers between  ++
38 C!!!                   zero (not included) and one (also not incl.). ++
39 C!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
40 C!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
41 C!!!               which is integer between zero and MAXLEV, or if   ++
42 C!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
43 C!!!               should be set to zero unless restarting at a break++ 
44 C!!!               point given by output of RLUXAT (see RLUXAT).     ++
45 C!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
46 C!!!               which can be used to restart the RANLUX generator ++
47 C!!!               at the current point by calling RLUXGO.  K1 and K2++
48 C!!!               specify how many numbers were generated since the ++
49 C!!!               initialization with LUX and INT.  The restarting  ++
50 C!!!               skips over  K1+K2*E9   numbers, so it can be long.++
51 C!!!   A more efficient but less convenient way of restarting is by: ++
52 C!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
53 C!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
54 C!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
55 C!!!                 32-bit integer seeds, to be used for restarting ++
56 C!!!      ISVEC must be dimensioned 25 in the calling program        ++
57 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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./
72 CCC  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
78 C                               default
79 C  Luxury Level   0     1     2   *3*    4
80       DATA NDSKIP/0,   24,   73,  199,  365 /
81 Corresponds to p=24    48    97   223   389
82 C     time factor 1     2     3     6    10   on slow workstation
83 C                 1    1.5    2     3     5   on fast mainframe
84 C
85 C  NOTYET is .TRUE. if no initialization has been performed yet.
86 C              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
119 C
120 C          The Generator proper: "Subtract-with-borrow",
121 C          as proposed by Marsaglia and Zaman,
122 C          Florida State University, March, 1989
123 C
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
136 C  small numbers (with less than 12 "significant" bits) are "padded".
137       IF (UNI .LT. TWOM12)  THEN
138          RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
139 C        and zero is forbidden in case someone takes a logarithm
140          IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24
141       ENDIF
142 C        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
167 C
168 C           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
206 C
207 C                    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
215 C
216 C                    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
223 C
224 C                    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
278 C        If restarting at a break point, skip K1 + IGIGA*K2
279 C        Note that this is the number of numbers delivered to
280 C        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
300 C         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
307 C       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