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