Re: ran2 from numerical recipes
- To: mathgroup at smc.vnet.net
- Subject: [mg42370] Re: ran2 from numerical recipes
- From: "Lawrence A. Walker Jr." <lwalker701 at earthlink.net>
- Date: Thu, 3 Jul 2003 06:10:28 -0400 (EDT)
- References: <bds0gc$626$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi Daniel,
I implemented ran1() from Numerical recipes.
First here is the C function that calls ran1() that is capable of
returning random numbers between a given range ...
/**********************************************************/
double random(double xmin, double xmax, long *idum)
{
double r;
r=ran1(idum);
return (double)(r*xmax + xmin - r*xmin);
}
/**********************************************************/
Next, some more C code that makes it easier to implement in Mathlink.
/**********************************************************/
double ml_random(double xmin, double xmax)
{
double r;
long idum;
idum=-time(NULL);
r=random(xmin, xmax, &idum);
return r;
}
And finally the template ....
:Begin:
:Function: ml_random
:Pattern: Random2[xmin_,xmax_]
:Arguments: {xmin, xmax}
:ArgumentTypes: {Real, Real}
:ReturnType: Real
:End:
Lawrence
Daniel Reeves wrote:
> I don't suppose anyone has implemented ran2() from Numerical Recipes in
> Mathematica? I know the built-in Random[] is good but I need a Mathematica
> implementation of ran2 in order to verify some C code. Come to think of
> it, a C implementation of Random[Real] would also do the trick.
>
> Here's ran2 in C:
>
> #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)
>
> float ran2(long *idum)
> {
> int j;
> long k;
> static long idum2=123456789;
> static long iy=0;
> static long iv[NTAB];
> float temp;
>
> if (*idum <= 0) {
> if (-(*idum) < 1) *idum=1;
> else *idum = -(*idum);
> idum2=(*idum);
> for (j=NTAB+7;j>=0;j--) {
> k=(*idum)/IQ1;
> *idum=IA1*(*idum-k*IQ1)-k*IR1;
> if (*idum < 0) *idum += IM1;
> if (j < NTAB) iv[j] = *idum;
> }
> iy=iv[0];
> }
> k=(*idum)/IQ1;
> *idum=IA1*(*idum-k*IQ1)-k*IR1;
> if (*idum < 0) *idum += IM1;
> k=idum2/IQ2;
> idum2=IA2*(idum2-k*IQ2)-k*IR2;
> if (idum2 < 0) idum2 += IM2;
> j=iy/NDIV;
> iy=iv[j]-idum2;
> iv[j] = *idum;
> if (iy < 1) iy += IMM1;
> if ((temp=AM*iy) > RNMX) return RNMX;
> else return temp;
> }
>
--
Lawrence A. Walker Jr.
http://www.kingshonor.com