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