1: /* another library of funtions to support the rand 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 <math.h>
28: #include <stdlib.h>
29: #define abs(x) ((x) >= 0 ? (x) : -(x))
30: #define min(a,b) ((a) <= (b) ? (a) : (b))
31: #define max(a,b) ((a) >= (b) ? (a) : (b))
32:
33: /*
34: **********************************************************************
35: float genbet(float aa,float bb)
36: GeNerate BETa random deviate
37: Function
38: Returns a single random deviate from the beta distribution with
39: parameters A and B. The density of the beta is
40: x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1
41: Arguments
42: aa --> First parameter of the beta distribution
43:
44: bb --> Second parameter of the beta distribution
45:
46: Method
47: R. C. H. Cheng
48: Generating Beta Variatew with Nonintegral Shape Parameters
49: Communications of the ACM, 21:317-322 (1978)
50: (Algorithms BB and BC)
51: **********************************************************************
52: */
53: float genbet(float aa,float bb)
54: {
55: #define expmax 89.0
56: #define infnty 1.0E38
57: static float olda = -1.0;
58: static float oldb = -1.0;
59: static float genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;
60: static long qsame;
61:
62: qsame = olda == aa && oldb == bb;
63: if(qsame) goto S20;
64: if(!(aa <= 0.0 || bb <= 0.0)) goto S10;
65: puts(" AA or BB <= 0 in GENBET - Abort!");
66: printf(" AA: %16.6E BB %16.6E\n",aa,bb);
67: exit(1);
68: S10:
69: olda = aa;
70: oldb = bb;
71: S20:
72: if(!(min(aa,bb) > 1.0)) goto S100;
73: /*
74: Algorithm BB
75: Initialize
76: */
77: if(qsame) goto S30;
78: a = min(aa,bb);
79: b = max(aa,bb);
80: alpha = a+b;
81: beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
82: gamma = a+1.0/beta;
83: S30:
84: S40:
85: u1 = ranf();
86: /*
87: Step 1
88: */
89: u2 = ranf();
90: v = beta*log(u1/(1.0-u1));
91: if(!(v > expmax)) goto S50;
92: w = infnty;
93: goto S60;
94: S50:
95: w = a*exp(v);
96: S60:
97: z = pow(u1,2.0)*u2;
98: r = gamma*v-1.3862944;
99: s = a+r-w;
100: /*
101: Step 2
102: */
103: if(s+2.609438 >= 5.0*z) goto S70;
104: /*
105: Step 3
106: */
107: t = log(z);
108: if(s > t) goto S70;
109: /*
110: Step 4
111: */
112: if(r+alpha*log(alpha/(b+w)) < t) goto S40;
113: S70:
114: /*
115: Step 5
116: */
117: if(!(aa == a)) goto S80;
118: genbet = w/(b+w);
119: goto S90;
120: S80:
121: genbet = b/(b+w);
122: S90:
123: goto S230;
124: S100:
125: /*
126: Algorithm BC
127: Initialize
128: */
129: if(qsame) goto S110;
130: a = max(aa,bb);
131: b = min(aa,bb);
132: alpha = a+b;
133: beta = 1.0/b;
134: delta = 1.0+a-b;
135: k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778);
136: k2 = 0.25+(0.5+0.25/delta)*b;
137: S110:
138: S120:
139: u1 = ranf();
140: /*
141: Step 1
142: */
143: u2 = ranf();
144: if(u1 >= 0.5) goto S130;
145: /*
146: Step 2
147: */
148: y = u1*u2;
149: z = u1*y;
150: if(0.25*u2+z-y >= k1) goto S120;
151: goto S170;
152: S130:
153: /*
154: Step 3
155: */
156: z = pow(u1,2.0)*u2;
157: if(!(z <= 0.25)) goto S160;
158: v = beta*log(u1/(1.0-u1));
159: if(!(v > expmax)) goto S140;
160: w = infnty;
161: goto S150;
162: S140:
163: w = a*exp(v);
164: S150:
165: goto S200;
166: S160:
167: if(z >= k2) goto S120;
168: S170:
169: /*
170: Step 4
171: Step 5
172: */
173: v = beta*log(u1/(1.0-u1));
174: if(!(v > expmax)) goto S180;
175: w = infnty;
176: goto S190;
177: S180:
178: w = a*exp(v);
179: S190:
180: if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120;
181: S200:
182: /*
183: Step 6
184: */
185: if(!(a == aa)) goto S210;
186: genbet = w/(b+w);
187: goto S220;
188: S210:
189: genbet = b/(b+w);
190: S230:
191: S220:
192: return genbet;
193: #undef expmax
194: #undef infnty
195: }
196: float genchi(float df)
197: /*
198: **********************************************************************
199: float genchi(float df)
200: Generate random value of CHIsquare variable
201: Function
202: Generates random deviate from the distribution of a chisquare
203: with DF degrees of freedom random variable.
204: Arguments
205: df --> Degrees of freedom of the chisquare
206: (Must be positive)
207:
208: Method
209: Uses relation between chisquare and gamma.
210: **********************************************************************
211: */
212: {
213: static float genchi;
214:
215: if(!(df <= 0.0)) goto S10;
216: puts("DF <= 0 in GENCHI - ABORT");
217: printf("Value of DF: %16.6E\n",df);
218: exit(1);
219: S10:
220: genchi = 2.0*gengam(1.0,df/2.0);
221: return genchi;
222: }
223: float genexp(float av)
224: /*
225: **********************************************************************
226: float genexp(float av)
227: GENerate EXPonential random deviate
228: Function
229: Generates a single random deviate from an exponential
230: distribution with mean AV.
231: Arguments
232: av --> The mean of the exponential distribution from which
233: a random deviate is to be generated.
234: Method
235: Renames SEXPO from TOMS as slightly modified by BWB to use RANF
236: instead of SUNIF.
237: For details see:
238: Ahrens, J.H. and Dieter, U.
239: Computer Methods for Sampling From the
240: Exponential and Normal Distributions.
241: Comm. ACM, 15,10 (Oct. 1972), 873 - 882.
242: **********************************************************************
243: */
244: {
245: static float genexp;
246:
247: genexp = sexpo()*av;
248: return genexp;
249: }
250: float genf(float dfn,float dfd)
251: /*
252: **********************************************************************
253: float genf(float dfn,float dfd)
254: GENerate random deviate from the F distribution
255: Function
256: Generates a random deviate from the F (variance ratio)
257: distribution with DFN degrees of freedom in the numerator
258: and DFD degrees of freedom in the denominator.
259: Arguments
260: dfn --> Numerator degrees of freedom
261: (Must be positive)
262: dfd --> Denominator degrees of freedom
263: (Must be positive)
264: Method
265: Directly generates ratio of chisquare variates
266: **********************************************************************
267: */
268: {
269: static float genf,xden,xnum;
270:
271: if(!(dfn <= 0.0 || dfd <= 0.0)) goto S10;
272: puts("Degrees of freedom nonpositive in GENF - abort!");
273: printf("DFN value: %16.6EDFD value: %16.6E\n",dfn,dfd);
274: exit(1);
275: S10:
276: xnum = genchi(dfn)/dfn;
277: /*
278: GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD )
279: */
280: xden = genchi(dfd)/dfd;
281: if(!(xden <= 9.999999999998E-39*xnum)) goto S20;
282: puts(" GENF - generated numbers would cause overflow");
283: printf(" Numerator %16.6E Denominator %16.6E\n",xnum,xden);
284: puts(" GENF returning 1.0E38");
285: genf = 1.0E38;
286: goto S30;
287: S20:
288: genf = xnum/xden;
289: S30:
290: return genf;
291: }
292: float gengam(float a,float r)
293: /*
294: **********************************************************************
295: float gengam(float a,float r)
296: GENerates random deviates from GAMma distribution
297: Function
298: Generates random deviates from the gamma distribution whose
299: density is
300: (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X)
301: Arguments
302: a --> Location parameter of Gamma distribution
303: r --> Shape parameter of Gamma distribution
304: Method
305: Renames SGAMMA from TOMS as slightly modified by BWB to use RANF
306: instead of SUNIF.
307: For details see:
308: (Case R >= 1.0)
309: Ahrens, J.H. and Dieter, U.
310: Generating Gamma Variates by a
311: Modified Rejection Technique.
312: Comm. ACM, 25,1 (Jan. 1982), 47 - 54.
313: Algorithm GD
314: (Case 0.0 <= R <= 1.0)
315: Ahrens, J.H. and Dieter, U.
316: Computer Methods for Sampling from Gamma,
317: Beta, Poisson and Binomial Distributions.
318: Computing, 12 (1974), 223-246/
319: Adapted algorithm GS.
320: **********************************************************************
321: */
322: {
323: static float gengam;
324:
325: gengam = sgamma(r);
326: gengam /= a;
327: return gengam;
328: }
329: void genmn(float *parm,float *x,float *work)
330: /*
331: **********************************************************************
332: void genmn(float *parm,float *x,float *work)
333: GENerate Multivariate Normal random deviate
334: Arguments
335: parm --> Parameters needed to generate multivariate normal
336: deviates (MEANV and Cholesky decomposition of
337: COVM). Set by a previous call to SETGMN.
338: 1 : 1 - size of deviate, P
339: 2 : P + 1 - mean vector
340: P+2 : P*(P+3)/2 + 1 - upper half of cholesky
341: decomposition of cov matrix
342: x <-- Vector deviate generated.
343: work <--> Scratch array
344: Method
345: 1) Generate P independent standard normal deviates - Ei ~ N(0,1)
346: 2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM
347: 3) trans(A)E + MEANV ~ N(MEANV,COVM)
348: **********************************************************************
349: */
350: {
351: static long i,icount,j,p,D1,D2,D3,D4;
352: static float ae;
353:
354: p = (long) (*parm);
355: /*
356: Generate P independent normal deviates - WORK ~ N(0,1)
357: */
358: for(i=1; i<=p; i++) *(work+i-1) = snorm();
359: for(i=1,D3=1,D4=(p-i+D3)/D3; D4>0; D4--,i+=D3) {
360: /*
361: PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky
362: decomposition of the desired covariance matrix.
363: trans(A)(1,1) = PARM(P+2)
364: trans(A)(2,1) = PARM(P+3)
365: trans(A)(2,2) = PARM(P+2+P)
366: trans(A)(3,1) = PARM(P+4)
367: trans(A)(3,2) = PARM(P+3+P)
368: trans(A)(3,3) = PARM(P+2-1+2P) ...
369: trans(A)*WORK + MEANV ~ N(MEANV,COVM)
370: */
371: icount = 0;
372: ae = 0.0;
373: for(j=1,D1=1,D2=(i-j+D1)/D1; D2>0; D2--,j+=D1) {
374: icount += (j-1);
375: ae += (*(parm+i+(j-1)*p-icount+p)**(work+j-1));
376: }
377: *(x+i-1) = ae+*(parm+i);
378: }
379: }
380: float gennch(float df,float xnonc)
381: /*
382: **********************************************************************
383: float gennch(float df,float xnonc)
384: Generate random value of Noncentral CHIsquare variable
385: Function
386: Generates random deviate from the distribution of a noncentral
387: chisquare with DF degrees of freedom and noncentrality parameter
388: xnonc.
389: Arguments
390: df --> Degrees of freedom of the chisquare
391: (Must be > 1.0)
392: xnonc --> Noncentrality parameter of the chisquare
393: (Must be >= 0.0)
394: Method
395: Uses fact that noncentral chisquare is the sum of a chisquare
396: deviate with DF-1 degrees of freedom plus the square of a normal
397: deviate with mean XNONC and standard deviation 1.
398: **********************************************************************
399: */
400: {
401: static float gennch;
402:
403: if(!(df <= 1.0 || xnonc < 0.0)) goto S10;
404: puts("DF <= 1 or XNONC < 0 in GENNCH - ABORT");
405: printf("Value of DF: %16.6E Value of XNONC%16.6E\n",df,xnonc);
406: exit(1);
407: S10:
408: gennch = genchi(df-1.0)+pow(gennor(sqrt(xnonc),1.0),2.0);
409: return gennch;
410: }
411: float gennf(float dfn,float dfd,float xnonc)
412: /*
413: **********************************************************************
414: float gennf(float dfn,float dfd,float xnonc)
415: GENerate random deviate from the Noncentral F distribution
416: Function
417: Generates a random deviate from the noncentral F (variance ratio)
418: distribution with DFN degrees of freedom in the numerator, and DFD
419: degrees of freedom in the denominator, and noncentrality parameter
420: XNONC.
421: Arguments
422: dfn --> Numerator degrees of freedom
423: (Must be >= 1.0)
424: dfd --> Denominator degrees of freedom
425: (Must be positive)
426: xnonc --> Noncentrality parameter
427: (Must be nonnegative)
428: Method
429: Directly generates ratio of noncentral numerator chisquare variate
430: to central denominator chisquare variate.
431: **********************************************************************
432: */
433: {
434: static float gennf,xden,xnum;
435: static long qcond;
436:
437: qcond = dfn <= 1.0 || dfd <= 0.0 || xnonc < 0.0;
438: if(!qcond) goto S10;
439: puts("In GENNF - Either (1) Numerator DF <= 1.0 or");
440: puts("(2) Denominator DF < 0.0 or ");
441: puts("(3) Noncentrality parameter < 0.0");
442: printf("DFN value: %16.6EDFD value: %16.6EXNONC value: \n%16.6E\n",dfn,dfd,
443: xnonc);
444: exit(1);
445: S10:
446: xnum = gennch(dfn,xnonc)/dfn;
447: /*
448: GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD )
449: */
450: xden = genchi(dfd)/dfd;
451: if(!(xden <= 9.999999999998E-39*xnum)) goto S20;
452: puts(" GENNF - generated numbers would cause overflow");
453: printf(" Numerator %16.6E Denominator %16.6E\n",xnum,xden);
454: puts(" GENNF returning 1.0E38");
455: gennf = 1.0E38;
456: goto S30;
457: S20:
458: gennf = xnum/xden;
459: S30:
460: return gennf;
461: }
462:
463: float gennor(float av,float sd)
464: /*
465: **********************************************************************
466: float gennor(float av,float sd)
467: GENerate random deviate from a NORmal distribution
468: Function
469: Generates a single random deviate from a normal distribution
470: with mean, AV, and standard deviation, SD.
471: Arguments
472: av --> Mean of the normal distribution.
473: sd --> Standard deviation of the normal distribution.
474: Method
475: Renames SNORM from TOMS as slightly modified by BWB to use RANF
476: instead of SUNIF.
477: For details see:
478: Ahrens, J.H. and Dieter, U.
479: Extensions of Forsythe's Method for Random
480: Sampling from the Normal Distribution.
481: Math. Comput., 27,124 (Oct. 1973), 927 - 937.
482: **********************************************************************
483: */
484: {
485: float gennor;
486: float tmp_f;
487:
488: tmp_f = snorm();
489:
490: gennor = sd*tmp_f+av;
491: return (gennor);
492: }
493:
494: float capa_gennor(double *num_d, float av,float sd)
495: /*
496: **********************************************************************
497: float gennor(float av,float sd)
498: GENerate random deviate from a NORmal distribution
499: Function
500: Generates a single random deviate from a normal distribution
501: with mean, AV, and standard deviation, SD.
502: Arguments
503: av --> Mean of the normal distribution.
504: sd --> Standard deviation of the normal distribution.
505: Method
506: Renames SNORM from TOMS as slightly modified by BWB to use RANF
507: instead of SUNIF.
508: For details see:
509: Ahrens, J.H. and Dieter, U.
510: Extensions of Forsythe's Method for Random
511: Sampling from the Normal Distribution.
512: Math. Comput., 27,124 (Oct. 1973), 927 - 937.
513: **********************************************************************
514: */
515: {
516: float gen_num;
517: float tmp_f;
518:
519: tmp_f = snorm();
520:
521: gen_num = sd*tmp_f+av;
522: /* printf("SNORM()=%f,GENNOR()=%f,%f*%f+%f\n",tmp_f,gen_num,sd,tmp_f,av); */
523: *num_d = (double)gen_num;
524:
525: gen_num = (float)37.358341;
526: return (gen_num);
527: }
528:
529:
530: void genprm(long *iarray,int larray)
531: /*
532: **********************************************************************
533: void genprm(long *iarray,int larray)
534: GENerate random PeRMutation of iarray
535: Arguments
536: iarray <--> On output IARRAY is a random permutation of its
537: value on input
538: larray <--> Length of IARRAY
539: **********************************************************************
540: */
541: {
542: static long i,itmp,iwhich,D1,D2;
543:
544: for(i=1,D1=1,D2=(larray-i+D1)/D1; D2>0; D2--,i+=D1) {
545: iwhich = ignuin(i,larray);
546: itmp = *(iarray+iwhich-1);
547: *(iarray+iwhich-1) = *(iarray+i-1);
548: *(iarray+i-1) = itmp;
549: }
550: }
551: float genunf(float low,float high)
552: /*
553: **********************************************************************
554: float genunf(float low,float high)
555: GeNerate Uniform Real between LOW and HIGH
556: Function
557: Generates a real uniformly distributed between LOW and HIGH.
558: Arguments
559: low --> Low bound (exclusive) on real value to be generated
560: high --> High bound (exclusive) on real value to be generated
561: **********************************************************************
562: */
563: {
564: static float genunf;
565:
566: if(!(low > high)) goto S10;
567: printf("LOW > HIGH in GENUNF: LOW %16.6E HIGH: %16.6E\n",low,high);
568: puts("Abort");
569: exit(1);
570: S10:
571: genunf = low+(high-low)*ranf();
572: return genunf;
573: }
574: void gscgn(long getset,long *g)
575: /*
576: **********************************************************************
577: void gscgn(long getset,long *g)
578: Get/Set GeNerator
579: Gets or returns in G the number of the current generator
580: Arguments
581: getset --> 0 Get
582: 1 Set
583: g <-- Number of the current random number generator (1..32)
584: **********************************************************************
585: */
586: {
587: #define numg 32L
588: static long curntg = 1;
589: if(getset == 0) *g = curntg;
590: else {
591: if(*g < 0 || *g > numg) {
592: puts(" Generator number out of range in GSCGN");
593: exit(0);
594: }
595: curntg = *g;
596: }
597: #undef numg
598: }
599: void gsrgs(long getset,long *qvalue)
600: /*
601: **********************************************************************
602: void gsrgs(long getset,long *qvalue)
603: Get/Set Random Generators Set
604: Gets or sets whether random generators set (initialized).
605: Initially (data statement) state is not set
606: If getset is 1 state is set to qvalue
607: If getset is 0 state returned in qvalue
608: **********************************************************************
609: */
610: {
611: static long qinit = 0;
612:
613: if(getset == 0) *qvalue = qinit;
614: else qinit = *qvalue;
615: }
616: void gssst(long getset,long *qset)
617: /*
618: **********************************************************************
619: void gssst(long getset,long *qset)
620: Get or Set whether Seed is Set
621: Initialize to Seed not Set
622: If getset is 1 sets state to Seed Set
623: If getset is 0 returns T in qset if Seed Set
624: Else returns F in qset
625: **********************************************************************
626: */
627: {
628: static long qstate = 0;
629: if(getset != 0) qstate = 1;
630: else *qset = qstate;
631: }
632: long ignbin(long n,float pp)
633: /*
634: **********************************************************************
635: long ignbin(long n,float pp)
636: GENerate BINomial random deviate
637: Function
638: Generates a single random deviate from a binomial
639: distribution whose number of trials is N and whose
640: probability of an event in each trial is P.
641: Arguments
642: n --> The number of trials in the binomial distribution
643: from which a random deviate is to be generated.
644: p --> The probability of an event in each trial of the
645: binomial distribution from which a random deviate
646: is to be generated.
647: ignbin <-- A random deviate yielding the number of events
648: from N independent trials, each of which has
649: a probability of event P.
650: Method
651: This is algorithm BTPE from:
652: Kachitvichyanukul, V. and Schmeiser, B. W.
653: Binomial Random Variate Generation.
654: Communications of the ACM, 31, 2
655: (February, 1988) 216.
656: **********************************************************************
657: SUBROUTINE BTPEC(N,PP,ISEED,JX)
658: BINOMIAL RANDOM VARIATE GENERATOR
659: MEAN .LT. 30 -- INVERSE CDF
660: MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA
661: FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE
662: (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
663: THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
664: BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
665: BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE
666: RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
667: USABLE ALGORITHM.
668: REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
669: "BINOMIAL RANDOM VARIATE GENERATION,"
670: COMMUNICATIONS OF THE ACM, FORTHCOMING
671: WRITTEN: SEPTEMBER 1980.
672: LAST REVISED: MAY 1985, JULY 1987
673: REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER
674: GENERATOR
675: ARGUMENTS
676: N : NUMBER OF BERNOULLI TRIALS (INPUT)
677: PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
678: ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT)
679: JX: RANDOMLY GENERATED OBSERVATION (OUTPUT)
680: VARIABLES
681: PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
682: NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
683: XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
684: P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
685: FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
686: M: INTEGER VALUE OF THE CURRENT MODE
687: FM: FLOATING POINT VALUE OF THE CURRENT MODE
688: XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
689: P1: AREA OF THE TRIANGLE
690: C: HEIGHT OF THE PARALLELOGRAMS
691: XM: CENTER OF THE TRIANGLE
692: XL: LEFT END OF THE TRIANGLE
693: XR: RIGHT END OF THE TRIANGLE
694: AL: TEMPORARY VARIABLE
695: XLL: RATE FOR THE LEFT EXPONENTIAL TAIL
696: XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL
697: P2: AREA OF THE PARALLELOGRAMS
698: P3: AREA OF THE LEFT EXPONENTIAL TAIL
699: P4: AREA OF THE RIGHT EXPONENTIAL TAIL
700: U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
701: FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
702: FROM THE REGION
703: V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
704: (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
705: REJECT THE CANDIDATE VALUE
706: IX: INTEGER CANDIDATE VALUE
707: X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
708: AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
709: K: ABSOLUTE VALUE OF (IX-M)
710: F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
711: ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
712: ALSO USED IN THE INVERSE TRANSFORMATION
713: R: THE RATIO P/Q
714: G: CONSTANT USED IN CALCULATION OF PROBABILITY
715: MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
716: OF F WHEN IX IS GREATER THAN M
717: IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
718: CALCULATION OF F WHEN IX IS LESS THAN M
719: I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
720: AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
721: YNORM: LOGARITHM OF NORMAL BOUND
722: ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
723: X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
724: USED IN THE FINAL ACCEPT/REJECT TEST
725: QN: PROBABILITY OF NO SUCCESS IN N TRIALS
726: REMARK
727: IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
728: SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME
729: COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
730: ARE NOT INVOLVED.
731: ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
732: GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
733: TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
734: **********************************************************************
735: *****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
736: */
737: {
738: static float psave = -1.0;
739: static long nsave = -1;
740: static long ignbin,i,ix,ix1,k,m,mp,T1;
741: static float al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,q,qn,r,u,v,w,w2,x,x1,
742: x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2;
743:
744: if(pp != psave) goto S10;
745: if(n != nsave) goto S20;
746: if(xnp < 30.0) goto S150;
747: goto S30;
748: S10:
749: /*
750: *****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
751: */
752: psave = pp;
753: p = min(psave,1.0-psave);
754: q = 1.0-p;
755: S20:
756: xnp = n*p;
757: nsave = n;
758: if(xnp < 30.0) goto S140;
759: ffm = xnp+p;
760: m = ffm;
761: fm = m;
762: xnpq = xnp*q;
763: p1 = (long) (2.195*sqrt(xnpq)-4.6*q)+0.5;
764: xm = fm+0.5;
765: xl = xm-p1;
766: xr = xm+p1;
767: c = 0.134+20.5/(15.3+fm);
768: al = (ffm-xl)/(ffm-xl*p);
769: xll = al*(1.0+0.5*al);
770: al = (xr-ffm)/(xr*q);
771: xlr = al*(1.0+0.5*al);
772: p2 = p1*(1.0+c+c);
773: p3 = p2+c/xll;
774: p4 = p3+c/xlr;
775: S30:
776: /*
777: *****GENERATE VARIATE
778: */
779: u = ranf()*p4;
780: v = ranf();
781: /*
782: TRIANGULAR REGION
783: */
784: if(u > p1) goto S40;
785: ix = xm-p1*v+u;
786: goto S170;
787: S40:
788: /*
789: PARALLELOGRAM REGION
790: */
791: if(u > p2) goto S50;
792: x = xl+(u-p1)/c;
793: v = v*c+1.0-abs(xm-x)/p1;
794: if(v > 1.0 || v <= 0.0) goto S30;
795: ix = x;
796: goto S70;
797: S50:
798: /*
799: LEFT TAIL
800: */
801: if(u > p3) goto S60;
802: ix = xl+log(v)/xll;
803: if(ix < 0) goto S30;
804: v *= ((u-p2)*xll);
805: goto S70;
806: S60:
807: /*
808: RIGHT TAIL
809: */
810: ix = xr-log(v)/xlr;
811: if(ix > n) goto S30;
812: v *= ((u-p3)*xlr);
813: S70:
814: /*
815: *****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
816: */
817: k = abs(ix-m);
818: if(k > 20 && k < xnpq/2-1) goto S130;
819: /*
820: EXPLICIT EVALUATION
821: */
822: f = 1.0;
823: r = p/q;
824: g = (n+1)*r;
825: T1 = m-ix;
826: if(T1 < 0) goto S80;
827: else if(T1 == 0) goto S120;
828: else goto S100;
829: S80:
830: mp = m+1;
831: for(i=mp; i<=ix; i++) f *= (g/i-r);
832: goto S120;
833: S100:
834: ix1 = ix+1;
835: for(i=ix1; i<=m; i++) f /= (g/i-r);
836: S120:
837: if(v <= f) goto S170;
838: goto S30;
839: S130:
840: /*
841: SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X))
842: */
843: amaxp = k/xnpq*((k*(k/3.0+0.625)+0.1666666666666)/xnpq+0.5);
844: ynorm = -(k*k/(2.0*xnpq));
845: alv = log(v);
846: if(alv < ynorm-amaxp) goto S170;
847: if(alv > ynorm+amaxp) goto S30;
848: /*
849: STIRLING'S FORMULA TO MACHINE ACCURACY FOR
850: THE FINAL ACCEPTANCE/REJECTION TEST
851: */
852: x1 = ix+1.0;
853: f1 = fm+1.0;
854: z = n+1.0-fm;
855: w = n-ix+1.0;
856: z2 = z*z;
857: x2 = x1*x1;
858: f2 = f1*f1;
859: w2 = w*w;
860: if(alv <= xm*log(f1/x1)+(n-m+0.5)*log(z/w)+(ix-m)*log(w*p/(x1*q))+(13860.0-
861: (462.0-(132.0-(99.0-140.0/f2)/f2)/f2)/f2)/f1/166320.0+(13860.0-(462.0-
862: (132.0-(99.0-140.0/z2)/z2)/z2)/z2)/z/166320.0+(13860.0-(462.0-(132.0-
863: (99.0-140.0/x2)/x2)/x2)/x2)/x1/166320.0+(13860.0-(462.0-(132.0-(99.0
864: -140.0/w2)/w2)/w2)/w2)/w/166320.0) goto S170;
865: goto S30;
866: S140:
867: /*
868: INVERSE CDF LOGIC FOR MEAN LESS THAN 30
869: */
870: qn = pow(q,(double)n);
871: r = p/q;
872: g = r*(n+1);
873: S150:
874: ix = 0;
875: f = qn;
876: u = ranf();
877: S160:
878: if(u < f) goto S170;
879: if(ix > 110) goto S150;
880: u -= f;
881: ix += 1;
882: f *= (g/ix-r);
883: goto S160;
884: S170:
885: if(psave > 0.5) ix = n-ix;
886: ignbin = ix;
887: return ignbin;
888: }
889: long ignpoi(float mu)
890: /*
891: **********************************************************************
892: long ignpoi(float mu)
893: GENerate POIsson random deviate
894: Function
895: Generates a single random deviate from a Poisson
896: distribution with mean AV.
897: Arguments
898: av --> The mean of the Poisson distribution from which
899: a random deviate is to be generated.
900: genexp <-- The random deviate.
901: Method
902: Renames KPOIS from TOMS as slightly modified by BWB to use RANF
903: instead of SUNIF.
904: For details see:
905: Ahrens, J.H. and Dieter, U.
906: Computer Generation of Poisson Deviates
907: From Modified Normal Distributions.
908: ACM Trans. Math. Software, 8, 2
909: (June 1982),163-179
910: **********************************************************************
911: **********************************************************************
912:
913:
914: P O I S S O N DISTRIBUTION
915:
916:
917: **********************************************************************
918: **********************************************************************
919:
920: FOR DETAILS SEE:
921:
922: AHRENS, J.H. AND DIETER, U.
923: COMPUTER GENERATION OF POISSON DEVIATES
924: FROM MODIFIED NORMAL DISTRIBUTIONS.
925: ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
926:
927: (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
928:
929: **********************************************************************
930: INTEGER FUNCTION IGNPOI(IR,MU)
931: INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
932: MU=MEAN MU OF THE POISSON DISTRIBUTION
933: OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
934: MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
935: TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
936: COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
937: SEPARATION OF CASES A AND B
938: */
939: {
940: extern float fsign( float num, float sign );
941: static float a0 = -0.5;
942: static float a1 = 0.3333333;
943: static float a2 = -0.2500068;
944: static float a3 = 0.2000118;
945: static float a4 = -0.1661269;
946: static float a5 = 0.1421878;
947: static float a6 = -0.1384794;
948: static float a7 = 0.125006;
949: static float muold = 0.0;
950: static float muprev = 0.0;
951: static float fact[10] = {
952: 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0
953: };
954: static long ignpoi,j,k,kflag,l,m;
955: static float b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,fk,fx,fy,g,omega,p,p0,px,py,q,s,
956: t,u,v,x,xx,pp[35];
957:
958: if(mu == muprev) goto S10;
959: if(mu < 10.0) goto S120;
960: /*
961: C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
962: */
963: muprev = mu;
964: s = sqrt(mu);
965: d = 6.0*mu*mu;
966: /*
967: THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
968: PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
969: IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
970: */
971: l = (long) (mu-1.1484);
972: S10:
973: /*
974: STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
975: */
976: g = mu+s*snorm();
977: if(g < 0.0) goto S20;
978: ignpoi = (long) (g);
979: /*
980: STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
981: */
982: if(ignpoi >= l) return ignpoi;
983: /*
984: STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
985: */
986: fk = (float)ignpoi;
987: difmuk = mu-fk;
988: u = ranf();
989: if(d*u >= difmuk*difmuk*difmuk) return ignpoi;
990: S20:
991: /*
992: STEP P. PREPARATIONS FOR STEPS Q AND H.
993: (RECALCULATIONS OF PARAMETERS IF NECESSARY)
994: .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
995: THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
996: APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
997: C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
998: */
999: if(mu == muold) goto S30;
1000: muold = mu;
1001: omega = 0.3989423/s;
1002: b1 = 4.166667E-2/mu;
1003: b2 = 0.3*b1*b1;
1004: c3 = 0.1428571*b1*b2;
1005: c2 = b2-15.0*c3;
1006: c1 = b1-6.0*b2+45.0*c3;
1007: c0 = 1.0-b1+3.0*b2-15.0*c3;
1008: c = 0.1069/mu;
1009: S30:
1010: if(g < 0.0) goto S50;
1011: /*
1012: 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
1013: */
1014: kflag = 0;
1015: goto S70;
1016: S40:
1017: /*
1018: STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
1019: */
1020: if(fy-u*fy <= py*exp(px-fx)) return ignpoi;
1021: S50:
1022: /*
1023: STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
1024: DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
1025: (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
1026: */
1027: e = sexpo();
1028: u = ranf();
1029: u += (u-1.0);
1030: t = 1.8+fsign(e,u);
1031: if(t <= -0.6744) goto S50;
1032: ignpoi = (long) (mu+s*t);
1033: fk = (float)ignpoi;
1034: difmuk = mu-fk;
1035: /*
1036: 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
1037: */
1038: kflag = 1;
1039: goto S70;
1040: S60:
1041: /*
1042: STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
1043: */
1044: if(c*fabs(u) > py*exp(px+e)-fy*exp(fx+e)) goto S50;
1045: return ignpoi;
1046: S70:
1047: /*
1048: STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
1049: CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
1050: */
1051: if(ignpoi >= 10) goto S80;
1052: px = -mu;
1053: py = pow(mu,(double)ignpoi)/ *(fact+ignpoi);
1054: goto S110;
1055: S80:
1056: /*
1057: CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
1058: A0-A7 FOR ACCURACY WHEN ADVISABLE
1059: .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
1060: */
1061: del = 8.333333E-2/fk;
1062: del -= (4.8*del*del*del);
1063: v = difmuk/fk;
1064: if(fabs(v) <= 0.25) goto S90;
1065: px = fk*log(1.0+v)-difmuk-del;
1066: goto S100;
1067: S90:
1068: px = fk*v*v*(((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0)-del;
1069: S100:
1070: py = 0.3989423/sqrt(fk);
1071: S110:
1072: x = (0.5-difmuk)/s;
1073: xx = x*x;
1074: fx = -0.5*xx;
1075: fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0);
1076: if(kflag <= 0) goto S40;
1077: goto S60;
1078: S120:
1079: /*
1080: C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
1081: */
1082: muprev = 0.0;
1083: if(mu == muold) goto S130;
1084: muold = mu;
1085: m = max(1L,(long) (mu));
1086: l = 0;
1087: p = exp(-mu);
1088: q = p0 = p;
1089: S130:
1090: /*
1091: STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
1092: */
1093: u = ranf();
1094: ignpoi = 0;
1095: if(u <= p0) return ignpoi;
1096: /*
1097: STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
1098: PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
1099: (0.458=PP(9) FOR MU=10)
1100: */
1101: if(l == 0) goto S150;
1102: j = 1;
1103: if(u > 0.458) j = min(l,m);
1104: for(k=j; k<=l; k++) {
1105: if(u <= *(pp+k-1)) goto S180;
1106: }
1107: if(l == 35) goto S130;
1108: S150:
1109: /*
1110: STEP C. CREATION OF NEW POISSON PROBABILITIES P
1111: AND THEIR CUMULATIVES Q=PP(K)
1112: */
1113: l += 1;
1114: for(k=l; k<=35; k++) {
1115: p = p*mu/(float)k;
1116: q += p;
1117: *(pp+k-1) = q;
1118: if(u <= q) goto S170;
1119: }
1120: l = 35;
1121: goto S130;
1122: S170:
1123: l = k;
1124: S180:
1125: ignpoi = k;
1126: return ignpoi;
1127: }
1128: long ignuin(long low,long high)
1129: /*
1130: **********************************************************************
1131: long ignuin(long low,long high)
1132: GeNerate Uniform INteger
1133: Function
1134: Generates an integer uniformly distributed between LOW and HIGH.
1135: Arguments
1136: low --> Low bound (inclusive) on integer value to be generated
1137: high --> High bound (inclusive) on integer value to be generated
1138: Note
1139: If (HIGH-LOW) > 2,147,483,561 prints error message on * unit and
1140: stops the program.
1141: **********************************************************************
1142: IGNLGI generates integers between 1 and 2147483562
1143: MAXNUM is 1 less than maximum generable value
1144: */
1145: {
1146: #define maxnum 2147483561L
1147: static long ignuin,ign,maxnow,range,ranp1;
1148:
1149: if(!(low > high)) goto S10;
1150: puts(" low > high in ignuin - ABORT");
1151: exit(1);
1152:
1153: S10:
1154: range = high-low;
1155: if(!(range > maxnum)) goto S20;
1156: puts(" high - low too large in ignuin - ABORT");
1157: exit(1);
1158:
1159: S20:
1160: if(!(low == high)) goto S30;
1161: ignuin = low;
1162: return ignuin;
1163:
1164: S30:
1165: /*
1166: Number to be generated should be in range 0..RANGE
1167: Set MAXNOW so that the number of integers in 0..MAXNOW is an
1168: integral multiple of the number in 0..RANGE
1169: */
1170: ranp1 = range+1;
1171: maxnow = maxnum/ranp1*ranp1;
1172: S40:
1173: ign = ignlgi()-1;
1174: if(!(ign <= maxnow)) goto S50;
1175: ignuin = low+ign%ranp1;
1176: return ignuin;
1177: S50:
1178: goto S40;
1179: #undef maxnum
1180: #undef err1
1181: #undef err2
1182: }
1183: long lennob( char *str )
1184: /*
1185: Returns the length of str ignoring trailing blanks but not
1186: other white space.
1187: */
1188: {
1189: long i, i_nb;
1190:
1191: for (i=0, i_nb= -1L; *(str+i); i++)
1192: if ( *(str+i) != ' ' ) i_nb = i;
1193: return (i_nb+1);
1194: }
1195: long mltmod(long a,long s,long m)
1196: /*
1197: **********************************************************************
1198: long mltmod(long a,long s,long m)
1199: Returns (A*S) MOD M
1200: This is a transcription from Pascal to Fortran of routine
1201: MULtMod_Decompos from the paper
1202: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
1203: with Splitting Facilities." ACM Transactions on Mathematical
1204: Software, 17:98-111 (1991)
1205: Arguments
1206: a, s, m -->
1207: **********************************************************************
1208: */
1209: {
1210: #define h 32768L
1211: static long mltmod,a0,a1,k,p,q,qh,rh;
1212: /*
1213: H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
1214: machine. On a different machine recompute H
1215: */
1216: if(!(a <= 0 || a >= m || s <= 0 || s >= m)) goto S10;
1217: puts(" a, m, s out of order in mltmod - ABORT!");
1218: printf(" a = %12ld s = %12ld m = %12ld\n",a,s,m);
1219: puts(" mltmod requires: 0 < a < m; 0 < s < m");
1220: exit(1);
1221: S10:
1222: if(!(a < h)) goto S20;
1223: a0 = a;
1224: p = 0;
1225: goto S120;
1226: S20:
1227: a1 = a/h;
1228: a0 = a-h*a1;
1229: qh = m/h;
1230: rh = m-h*qh;
1231: if(!(a1 >= h)) goto S50;
1232: a1 -= h;
1233: k = s/qh;
1234: p = h*(s-k*qh)-k*rh;
1235: S30:
1236: if(!(p < 0)) goto S40;
1237: p += m;
1238: goto S30;
1239: S40:
1240: goto S60;
1241: S50:
1242: p = 0;
1243: S60:
1244: /*
1245: P = (A2*S*H)MOD M
1246: */
1247: if(!(a1 != 0)) goto S90;
1248: q = m/a1;
1249: k = s/q;
1250: p -= (k*(m-a1*q));
1251: if(p > 0) p -= m;
1252: p += (a1*(s-k*q));
1253: S70:
1254: if(!(p < 0)) goto S80;
1255: p += m;
1256: goto S70;
1257: S90:
1258: S80:
1259: k = p/qh;
1260: /*
1261: P = ((A2*H + A1)*S)MOD M
1262: */
1263: p = h*(p-k*qh)-k*rh;
1264: S100:
1265: if(!(p < 0)) goto S110;
1266: p += m;
1267: goto S100;
1268: S120:
1269: S110:
1270: if(!(a0 != 0)) goto S150;
1271: /*
1272: P = ((A2*H + A1)*H*S)MOD M
1273: */
1274: q = m/a0;
1275: k = s/q;
1276: p -= (k*(m-a0*q));
1277: if(p > 0) p -= m;
1278: p += (a0*(s-k*q));
1279: S130:
1280: if(!(p < 0)) goto S140;
1281: p += m;
1282: goto S130;
1283: S150:
1284: S140:
1285: mltmod = p;
1286: return mltmod;
1287: #undef h
1288: }
1289: void phrtsd(char* phrase,long *seed1,long *seed2)
1290: /*
1291: **********************************************************************
1292: void phrtsd(char* phrase,long *seed1,long *seed2)
1293: PHRase To SeeDs
1294:
1295: Function
1296:
1297: Uses a phrase (character string) to generate two seeds for the RGN
1298: random number generator.
1299: Arguments
1300: phrase --> Phrase to be used for random number generation
1301:
1302: seed1 <-- First seed for generator
1303:
1304: seed2 <-- Second seed for generator
1305:
1306: Note
1307:
1308: Trailing blanks are eliminated before the seeds are generated.
1309: Generated seed values will fall in the range 1..2^30
1310: (1..1,073,741,824)
1311: **********************************************************************
1312: */
1313: {
1314:
1315: static char table[] =
1316: "abcdefghijklmnopqrstuvwxyz\
1317: ABCDEFGHIJKLMNOPQRSTUVWXYZ\
1318: 0123456789\
1319: !@#$%^&*()_+[];:'\\\"<>?,./";
1320:
1321: long ix;
1322:
1323: static long twop30 = 1073741824L;
1324: static long shift[5] = {
1325: 1L,64L,4096L,262144L,16777216L
1326: };
1327: static long i,ichr,j,lphr,values[5];
1328: extern long lennob(char *str);
1329:
1330: *seed1 = 1234567890L;
1331: *seed2 = 123456789L;
1332: lphr = lennob(phrase);
1333: if(lphr < 1) return;
1334: for(i=0; i<=(lphr-1); i++) {
1335: for (ix=0; table[ix]; ix++) if (*(phrase+i) == table[ix]) break;
1336: if (!table[ix]) ix = 0;
1337: ichr = ix % 64;
1338: if(ichr == 0) ichr = 63;
1339: for(j=1; j<=5; j++) {
1340: *(values+j-1) = ichr-j;
1341: if(*(values+j-1) < 1) *(values+j-1) += 63;
1342: }
1343: for(j=1; j<=5; j++) {
1344: *seed1 = ( *seed1+*(shift+j-1)**(values+j-1) ) % twop30;
1345: *seed2 = ( *seed2+*(shift+j-1)**(values+6-j-1) ) % twop30;
1346: }
1347: }
1348: #undef twop30
1349: }
1350: float ranf(void)
1351: /*
1352: **********************************************************************
1353: float ranf(void)
1354: RANDom number generator as a Function
1355: Returns a random floating point number from a uniform distribution
1356: over 0 - 1 (endpoints of this interval are not returned) using the
1357: current generator
1358: This is a transcription from Pascal to Fortran of routine
1359: Uniform_01 from the paper
1360: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
1361: with Splitting Facilities." ACM Transactions on Mathematical
1362: Software, 17:98-111 (1991)
1363: **********************************************************************
1364: */
1365: {
1366: static float ranf;
1367: long tmp_l;
1368: double tmp_d;
1369: /*
1370: 4.656613057E-10 is 1/M1 M1 is set in a data statement in IGNLGI
1371: and is currently 2147483563. If M1 changes, change this also.
1372: */
1373: tmp_l = ignlgi();
1374: tmp_d = (double)tmp_l * (double)4.656613057E-10;
1375: ranf = (float)tmp_d;
1376: /* printf("RANF()=%f\n",ranf); */
1377: return ranf;
1378: }
1379: float capa_ranf(void)
1380: /*
1381: **********************************************************************
1382: float ranf(void)
1383: RANDom number generator as a Function
1384: Returns a random floating point number from a uniform distribution
1385: over 0 - 1 (endpoints of this interval are not returned) using the
1386: current generator
1387: This is a transcription from Pascal to Fortran of routine
1388: Uniform_01 from the paper
1389: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
1390: with Splitting Facilities." ACM Transactions on Mathematical
1391: Software, 17:98-111 (1991)
1392: **********************************************************************
1393: */
1394: {
1395: float ran_f;
1396: long my_ran;
1397: double my_doub;
1398: /*
1399: 4.656613057E-10 is 1/M1 M1 is set in a data statement in IGNLGI
1400: and is currently 2147483563. If M1 changes, change this also.
1401: */
1402: my_ran = ignlgi();
1403: /* printf("MY_ignlgi=%ld -- first time\n",my_ran); */
1404: /* ran_f = my_ran * 4.656613057E-10; */
1405:
1406: my_doub = (double)my_ran * (double)4.656613057E-10;
1407: printf("MY_ranf in double=%.15g -- first time\n",my_doub);
1408: ran_f = (float)my_doub;
1409: return (ran_f);
1410: }
1411:
1412: void setgmn(float *meanv,float *covm,long p,float *parm)
1413: /*
1414: **********************************************************************
1415: void setgmn(float *meanv,float *covm,long p,float *parm)
1416: SET Generate Multivariate Normal random deviate
1417: Function
1418: Places P, MEANV, and the Cholesky factoriztion of COVM
1419: in GENMN.
1420: Arguments
1421: meanv --> Mean vector of multivariate normal distribution.
1422: covm <--> (Input) Covariance matrix of the multivariate
1423: normal distribution
1424: (Output) Destroyed on output
1425: p --> Dimension of the normal, or length of MEANV.
1426: parm <-- Array of parameters needed to generate multivariate norma
1427: deviates (P, MEANV and Cholesky decomposition of
1428: COVM).
1429: 1 : 1 - P
1430: 2 : P + 1 - MEANV
1431: P+2 : P*(P+3)/2 + 1 - Cholesky decomposition of COVM
1432: Needed dimension is (p*(p+3)/2 + 1)
1433: **********************************************************************
1434: */
1435: {
1436: extern void spofa(float *a,long lda,long n,long *info);
1437: static long T1;
1438: static long i,icount,info,j,D2,D3,D4,D5;
1439: T1 = p*(p+3)/2+1;
1440: /*
1441: TEST THE INPUT
1442: */
1443: if(!(p <= 0)) goto S10;
1444: puts("P nonpositive in SETGMN");
1445: printf("Value of P: %12ld\n",p);
1446: exit(1);
1447: S10:
1448: *parm = p;
1449: /*
1450: PUT P AND MEANV INTO PARM
1451: */
1452: for(i=2,D2=1,D3=(p+1-i+D2)/D2; D3>0; D3--,i+=D2) *(parm+i-1) = *(meanv+i-2);
1453: /*
1454: Cholesky decomposition to find A s.t. trans(A)*(A) = COVM
1455: */
1456: spofa(covm,p,p,&info);
1457: if(!(info != 0)) goto S30;
1458: puts(" COVM not positive definite in SETGMN");
1459: exit(1);
1460: S30:
1461: icount = p+1;
1462: /*
1463: PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM
1464: COVM(1,1) = PARM(P+2)
1465: COVM(1,2) = PARM(P+3)
1466: :
1467: COVM(1,P) = PARM(2P+1)
1468: COVM(2,2) = PARM(2P+2) ...
1469: */
1470: for(i=1,D4=1,D5=(p-i+D4)/D4; D5>0; D5--,i+=D4) {
1471: for(j=i-1; j<p; j++) {
1472: icount += 1;
1473: *(parm+icount-1) = *(covm+i-1+j*p);
1474: }
1475: }
1476: }
1477: float sexpo(void)
1478: /*
1479: **********************************************************************
1480:
1481:
1482: (STANDARD-) E X P O N E N T I A L DISTRIBUTION
1483:
1484:
1485: **********************************************************************
1486: **********************************************************************
1487:
1488: FOR DETAILS SEE:
1489:
1490: AHRENS, J.H. AND DIETER, U.
1491: COMPUTER METHODS FOR SAMPLING FROM THE
1492: EXPONENTIAL AND NORMAL DISTRIBUTIONS.
1493: COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
1494:
1495: ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
1496: 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1497:
1498: Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1499: SUNIF. The argument IR thus goes away.
1500:
1501: **********************************************************************
1502: Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N
1503: (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
1504: */
1505: {
1506: static float q[8] = {
1507: 0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,1.0
1508: };
1509: static long i;
1510: static float sexpo,a,u,ustar,umin;
1511: static float *q1 = q;
1512: a = 0.0;
1513: u = ranf();
1514: goto S30;
1515: S20:
1516: a += *q1;
1517: S30:
1518: u += u;
1519: if(u <= 1.0) goto S20;
1520: u -= 1.0;
1521: if(u > *q1) goto S60;
1522: sexpo = a+u;
1523: return sexpo;
1524: S60:
1525: i = 1;
1526: ustar = ranf();
1527: umin = ustar;
1528: S70:
1529: ustar = ranf();
1530: if(ustar < umin) umin = ustar;
1531: i += 1;
1532: if(u > *(q+i-1)) goto S70;
1533: sexpo = a+umin**q1;
1534: return sexpo;
1535: }
1536: float sgamma(float a)
1537: /*
1538: **********************************************************************
1539:
1540:
1541: (STANDARD-) G A M M A DISTRIBUTION
1542:
1543:
1544: **********************************************************************
1545: **********************************************************************
1546:
1547: PARAMETER A >= 1.0 !
1548:
1549: **********************************************************************
1550:
1551: FOR DETAILS SEE:
1552:
1553: AHRENS, J.H. AND DIETER, U.
1554: GENERATING GAMMA VARIATES BY A
1555: MODIFIED REJECTION TECHNIQUE.
1556: COMM. ACM, 25,1 (JAN. 1982), 47 - 54.
1557:
1558: STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER
1559: (STRAIGHTFORWARD IMPLEMENTATION)
1560:
1561: Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1562: SUNIF. The argument IR thus goes away.
1563:
1564: **********************************************************************
1565:
1566: PARAMETER 0.0 < A < 1.0 !
1567:
1568: **********************************************************************
1569:
1570: FOR DETAILS SEE:
1571:
1572: AHRENS, J.H. AND DIETER, U.
1573: COMPUTER METHODS FOR SAMPLING FROM GAMMA,
1574: BETA, POISSON AND BINOMIAL DISTRIBUTIONS.
1575: COMPUTING, 12 (1974), 223 - 246.
1576:
1577: (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)
1578:
1579: **********************************************************************
1580: INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
1581: OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
1582: COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
1583: COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
1584: COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
1585: PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
1586: SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
1587: */
1588: {
1589: extern float fsign( float num, float sign );
1590: static float q1 = 4.166669E-2;
1591: static float q2 = 2.083148E-2;
1592: static float q3 = 8.01191E-3;
1593: static float q4 = 1.44121E-3;
1594: static float q5 = -7.388E-5;
1595: static float q6 = 2.4511E-4;
1596: static float q7 = 2.424E-4;
1597: static float a1 = 0.3333333;
1598: static float a2 = -0.250003;
1599: static float a3 = 0.2000062;
1600: static float a4 = -0.1662921;
1601: static float a5 = 0.1423657;
1602: static float a6 = -0.1367177;
1603: static float a7 = 0.1233795;
1604: static float e1 = 1.0;
1605: static float e2 = 0.4999897;
1606: static float e3 = 0.166829;
1607: static float e4 = 4.07753E-2;
1608: static float e5 = 1.0293E-2;
1609: static float aa = 0.0;
1610: static float aaa = 0.0;
1611: static float sqrt32 = 5.656854;
1612: static float sgamma,s2,s,d,t,x,u,r,q0,b,si,c,v,q,e,w,p;
1613: if(a == aa) goto S10;
1614: if(a < 1.0) goto S120;
1615: /*
1616: STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED
1617: */
1618: aa = a;
1619: s2 = a-0.5;
1620: s = sqrt(s2);
1621: d = sqrt32-12.0*s;
1622: S10:
1623: /*
1624: STEP 2: T=STANDARD NORMAL DEVIATE,
1625: X=(S,1/2)-NORMAL DEVIATE.
1626: IMMEDIATE ACCEPTANCE (I)
1627: */
1628: t = snorm();
1629: x = s+0.5*t;
1630: sgamma = x*x;
1631: if(t >= 0.0) return sgamma;
1632: /*
1633: STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
1634: */
1635: u = ranf();
1636: if(d*u <= t*t*t) return sgamma;
1637: /*
1638: STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
1639: */
1640: if(a == aaa) goto S40;
1641: aaa = a;
1642: r = 1.0/ a;
1643: q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
1644: /*
1645: APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
1646: THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
1647: C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
1648: */
1649: if(a <= 3.686) goto S30;
1650: if(a <= 13.022) goto S20;
1651: /*
1652: CASE 3: A .GT. 13.022
1653: */
1654: b = 1.77;
1655: si = 0.75;
1656: c = 0.1515/s;
1657: goto S40;
1658: S20:
1659: /*
1660: CASE 2: 3.686 .LT. A .LE. 13.022
1661: */
1662: b = 1.654+7.6E-3*s2;
1663: si = 1.68/s+0.275;
1664: c = 6.2E-2/s+2.4E-2;
1665: goto S40;
1666: S30:
1667: /*
1668: CASE 1: A .LE. 3.686
1669: */
1670: b = 0.463+s-0.178*s2;
1671: si = 1.235;
1672: c = 0.195/s-7.9E-2+1.6E-2*s;
1673: S40:
1674: /*
1675: STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE
1676: */
1677: if(x <= 0.0) goto S70;
1678: /*
1679: STEP 6: CALCULATION OF V AND QUOTIENT Q
1680: */
1681: v = t/(s+s);
1682: if(fabs(v) <= 0.25) goto S50;
1683: q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
1684: goto S60;
1685: S50:
1686: q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
1687: S60:
1688: /*
1689: STEP 7: QUOTIENT ACCEPTANCE (Q)
1690: */
1691: if(log(1.0-u) <= q) return sgamma;
1692: S70:
1693: /*
1694: STEP 8: E=STANDARD EXPONENTIAL DEVIATE
1695: U= 0,1 -UNIFORM DEVIATE
1696: T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
1697: */
1698: e = sexpo();
1699: u = ranf();
1700: u += (u-1.0);
1701: t = b+fsign(si*e,u);
1702: /*
1703: STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719
1704: */
1705: if(t < -0.7187449) goto S70;
1706: /*
1707: STEP 10: CALCULATION OF V AND QUOTIENT Q
1708: */
1709: v = t/(s+s);
1710: if(fabs(v) <= 0.25) goto S80;
1711: q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
1712: goto S90;
1713: S80:
1714: q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
1715: S90:
1716: /*
1717: STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
1718: */
1719: if(q <= 0.0) goto S70;
1720: if(q <= 0.5) goto S100;
1721: w = exp(q)-1.0;
1722: goto S110;
1723: S100:
1724: w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
1725: S110:
1726: /*
1727: IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
1728: */
1729: if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
1730: x = s+0.5*t;
1731: sgamma = x*x;
1732: return sgamma;
1733: S120:
1734: /*
1735: ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.))
1736: */
1737: aa = 0.0;
1738: b = 1.0+0.3678794*a;
1739: S130:
1740: p = b*ranf();
1741: if(p >= 1.0) goto S140;
1742: sgamma = exp(log(p)/ a);
1743: if(sexpo() < sgamma) goto S130;
1744: return sgamma;
1745: S140:
1746: sgamma = -log((b-p)/ a);
1747: if(sexpo() < (1.0-a)*log(sgamma)) goto S130;
1748: return sgamma;
1749: }
1750: float snorm(void)
1751: /*
1752: **********************************************************************
1753:
1754:
1755: (STANDARD-) N O R M A L DISTRIBUTION
1756:
1757:
1758: **********************************************************************
1759: **********************************************************************
1760:
1761: FOR DETAILS SEE:
1762:
1763: AHRENS, J.H. AND DIETER, U.
1764: EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM
1765: SAMPLING FROM THE NORMAL DISTRIBUTION.
1766: MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.
1767:
1768: ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'
1769: (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1770:
1771: Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1772: SUNIF. The argument IR thus goes away.
1773:
1774: **********************************************************************
1775: THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
1776: H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
1777: */
1778: {
1779: static float a[32] = {
1780: 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
1781: 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
1782: 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
1783: 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
1784: 1.862732,2.153875
1785: };
1786: static float d[31] = {
1787: 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
1788: 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
1789: 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
1790: 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
1791: };
1792: static float t[31] = {
1793: 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
1794: 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
1795: 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
1796: 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
1797: 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
1798: };
1799: static float h[31] = {
1800: 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
1801: 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
1802: 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
1803: 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
1804: 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
1805: };
1806: static long i;
1807: static float snorm,u,s,ustar,aa,w,y,tt;
1808: u = ranf();
1809: s = 0.0;
1810: if(u > 0.5) s = 1.0;
1811: u += (u-s);
1812: u = 32.0*u;
1813: i = (long) (u);
1814: if(i == 32) i = 31;
1815: if(i == 0) goto S100;
1816: /*
1817: START CENTER
1818: */
1819: ustar = u-(float)i;
1820: aa = *(a+i-1);
1821: S40:
1822: if(ustar <= *(t+i-1)) goto S60;
1823: w = (ustar-*(t+i-1))**(h+i-1);
1824: S50:
1825: /*
1826: EXIT (BOTH CASES)
1827: */
1828: y = aa+w;
1829: snorm = y;
1830: if(s == 1.0) snorm = -y;
1831: return snorm;
1832: S60:
1833: /*
1834: CENTER CONTINUED
1835: */
1836: u = ranf();
1837: w = u*(*(a+i)-aa);
1838: tt = (0.5*w+aa)*w;
1839: goto S80;
1840: S70:
1841: tt = u;
1842: ustar = ranf();
1843: S80:
1844: if(ustar > tt) goto S50;
1845: u = ranf();
1846: if(ustar >= u) goto S70;
1847: ustar = ranf();
1848: goto S40;
1849: S100:
1850: /*
1851: START TAIL
1852: */
1853: i = 6;
1854: aa = *(a+31);
1855: goto S120;
1856: S110:
1857: aa += *(d+i-1);
1858: i += 1;
1859: S120:
1860: u += u;
1861: if(u < 1.0) goto S110;
1862: u -= 1.0;
1863: S140:
1864: w = u**(d+i-1);
1865: tt = (0.5*w+aa)*w;
1866: goto S160;
1867: S150:
1868: tt = u;
1869: S160:
1870: ustar = ranf();
1871: if(ustar > tt) goto S50;
1872: u = ranf();
1873: if(ustar >= u) goto S150;
1874: u = ranf();
1875: goto S140;
1876: }
1877: float fsign( float num, float sign )
1878: /* Transfers sign of argument sign to argument num */
1879: {
1880: if ( ( sign>0.0f && num<0.0f ) || ( sign<0.0f && num>0.0f ) )
1881: return -num;
1882: else return num;
1883: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/back.gif){kind=icon}
Return to ranlib.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [LON-CAPA] / capa / capa51 / pProj