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