random.c

Go to the documentation of this file.

/// \file
///  Simple random number generators for uniform and Gaussian
///  distributions.  The generator in lcg61 and the hash in mkSeed aren't
///  really industrial strength, but they're more than good enough for
///  present purposes.

#include "random.h"

#include <math.h>

/// \details
///  Use the Box-Muller method to sample a Gaussian distribution with
///  zero mean and unit variance.  To ensure the same input seed always
///  generates the same returned value we do not use the standard
///  technique of saving one of the two generated randoms for the next
///  call.
///
///  \param [in,out] seed  Seed for generator.
///
///  \return A pseudo-random number in the interval (-infinity, infinity).
real_t gasdev(uint64_t* seed)
{
   real_t rsq,v1,v2;
   do 
   {
      v1 = 2.0*lcg61(seed)-1.0;
      v2 = 2.0*lcg61(seed)-1.0;
      rsq = v1*v1+v2*v2;
   } while (rsq >= 1.0 || rsq == 0.0);

   return v2 * sqrt(-2.0*log(rsq)/rsq);
} 

/// \details
///  A 61-bit prime modulus linear congruential generator with
/// modulus = 2^61 -1.
///
///  \param [in,out] seed  Seed for generator.
///
///  \return A pseudo-random number in the interval [0, 1].
double lcg61(uint64_t* seed)
{
   static const double convertToDouble = 1.0/UINT64_C(2305843009213693951);

   *seed *= UINT64_C(437799614237992725);
   *seed %= UINT64_C(2305843009213693951);

   return *seed*convertToDouble;
}

/// \details
///  Forms a 64-bit seed for lcg61 from the combination of 2 32-bit Knuth
///  multiplicative hashes, then runs off 10 values to pass up the worst
///  of the early low-bit correlations.
///
///  \param [in] id An id number such as an atom gid that is unique to
///  each entity that requires random numbers.
///
///  \param [in] callSite A unique number for each call site in the code
///  that needs to generate random seeds.  Using a different value for
///  callSite allows different parts of the code to obtain different 
///  random streams for the same id.
/// 
///  \return A 64-bit seed that is unique to the id and call site.
uint64_t mkSeed(uint32_t id, uint32_t callSite)
{
   uint32_t s1 = id * UINT32_C(2654435761);
   uint32_t s2 = (id+callSite) * UINT32_C(2654435761);

   uint64_t iSeed = (UINT64_C(0x100000000) * s1) + s2;
   for (unsigned jj=0; jj<10; ++jj)
      lcg61(&iSeed);
      
   return iSeed;
}