Eureka.im

Simply the best answers.

Get 100$ bonus

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);
}





Show Form
No comments yet. Be the first to add a comment!