generate-random-numbers.c

double
uniform (void)
// Return a uniformly distributed random variable between 0 and 1.
{
#define ENTROPY_PER_RAND 31
#if RAND_MAX != (1<<ENTROPY_PER_RAND)-1
#error Please fix ENTROPY_PER_RAND
#endif
  int entropy = 0;
  double val = 0.5;
  while ( 1 )
    {
      int exp;
      frexp (val, &exp);
#if 0
      fprintf (stderr, "DEBUG: val=%a (exp=%d, have %d bits of entropy, seek %d)\n",
               val, exp, entropy, 53-exp);
#endif
      if ( entropy >= 53 - exp )
        break;
      val = (val + rand())/(RAND_MAX+1.);
      entropy += ENTROPY_PER_RAND;
    }
  return val;
}

double
gaussian (void)
// Return a gaussian variable with standard deviation 1, using Box-Muller algorithm.
{
  static char gaussian_store_d;
  static double gaussian_store;
  if ( gaussian_store_d )
    {
      gaussian_store_d = 0;
      return gaussian_store;
    }
  double u0 = uniform ();
  double u1 = uniform ();
  gaussian_store_d = 1;
  gaussian_store = sqrt(-2*log(u0)) * sin(2*M_PI*u1);
  return sqrt(-2*log(u0)) * cos(2*M_PI*u1);
}