[Bio] / Numeric / Src / ranlib.c Repository:
ViewVC logotype

Annotation of /Numeric/Src/ranlib.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3