ran2 from numerical recipes

[Daniel Reeves]

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

--
Daniel Reeves -- http://ai.eecs.umich.edu/people/dreeves/

"The complete lack of evidence is the surest sign that the conspiracy
is working."

[Lawrence A. Walker Jr.]

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

--
Lawrence A. Walker Jr.
http://www.kingshonor.com

[Ezine]

There are 2 approaches you can take:

1) Add a MathLink template to call the C function from Mathematica. This is
discussed in the current issue of The Mathematica Ezine

http://omegaconsultinggroup.com/Services/ezv2i06.html

2) Convert the C code to Mathematica code. I don't see any functions that
don't translate directly. Simply convert the syntax.

if (test) command; becomes If[test, command].
if (test) cmd1 else cmd2; becomes If[test, cmd1, cmd2].
for(init, test, mod) cmd; becomes For[init, test, mod, cmd].
return value; becomes Return[value].

It should be pretty straight-forward.

>Daniel Reeves -- http://ai.eecs.umich.edu/people/dreeves/
>
>"The complete lack of evidence is the surest sign that the conspiracy
>is working."

--------------------------------------------------------------
Omega Consulting
"The final answer to your Mathematica needs"
http://omegaconsultinggroup.com