294 Chapter 7.Random Numbers 0 reject http://www.nr Permission is read able files .com or call 1-800-872- (including this one) granted fori internet 7423 1988-1992 by Cambridge University Press. from NUMERICAL RECIPES IN Figure 7.3.2.Rejection method as applied to an integer-valued distribution.The method is performed (North America to any server computer. t users to make one paper THE on the step function shown as a dashed line,yielding a real-valued deviate.This deviate is rounded down to the next lower integer,which is output. 是 ART #include <math.h> #define PI3.141592654 st st Programs for their float poidev(float xm,long *idum) Returns as a floating-point number an integer value that is a random deviate drawn from a Poisson distribution of mean xm,using ran1(idum)as a source of uniform random deviates. to dir float gammln(float xx); float ran1(long *idum); ectcustser static float sq,alxm,g,oldm=(-1.0); oldm is a flag for whether xm has changed 18881920 OF SCIENTIFIC COMPUTING(ISBN float em,t,yi since last call. 1f(xm<12.0){ Use direct method. Further reproduction. Numerical Recipes 12-521 if(xm！=oldm){ oldm=xm; If xm is new,compute the exponential. 43108 g=exp(-xm) 2 em=-1: t=1.0; (outside do Instead of adding exponential deviates it is equiv- Software. ++em; alent to multiply uniform deviates.We never North t *=ran1(idum); actually have to take the log.merely com- while (t g); pare to the pre-computed exponential else Use rejection method. 1f(xm！=oldm)[ If xm has changed since the last call,then pre- visit website machine oldm=xm; compute some functions that occur below. sq=sqrt(2.0*xm); alxm=log(xm); g=xm*alxm-gammln(xm+1.0); The function gammin is the natural log of the gamma function,as given in $6.1. do do y is a deviate from a Lorentzian comparison func- y=tan(PI*ran1(idum)); tion.294 Chapter 7. Random Numbers Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 0 12345 accept reject in 1 Figure 7.3.2. Rejection method as applied to an integer-valued distribution. The method is performed on the step function shown as a dashed line, yielding a real-valued deviate. This deviate is rounded down to the next lower integer, which is output. #include <math.h> #define PI 3.141592654 float poidev(float xm, long *idum) Returns as a floating-point number an integer value that is a random deviate drawn from a Poisson distribution of mean xm, using ran1(idum) as a source of uniform random deviates. { float gammln(float xx); float ran1(long *idum); static float sq,alxm,g,oldm=(-1.0); oldm is a flag for whether xm has changed float em,t,y; since last call. if (xm < 12.0) { Use direct method. if (xm != oldm) { oldm=xm; g=exp(-xm); If xm is new, compute the exponential. } em = -1; t=1.0; do { Instead of adding exponential deviates it is equiv￾alent to multiply uniform deviates. We never actually have to take the log, merely com￾pare to the pre-computed exponential. ++em; t *= ran1(idum); } while (t > g); } else { Use rejection method. if (xm != oldm) { If xm has changed since the last call, then pre￾oldm=xm; compute some functions that occur below. sq=sqrt(2.0*xm); alxm=log(xm); g=xm*alxm-gammln(xm+1.0); The function gammln is the natural log of the gamma function, as given in §6.1. } do { do { y is a deviate from a Lorentzian comparison func￾y=tan(PI*ran1(idum)); tion