trapping/docs/csource/random.c

185 lines
3.9 KiB
C

/* $Id$ */
/*
* random.c
*
* written by Sebastian Voecking <seb.voeck@uni-muenster.de>
*
* For details see random.h
*/
#include "random.h"
#include <stdlib.h>
static int method = RANDOM_STDLIB;
static int rand_seed=-1; // for random calculation
static double get_random_1_cw();
static void subrn(double* u, int len);
static double random_james();
void random_set_method(RandomMethode meth)
{
method = meth;
}
double random_get()
{
switch(method) {
case RANDOM_STDLIB:
return ((double)rand())/RAND_MAX;
case RANDOM_CW:
return get_random_1_cw();
case RANDOM_JAMES:
return random_james();
default:
return 0;
}
}
void random_seed(int seed)
{
switch(method) {
case RANDOM_STDLIB:
srand(seed);
break;
case RANDOM_CW:
rand_seed = seed;
break;
case RANDOM_JAMES:
break;
}
}
double get_random_1_cw()
// function gives a random value between 0 and 1.
// NO imput value needed.
// numerical rec. version provided by Ch. Weinheimer
{
#define IA 16807
#define IM 2147483647
#define AM (1./IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define RNMX (1.-1.2e-7)
int j, k;
static int iy=0, iv[NTAB];
double temp;
if (rand_seed<=0 || !iy){
if (-rand_seed<1) rand_seed=1;
else rand_seed=-rand_seed;
for (j=NTAB+7; j>=0; j--){
k=rand_seed/IQ;
rand_seed=IA*(rand_seed-k*IQ)-IR*k;
if (rand_seed<0) rand_seed+=IM;
if (j<NTAB) iv[j]=rand_seed;
}
iy=iv[0];
}
k=rand_seed/IQ;
rand_seed=IA*(rand_seed-k*IQ)-IR*k;
if (rand_seed<0) rand_seed+=IM;
j=iy/NDIV;
iy=iv[j];
iv[j]=rand_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 RNMX
}
void subrn(double *u,int len)
{
// This subroutine computes random numbers u[1],...,u[len]
// in the (0,1) interval. It uses the 0<IJKLRANDOM<900000000
// integer as initialization seed.
// In the calling program the dimension
// of the u[] vector should be larger than len (the u[0] value is
// not used).
// For each IJKLRANDOM
// numbers the program computes completely independent random number
// sequences (see: F. James, Comp. Phys. Comm. 60 (1990) 329, sec. 3.3).
//
// remark by T. Thuemmler:
// same random numbers appear each time one restarts the whole program
//
static long IJKLRANDOM=100;
static int iff=0;
static long ijkl,ij,kl,i,j,k,l,ii,jj,m,i97,j97,ivec;
static float s,t,uu[98],c,cd,cm,uni;
if(iff==0)
{ ijkl=IJKLRANDOM;
if(ijkl<1 || ijkl>=900000000) ijkl=1;
ij=ijkl/30082;
kl=ijkl-30082*ij;
i=((ij/177)%177)+2;
j=(ij%177)+2;
k=((kl/169)%178)+1;
l=kl%169;
for(ii=1;ii<=97;ii++)
{ s=0; t=0.5;
for(jj=1;jj<=24;jj++)
{ m=(((i*j)%179)*k)%179;
i=j; j=k; k=m;
l=(53*l+1)%169;
if((l*m)%64 >= 32) s=s+t;
t=0.5*t;
}
uu[ii]=s;
}
c=362436./16777216.;
cd=7654321./16777216.;
cm=16777213./16777216.;
i97=97;
j97=33;
iff=1;
}
for(ivec=1;ivec<=len;ivec++)
{ uni=uu[i97]-uu[j97];
if(uni<0.) uni=uni+1.;
uu[i97]=uni;
i97=i97-1;
if(i97==0) i97=97;
j97=j97-1;
if(j97==0) j97=97;
c=c-cd;
if(c<0.) c=c+cm;
uni=uni-c;
if(uni<0.) uni=uni+1.;
if(uni==0.)
{ uni=uu[j97]*0.59604644775391e-07;
if(uni==0.) uni=0.35527136788005e-14;
}
u[ivec]=uni;
}
// cout << endl<< "random: " << u[1] << endl << flush;
return;
}
double random_james()
{
// This function computes 1 random number in the (0,1) interval,
// using the subrn subroutine.
double u[2];
subrn(u,1);
return u[1];
}
///////////////////////////////////////////////////////////