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