[tor-commits] [stegotorus/master] Add geometric distribution generator to rng.

[zwol at torproject.org]

commit b6de45bcd5c535a16c8bf18ae8f5e7e1ea27e5aa
Author: Zack Weinberg <zackw at 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