Dummy implementation of Alarms
[u/mrichter/AliRoot.git] / MEVSIM / ranlux2.F
1 C*
2 C* $Id$
3 C*
4 C#include "gen/pilot.h"
5       SUBROUTINE RANLUX2(RVEC,LENV,Input_seed)
6 C         Subtract-and-borrow random number generator proposed by
7 C         Marsaglia and Zaman, implemented by F. James with the name
8 C         RCARRY in 1991, and later improved by Martin Luescher
9 C         in 1993 to produce "Luxury Pseudorandom Numbers".
10 C     Fortran 77 coded by F. James, 1993
11 C
12 C   LUXURY LEVELS.
13 C   ------ ------      The available luxury levels are:
14 C
15 C  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
16 C           and Zaman, very long period, but fails many tests.
17 C  level 1  (p=48): considerable improvement in quality over level 0,
18 C           now passes the gap test, but still fails spectral test.
19 C  level 2  (p=97): passes all known tests, but theoretically still
20 C           defective.
21 C  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
22 C           correlations have very small chance of being observed.
23 C  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
24 C
25 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
26 C!!!  Calling sequences for RANLUX:                                  ++
27 C!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
28 C!!!                   32-bit random floating point numbers between  ++
29 C!!!                   zero (not included) and one (also not incl.). ++
30 C!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
31 C!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
32 C!!!               which is integer between zero and MAXLEV, or if   ++
33 C!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
34 C!!!               should be set to zero unless restarting at a break++ 
35 C!!!               point given by output of RLUXAT (see RLUXAT).     ++
36 C!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
37 C!!!               which can be used to restart the RANLUX generator ++
38 C!!!               at the current point by calling RLUXGO.  K1 and K2++
39 C!!!               specify how many numbers were generated since the ++
40 C!!!               initialization with LUX and INT.  The restarting  ++
41 C!!!               skips over  K1+K2*E9   numbers, so it can be long.++
42 C!!!   A more efficient but less convenient way of restarting is by: ++
43 C!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
44 C!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
45 C!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
46 C!!!                 32-bit integer seeds, to be used for restarting ++
47 C!!!      ISVEC must be dimensioned 25 in the calling program        ++
48 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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./
63 CCC  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
69 C                               default
70 C  Luxury Level   0     1     2   *3*    4
71       DATA NDSKIP/0,   24,   73,  199,  365 /
72 Corresponds to p=24    48    97   223   389
73 C     time factor 1     2     3     6    10   on slow workstation
74 C                 1    1.5    2     3     5   on fast mainframe
75 C
76 C  NOTYET is .TRUE. if no initialization has been performed yet.
77 C              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
110 C
111 C          The Generator proper: "Subtract-with-borrow",
112 C          as proposed by Marsaglia and Zaman,
113 C          Florida State University, March, 1989
114 C
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
127 C  small numbers (with less than 12 "significant" bits) are "padded".
128       IF (UNI .LT. TWOM12)  THEN
129          RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
130 C        and zero is forbidden in case someone takes a logarithm
131          IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24
132       ENDIF
133 C        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
158 C
159 C           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
197 C
198 C                    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
206 C
207 C                    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
214 C
215 C                    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
269 C        If restarting at a break point, skip K1 + IGIGA*K2
270 C        Note that this is the number of numbers delivered to
271 C        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
291 C         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
298 C       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