nps_random.c

Go to the documentation of this file.

#include "nps_random.h"


#include <math.h>
#include <limits.h>


#if 0
/*
 * R250
 * Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
 * Shift-Register Sequence Random Number Generator",
 * Journal of Computational Physics, V.40
 *
 */

static void         r250_init(int seed);
//static unsigned int r250( void );
static double       dr250(void);

/*
 * randclg
 * Linear Congruential Method, the "minimal standard generator"
 * Park & Miller, 1988, Comm of the ACM, 31(10), pp. 1192-1201
 *
 */

static long              set_seed(long);
//static long              get_seed(void);
static unsigned long int randlcg(void);
#endif


void double_vect3_add_gaussian_noise(struct DoubleVect3 *vect, struct DoubleVect3 *std_dev)
{
  vect->x += get_gaussian_noise() * std_dev->x;
  vect->y += get_gaussian_noise() * std_dev->y;
  vect->z += get_gaussian_noise() * std_dev->z;
}

void float_vect3_add_gaussian_noise(struct FloatVect3 *vect, struct FloatVect3 *std_dev)
{
  vect->x += get_gaussian_noise() * std_dev->x;
  vect->y += get_gaussian_noise() * std_dev->y;
  vect->z += get_gaussian_noise() * std_dev->z;
}

void float_rates_add_gaussian_noise(struct FloatRates *vect, struct FloatRates *std_dev)
{
  vect->p += get_gaussian_noise() * std_dev->p;
  vect->q += get_gaussian_noise() * std_dev->q;
  vect->r += get_gaussian_noise() * std_dev->r;
}



void double_vect3_get_gaussian_noise(struct DoubleVect3 *vect, struct DoubleVect3 *std_dev)
{
  vect->x = get_gaussian_noise() * std_dev->x;
  vect->y = get_gaussian_noise() * std_dev->y;
  vect->z = get_gaussian_noise() * std_dev->z;
}


void double_vect3_update_random_walk(struct DoubleVect3 *rw, struct DoubleVect3 *std_dev, double dt, double thau)
{
  struct DoubleVect3 drw;
  double_vect3_get_gaussian_noise(&drw, std_dev);
  struct DoubleVect3 tmp;
  VECT3_SMUL(tmp, *rw, (-1. / thau));
  VECT3_ADD(drw, tmp);
  VECT3_SMUL(drw, drw, dt);
  VECT3_ADD(*rw, drw);
}





#if 0
/*
   http://www.taygeta.com/random/gaussian.html
*/
double get_gaussian_noise(void)
{

  double x1;
  static int nb_call = 0;
  static double x2, w;
  if (nb_call == 0) { r250_init(0); }
  nb_call++;
  if (nb_call % 2) {
    do {
      x1 = 2.0 * dr250() - 1.0;
      x2 = 2.0 * dr250() - 1.0;
      w = x1 * x1 + x2 * x2;
    } while (w >= 1.0);

    w = sqrt((-2.0 * log(w)) / w);
    return x1 * w;
  } else {
    return x2 * w;
  }
}
#else
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdlib.h>
double get_gaussian_noise(void)
{
  static gsl_rng *r = NULL;
  // select random number generator
  if (!r) { r = gsl_rng_alloc(gsl_rng_mt19937); }
  return gsl_ran_gaussian(r, 1.);
}
#endif


#if 0
/*
 * R250
 * Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
 * Shift-Register Sequence Random Number Generator",
 * Journal of Computational Physics, V.40
 *
 */

/* defines to allow for 16 or 32 bit integers */
#define BITS 31


#if WORD_BIT == 32
#ifndef BITS
#define BITS  32
#endif
#else
#ifndef BITS
#define BITS    16
#endif
#endif

#if BITS == 31
#define MSB          0x40000000L
#define ALL_BITS     0x7fffffffL
#define HALF_RANGE   0x20000000L
#define STEP         7
#endif

#if BITS == 32
#define MSB          0x80000000L
#define ALL_BITS     0xffffffffL
#define HALF_RANGE   0x40000000L
#define STEP         7
#endif

#if BITS == 16
#define MSB         0x8000
#define ALL_BITS    0xffff
#define HALF_RANGE  0x4000
#define STEP        11
#endif

static unsigned int r250_buffer[ 250 ];
static int r250_index;

static void r250_init(int sd)
{
  int j, k;
  unsigned int mask, msb;
  set_seed(sd);

  r250_index = 0;
  for (j = 0; j < 250; j++) {    /* fill r250 buffer with BITS-1 bit values */
    r250_buffer[j] = randlcg();
  }

  for (j = 0; j < 250; j++) /* set some MSBs to 1 */
    if (randlcg() > HALF_RANGE) {
      r250_buffer[j] |= MSB;
    }

  msb = MSB;          /* turn on diagonal bit */
  mask = ALL_BITS;  /* turn off the leftmost bits */

  for (j = 0; j < BITS; j++)  {
    k = STEP * j + 3; /* select a word to operate on */
    r250_buffer[k] &= mask; /* turn off bits left of the diagonal */
    r250_buffer[k] |= msb;  /* turn on the diagonal bit */
    mask >>= 1;
    msb  >>= 1;
  }

}

#if 0
/* returns a random unsigned integer */
static unsigned int r250(void)
{
  register int  j;
  register unsigned int new_rand;

  if (r250_index >= 147) {
    j = r250_index - 147;  /* wrap pointer around */
  } else {
    j = r250_index + 103;
  }

  new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
  r250_buffer[ r250_index ] = new_rand;

  if (r250_index >= 249) { /* increment pointer for next time */
    r250_index = 0;
  } else {
    r250_index++;
  }

  return new_rand;

}
#endif

/* returns a random double in range 0..1 */
static double dr250(void)
{
  register int  j;
  register unsigned int new_rand;

  if (r250_index >= 147) {
    j = r250_index - 147;  /* wrap pointer around */
  } else {
    j = r250_index + 103;
  }

  new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
  r250_buffer[ r250_index ] = new_rand;

  if (r250_index >= 249) { /* increment pointer for next time */
    r250_index = 0;
  } else {
    r250_index++;
  }

  return (double)new_rand / ALL_BITS;

}

/*
 * randclg
 * Linear Congruential Method, the "minimal standard generator"
 * Park & Miller, 1988, Comm of the ACM, 31(10), pp. 1192-1201
 *
 */

static long int quotient  = LONG_MAX / 16807L;
static long int my_remainder = LONG_MAX % 16807L;

static long int seed_val = 1L;

static long set_seed(long int sd)
{
  return seed_val = sd;
}

//static long get_seed(void) {
//  return seed_val;
//}

/* returns a random unsigned integer */
unsigned long int randlcg()
{
  if (seed_val <= quotient) {
    seed_val = (seed_val * 16807L) % LONG_MAX;
  } else {
    long int high_part = seed_val / quotient;
    long int low_part  = seed_val % quotient;

    long int test = 16807L * low_part - my_remainder * high_part;

    if (test > 0) {
      seed_val = test;
    } else {
      seed_val = test + LONG_MAX;
    }

  }

  return seed_val;
}

#endif /* 0*/