Histogram ranges changed to cut off saturation peak and noise
[u/mrichter/AliRoot.git] / HBTP / ranlux.f
1 C*
2 C* $Id$
3 C*
4 C#include "gen/pilot.h"
5       SUBROUTINE RANLUX(RVEC,LENV)
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       LOGICAL NOTYET
60       DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
61       DATA I24,J24,CARRY/24,10,0./
62 C                               default
63 C  Luxury Level   0     1     2   *3*    4
64       DATA NDSKIP/0,   24,   73,  199,  365 /
65 Corresponds to p=24    48    97   223   389
66 C     time factor 1     2     3     6    10   on slow workstation
67 C                 1    1.5    2     3     5   on fast mainframe
68 C
69 C  NOTYET is .TRUE. if no initialization has been performed yet.
70 C              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
103 C
104 C          The Generator proper: "Subtract-with-borrow",
105 C          as proposed by Marsaglia and Zaman,
106 C          Florida State University, March, 1989
107 C
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
120 C  small numbers (with less than 12 "significant" bits) are "padded".
121       IF (UNI .LT. TWOM12)  THEN
122          RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
123 C        and zero is forbidden in case someone takes a logarithm
124          IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24
125       ENDIF
126 C        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
151 C
152 C           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
190 C
191 C                    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
199 C
200 C                    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
207 C
208 C                    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
262 C        If restarting at a break point, skip K1 + IGIGA*K2
263 C        Note that this is the number of numbers delivered to
264 C        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
284 C         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
291 C       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