Re: ran2 from numerical recipes
- To: mathgroup at smc.vnet.net
- Subject: [mg42444] Re: [mg42333] ran2 from numerical recipes
- From: Ezine <ezine at omegaconsultinggroup.com>
- Date: Tue, 8 Jul 2003 04:37:25 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
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. At 07:45 AM 7/1/2003, 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; >} > >-- >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