generate random numbers
UDF to generate many types of random number distributions: uniform, gaussian, poisson, binomial and exponential. A UDF has been written as a DEFINE_ON_DEMAND routine. The mainf unction is used to call a function called udf_random(int distribution, long* seed, real (*fcn)(long* theSeed, real r1, int i1) Arguments are as follows distribution = type of random number distribution: random_uniform = uniformly distributed random_gaussian = gaussian distribution mean=0 std dev = 1.0 random_binomial = binomial distributed random_poisson = poisson distributed random_exponential =exponential distribution NOTE: See NUMERICAL RECIPES for background material NOTE: All algorithms are very slightly modified from the original NUMERICAL RECIPES routines. /* * * UDF to generate random numbers from ran2 of Numerical Recipes in C * * by * * William Wangard, Ph.D. * ww@fluent.com * Fluent, Inc. * * All rights reserved, 2003. * * Fluent, Inc. does not guarantee the rigor of this User defined function. * * ran1 = function for quick random numbers * ran2 = function for slower but better random numbers * * r1 = real parameter used in binomial and poisson distributions * i1 = int parameter used in binomail distribution * * seed = long used for seeding computations * * */ #if _STANDALONE #include "math.h" #include "stdio.h" #define real double #define M_PI 3.141592654 #else #include "udf.h" #endif real udf_random(int method, long *seed, real (*fcn)(long* theSeed), real r1, int i1); real ran1(long *seed); real ran2(long *seed); real gasdev(long *seed, real (*fcn)(long* theSeed)); real bnldev(real pp, int n, long *seed, real (*fcn)(long* theSeed)); real poidev(real xm, long *seed, real (*fcn)(long* theSeed)); real expdev(long *seed, real (*fcn)(long* theSeed)); real gammln(real xx); static enum { random_uniform=0, random_gaussian, random_binomial, random_poisson, random_exponential } random_types; #if _STANDALONE int main() { #else DEFINE_ON_DEMAND(generate_random_number) { #endif long seed=0; real x; int i; int method = random_gaussian; int i1 = 16; reaol r1 = 0.1; /* Generate 100 uniformly distributed random numbers with ran1 */ for(i=0; i<100; i++) { x = udf_random(random_uniform,&seed,&ran1,r1,i1); printf("%fn",x); } /* Generate 100 gaussian random numbers with ran2 */ for(i=0; i<100; i++) { x = udf_random(random_gaussian,&seed,&ran2,r1,i1); printf("%fn",x); } /* Generate 100 poisson random numbers with ran2 */ for(i=0; i<100; i++) { x = udf_random(random_poisson,&seed,&ran2,r1,i1); printf("%fn",x); } return 0; } real udf_random(int method, long *seed, real (*fcn)(long* theSeed), real r1, int i1) { real x = 0.; switch(method){ case random_uniform: x = (*fcn)(seed); break; case random_gaussian: x = gasdev(seed, fcn); break; case random_binomial: x = bnldev(r1, i1, seed, fcn); break; case random_poisson: x = poidev(r1, seed, fcn); break; case random_exponential: x = expdev(seed, fcn); break; default: break; } return x; } real bnldev(real pp, int n, long *seed, real (*fcn)(long* theSeed)) { int j; static int nold=(-1); real am,em,g,angle,p,bnl,sq,t,y; static real pold=(-1.0),pc,plog,pclog,en,oldg; p=(pp <= 0.5 ? pp : 1.0-pp); am=n*p; if (n < 25) { bnl=0.0; for (j=1;j<=n;j++) if ((*fcn)(seed) < p) ++bnl; } else if (am < 1.0) { g=exp(-am); t=1.0; for (j=0;j<=n;j++) { t *= (*fcn)(seed); if (t < g) break; } bnl=(j <= n ? j : n); } else { if (n != nold) { en=n; oldg=gammln(en+1.0); nold=n; } if (p != pold) { pc=1.0-p; plog=log(p); pclog=log(pc); pold=p; } sq=sqrt(2.0*am*pc); do { do { angle=M_PI*(*fcn)(seed); y=tan(angle); em=sq*y+am; } while (em < 0.0 || em >= (en+1.0)); em=floor(em); t=1.2*sq*(1.0+y*y)*exp(oldg-gammln(em+1.0) -gammln(en-em+1.0)+em*plog+(en-em)*pclog); } while ((*fcn)(seed) > t); bnl=em; } if (p != pp) bnl=n-bnl; return bnl; } real expdev(long *seed, real (*fcn)(long* theSeed)) { real dum; do dum=(*fcn)(seed); while (dum == 0.0); return -log(dum); } real gasdev(long *seed, real (*fcn)(long* theSeed)) { static int iset=0; static real gset; real fac,rsq,v1,v2; if (iset == 0) { do { v1=2.0*(*fcn)(seed)-1.0; v2=2.0*(*fcn)(seed)-1.0; rsq=v1*v1+v2*v2; } while (rsq >= 1.0 || rsq == 0.0); fac=sqrt(-2.0*log(rsq)/rsq); gset=v1*fac; iset=1; return v2*fac; } else { iset=0; return gset; } } real poidev(real xm, long *seed, real (*fcn)(long* theSeed)) { static real sq,alxm,g,oldm=(-1.0); real em,t,y; if (xm < 12.0) { if (xm != oldm) { oldm=xm; g=exp(-xm); } em = -1; t=1.0; do { ++em; t *= (*fcn)(seed); } while (t > g); } else { if (xm != oldm) { oldm=xm; sq=sqrt(2.0*xm); alxm=log(xm); g=xm*alxm-gammln(xm+1.0); } do { do { y=tan(M_PI*(*fcn)(seed)); em=sq*y+xm; } while (em < 0.0); em=floor(em); t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g); } while ((*fcn)(seed) > t); } return em; } #define IM1 2147483563 #define IM2 2147483399 #define AM (1.0/IM1) #define IMM1 (IM1-1) #define IA1 40014 #define IA2 40692 #define IQ1 53668 #define IQ2 52774 #define IR1 12211 #define IR2 3791 #define NTAB 32 #define NDIV (1+IMM1/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) real ran2(long *seed) { int j; long k; static long seed2=123456789; static long iy=0; static long iv[NTAB]; real temp; if (*seed <= 0) { if (-(*seed) < 1) *seed=1; else *seed = -(*seed); seed2=(*seed); for (j=NTAB+7;j>=0;j--) { k=(*seed)/IQ1; *seed=IA1*(*seed-k*IQ1)-k*IR1; if (*seed < 0) *seed += IM1; if (j < NTAB) iv[j] = *seed; } iy=iv[0]; } k=(*seed)/IQ1; *seed=IA1*(*seed-k*IQ1)-k*IR1; if (*seed < 0) *seed += IM1; k=seed2/IQ2; seed2=IA2*(seed2-k*IQ2)-k*IR2; if (seed2 < 0) seed2 += IM2; j=iy/NDIV; iy=iv[j]-seed2; iv[j] = *seed; if (iy < 1) iy += IMM1; if ((temp=AM*iy) > RNMX) { return RNMX;} else { return temp;} } #undef IM1 #undef IM2 #undef AM #undef IMM1 #undef IA1 #undef IA2 #undef IQ1 #undef IQ2 #undef IR1 #undef IR2 #undef NTAB #undef NDIV #undef EPS #undef RNMX #define IA 16807 #define IM 2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define NTAB 32 #define NDIV (1+(IM-1)/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) real ran1(long *seed) { int j; long k; static long iy=0; static long iv[NTAB]; real temp; if (*seed <= 0 || !iy) { if (-(*seed) < 1) *seed=1; else *seed = -(*seed); for (j=NTAB+7;j>=0;j--) { k=(*seed)/IQ; *seed=IA*(*seed-k*IQ)-IR*k; if (*seed < 0) *seed += IM; if (j < NTAB) iv[j] = *seed; } iy=iv[0]; } k=(*seed)/IQ; *seed=IA*(*seed-k*IQ)-IR*k; if (*seed < 0) *seed += IM; j=iy/NDIV; iy=iv[j]; iv[j] = *seed; if ((temp=AM*iy) > RNMX) return RNMX; else return temp; } #undef IA #undef IM #undef AM #undef IQ #undef IR #undef NTAB #undef NDIV #undef EPS #undef RNMX real gammln(real xx) { double x,y,tmp,ser; static double cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5}; int j; y=x=xx; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.000000000190015; for (j=0;j<=5;j++) ser += cof[j]/++y; return -tmp+log(2.5066282746310005*ser/x); } |
||
![]()
|