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