1 #include "randlib.h"
2 #include <stdio.h>
3 #include <stdlib.h>
4 static long Xm1,Xm2,Xa1,Xa2,Xcg1[32],Xcg2[32],Xa1w,Xa2w,Xig1[32],Xig2[32],
5 Xlg1[32],Xlg2[32],Xa1vw,Xa2vw;
6 static 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 static long g,i,ib1,ib2;
29 static 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 static long g;
74 static 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 static long ignlgi,curntg,k,s1,s2,z;
112 static 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 C 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 WGR, 12/19/00: replaced S10, S20, etc. with C blocks {} per
162 original paper.
163 **********************************************************************
164 */
165 {
166 #define numg 32L
167 extern void gsrgs(long getset,long *qvalue);
168 extern void gscgn(long getset,long *g);
169 extern long Xm1,Xm2,Xa1w,Xa2w,Xig1[],Xig2[],Xlg1[],Xlg2[],Xcg1[],Xcg2[];
170 static long g;
171 static long qrgnin;
172 /*
173 Abort unless random number generator initialized
174 */
175 gsrgs(0L,&qrgnin);
176 if (! qrgnin) {
177 fprintf(stderr,"%s\n",
178 " INITGN called before random number generator initialized -- abort!");
179 exit(1);
180 }
181
182 gscgn(0L,&g);
183 if(isdtyp == -1) { /* Initial seed */
184 *(Xlg1+g-1) = *(Xig1+g-1);
185 *(Xlg2+g-1) = *(Xig2+g-1);
186 }
187 else if (isdtyp == 0) { ; } /* Last seed */
188 else if (isdtyp == 1) { /* New seed */
189 *(Xlg1+g-1) = mltmod(Xa1w,*(Xlg1+g-1),Xm1);
190 *(Xlg2+g-1) = mltmod(Xa2w,*(Xlg2+g-1),Xm2);
191 } else {
192 fprintf(stderr,"%s\n","isdtyp not in range in INITGN");
193 exit(1);
194 }
195
196 *(Xcg1+g-1) = *(Xlg1+g-1);
197 *(Xcg2+g-1) = *(Xlg2+g-1);
198 #undef numg
199 }
200
inrgcm(void)201 void inrgcm(void)
202 /*
203 **********************************************************************
204 void inrgcm(void)
205 INitialize Random number Generator CoMmon
206 Function
207 Initializes common area for random number generator. This saves
208 the nuisance of a BLOCK DATA routine and the difficulty of
209 assuring that the routine is loaded with the other routines.
210 **********************************************************************
211 */
212 {
213 #define numg 32L
214 extern void gsrgs(long getset,long *qvalue);
215 extern long Xm1,Xm2,Xa1,Xa2,Xa1w,Xa2w,Xa1vw,Xa2vw;
216 extern long Xqanti[];
217 static long T1;
218 static long i;
219 /*
220 V=20; W=30;
221 A1W = MOD(A1**(2**W),M1) A2W = MOD(A2**(2**W),M2)
222 A1VW = MOD(A1**(2**(V+W)),M1) A2VW = MOD(A2**(2**(V+W)),M2)
223 If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed.
224 An efficient way to precompute a**(2*j) MOD m is to start with
225 a and square it j times modulo m using the function MLTMOD.
226 */
227 Xm1 = 2147483563L;
228 Xm2 = 2147483399L;
229 Xa1 = 40014L;
230 Xa2 = 40692L;
231 Xa1w = 1033780774L;
232 Xa2w = 1494757890L;
233 Xa1vw = 2082007225L;
234 Xa2vw = 784306273L;
235 for(i=0; i<numg; i++) *(Xqanti+i) = 0;
236 T1 = 1;
237 /*
238 Tell the world that common has been initialized
239 */
240 gsrgs(1L,&T1);
241 #undef numg
242 }
setall(long iseed1,long iseed2)243 void setall(long iseed1,long iseed2)
244 /*
245 **********************************************************************
246 void setall(long iseed1,long iseed2)
247 SET ALL random number generators
248 Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
249 initial seeds of the other generators are set accordingly, and
250 all generators states are set to these seeds.
251 This is a transcription from Pascal to Fortran of routine
252 Set_Initial_Seed from the paper
253 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
254 with Splitting Facilities." ACM Transactions on Mathematical
255 Software, 17:98-111 (1991)
256 Arguments
257 iseed1 -> First of two integer seeds
258 iseed2 -> Second of two integer seeds
259 **********************************************************************
260 */
261 {
262 #define numg 32L
263 extern void gsrgs(long getset,long *qvalue);
264 extern void gssst(long getset,long *qset);
265 extern void gscgn(long getset,long *g);
266 extern long Xm1,Xm2,Xa1vw,Xa2vw,Xig1[],Xig2[];
267 static long T1;
268 static long g,ocgn;
269 static long qrgnin;
270 T1 = 1;
271 /*
272 TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE
273 HAS BEEN CALLED.
274 */
275 gssst(1,&T1);
276 gscgn(0L,&ocgn);
277 /*
278 Initialize Common Block if Necessary
279 */
280 gsrgs(0L,&qrgnin);
281 if(!qrgnin) inrgcm();
282 *Xig1 = iseed1;
283 *Xig2 = iseed2;
284 initgn(-1L);
285 for(g=2; g<=numg; g++) {
286 *(Xig1+g-1) = mltmod(Xa1vw,*(Xig1+g-2),Xm1);
287 *(Xig2+g-1) = mltmod(Xa2vw,*(Xig2+g-2),Xm2);
288 gscgn(1L,&g);
289 initgn(-1L);
290 }
291 gscgn(1L,&ocgn);
292 #undef numg
293 }
setant(long qvalue)294 void setant(long qvalue)
295 /*
296 **********************************************************************
297 void setant(long qvalue)
298 SET ANTithetic
299 Sets whether the current generator produces antithetic values. If
300 X is the value normally returned from a uniform [0,1] random
301 number generator then 1 - X is the antithetic value. If X is the
302 value normally returned from a uniform [0,N] random number
303 generator then N - 1 - X is the antithetic value.
304 All generators are initialized to NOT generate antithetic values.
305 This is a transcription from Pascal to Fortran of routine
306 Set_Antithetic from the paper
307 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
308 with Splitting Facilities." ACM Transactions on Mathematical
309 Software, 17:98-111 (1991)
310 Arguments
311 qvalue -> nonzero if generator G is to generating antithetic
312 values, otherwise zero
313 **********************************************************************
314 */
315 {
316 #define numg 32L
317 extern void gsrgs(long getset,long *qvalue);
318 extern void gscgn(long getset,long *g);
319 extern long Xqanti[];
320 static long g;
321 static long qrgnin;
322 /*
323 Abort unless random number generator initialized
324 */
325 gsrgs(0L,&qrgnin);
326 if(qrgnin) goto S10;
327 fprintf(stderr,"%s\n",
328 " SETANT called before random number generator initialized -- abort!");
329 exit(1);
330 S10:
331 gscgn(0L,&g);
332 Xqanti[g-1] = qvalue;
333 #undef numg
334 }
setsd(long iseed1,long iseed2)335 void setsd(long iseed1,long iseed2)
336 /*
337 **********************************************************************
338 void setsd(long iseed1,long iseed2)
339 SET S-ee-D of current generator
340 Resets the initial seed of the current generator to ISEED1 and
341 ISEED2. The seeds of the other generators remain unchanged.
342 This is a transcription from Pascal to Fortran of routine
343 Set_Seed from the paper
344 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
345 with Splitting Facilities." ACM Transactions on Mathematical
346 Software, 17:98-111 (1991)
347 Arguments
348 iseed1 -> First integer seed
349 iseed2 -> Second integer seed
350 **********************************************************************
351 */
352 {
353 #define numg 32L
354 extern void gsrgs(long getset,long *qvalue);
355 extern void gscgn(long getset,long *g);
356 extern long Xig1[],Xig2[];
357 static long g;
358 static long qrgnin;
359 /*
360 Abort unless random number generator initialized
361 */
362 gsrgs(0L,&qrgnin);
363 if(qrgnin) goto S10;
364 fprintf(stderr,"%s\n",
365 " SETSD called before random number generator initialized -- abort!");
366 exit(1);
367 S10:
368 gscgn(0L,&g);
369 *(Xig1+g-1) = iseed1;
370 *(Xig2+g-1) = iseed2;
371 initgn(-1L);
372 #undef numg
373 }
374