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