mirror of
https://bitbucket.org/cosmicvoids/vide_public.git
synced 2025-07-05 15:51:12 +00:00
126 lines
3 KiB
C
126 lines
3 KiB
C
#include "error.h"
|
|
#include "randoms.h"
|
|
|
|
/* ran2 from Numerical Recipes 2nd ed. modified to conform to our interface */
|
|
/* You win $1000 from Press et al. if it fails non-trivially */
|
|
|
|
#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 NDIV (1+IMM1/NTAB)
|
|
#define EPS 1.2e-7
|
|
#define RNMX (1.0-EPS)
|
|
|
|
void
|
|
ran_init(int seed, ran_state *s)
|
|
{
|
|
int j;
|
|
long k;
|
|
|
|
if (seed < 1)
|
|
Error("Bad seed in ran_init2 (%d)\n", seed);
|
|
s->idum = s->idum2 = seed;
|
|
s->next_norml_ok = 0;
|
|
for (j=NTAB+7;j>=0;j--) {
|
|
k=(s->idum)/IQ1;
|
|
s->idum=IA1*(s->idum-k*IQ1)-k*IR1;
|
|
if (s->idum < 0) s->idum += IM1;
|
|
if (j < NTAB) s->iv[j] = s->idum;
|
|
}
|
|
s->iy=s->iv[0];
|
|
s->did_init = IQ1;
|
|
}
|
|
|
|
float
|
|
uniform_rand(ran_state *s)
|
|
{
|
|
int j;
|
|
long k;
|
|
float temp;
|
|
|
|
if (s->did_init != IQ1)
|
|
Error("You forgot to call ran_init2\n");
|
|
k=(s->idum)/IQ1;
|
|
s->idum=IA1*(s->idum-k*IQ1)-k*IR1;
|
|
if (s->idum < 0) s->idum += IM1;
|
|
k=s->idum2/IQ2;
|
|
s->idum2=IA2*(s->idum2-k*IQ2)-k*IR2;
|
|
if (s->idum2 < 0) s->idum2 += IM2;
|
|
j=s->iy/NDIV;
|
|
s->iy=s->iv[j]-s->idum2;
|
|
s->iv[j] = s->idum;
|
|
if (s->iy < 1) s->iy += IMM1;
|
|
if ((temp=AM*s->iy) > RNMX) return RNMX;
|
|
else return temp;
|
|
}
|
|
|
|
float normal_rand(ran_state *st)
|
|
/*
|
|
This is the Polar method for normal distributions, as described on or near
|
|
page 104 of Knuth, Semi-numerical Algorithms. To quote Knuth, "The polar
|
|
method is quite slow, but it has essentially perfect accuracy, and it is very
|
|
easy to write a program for the polar method..." 'nuf said. Algorithm due
|
|
to Box, Muller and Marsaglia.
|
|
*/
|
|
{
|
|
float v1, v2; /* uniformly distributed on [-1, 1) */
|
|
float s; /* radius of a point pulled from a uniform circle */
|
|
float foo; /* A useful intermediate value. */
|
|
double log(double), sqrt(double);
|
|
|
|
if(st->next_norml_ok){
|
|
st->next_norml_ok = 0;
|
|
return st->next_norml;
|
|
}
|
|
|
|
do{
|
|
v1 = 2.0F * uniform_rand(st) - 1.0F;
|
|
v2 = 2.0F * uniform_rand(st) - 1.0F;
|
|
s = v1*v1 + v2*v2;
|
|
} while(s >= 1.0F);
|
|
foo = sqrt( -2.0F * log(s)/s);
|
|
st->next_norml_ok = 1;
|
|
st->next_norml = v1*foo;
|
|
return v2*foo;
|
|
}
|
|
|
|
/* Return a uniform point in a ndim-sphere by rejection. */
|
|
/* If ndim is large (bigger than 4 or so), this becomes very inefficient */
|
|
/* Return the radius-squard of the result. */
|
|
float sphere_rand(ran_state *st, int ndim, float *x)
|
|
{
|
|
int k;
|
|
float rsqx;
|
|
|
|
do {
|
|
rsqx = (float)0.0;
|
|
for (k = 0; k < ndim; k++) {
|
|
x[k] = uniform_rand(st)*2.0F - 1.0F; /* a pt in (-1,1) */
|
|
rsqx += x[k] * x[k];
|
|
}
|
|
} while (rsqx > (float)1.0);
|
|
return rsqx;
|
|
}
|
|
|
|
/* Return a uniform point in a ndim-cube */
|
|
/* Return the radius-squard of the result. */
|
|
float cube_rand(ran_state *st, int ndim, float *x)
|
|
{
|
|
int k;
|
|
float rsqx;
|
|
|
|
rsqx = (float)0.0;
|
|
for (k = 0; k < ndim; k++) {
|
|
x[k] = uniform_rand(st)*2.0F - 1.0F; /* a pt in (-1,1) */
|
|
rsqx += x[k] * x[k];
|
|
}
|
|
return rsqx;
|
|
}
|
|
|