commit b6de45bcd5c535a16c8bf18ae8f5e7e1ea27e5aa Author: Zack Weinberg <zackw@cmu.edu> Date: Tue Jan 10 01:09:59 2012 +0000 Add geometric distribution generator to rng. git-svn-id: svn+ssh://spartan.csl.sri.com/svn/private/DEFIANCE@212 a58ff0ac-194c-e011-a152-003048836090 --- src/rng.cc | 59 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/rng.h | 9 +++++++++ 2 files changed, 68 insertions(+), 0 deletions(-) diff --git a/src/rng.cc b/src/rng.cc index 9bde6cc..95bea3a 100644 --- a/src/rng.cc +++ b/src/rng.cc @@ -5,6 +5,8 @@ #include "util.h" #include "rng.h" +#include <limits> +#include <math.h> #include <cryptopp/osrng.h> /* Note: this file wraps a C++ library into a C-style program and must @@ -86,3 +88,60 @@ rng_range(unsigned int min, unsigned int max) } CATCH_ALL_EXCEPTIONS(-1); } + +/** Internal use only (can be externalized if someone has a good use + * for it): generate a random double-precision floating-point number + * in the range [0.0, 1.0). Implementation tactic from "Common Lisp + * the Language, 2nd Edition", section 12.9. Assumes IEEE754. + */ +static double +rng_double() +{ + union ieee754_double { + double d; + uint64_t i; + }; + + union ieee754_double n; + + /* This may waste up to 12 bits of randomness on each call, + depending on how clever GenerateWord32 is internally; but the + implementation is much simpler than if we used GenerateBlock. */ + try { + rng_init(); + n.i = (0x3FF0000000000000ULL | + (uint64_t(rng->GenerateWord32(0, 0x000FFFFFu)) << 32) | + uint64_t(rng->GenerateWord32())); + } CATCH_ALL_EXCEPTIONS(std::numeric_limits<double>::quiet_NaN()); + + return n.d - 1.0; +} + +/** Return a random integer in the range [0, hi), geometrically + * distributed over that range, with expected value 'xv'. + * (The rate parameter 'lambda' that's usually used to characterize + * the geometric/exponential distribution is equal to 1/xv.) + * 'hi' must be no more than INT_MAX+1, as for 'rng_range'. + * 'xv' must be greater than 0 and less than 'hi'. + */ +int +rng_range_geom(unsigned int hi, unsigned int xv) +{ + log_assert(hi <= ((unsigned int)INT_MAX)+1); + log_assert(0 < xv && xv < hi); + + double U = rng_double(); + if (isnan(U)) + return -1; + + /* Inverse transform sampling: + T = (-ln U)/lambda; lambda=1/(xv-lo); therefore T = (xv-lo) * -ln(U). + Minor wrinkle: rng_double() produces [0, 1) but we want (0, 1] to + avoid hitting the undefined log(0). This is what nextafter() is for. */ + + double T = -log(nextafter(U, 2.0)) * xv; + + /* Technically we should rejection-sample here instead of clamping, but + that would make this not a constant-time operation. */ + return std::min(hi-1, std::max(0U, (unsigned int)floor(T))); +} diff --git a/src/rng.h b/src/rng.h index 43955bc..02b0947 100644 --- a/src/rng.h +++ b/src/rng.h @@ -19,4 +19,13 @@ int rng_int(unsigned int max); */ int rng_range(unsigned int min, unsigned int max); +/** Return a random integer in the range [0, hi), geometrically + * distributed over that range, with expected value 'xv'. + * (The rate parameter 'lambda' that's usually used to characterize + * the geometric/exponential distribution is equal to 1/xv.) + * 'hi' must be no more than INT_MAX+1, as for 'rng_range'. + * 'xv' must be greater than 0 and less than 'hi'. + */ +int rng_range_geom(unsigned int hi, unsigned int xv); + #endif