c28b8f8d |
1 | C* |
2 | C* $Id$ |
3 | C* |
c28b8f8d |
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 |