gaussdev.c

Go to the documentation of this file.

/*****************************************************************************/
long int GAUSSDEV_SEED;
/*****************************************************************************/
#include "libtpcmodel.h"
/*****************************************************************************/
#ifndef RAND_MAX
#define RAND_MAX 32767
#endif
#define RS_SCALE (1.0 / (1.0 + RAND_MAX))
/*****************************************************************************/
 
/*****************************************************************************/
double gaussdev()
{
  static int ready=0, first=1;
  static double dev;
  double fac, rsq, a, b;
  if(first) {first=0; init_gaussdev();}
 
  /* If we don't have deviate already, then we'll have to make one */
  if(!ready) {
    do {
      //a = 2.*(double)rand()/(double)RAND_MAX - 1.0; 
      //b = 2.*(double)rand()/(double)RAND_MAX - 1.0; 
      a = 2.*drand() - 1.0; 
      b = 2.*drand() - 1.0; 
      rsq = a*a + b*b; 
    } while (rsq>=1.0 || rsq==0.0);
   
    fac = sqrt(-2.0*log(rsq)/rsq);
    dev=a*fac; ready=1; 
    return(b*fac); 
  } else { /* dev is ready so return it */
    ready=0;
    return(dev);
  }
}
/*****************************************************************************/
 
/*****************************************************************************/ 
unsigned int drandSeed(
  short int seed
) {
  unsigned int li;
#ifdef HAVE_TIMESPEC_GET
  struct timespec ts;
  timespec_get(&ts, TIME_UTC);
  li=((ts.tv_sec % 10000)*523 ^ ts.tv_nsec*10) ^ ((getpid() % 1000)*983);
#elif defined HAVE_CLOCK_GETTIME
  struct timespec ts;
  clock_gettime(CLOCK_REALTIME, &ts);
  li=((ts.tv_sec % 10000)*523 ^ ts.tv_nsec*10) ^ ((getpid() % 1000)*983);
#elif defined HAVE_GETTIMEOFDAY
  struct timeval tv;
  gettimeofday(&tv, 0);
  li=((tv.tv_sec % 10000)*523 ^ tv.tv_usec*13) ^ ((getpid() % 1000)*983);
#else
  li=(unsigned int)time(NULL)+(unsigned int)getpid();
#endif
  li+=(unsigned int)rand();
  if(seed) srand(li);
//printf("seed := %u\n", li); printf("RAND_MAX := %u\n", RAND_MAX);
  return(li);
}
/*****************************************************************************/
 
/*****************************************************************************/ 
void init_gaussdev()
{
  if(GAUSSDEV_SEED<1L) GAUSSDEV_SEED=893165470L;
  srand(GAUSSDEV_SEED);
}
/*****************************************************************************/
 
/*****************************************************************************/
double gaussdev2()
{
  static int ready=0;
  static double dev;
  double fac, rsq, a, b;
 
  /* If we don't have deviate already, then we'll have to make one */
  if(!ready) {
    do {
      a = 2.*drand() - 1.0; 
      b = 2.*drand() - 1.0; 
      rsq = a*a + b*b; 
    } while (rsq>=1.0 || rsq==0.0);
   
    fac = sqrt(-2.0*log(rsq)/rsq);
    dev=a*fac; ready=1; 
    return(b*fac); 
  } else { /* dev is ready so return it */
    ready=0;
    return(dev);
  }
}
/*****************************************************************************/
 
/*****************************************************************************/
double drand()
{
  double d, s;
  s=1.0/(1.0+RAND_MAX);
  do {
    d = ( ( s*rand() + rand() )*s + rand() ) * s;
  } while(d>=1.0);
  return d;
}
/*****************************************************************************/
 
/*****************************************************************************/
int rand_range(
  int nr,
  double *d,
  double low,
  double up,
  int type
) {
  int i;
  double dif, v, stl, stu;
 
  if(nr<1) return 0;
  if(d==NULL || type<0 || type>1) return 1;
 
  dif=up-low; if(dif<0.0) return 2;
  if(dif==0.0) {
    for(i=0; i<nr; i++) d[i]=low;
    return 0;
  }
 
  if(type==0) {
    for(i=0; i<nr; i++) d[i] = drand()*dif + low;
  } else if(type==1) {
    stl=copysign(sqrt(fabs(low)),low); if(!isnormal(stl)) stl=0.0;
    stu=copysign(sqrt(fabs(up)), up); if(!isnormal(stu)) stu=0.0;
    dif=stu-stl;
    for(i=0; i<nr; i++) {v=drand()*dif+stl; d[i]=copysign(v*v, v);}
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/