1 #include "randlib.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4 long Xm1,Xm2,Xa1,Xa2,Xcg1[32],Xcg2[32],Xa1w,Xa2w,Xig1[32],Xig2[32],
5 Xlg1[32],Xlg2[32],Xa1vw,Xa2vw;
6 long Xqanti[32];
advnst(long k)7 void advnst(long k)
8 /*
9 **********************************************************************
10 void advnst(long k)
11 ADV-a-N-ce ST-ate
12 Advances the state of the current generator by 2^K values and
13 resets the initial seed to that value.
14 This is a transcription from Pascal to Fortran of routine
15 Advance_State from the paper
16 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
17 with Splitting Facilities." ACM Transactions on Mathematical
18 Software, 17:98-111 (1991)
19 Arguments
20 k -> The generator is advanced by2^K values
21 **********************************************************************
22 */
23 {
24 #define numg 32L
25 extern void gsrgs(long getset,long *qvalue);
26 extern void gscgn(long getset,long *g);
27 extern long Xm1,Xm2,Xa1,Xa2,Xcg1[],Xcg2[];
28 long g,i,ib1,ib2;
29 long qrgnin;
30 /*
31 Abort unless random number generator initialized
32 */
33 gsrgs(0L,&qrgnin);
34 if(qrgnin) goto S10;
35 fputs(" ADVNST called before random generator initialized - ABORT\n",
36 stderr);
37 exit(1);
38 S10:
39 gscgn(0L,&g);
40 ib1 = Xa1;
41 ib2 = Xa2;
42 for(i=1; i<=k; i++) {
43 ib1 = mltmod(ib1,ib1,Xm1);
44 ib2 = mltmod(ib2,ib2,Xm2);
45 }
46 setsd(mltmod(ib1,*(Xcg1+g-1),Xm1),mltmod(ib2,*(Xcg2+g-1),Xm2));
47 /*
48 NOW, IB1 = A1**K AND IB2 = A2**K
49 */
50 #undef numg
51 }
getsd(long * iseed1,long * iseed2)52 void getsd(long *iseed1,long *iseed2)
53 /*
54 **********************************************************************
55 void getsd(long *iseed1,long *iseed2)
56 GET SeeD
57 Returns the value of two integer seeds of the current generator
58 This is a transcription from Pascal to Fortran of routine
59 Get_State from the paper
60 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
61 with Splitting Facilities." ACM Transactions on Mathematical
62 Software, 17:98-111 (1991)
63 Arguments
64 iseed1 <- First integer seed of generator G
65 iseed2 <- Second integer seed of generator G
66 **********************************************************************
67 */
68 {
69 #define numg 32L
70 extern void gsrgs(long getset,long *qvalue);
71 extern void gscgn(long getset,long *g);
72 extern long Xcg1[],Xcg2[];
73 long g;
74 long qrgnin;
75 /*
76 Abort unless random number generator initialized
77 */
78 gsrgs(0L,&qrgnin);
79 if(qrgnin) goto S10;
80 fprintf(stderr,"%s\n",
81 " GETSD called before random number generator initialized -- abort!");
82 exit(0);
83 S10:
84 gscgn(0L,&g);
85 *iseed1 = *(Xcg1+g-1);
86 *iseed2 = *(Xcg2+g-1);
87 #undef numg
88 }
ignlgi(void)89 long ignlgi(void)
90 /*
91 **********************************************************************
92 long ignlgi(void)
93 GeNerate LarGe Integer
94 Returns a random integer following a uniform distribution over
95 (1, 2147483562) using the current generator.
96 This is a transcription from Pascal to Fortran of routine
97 Random from the paper
98 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
99 with Splitting Facilities." ACM Transactions on Mathematical
100 Software, 17:98-111 (1991)
101 **********************************************************************
102 */
103 {
104 #define numg 32L
105 extern void gsrgs(long getset,long *qvalue);
106 extern void gssst(long getset,long *qset);
107 extern void gscgn(long getset,long *g);
108 extern void inrgcm(void);
109 extern long Xm1,Xm2,Xa1,Xa2,Xcg1[],Xcg2[];
110 extern long Xqanti[];
111 long ignlgi,curntg,k,s1,s2,z;
112 long qqssd,qrgnin;
113 /*
114 IF THE RANDOM NUMBER PACKAGE HAS NOT BEEN INITIALIZED YET, DO SO.
115 IT CAN BE INITIALIZED IN ONE OF TWO WAYS : 1) THE FIRST CALL TO
116 THIS ROUTINE 2) A CALL TO SETALL.
117 */
118 gsrgs(0L,&qrgnin);
119 if(!qrgnin) inrgcm();
120 gssst(0,&qqssd);
121 if(!qqssd) setall(1234567890L,123456789L);
122 /*
123 Get Current Generator
124 */
125 gscgn(0L,&curntg);
126 s1 = *(Xcg1+curntg-1);
127 s2 = *(Xcg2+curntg-1);
128 k = s1/53668L;
129 s1 = Xa1*(s1-k*53668L)-k*12211;
130 if(s1 < 0) s1 += Xm1;
131 k = s2/52774L;
132 s2 = Xa2*(s2-k*52774L)-k*3791;
133 if(s2 < 0) s2 += Xm2;
134 *(Xcg1+curntg-1) = s1;
135 *(Xcg2+curntg-1) = s2;
136 z = s1-s2;
137 if(z < 1) z += (Xm1-1);
138 if(*(Xqanti+curntg-1)) z = Xm1-z;
139 ignlgi = z;
140 return ignlgi;
141 #undef numg
142 }
initgn(long isdtyp)143 void initgn(long isdtyp)
144 /*
145 **********************************************************************
146 void initgn(long isdtyp)
147 INIT-ialize current G-e-N-erator
148 Reinitializes the state of the current generator
149 This is a transcription from Pascal to Fortran of routine
150 Init_Generator from the paper
151 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
152 with Splitting Facilities." ACM Transactions on Mathematical
153 Software, 17:98-111 (1991)
154 Arguments
155 isdtyp -> The state to which the generator is to be set
156 isdtyp = -1 => sets the seeds to their initial value
157 isdtyp = 0 => sets the seeds to the first value of
158 the current block
159 isdtyp = 1 => sets the seeds to the first value of
160 the next block
161 **********************************************************************
162 */
163 {
164 #define numg 32L
165 extern void gsrgs(long getset,long *qvalue);
166 extern void gscgn(long getset,long *g);
167 extern long Xm1,Xm2,Xa1w,Xa2w,Xig1[],Xig2[],Xlg1[],Xlg2[],Xcg1[],Xcg2[];
168 long g;
169 long qrgnin;
170 /*
171 Abort unless random number generator initialized
172 */
173 gsrgs(0L,&qrgnin);
174 if(qrgnin) goto S10;
175 fprintf(stderr,"%s\n",
176 " INITGN called before random number generator initialized -- abort!");
177 exit(1);
178 S10:
179 gscgn(0L,&g);
180 if(-1 != isdtyp) goto S20;
181 *(Xlg1+g-1) = *(Xig1+g-1);
182 *(Xlg2+g-1) = *(Xig2+g-1);
183 goto S50;
184 S20:
185 if(0 != isdtyp) goto S30;
186 goto S50;
187 S30:
188 /*
189 do nothing
190 */
191 if(1 != isdtyp) goto S40;
192 *(Xlg1+g-1) = mltmod(Xa1w,*(Xlg1+g-1),Xm1);
193 *(Xlg2+g-1) = mltmod(Xa2w,*(Xlg2+g-1),Xm2);
194 goto S50;
195 S40:
196 fprintf(stderr,"%s\n","isdtyp not in range in INITGN");
197 exit(1);
198 S50:
199 *(Xcg1+g-1) = *(Xlg1+g-1);
200 *(Xcg2+g-1) = *(Xlg2+g-1);
201 #undef numg
202 }
inrgcm(void)203 void inrgcm(void)
204 /*
205 **********************************************************************
206 void inrgcm(void)
207 INitialize Random number Generator CoMmon
208 Function
209 Initializes common area for random number generator. This saves
210 the nuisance of a BLOCK DATA routine and the difficulty of
211 assuring that the routine is loaded with the other routines.
212 **********************************************************************
213 */
214 {
215 #define numg 32L
216 extern void gsrgs(long getset,long *qvalue);
217 extern long Xm1,Xm2,Xa1,Xa2,Xa1w,Xa2w,Xa1vw,Xa2vw;
218 extern long Xqanti[];
219 long T1;
220 long i;
221 /*
222 V=20; W=30;
223 A1W = MOD(A1**(2**W),M1) A2W = MOD(A2**(2**W),M2)
224 A1VW = MOD(A1**(2**(V+W)),M1) A2VW = MOD(A2**(2**(V+W)),M2)
225 If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed.
226 An efficient way to precompute a**(2*j) MOD m is to start with
227 a and square it j times modulo m using the function MLTMOD.
228 */
229 Xm1 = 2147483563L;
230 Xm2 = 2147483399L;
231 Xa1 = 40014L;
232 Xa2 = 40692L;
233 Xa1w = 1033780774L;
234 Xa2w = 1494757890L;
235 Xa1vw = 2082007225L;
236 Xa2vw = 784306273L;
237 for(i=0; i<numg; i++) *(Xqanti+i) = 0;
238 T1 = 1;
239 /*
240 Tell the world that common has been initialized
241 */
242 gsrgs(1L,&T1);
243 #undef numg
244 }
setall(long iseed1,long iseed2)245 void setall(long iseed1,long iseed2)
246 /*
247 **********************************************************************
248 void setall(long iseed1,long iseed2)
249 SET ALL random number generators
250 Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
251 initial seeds of the other generators are set accordingly, and
252 all generators states are set to these seeds.
253 This is a transcription from Pascal to Fortran of routine
254 Set_Initial_Seed from the paper
255 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
256 with Splitting Facilities." ACM Transactions on Mathematical
257 Software, 17:98-111 (1991)
258 Arguments
259 iseed1 -> First of two integer seeds
260 iseed2 -> Second of two integer seeds
261 **********************************************************************
262 */
263 {
264 #define numg 32L
265 extern void gsrgs(long getset,long *qvalue);
266 extern void gssst(long getset,long *qset);
267 extern void gscgn(long getset,long *g);
268 extern long Xm1,Xm2,Xa1vw,Xa2vw,Xig1[],Xig2[];
269 long T1;
270 long g,ocgn;
271 long qrgnin;
272 T1 = 1;
273 /*
274 TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE
275 HAS BEEN CALLED.
276 */
277 gssst(1,&T1);
278 gscgn(0L,&ocgn);
279 /*
280 Initialize Common Block if Necessary
281 */
282 gsrgs(0L,&qrgnin);
283 if(!qrgnin) inrgcm();
284 *Xig1 = iseed1;
285 *Xig2 = iseed2;
286 initgn(-1L);
287 for(g=2; g<=numg; g++) {
288 *(Xig1+g-1) = mltmod(Xa1vw,*(Xig1+g-2),Xm1);
289 *(Xig2+g-1) = mltmod(Xa2vw,*(Xig2+g-2),Xm2);
290 gscgn(1L,&g);
291 initgn(-1L);
292 }
293 gscgn(1L,&ocgn);
294 #undef numg
295 }
setant(long qvalue)296 void setant(long qvalue)
297 /*
298 **********************************************************************
299 void setant(long qvalue)
300 SET ANTithetic
301 Sets whether the current generator produces antithetic values. If
302 X is the value normally returned from a uniform [0,1] random
303 number generator then 1 - X is the antithetic value. If X is the
304 value normally returned from a uniform [0,N] random number
305 generator then N - 1 - X is the antithetic value.
306 All generators are initialized to NOT generate antithetic values.
307 This is a transcription from Pascal to Fortran of routine
308 Set_Antithetic from the paper
309 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
310 with Splitting Facilities." ACM Transactions on Mathematical
311 Software, 17:98-111 (1991)
312 Arguments
313 qvalue -> nonzero if generator G is to generating antithetic
314 values, otherwise zero
315 **********************************************************************
316 */
317 {
318 #define numg 32L
319 extern void gsrgs(long getset,long *qvalue);
320 extern void gscgn(long getset,long *g);
321 extern long Xqanti[];
322 long g;
323 long qrgnin;
324 /*
325 Abort unless random number generator initialized
326 */
327 gsrgs(0L,&qrgnin);
328 if(qrgnin) goto S10;
329 fprintf(stderr,"%s\n",
330 " SETANT called before random number generator initialized -- abort!");
331 exit(1);
332 S10:
333 gscgn(0L,&g);
334 Xqanti[g-1] = qvalue;
335 #undef numg
336 }
setsd(long iseed1,long iseed2)337 void setsd(long iseed1,long iseed2)
338 /*
339 **********************************************************************
340 void setsd(long iseed1,long iseed2)
341 SET S-ee-D of current generator
342 Resets the initial seed of the current generator to ISEED1 and
343 ISEED2. The seeds of the other generators remain unchanged.
344 This is a transcription from Pascal to Fortran of routine
345 Set_Seed from the paper
346 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
347 with Splitting Facilities." ACM Transactions on Mathematical
348 Software, 17:98-111 (1991)
349 Arguments
350 iseed1 -> First integer seed
351 iseed2 -> Second integer seed
352 **********************************************************************
353 */
354 {
355 #define numg 32L
356 extern void gsrgs(long getset,long *qvalue);
357 extern void gscgn(long getset,long *g);
358 extern long Xig1[],Xig2[];
359 long g;
360 long qrgnin;
361 /*
362 Abort unless random number generator initialized
363 */
364 gsrgs(0L,&qrgnin);
365 if(qrgnin) goto S10;
366 fprintf(stderr,"%s\n",
367 " SETSD called before random number generator initialized -- abort!");
368 exit(1);
369 S10:
370 gscgn(0L,&g);
371 *(Xig1+g-1) = iseed1;
372 *(Xig2+g-1) = iseed2;
373 initgn(-1L);
374 #undef numg
375 }
376