gaussdev.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <sys/time.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
/*****************************************************************************/
#include "tpcrand.h"
/*****************************************************************************/
 
/*****************************************************************************/ 
unsigned int drandSeed(
  short int seed
) {
  unsigned int li;
#if defined HAVE_TIMESPEC_GET
  struct timespec ts;
  timespec_get(&ts, TIME_UTC);
  //li=(unsigned int)ts.tv_sec+(unsigned int)ts.tv_nsec+(unsigned int)getpid();
  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=(unsigned int)ts.tv_sec+(unsigned int)ts.tv_nsec+(unsigned int)getpid();
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
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 drandRange(
  unsigned int nr,
  double *d,
  double low,
  double up,
  int type
) {
  unsigned 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;
}
/*****************************************************************************/
 
/*****************************************************************************/
double drandGaussian()
{
  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 drandExponential(
  double mean
) {
  double r;
  do {
    r=drand(); 
  } while(r==0.0);
  return(-mean*log(r)); 
}
/*****************************************************************************/
 
/*****************************************************************************/