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