]>
Commit | Line | Data |
---|---|---|
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 |