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 followsdistribution = type of random number distribution:random_uniform = uniformly distributedrandom_gaussian = gaussian distribution mean=0 std dev = 1.0random_binomial = binomial distributedrandom_poisson = poisson distributedrandom_exponential =exponential distributionNOTE: See NUMERICAL RECIPES for background materialNOTE: 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"#endifreal 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 _STANDALONEint main(){#elseDEFINE_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 RNMXreal 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);}