[Author Prev][Author Next][Thread Prev][Thread Next][Author Index][Thread Index]

[tor-commits] [tor/master] Implement and test probability distributions used by WTF-PAD.



commit 2ccf3268375cd46e8c948e94ba58e0d2f03fe722
Author: George Kadianakis <desnacked@xxxxxxxxxx>
Date:   Tue Nov 27 01:56:23 2018 +0200

    Implement and test probability distributions used by WTF-PAD.
    
    This project introduces the prob_distr.c subsystem which implements all the
    probability distributions that WTF-PAD needs. It also adds unittests for all of
    them.
    
    Code and tests courtesy of Riastradh.
    
    Co-authored-by: Taylor R Campbell <campbell+tor@xxxxxxxxxx>
    Co-authored-by: Mike Perry <mikeperry-git@xxxxxxxxxxxxxx>
---
 src/core/or/circuitpadding.c    |    8 +-
 src/lib/crypt_ops/crypto_rand.c |   11 +
 src/lib/crypt_ops/crypto_rand.h |    1 +
 src/lib/math/.may_include       |    2 +
 src/lib/math/fp.c               |   25 +
 src/lib/math/fp.h               |    1 +
 src/lib/math/include.am         |    6 +-
 src/lib/math/prob_distr.c       | 1628 +++++++++++++++++++++++++++++++++++++++
 src/lib/math/prob_distr.h       |  156 ++++
 src/test/include.am             |    2 +
 src/test/prob_distr_mpfr_ref.c  |   64 ++
 src/test/test.c                 |    1 +
 src/test/test.h                 |    2 +
 src/test/test_prob_distr.c      | 1414 ++++++++++++++++++++++++++++++++++
 src/test/test_slow.c            |    1 +
 15 files changed, 3316 insertions(+), 6 deletions(-)

diff --git a/src/core/or/circuitpadding.c b/src/core/or/circuitpadding.c
index 6cfbf4ba5..a9d927619 100644
--- a/src/core/or/circuitpadding.c
+++ b/src/core/or/circuitpadding.c
@@ -516,10 +516,10 @@ circpad_distribution_sample(circpad_distribution_t dist)
        * param1 is Alpha, param2 is Beta */
       return dist.param1 * pow(p/(1.0-p), 1.0/dist.param2);
     case CIRCPAD_DIST_GEOMETRIC:
-      p = crypto_rand_double();
-      /* https://github.com/distributions-io/geometric-quantile/
-       * param1 is 'p' (success probability) */
-      return ceil(tor_mathlog(1.0-p)/tor_mathlog(1.0-dist.param1));
+      {
+        /* param1 is 'p' (success probability) */
+        return geometric_sample(dist.param1);
+      }
     case CIRCPAD_DIST_WEIBULL:
       p = crypto_rand_double();
       /* https://en.wikipedia.org/wiki/Weibull_distribution \
diff --git a/src/lib/crypt_ops/crypto_rand.c b/src/lib/crypt_ops/crypto_rand.c
index cffd0610f..7a2c417e5 100644
--- a/src/lib/crypt_ops/crypto_rand.c
+++ b/src/lib/crypt_ops/crypto_rand.c
@@ -529,6 +529,17 @@ crypto_rand_unmocked(char *to, size_t n)
 }
 
 /**
+ * Draw an unsigned 32-bit integer uniformly at random.
+ */
+uint32_t
+crypto_rand_uint32(void)
+{
+  uint32_t rand;
+  crypto_rand((void*)&rand, sizeof(rand));
+  return rand;
+}
+
+/**
  * Return a pseudorandom integer, chosen uniformly from the values
  * between 0 and <b>max</b>-1 inclusive.  <b>max</b> must be between 1 and
  * INT_MAX+1, inclusive.
diff --git a/src/lib/crypt_ops/crypto_rand.h b/src/lib/crypt_ops/crypto_rand.h
index 0c538d81a..61fd82c80 100644
--- a/src/lib/crypt_ops/crypto_rand.h
+++ b/src/lib/crypt_ops/crypto_rand.h
@@ -27,6 +27,7 @@ int crypto_rand_int(unsigned int max);
 int crypto_rand_int_range(unsigned int min, unsigned int max);
 uint64_t crypto_rand_uint64_range(uint64_t min, uint64_t max);
 time_t crypto_rand_time_range(time_t min, time_t max);
+uint32_t crypto_rand_uint32(void);
 uint64_t crypto_rand_uint64(uint64_t max);
 double crypto_rand_double(void);
 struct tor_weak_rng_t;
diff --git a/src/lib/math/.may_include b/src/lib/math/.may_include
index 1fd26864d..f8bc264a5 100644
--- a/src/lib/math/.may_include
+++ b/src/lib/math/.may_include
@@ -3,3 +3,5 @@ orconfig.h
 lib/cc/*.h
 lib/log/*.h
 lib/math/*.h
+lib/testsupport/*.h
+lib/crypt_ops/*.h
diff --git a/src/lib/math/fp.c b/src/lib/math/fp.c
index d5989db63..57082fa46 100644
--- a/src/lib/math/fp.c
+++ b/src/lib/math/fp.c
@@ -117,3 +117,28 @@ ENABLE_GCC_WARNING(double-promotion)
 ENABLE_GCC_WARNING(float-conversion)
 #endif
 }
+
+/* isinf() wrapper for tor */
+int
+tor_isinf(double x)
+{
+  /* Same as above, work around the "double promotion" warnings */
+#if defined(MINGW_ANY) && GCC_VERSION >= 409
+#define PROBLEMATIC_FLOAT_CONVERSION_WARNING
+DISABLE_GCC_WARNING(float-conversion)
+#endif /* defined(MINGW_ANY) && GCC_VERSION >= 409 */
+#if defined(__clang__)
+#if __has_warning("-Wdouble-promotion")
+#define PROBLEMATIC_DOUBLE_PROMOTION_WARNING
+DISABLE_GCC_WARNING(double-promotion)
+#endif
+#endif /* defined(__clang__) */
+  return isinf(x);
+#ifdef PROBLEMATIC_DOUBLE_PROMOTION_WARNING
+ENABLE_GCC_WARNING(double-promotion)
+#endif
+#ifdef PROBLEMATIC_FLOAT_CONVERSION_WARNING
+ENABLE_GCC_WARNING(float-conversion)
+#endif
+}
+
diff --git a/src/lib/math/fp.h b/src/lib/math/fp.h
index e27b8f8d8..ddf3ed24d 100644
--- a/src/lib/math/fp.h
+++ b/src/lib/math/fp.h
@@ -19,5 +19,6 @@ double tor_mathlog(double d) ATTR_CONST;
 long tor_lround(double d) ATTR_CONST;
 int64_t tor_llround(double d) ATTR_CONST;
 int64_t clamp_double_to_int64(double number);
+int tor_isinf(double x);
 
 #endif
diff --git a/src/lib/math/include.am b/src/lib/math/include.am
index b088b3f3c..6d65ce90a 100644
--- a/src/lib/math/include.am
+++ b/src/lib/math/include.am
@@ -7,7 +7,8 @@ endif
 
 src_lib_libtor_math_a_SOURCES =	\
 		src/lib/math/fp.c		\
-		src/lib/math/laplace.c
+		src/lib/math/laplace.c 	\
+		src/lib/math/prob_distr.c
 
 
 src_lib_libtor_math_testing_a_SOURCES = \
@@ -17,4 +18,5 @@ src_lib_libtor_math_testing_a_CFLAGS = $(AM_CFLAGS) $(TEST_CFLAGS)
 
 noinst_HEADERS +=				\
 		src/lib/math/fp.h		\
-		src/lib/math/laplace.h
+		src/lib/math/laplace.h  \
+		src/lib/math/prob_distr.h
diff --git a/src/lib/math/prob_distr.c b/src/lib/math/prob_distr.c
new file mode 100644
index 000000000..832d3b4d9
--- /dev/null
+++ b/src/lib/math/prob_distr.c
@@ -0,0 +1,1628 @@
+/* Copyright (c) 2018, The Tor Project, Inc. */
+/* See LICENSE for licensing information */
+
+/**
+ * \file prob_distr.c
+ *
+ * \brief
+ *  Implements various probability distributions.
+ *  Almost all code is courtesy of Riastradh.
+ *
+ * \details
+ * Here are some details that might help you understand this file:
+ *
+ * - Throughout this file, `eps' means the largest relative error of a
+ *   correctly rounded floating-point operation, which in binary64
+ *   floating-point arithmetic is 2^-53.  Here the relative error of a
+ *   true value x from a computed value y is |x - y|/|x|.  This
+ *   definition of epsilon is conventional for numerical analysts when
+ *   writing error analyses.  (If your libm doesn't provide correctly
+ *   rounded exp and log, their relative error is usually below 2*2^-53
+ *   and probably closer to 1.1*2^-53 instead.)
+ *
+ *   The C constant DBL_EPSILON is actually twice this, and should
+ *   perhaps rather be named ulp(1) -- that is, it is the distance from
+ *   1 to the next greater floating-point number, which is usually of
+ *   more interest to programmers and hardware engineers.
+ *
+ *   Since this file is concerned mainly with error bounds rather than
+ *   with low-level bit-hacking of floating-point numbers, we adopt the
+ *   numerical analysts' definition in the comments, though we do use
+ *   DBL_EPSILON in a handful of places where it is convenient to use
+ *   some function of eps = DBL_EPSILON/2 in a case analysis.
+ *
+ * - In various functions (e.g. sample_log_logistic()) we jump through hoops so
+ *   that we can use reals closer to 0 than closer to 1, since we achieve much
+ *   greater accuracy for floating point numbers near 0. In particular, we can
+ *   represent differences as small as 10^-300 for numbers near 0, but of no
+ *   less than 10^-16 for numbers near 1.
+ **/
+
+#define PROB_DISTR_PRIVATE
+
+#include "orconfig.h"
+
+#include "lib/math/prob_distr.h"
+
+#include "lib/crypt_ops/crypto_rand.h"
+#include "lib/cc/ctassert.h"
+
+#include <float.h>
+#include <math.h>
+#include <stddef.h>
+
+/** Validators for downcasting macros below */
+#define validate_container_of(PTR, TYPE, FIELD)                         \
+  (0 * sizeof((PTR) - &((TYPE *)(((char *)(PTR)) -                      \
+      offsetof(TYPE, FIELD)))->FIELD))
+#define validate_const_container_of(PTR, TYPE, FIELD)                   \
+  (0 * sizeof((PTR) - &((const TYPE *)(((const char *)(PTR)) -          \
+      offsetof(TYPE, FIELD)))->FIELD))
+/** Downcasting macro */
+#define container_of(PTR, TYPE, FIELD)                                  \
+  ((TYPE *)(((char *)(PTR)) - offsetof(TYPE, FIELD))                    \
+    + validate_container_of(PTR, TYPE, FIELD))
+/** Constified downcasting macro */
+#define const_container_of(PTR, TYPE, FIELD)                            \
+  ((const TYPE *)(((const char *)(PTR)) - offsetof(TYPE, FIELD))        \
+    + validate_const_container_of(PTR, TYPE, FIELD))
+
+/**
+ * Count number of one bits in 32-bit word.
+ */
+static unsigned
+bitcount32(uint32_t x)
+{
+
+  /* Count two-bit groups.  */
+  x -= (x >> 1) & UINT32_C(0x55555555);
+
+  /* Count four-bit groups.  */
+  x = ((x >> 2) & UINT32_C(0x33333333)) + (x & UINT32_C(0x33333333));
+
+  /* Count eight-bit groups.  */
+  x = (x + (x >> 4)) & UINT32_C(0x0f0f0f0f);
+
+  /* Sum all eight-bit groups, and extract the sum.  */
+  return (x * UINT32_C(0x01010101)) >> 24;
+}
+
+/**
+ * Count leading zeros in 32-bit word.
+ */
+static unsigned
+clz32(uint32_t x)
+{
+
+  /* Round up to a power of two.  */
+  x |= x >> 1;
+  x |= x >> 2;
+  x |= x >> 4;
+  x |= x >> 8;
+  x |= x >> 16;
+
+  /* Subtract count of one bits from 32.  */
+  return (32 - bitcount32(x));
+}
+
+/*
+ * Some lemmas that will be used throughout this file to prove various error
+ * bounds:
+ *
+ * Lemma 1.  If |d| <= 1/2, then 1/(1 + d) <= 2.
+ *
+ * Proof.  If 0 <= d <= 1/2, then 1 + d >= 1, so that 1/(1 + d) <= 1.
+ * If -1/2 <= d <= 0, then 1 + d >= 1/2, so that 1/(1 + d) <= 2.  QED.
+ *
+ * Lemma 2. If b = a*(1 + d)/(1 + d') for |d'| < 1/2 and nonzero a, b,
+ * then b = a*(1 + e) for |e| <= 2|d' - d|.
+ *
+ * Proof.  |a - b|/|a|
+ *             = |a - a*(1 + d)/(1 + d')|/|a|
+ *             = |1 - (1 + d)/(1 + d')|
+ *             = |(1 + d' - 1 - d)/(1 + d')|
+ *             = |(d' - d)/(1 + d')|
+ *            <= 2|d' - d|, by Lemma 1,
+ *
+ * QED.
+ *
+ * Lemma 3.  For |d|, |d'| < 1/4,
+ *
+ *     |log((1 + d)/(1 + d'))| <= 4|d - d'|.
+ *
+ * Proof.  Write
+ *
+ *     log((1 + d)/(1 + d'))
+ *      = log(1 + (1 + d)/(1 + d') - 1)
+ *      = log(1 + (1 + d - 1 - d')/(1 + d')
+ *      = log(1 + (d - d')/(1 + d')).
+ *
+ * By Lemma 1, |(d - d')/(1 + d')| < 2|d' - d| < 1, so the Taylor
+ * series of log(1 + x) converges absolutely for (d - d')/(1 + d'),
+ * and thus we have
+ *
+ *     |log(1 + (d - d')/(1 + d'))|
+ *      = |\sum_{n=1}^\infty ((d - d')/(1 + d'))^n/n|
+ *     <= \sum_{n=1}^\infty |(d - d')/(1 + d')|^n/n
+ *     <= \sum_{n=1}^\infty |2(d' - d)|^n/n
+ *     <= \sum_{n=1}^\infty |2(d' - d)|^n
+ *      = 1/(1 - |2(d' - d)|)
+ *     <= 4|d' - d|,
+ *
+ * QED.
+ *
+ * Lemma 4.  If 1/e <= 1 + x <= e, then
+ *
+ *     log(1 + (1 + d) x) = (1 + d') log(1 + x)
+ *
+ * for |d'| < 8|d|.
+ *
+ * Proof.  Write
+ *
+ *     log(1 + (1 + d) x)
+ *     = log(1 + x + x*d)
+ *     = log((1 + x) (1 + x + x*d)/(1 + x))
+ *     = log(1 + x) + log((1 + x + x*d)/(1 + x))
+ *     = log(1 + x) (1 + log((1 + x + x*d)/(1 + x))/log(1 + x)).
+ *
+ * The relative error is bounded by
+ *
+ *     |log((1 + x + x*d)/(1 + x))/log(1 + x)|
+ *     <= 4|x + x*d - x|/|log(1 + x)|, by Lemma 3,
+ *      = 4|x*d|/|log(1 + x)|
+ *      < 8|d|,
+ *
+ * since in this range 0 < 1 - 1/e < x/log(1 + x) <= e - 1 < 2.  QED.
+ */
+
+/**
+ * Compute the logistic function: f(x) = 1/(1 + e^{-x}) = e^x/(1 + e^x).
+ * Maps a log-odds-space probability in [-\infty, +\infty] into a direct-space
+ * probability in [0,1].  Inverse of logit.
+ *
+ * Ill-conditioned for large x; the identity logistic(-x) = 1 -
+ * logistic(x) and the function logistichalf(x) = logistic(x) - 1/2 may
+ * help to rearrange a computation.
+ *
+ * This implementation gives relative error bounded by 7 eps.
+ */
+STATIC double
+logistic(double x)
+{
+  if (x <= log(DBL_EPSILON/2)) {
+    /*
+     * If x <= log(DBL_EPSILON/2) = log(eps), then e^x <= eps. In this case
+     * we will approximate the logistic() function with e^x because the
+     * relative error is less than eps. Here is a calculation of the
+     * relative error between the logistic() function and e^x and a proof
+     * that it's less than eps:
+     *
+     *     |e^x - e^x/(1 + e^x)|/|e^x/(1 + e^x)|
+     *     <= |1 - 1/(1 + e^x)|*|1 + e^x|
+     *      = |e^x/(1 + e^x)|*|1 + e^x|
+     *      = |e^x|
+     *     <= eps.
+     */
+    return exp(x); /* return e^x */
+  } else if (x <= -log(DBL_EPSILON/2)) {
+    /*
+     * e^{-x} > 0, so 1 + e^{-x} > 1, and 0 < 1/(1 +
+     * e^{-x}) < 1; further, since e^{-x} < 1 + e^{-x}, we
+     * also have 0 < 1/(1 + e^{-x}) < 1.  Thus, if exp has
+     * relative error d0, + has relative error d1, and /
+     * has relative error d2, then we get
+     *
+     *     (1 + d2)/[(1 + (1 + d0) e^{-x})(1 + d1)]
+     *     = (1 + d0)/[1 + e^{-x} + d0 e^{-x}
+     *                     + d1 + d1 e^{-x} + d0 d1 e^{-x}]
+     *     = (1 + d0)/[(1 + e^{-x})
+     *                 * (1 + d0 e^{-x}/(1 + e^{-x})
+     *                      + d1/(1 + e^{-x})
+     *                      + d0 d1 e^{-x}/(1 + e^{-x}))].
+     *     = (1 + d0)/[(1 + e^{-x})(1 + d')]
+     *     = [1/(1 + e^{-x})] (1 + d0)/(1 + d')
+     *
+     * where
+     *
+     *     d' = d0 e^{-x}/(1 + e^{-x})
+     *          + d1/(1 + e^{-x})
+     *          + d0 d1 e^{-x}/(1 + e^{-x}).
+     *
+     * By Lemma 2 this relative error is bounded by
+     *
+     *     2|d0 - d'|
+     *      = 2|d0 - d0 e^{-x}/(1 + e^{-x})
+     *             - d1/(1 + e^{-x})
+     *             - d0 d1 e^{-x}/(1 + e^{-x})|
+     *     <= 2|d0| + 2|d0 e^{-x}/(1 + e^{-x})|
+     *             + 2|d1/(1 + e^{-x})|
+     *             + 2|d0 d1 e^{-x}/(1 + e^{-x})|
+     *     <= 2|d0| + 2|d0| + 2|d1| + 2|d0 d1|
+     *     <= 4|d0| + 2|d1| + 2|d0 d1|
+     *     <= 6 eps + 2 eps^2.
+     */
+    return 1/(1 + exp(-x));
+  } else {
+    /*
+     * e^{-x} <= eps, so the relative error of 1 from 1/(1
+     * + e^{-x}) is
+     *
+     *     |1/(1 + e^{-x}) - 1|/|1/(1 + e^{-x})|
+     *      = |e^{-x}/(1 + e^{-x})|/|1/(1 + e^{-x})|
+     *      = |e^{-x}|
+     *     <= eps.
+     *
+     * This computation avoids an intermediate overflow
+     * exception, although the effect on the result is
+     * harmless.
+     *
+     * XXX Should maybe raise inexact here.
+     */
+    return 1;
+  }
+}
+
+/**
+ * Compute the logit function: log p/(1 - p).  Defined on [0,1].  Maps
+ * a direct-space probability in [0,1] to a log-odds-space probability
+ * in [-\infty, +\infty].  Inverse of logistic.
+ *
+ * Ill-conditioned near 1/2 and 1; the identity logit(1 - p) =
+ * -logit(p) and the function logithalf(p0) = logit(1/2 + p0) may help
+ * to rearrange a computation for p in [1/(1 + e), 1 - 1/(1 + e)].
+ *
+ * This implementation gives relative error bounded by 10 eps.
+ */
+STATIC double
+logit(double p)
+{
+
+  /* logistic(-1) <= p <= logistic(+1) */
+  if (1/(1 + exp(1)) <= p && p <= 1/(1 + exp(-1))) {
+    /*
+     * For inputs near 1/2, we want to compute log1p(near
+     * 0) rather than log(near 1), so write this as:
+     *
+     * log(p/(1 - p)) = -log((1 - p)/p)
+     * = -log(1 + (1 - p)/p - 1)
+     * = -log(1 + (1 - p - p)/p)
+     * = -log(1 + (1 - 2p)/p).
+     *
+     * Since p = 2p/2 <= 1 <= 2*2p = 4p, the floating-point
+     * evaluation of 1 - 2p is exact; the only error arises
+     * from division and log1p.  First, note that if
+     * logistic(-1) <= p <= logistic(+1), (1 - 2p)/p lies
+     * in the bounds of Lemma 4.
+     *
+     * If division has relative error d0 and log1p has
+     * relative error d1, the outcome is
+     *
+     *     -(1 + d1) log(1 + (1 - 2p) (1 + d0)/p)
+     *     = -(1 + d1) (1 + d') log(1 + (1 - 2p)/p)
+     *     = -(1 + d1 + d' + d1 d') log(1 + (1 - 2p)/p).
+     *
+     * where |d'| < 8|d0| by Lemma 4.  The relative error
+     * is then bounded by
+     *
+     *     |d1 + d' + d1 d'|
+     *     <= |d1| + 8|d0| + 8|d1 d0|
+     *     <= 9 eps + 8 eps^2.
+     */
+    return -log1p((1 - 2*p)/p);
+  } else {
+    /*
+     * For inputs near 0, although 1 - p may be rounded to
+     * 1, it doesn't matter much because the magnitude of
+     * the result is so much larger.  For inputs near 1, we
+     * can compute 1 - p exactly, although the precision on
+     * the input is limited so we won't ever get more than
+     * about 700 for the output.
+     *
+     * If - has relative error d0, / has relative error d1,
+     * and log has relative error d2, then
+     *
+     *     (1 + d2) log((1 + d0) p/[(1 - p)(1 + d1)])
+     *     = (1 + d2) [log(p/(1 - p)) + log((1 + d0)/(1 + d1))]
+     *     = log(p/(1 - p)) + d2 log(p/(1 - p))
+     *       + (1 + d2) log((1 + d0)/(1 + d1))
+     *     = log(p/(1 - p))*[1 + d2 +
+     *         + (1 + d2) log((1 + d0)/(1 + d1))/log(p/(1 - p))]
+     *
+     * Since 0 <= p < logistic(-1) or logistic(+1) < p <=
+     * 1, we have |log(p/(1 - p))| > 1.  Hence this error
+     * is bounded by
+     *
+     *     |d2 + (1 + d2) log((1 + d0)/(1 + d1))/log(p/(1 - p))|
+     *     <= |d2| + |(1 + d2) log((1 + d0)/(1 + d1))
+     *                      / log(p/(1 - p))|
+     *     <= |d2| + |(1 + d2) log((1 + d0)/(1 + d1))|
+     *     <= |d2| + 4|(1 + d2) (d0 - d1)|, by Lemma 3,
+     *     <= |d2| + 4|d0 - d1 + d2 d0 - d1 d0|
+     *     <= |d2| + 4|d0| + 4|d1| + 4|d2 d0| + 4|d1 d0|
+     *     <= 9 eps + 8 eps^2.
+     */
+    return log(p/(1 - p));
+  }
+}
+
+/**
+ * Compute the logit function, translated in input by 1/2: logithalf(p)
+ * = logit(1/2 + p).  Defined on [-1/2, 1/2].  Inverse of logistichalf.
+ *
+ * Ill-conditioned near +/-1/2.  If |p0| > 1/2 - 1/(1 + e), it may be
+ * better to compute 1/2 + p0 or -1/2 - p0 and to use logit instead.
+ * This implementation gives relative error bounded by 34 eps.
+ */
+STATIC double
+logithalf(double p0)
+{
+
+  if (fabs(p0) <= 0.5 - 1/(1 + exp(1))) {
+    /*
+     * logit(1/2 + p0)
+     * = log((1/2 + p0)/(1 - (1/2 + p0)))
+     * = log((1/2 + p0)/(1/2 - p0))
+     * = log(1 + (1/2 + p0)/(1/2 - p0) - 1)
+     * = log(1 + (1/2 + p0 - (1/2 - p0))/(1/2 - p0))
+     * = log(1 + (1/2 + p0 - 1/2 + p0)/(1/2 - p0))
+     * = log(1 + 2 p0/(1/2 - p0))
+     *
+     * If the error of subtraction is d0, the error of
+     * division is d1, and the error of log1p is d2, then
+     * what we compute is
+     *
+     *  (1 + d2) log(1 + (1 + d1) 2 p0/[(1 + d0) (1/2 - p0)])
+     *  = (1 + d2) log(1 + (1 + d') 2 p0/(1/2 - p0))
+     *  = (1 + d2) (1 + d'') log(1 + 2 p0/(1/2 - p0))
+     *  = (1 + d2 + d'' + d2 d'') log(1 + 2 p0/(1/2 - p0)),
+     *
+     * where |d'| < 2|d0 - d1| <= 4 eps by Lemma 2, and
+     * |d''| < 8|d'| < 32 eps by Lemma 4 since
+     *
+     *  1/e <= 1 + 2*p0/(1/2 - p0) <= e
+     *
+     * when |p0| <= 1/2 - 1/(1 + e).  Hence the relative
+     * error is bounded by
+     *
+     *  |d2 + d'' + d2 d''|
+     *  <= |d2| + |d''| + |d2 d''|
+     *  <= |d1| + 32 |d0| + 32 |d1 d0|
+     *  <= 33 eps + 32 eps^2.
+     */
+    return log1p(2*p0/(0.5 - p0));
+  } else {
+    /*
+     * We have a choice of computing logit(1/2 + p0) or
+     * -logit(1 - (1/2 + p0)) = -logit(1/2 - p0).  It
+     * doesn't matter which way we do this: either way,
+     * since 1/2 p0 <= 1/2 <= 2 p0, the sum and difference
+     * are computed exactly.  So let's do the one that
+     * skips the final negation.
+     *
+     * The result is
+     *
+     *  (1 + d1) log((1 + d0) (1/2 + p0)/[(1 + d2) (1/2 - p0)])
+     *  = (1 + d1) (1 + log((1 + d0)/(1 + d2))
+     *                  / log((1/2 + p0)/(1/2 - p0)))
+     *    * log((1/2 + p0)/(1/2 - p0))
+     *  = (1 + d') log((1/2 + p0)/(1/2 - p0))
+     *  = (1 + d') logit(1/2 + p0)
+     *
+     * where
+     *
+     *  d' = d1 + log((1 + d0)/(1 + d2))/logit(1/2 + p0)
+     *       + d1 log((1 + d0)/(1 + d2))/logit(1/2 + p0).
+     *
+     * For |p| > 1/2 - 1/(1 + e), logit(1/2 + p0) > 1.
+     * Provided |d0|, |d2| < 1/4, by Lemma 3 we have
+     *
+     *  |log((1 + d0)/(1 + d2))| <= 4|d0 - d2|.
+     *
+     * Hence the relative error is bounded by
+     *
+     *  |d'| <= |d1| + 4|d0 - d2| + 4|d1| |d0 - d2|
+     *       <= |d1| + 4|d0| + 4|d2| + 4|d1 d0| + 4|d1 d2|
+     *       <= 9 eps + 8 eps^2.
+     */
+    return log((0.5 + p0)/(0.5 - p0));
+  }
+}
+
+/*
+ * The following random_uniform_01 is tailored for IEEE 754 binary64
+ * floating-point or smaller.  It can be adapted to larger
+ * floating-point formats like i387 80-bit or IEEE 754 binary128, but
+ * it may require sampling more bits.
+ */
+CTASSERT(FLT_RADIX == 2);
+CTASSERT(-DBL_MIN_EXP <= 1021);
+CTASSERT(DBL_MANT_DIG <= 53);
+
+/**
+ * Draw a floating-point number in [0, 1] with uniform distribution.
+ *
+ * Note that the probability of returning 0 is less than 2^-1074, so
+ * callers need not check for it.  However, callers that cannot handle
+ * rounding to 1 must deal with that, because it occurs with
+ * probability 2^-54, which is small but nonnegligible.
+ */
+STATIC double
+random_uniform_01(void)
+{
+  uint32_t z, x, hi, lo;
+  double s;
+
+  /*
+   * Draw an exponent, geometrically distributed, but give up if
+   * we get a run of more than 1088 zeros, which really means the
+   * system is broken.
+   */
+  z = 0;
+  while ((x = crypto_rand_uint32()) == 0) {
+    if (z >= 1088)
+      /* Your bit sampler is broken.  Go home.  */
+      return 0;
+    z += 32;
+  }
+  z += clz32(x);
+
+  /*
+   * Pick 32-bit halves of an odd normalized significand.
+   * Picking it odd breaks ties in the subsequent rounding, which
+   * occur only with measure zero in the uniform distribution on
+   * [0, 1].
+   */
+  hi = crypto_rand_uint32() | UINT32_C(0x80000000);
+  lo = crypto_rand_uint32() | UINT32_C(0x00000001);
+
+  /* Round to nearest scaled significand in [2^63, 2^64].  */
+  s = hi*(double)4294967296 + lo;
+
+  /* Rescale into [1/2, 1] and apply exponent in one swell foop.  */
+  return s * ldexp(1, -(64 + z));
+}
+
+/*******************************************************************/
+
+/* Functions for specific probability distributions start here: */
+
+/*
+ * Logistic(mu, sigma) distribution, supported on (-\infty,+\infty)
+ *
+ * This is the uniform distribution on [0,1] mapped into log-odds
+ * space, scaled by sigma and translated by mu.
+ *
+ * pdf(x) = e^{-(x - mu)/sigma} sigma (1 + e^{-(x - mu)/sigma})^2
+ * cdf(x) = 1/(1 + e^{-(x - mu)/sigma}) = logistic((x - mu)/sigma)
+ * sf(x) = 1 - cdf(x) = 1 - logistic((x - mu)/sigma = logistic(-(x - mu)/sigma)
+ * icdf(p) = mu + sigma log p/(1 - p) = mu + sigma logit(p)
+ * isf(p) = mu + sigma log (1 - p)/p = mu - sigma logit(p)
+ */
+
+/**
+ * Compute the CDF of the Logistic(mu, sigma) distribution: the
+ * logistic function.  Well-conditioned for negative inputs and small
+ * positive inputs; ill-conditioned for large positive inputs.
+ */
+STATIC double
+cdf_logistic(double x, double mu, double sigma)
+{
+  return logistic((x - mu)/sigma);
+}
+
+/**
+ * Compute the SF of the Logistic(mu, sigma) distribution: the logistic
+ * function reflected over the y axis.  Well-conditioned for positive
+ * inputs and small negative inputs; ill-conditioned for large negative
+ * inputs.
+ */
+STATIC double
+sf_logistic(double x, double mu, double sigma)
+{
+  return logistic(-(x - mu)/sigma);
+}
+
+/**
+ * Compute the inverse of the CDF of the Logistic(mu, sigma)
+ * distribution: the logit function.  Well-conditioned near 0;
+ * ill-conditioned near 1/2 and 1.
+ */
+STATIC double
+icdf_logistic(double p, double mu, double sigma)
+{
+  return mu + sigma*logit(p);
+}
+
+/**
+ * Compute the inverse of the SF of the Logistic(mu, sigma)
+ * distribution: the -logit function.  Well-conditioned near 0;
+ * ill-conditioned near 1/2 and 1.
+ */
+STATIC double
+isf_logistic(double p, double mu, double sigma)
+{
+  return mu - sigma*logit(p);
+}
+
+/*
+ * LogLogistic(alpha, beta) distribution, supported on (0, +\infty).
+ *
+ * This is the uniform distribution on [0,1] mapped into odds space,
+ * scaled by positive alpha and shaped by positive beta.
+ *
+ * Equivalent to computing exp of a Logistic(log alpha, 1/beta) sample.
+ * (Name arises because the pdf has LogLogistic(x; alpha, beta) =
+ * Logistic(log x; log alpha, 1/beta) and mathematicians got their
+ * covariance contravariant.)
+ *
+ * pdf(x) = (beta/alpha) (x/alpha)^{beta - 1}/(1 + (x/alpha)^beta)^2
+ *        = (1/e^mu sigma) (x/e^mu)^{1/sigma - 1} /
+ *              (1 + (x/e^mu)^{1/sigma})^2
+ * cdf(x) = 1/(1 + (x/alpha)^-beta) = 1/(1 + (x/e^mu)^{-1/sigma})
+ *        = 1/(1 + (e^{log x}/e^mu)^{-1/sigma})
+ *        = 1/(1 + (e^{log x - mu})^{-1/sigma})
+ *        = 1/(1 + e^{-(log x - mu)/sigma})
+ *        = logistic((log x - mu)/sigma)
+ *        = logistic((log x - log alpha)/(1/beta))
+ * sf(x) = 1 - 1/(1 + (x/alpha)^-beta)
+ *       = (x/alpha)^-beta/(1 + (x/alpha)^-beta)
+ *       = 1/((x/alpha)^beta + 1)
+ *       = 1/(1 + (x/alpha)^beta)
+ * icdf(p) = alpha (p/(1 - p))^{1/beta}
+ *         = alpha e^{logit(p)/beta}
+ *         = e^{mu + sigma logit(p)}
+ * isf(p) = alpha ((1 - p)/p)^{1/beta}
+ *        = alpha e^{-logit(p)/beta}
+ *        = e^{mu - sigma logit(p)}
+ */
+
+/**
+ * Compute the CDF of the LogLogistic(alpha, beta) distribution.
+ * Well-conditioned for all x and alpha, and the condition number
+ *
+ *      -beta/[1 + (x/alpha)^{-beta}]
+ *
+ * grows linearly with beta.
+ *
+ * Loosely, the relative error of this implementation is bounded by
+ *
+ *      4 eps + 2 eps^2 + O(beta eps),
+ *
+ * so don't bother trying this for beta anywhere near as large as
+ * 1/eps, around which point it levels off at 1.
+ */
+STATIC double
+cdf_log_logistic(double x, double alpha, double beta)
+{
+  /*
+   * Let d0 be the error of x/alpha; d1, of pow; d2, of +; and
+   * d3, of the final quotient.  The exponentiation gives
+   *
+   *    ((1 + d0) x/alpha)^{-beta}
+   *    = (x/alpha)^{-beta} (1 + d0)^{-beta}
+   *    = (x/alpha)^{-beta} (1 + (1 + d0)^{-beta} - 1)
+   *    = (x/alpha)^{-beta} (1 + d')
+   *
+   * where d' = (1 + d0)^{-beta} - 1.  If y = (x/alpha)^{-beta},
+   * the denominator is
+   *
+   *    (1 + d2) (1 + (1 + d1) (1 + d') y)
+   *    = (1 + d2) (1 + y + (d1 + d' + d1 d') y)
+   *    = 1 + y + (1 + d2) (d1 + d' + d1 d') y
+   *    = (1 + y) (1 + (1 + d2) (d1 + d' + d1 d') y/(1 + y))
+   *    = (1 + y) (1 + d''),
+   *
+   * where d'' = (1 + d2) (d1 + d' + d1 d') y/(1 + y).  The
+   * final result is
+   *
+   *    (1 + d3) / [(1 + d2) (1 + d'') (1 + y)]
+   *    = (1 + d''') / (1 + y)
+   *
+   * for |d'''| <= 2|d3 - d''| by Lemma 2 as long as |d''| < 1/2
+   * (which may not be the case for very large beta).  This
+   * relative error is therefore bounded by
+   *
+   *    |d'''|
+   *    <= 2|d3 - d''|
+   *    <= 2|d3| + 2|(1 + d2) (d1 + d' + d1 d') y/(1 + y)|
+   *    <= 2|d3| + 2|(1 + d2) (d1 + d' + d1 d')|
+   *     = 2|d3| + 2|d1 + d' + d1 d' + d2 d1 + d2 d' + d2 d1 d'|
+   *      <= 2|d3| + 2|d1| + 2|d'| + 2|d1 d'| + 2|d2 d1| + 2|d2 d'|
+   *         + 2|d2 d1 d'|
+   *      <= 4 eps + 2 eps^2 + (2 + 2 eps + 2 eps^2) |d'|.
+   *
+   * Roughly, |d'| = |(1 + d0)^{-beta} - 1| grows like beta eps,
+   * until it levels off at 1.
+   */
+  return 1/(1 + pow(x/alpha, -beta));
+}
+
+/**
+ * Compute the SF of the LogLogistic(alpha, beta) distribution.
+ * Well-conditioned for all x and alpha, and the condition number
+ *
+ *      beta/[1 + (x/alpha)^beta]
+ *
+ * grows linearly with beta.
+ *
+ * Loosely, the relative error of this implementation is bounded by
+ *
+ *      4 eps + 2 eps^2 + O(beta eps)
+ *
+ * so don't bother trying this for beta anywhere near as large as
+ * 1/eps, beyond which point it grows unbounded.
+ */
+STATIC double
+sf_log_logistic(double x, double alpha, double beta)
+{
+  /*
+   * The error analysis here is essentially the same as in
+   * cdf_log_logistic, except that rather than levelling off at
+   * 1, |(1 + d0)^beta - 1| grows unbounded.
+   */
+  return 1/(1 + pow(x/alpha, beta));
+}
+
+/**
+ * Compute the inverse of the CDF of the LogLogistic(alpha, beta)
+ * distribution.  Ill-conditioned for p near 1 and beta near 0 with
+ * condition number 1/[beta (1 - p)].
+ */
+STATIC double
+icdf_log_logistic(double p, double alpha, double beta)
+{
+  return alpha*pow(p/(1 - p), 1/beta);
+}
+
+/**
+ * Compute the inverse of the SF of the LogLogistic(alpha, beta)
+ * distribution.  Ill-conditioned for p near 1 and for large beta, with
+ * condition number -1/[beta (1 - p)].
+ */
+STATIC double
+isf_log_logistic(double p, double alpha, double beta)
+{
+  return alpha*pow((1 - p)/p, 1/beta);
+}
+
+/*
+ * Weibull(lambda, k) distribution, supported on (0, +\infty).
+ *
+ * pdf(x) = (k/lambda) (x/lambda)^{k - 1} e^{-(x/lambda)^k}
+ * cdf(x) = 1 - e^{-(x/lambda)^k}
+ * icdf(p) = lambda * (-log (1 - p))^{1/k}
+ * sf(x) = e^{-(x/lambda)^k}
+ * isf(p) = lambda * (-log p)^{1/k}
+ */
+
+/**
+ * Compute the CDF of the Weibull(lambda, k) distribution.
+ * Well-conditioned for small x and k, and for large lambda --
+ * condition number
+ *
+ *      -k (x/lambda)^k exp(-(x/lambda)^k)/[exp(-(x/lambda)^k) - 1]
+ *
+ * grows linearly with k, x^k, and lambda^{-k}.
+ */
+STATIC double
+cdf_weibull(double x, double lambda, double k)
+{
+  return -expm1(-pow(x/lambda, k));
+}
+
+/**
+ * Compute the SF of the Weibull(lambda, k) distribution.
+ * Well-conditioned for small x and k, and for large lambda --
+ * condition number
+ *
+ *      -k (x/lambda)^k
+ *
+ * grows linearly with k, x^k, and lambda^{-k}.
+ */
+STATIC double
+sf_weibull(double x, double lambda, double k)
+{
+  return exp(-pow(x/lambda, k));
+}
+
+/**
+ * Compute the inverse of the CDF of the Weibull(lambda, k)
+ * distribution.  Ill-conditioned for p near 1, and for k near 0;
+ * condition number is
+ *
+ *      (p/(1 - p))/(k log(1 - p)).
+ */
+STATIC double
+icdf_weibull(double p, double lambda, double k)
+{
+  return lambda*pow(-log1p(-p), 1/k);
+}
+
+/**
+ * Compute the inverse of the SF of the Weibull(lambda, k)
+ * distribution.  Ill-conditioned for p near 0, and for k near 0;
+ * condition number is
+ *
+ *      1/(k log(p)).
+ */
+STATIC double
+isf_weibull(double p, double lambda, double k)
+{
+  return lambda*pow(-log(p), 1/k);
+}
+
+/*
+ * GeneralizedPareto(mu, sigma, xi), supported on (mu, +\infty) for
+ * nonnegative xi, or (mu, mu - sigma/xi) for negative xi.
+ *
+ * Samples:
+ * = mu - sigma log U, if xi = 0;
+ * = mu + sigma (U^{-xi} - 1)/xi = mu + sigma*expm1(-xi log U)/xi, if xi =/= 0,
+ * where U is uniform on (0,1].
+ * = mu + sigma (e^{xi X} - 1)/xi,
+ * where X has standard exponential distribution.
+ *
+ * pdf(x) = sigma^{-1} (1 + xi (x - mu)/sigma)^{-(1 + 1/xi)}
+ * cdf(x) = 1 - (1 + xi (x - mu)/sigma)^{-1/xi}
+ *        = 1 - e^{-log(1 + xi (x - mu)/sigma)/xi}
+ *        --> 1 - e^{-(x - mu)/sigma}  as  xi --> 0
+ * sf(x) = (1 + xi (x - mu)/sigma)^{-1/xi}
+ *       --> e^{-(x - mu)/sigma}  as  xi --> 0
+ * icdf(p) = mu + sigma*(p^{-xi} - 1)/xi
+ *         = mu + sigma*expm1(-xi log p)/xi
+ *         --> mu + sigma*log p  as  xi --> 0
+ * isf(p) = mu + sigma*((1 - p)^{xi} - 1)/xi
+ *        = mu + sigma*expm1(-xi log1p(-p))/xi
+ *        --> mu + sigma*log1p(-p)  as  xi --> 0
+ */
+
+/**
+ * Compute the CDF of the GeneralizedPareto(mu, sigma, xi)
+ * distribution.  Well-conditioned everywhere.  For standard
+ * distribution (mu=0, sigma=1), condition number
+ *
+ *      (x/(1 + x xi)) / ((1 + x xi)^{1/xi} - 1)
+ *
+ * is bounded by 1, attained only at x = 0.
+ */
+STATIC double
+cdf_genpareto(double x, double mu, double sigma, double xi)
+{
+  double x_0 = (x - mu)/sigma;
+
+  /*
+   * log(1 + xi x_0)/xi
+   * = (-1/xi) \sum_{n=1}^\infty (-xi x_0)^n/n
+   * = (-1/xi) (-xi x_0 + \sum_{n=2}^\infty (-xi x_0)^n/n)
+   * = x_0 - (1/xi) \sum_{n=2}^\infty (-xi x_0)^n/n
+   * = x_0 - x_0 \sum_{n=2}^\infty (-xi x_0)^{n-1}/n
+   * = x_0 (1 - d),
+   *
+   * where d = \sum_{n=2}^\infty (-xi x_0)^{n-1}/n.  If |xi| <
+   * eps/4|x_0|, then
+   *
+   * |d| <= \sum_{n=2}^\infty (eps/4)^{n-1}/n
+   *     <= \sum_{n=2}^\infty (eps/4)^{n-1}
+   *      = \sum_{n=1}^\infty (eps/4)^n
+   *      = (eps/4) \sum_{n=0}^\infty (eps/4)^n
+   *      = (eps/4)/(1 - eps/4)
+   *      < eps/2
+   *
+   * for any 0 < eps < 2.  Thus, the relative error of x_0 from
+   * log(1 + xi x_0)/xi is bounded by eps.
+   */
+  if (fabs(xi) < 1e-17/x_0)
+    return -expm1(-x_0);
+  else
+    return -expm1(-log1p(xi*x_0)/xi);
+}
+
+/**
+ * Compute the SF of the GeneralizedPareto(mu, sigma, xi) distribution.
+ * For standard distribution (mu=0, sigma=1), ill-conditioned for xi
+ * near 0; condition number
+ *
+ *      -x (1 + x xi)^{(-1 - xi)/xi}/(1 + x xi)^{-1/xi}
+ *      = -x (1 + x xi)^{-1/xi - 1}/(1 + x xi)^{-1/xi}
+ *      = -(x/(1 + x xi)) (1 + x xi)^{-1/xi}/(1 + x xi)^{-1/xi}
+ *      = -x/(1 + x xi)
+ *
+ * is bounded by 1/xi.
+ */
+STATIC double
+sf_genpareto(double x, double mu, double sigma, double xi)
+{
+  double x_0 = (x - mu)/sigma;
+
+  if (fabs(xi) < 1e-17/x_0)
+    return exp(-x_0);
+  else
+    return exp(-log1p(xi*x_0)/xi);
+}
+
+/**
+ * Compute the inverse of the CDF of the GeneralizedPareto(mu, sigma,
+ * xi) distribution.  Ill-conditioned for p near 1; condition number is
+ *
+ *      xi (p/(1 - p))/(1 - (1 - p)^xi)
+ */
+STATIC double
+icdf_genpareto(double p, double mu, double sigma, double xi)
+{
+  /*
+   * To compute f(xi) = (U^{-xi} - 1)/xi = (e^{-xi log U} - 1)/xi
+   * for xi near zero (note f(xi) --> -log U as xi --> 0), write
+   * the absolutely convergent Taylor expansion
+   *
+   * f(xi) = (1/xi)*(-xi log U + \sum_{n=2}^\infty (-xi log U)^n/n!
+   *       = -log U + (1/xi)*\sum_{n=2}^\infty (-xi log U)^n/n!
+   *       = -log U + \sum_{n=2}^\infty xi^{n-1} (-log U)^n/n!
+   *       = -log U - log U \sum_{n=2}^\infty (-xi log U)^{n-1}/n!
+   *       = -log U (1 + \sum_{n=2}^\infty (-xi log U)^{n-1}/n!).
+   *
+   * Let d = \sum_{n=2}^\infty (-xi log U)^{n-1}/n!.  What do we
+   * lose if we discard it and use -log U as an approximation to
+   * f(xi)?  If |xi| < eps/-4log U, then
+   *
+   * |d| <= \sum_{n=2}^\infty |xi log U|^{n-1}/n!
+   *     <= \sum_{n=2}^\infty (eps/4)^{n-1}/n!
+   *     <= \sum_{n=1}^\infty (eps/4)^n
+   *      = (eps/4) \sum_{n=0}^\infty (eps/4)^n
+   *      = (eps/4)/(1 - eps/4)
+   *      < eps/2,
+   *
+   * for any 0 < eps < 2.  Hence, as long as |xi| < eps/-2log U,
+   * f(xi) = -log U (1 + d) for |d| <= eps/2.  |d| is the
+   * relative error of f(xi) from -log U; from this bound, the
+   * relative error of -log U from f(xi) is at most (eps/2)/(1 -
+   * eps/2) = eps/2 + (eps/2)^2 + (eps/2)^3 + ... < eps for 0 <
+   * eps < 1.  Since -log U < 1000 for all U in (0, 1] in
+   * binary64 floating-point, we can safely cut xi off at 1e-20 <
+   * eps/4000 and attain <1ulp error from series truncation.
+   */
+  if (fabs(xi) <= 1e-20)
+    return mu - sigma*log1p(-p);
+  else
+    return mu + sigma*expm1(-xi*log1p(-p))/xi;
+}
+
+/**
+ * Compute the inverse of the SF of the GeneralizedPareto(mu, sigma,
+ * xi) distribution.  Ill-conditioned for p near 1; conditon number is
+ *
+ *      -xi/(1 - p^{-xi})
+ */
+STATIC double
+isf_genpareto(double p, double mu, double sigma, double xi)
+{
+  if (fabs(xi) <= 1e-20)
+    return mu - sigma*log(p);
+  else
+    return mu + sigma*expm1(-xi*log(p))/xi;
+}
+
+/*******************************************************************/
+
+/**
+ * Deterministic samplers, parametrized by uniform integer and (0,1]
+ * samples.  No guarantees are made about _which_ mapping from the
+ * integer and (0,1] samples these use; all that is guaranteed is the
+ * distribution of the outputs conditioned on a uniform distribution on
+ * the inputs.  The automatic tests in test_prob_distr.c double-check
+ * the particular mappings we use.
+ *
+ * Beware: Unlike random_uniform_01(), these are not guaranteed to be
+ * supported on all possible outputs.  See Ilya Mironov, `On the
+ * Significance of the Least Significant Bits for Differential
+ * Privacy', for an example of what can go wrong if you try to use
+ * these to conceal information from an adversary but you expose the
+ * specific full-precision floating-point values.
+ *
+ * Note: None of these samplers use rejection sampling; they are all
+ * essentially inverse-CDF transforms with tweaks.  If you were to add,
+ * say, a Gamma sampler with the Marsaglia-Tsang method, you would have
+ * to parametrize it by a potentially infinite stream of uniform (and
+ * perhaps normal) samples rather than a fixed number, which doesn't
+ * make for quite as nice automatic testing as for these.
+ */
+
+/**
+ * Deterministically sample from the interval [a, b], indexed by a
+ * uniform random floating-point number p0 in (0, 1].
+ *
+ * Note that even if p0 is nonzero, the result may be equal to a, if
+ * ulp(a)/2 is nonnegligible, e.g. if a = 1.  For maximum resolution,
+ * arrange |a| <= |b|.
+ */
+STATIC double
+sample_uniform_interval(double p0, double a, double b)
+{
+  /*
+   * XXX Prove that the distribution is, in fact, uniform on
+   * [a,b], particularly around p0 = 1, or at least has very
+   * small deviation from uniform, quantified appropriately
+   * (e.g., like in Monahan 1984, or by KL divergence).  It
+   * almost certainly does but it would be nice to quantify the
+   * error.
+   */
+  if ((a <= 0 && 0 <= b) || (b <= 0 && 0 <= a)) {
+    /*
+     * When ab < 0, (1 - t) a + t b is monotonic, since for
+     * a <= b it is a sum of nondecreasing functions of t,
+     * and for b <= a, of nonincreasing functions of t.
+     * Further, clearly at 0 and 1 it attains a and b,
+     * respectively.  Hence it is bounded within [a, b].
+     */
+    return (1 - p0)*a + p0*b;
+  } else {
+    /*
+     * a + (b - a) t is monotonic -- it is obviously a
+     * nondecreasing function of t for a <= b.  Further, it
+     * attains a at 0, and while it may overshoot b at 1,
+     * we have a
+     *
+     * Theorem.  If 0 <= t < 1, then the floating-point
+     * evaluation of a + (b - a) t is bounded in [a, b].
+     *
+     * Lemma 1.  If 0 <= t < 1 is a floating-point number,
+     * then for any normal floating-point number x except
+     * the smallest in magnitude, |round(x*t)| < |x|.
+     *
+     * Proof.  WLOG, assume x >= 0.  Since the rounding
+     * function and t |---> x*t are nondecreasing, their
+     * composition t |---> round(x*t) is also
+     * nondecreasing, so it suffices to consider the
+     * largest floating-point number below 1, in particular
+     * t = 1 - ulp(1)/2.
+     *
+     * Case I: If x is a power of two, then the next
+     * floating-point number below x is x - ulp(x)/2 = x -
+     * x*ulp(1)/2 = x*(1 - ulp(1)/2) = x*t, so, since x*t
+     * is a floating-point number, multiplication is exact,
+     * and thus round(x*t) = x*t < x.
+     *
+     * Case II: If x is not a power of two, then the
+     * greatest lower bound of real numbers rounded to x is
+     * x - ulp(x)/2 = x - ulp(T(x))/2 = x - T(x)*ulp(1)/2,
+     * where T(X) is the largest power of two below x.
+     * Anything below this bound is rounded to a
+     * floating-point number smaller than x, and x*t = x*(1
+     * - ulp(1)/2) = x - x*ulp(1)/2 < x - T(x)*ulp(1)/2
+     * since T(x) < x, so round(x*t) < x*t < x.  QED.
+     *
+     * Lemma 2.  If x and y are subnormal, then round(x +
+     * y) = x + y.
+     *
+     * Proof.  It is a matter of adding the significands,
+     * since if we treat subnormals as having an implicit
+     * zero bit before the `binary' point, their exponents
+     * are all the same.  There is at most one carry/borrow
+     * bit, which can always be acommodated either in a
+     * subnormal, or, at largest, in the implicit one bit
+     * of a normal.
+     *
+     * Lemma 3.  Let x and y be floating-point numbers.  If
+     * round(x - y) is subnormal or zero, then it is equal
+     * to x - y.
+     *
+     * Proof.  Case I (equal): round(x - y) = 0 iff x = y;
+     * hence if round(x - y) = 0, then round(x - y) = 0 = x
+     * - y.
+     *
+     * Case II (subnormal/subnormal): If x and y are both
+     * subnormal, this follows directly from Lemma 2.
+     *
+     * Case IIIa (normal/subnormal): If x is normal and y
+     * is subnormal, then x and y must share sign, or else
+     * x - y would be larger than x and thus rounded to
+     * normal.  If s is the smallest normal positive
+     * floating-point number, |x| < 2s since by
+     * construction 2s - |y| is normal for all subnormal y.
+     * This means that x and y must have the same exponent,
+     * so the difference is the difference of significands,
+     * which is exact.
+     *
+     * Case IIIb (subnormal/normal): Same as case IIIa for
+     * -(y - x).
+     *
+     * Case IV (normal/normal): If x and y are both normal,
+     * then they must share sign, or else x - y would be
+     * larger than x and thus rounded to normal.  Note that
+     * |y| < 2|x|, for if |y| >= 2|x|, then |x| - |y| <=
+     * -|x| but -|x| is normal like x.  Also, |x|/2 < |y|:
+     * if |x|/2 is subnormal, it must hold because y is
+     * normal; if |x|/2 is normal, then |x|/2 >= s, so
+     * since |x| - |y| < s,
+     *
+     *  |x|/2 = |x| - |x|/2 <= |x| - s <= |y|;
+     *
+     * that is, |x|/2 < |y| < 2|x|, so by the Sterbenz
+     * lemma, round(x - y) = x - y.  QED.
+     *
+     * Proof of theorem.  WLOG, assume 0 <= a <= b.  Since
+     * round(a + round(round(b - a)*t) is nondecreasing in
+     * t and attains a at 0, the lower end of the bound is
+     * trivial; we must show the upper end of the bound
+     * strictly.  It suffices to show this for the largest
+     * floating-point number below 1, namely 1 - ulp(1)/2.
+     *
+     * Case I: round(b - a) is normal.  Then it is at most
+     * the smallest floating-point number above b - a.  By
+     * Lemma 1, round(round(b - a)*t) < round(b - a).
+     * Since the inequality is strict, and since
+     * round(round(b - a)*t) is a floating-point number
+     * below round(b - a), and since there are no
+     * floating-point numbers between b - a and round(b -
+     * a), we must have round(round(b - a)*t) < b - a.
+     * Then since y |---> round(a + y) is nondecreasing, we
+     * must have
+     *
+     *  round(a + round(round(b - a)*t))
+     *  <= round(a + (b - a))
+     *   = round(b) = b.
+     *
+     * Case II: round(b - a) is subnormal.  In this case,
+     * Lemma 1 falls apart -- we are not guaranteed the
+     * strict inequality.  However, by Lemma 3, the
+     * difference is exact: round(b - a) = b - a.  Thus,
+     *
+     *  round(a + round(round(b - a)*t))
+     *  <= round(a + round((b - a)*t))
+     *  <= round(a + (b - a))
+     *   = round(b)
+     *   = b,
+     *
+     * QED.
+     */
+
+    /* p0 is restricted to [0,1], but we use >= to silence -Wfloat-equal.  */
+    if (p0 >= 1)
+      return b;
+    return a + (b - a)*p0;
+  }
+}
+
+/**
+ * Deterministically sample from the standard logistic distribution,
+ * indexed by a uniform random 32-bit integer s and uniform random
+ * floating-point numbers t and p0 in (0, 1].
+ */
+STATIC double
+sample_logistic(uint32_t s, double t, double p0)
+{
+  double sign = (s & 1) ? -1 : +1;
+  double r;
+
+  /*
+   * We carve up the interval (0, 1) into subregions to compute
+   * the inverse CDF precisely:
+   *
+   * A = (0, 1/(1 + e)] ---> (-\infty, -1]
+   * B = [1/(1 + e), 1/2] ---> [-1, 0]
+   * C = [1/2, 1 - 1/(1 + e)] ---> [0, 1]
+   * D = [1 - 1/(1 + e), 1) ---> [1, +\infty)
+   *
+   * Cases D and C are mirror images of cases A and B,
+   * respectively, so we choose between them by the sign chosen
+   * by a fair coin toss.  We choose between cases A and B by a
+   * coin toss weighted by
+   *
+   *    2/(1 + e) = 1 - [1/2 - 1/(1 + e)]/(1/2):
+   *
+   * if it comes up heads, scale p0 into a uniform (0, 1/(1 + e)]
+   * sample p; if it comes up tails, scale p0 into a uniform (0,
+   * 1/2 - 1/(1 + e)] sample and compute the inverse CDF of p =
+   * 1/2 - p0.
+   */
+  if (t <= 2/(1 + exp(1))) {
+    /* p uniform in (0, 1/(1 + e)], represented by p.  */
+    p0 /= 1 + exp(1);
+    r = logit(p0);
+  } else {
+    /*
+     * p uniform in [1/(1 + e), 1/2), actually represented
+     * by p0 = 1/2 - p uniform in (0, 1/2 - 1/(1 + e)], so
+     * that p = 1/2 - p.
+     */
+    p0 *= 0.5 - 1/(1 + exp(1));
+    r = logithalf(p0);
+  }
+
+  /*
+   * We have chosen from the negative half of the standard
+   * logistic distribution, which is symmetric with the positive
+   * half.  Now use the sign to choose uniformly between them.
+   */
+  return sign*r;
+}
+
+/**
+ * Deterministically sample from the logistic distribution scaled by
+ * sigma and translated by mu.
+ */
+static double
+sample_logistic_locscale(uint32_t s, double t, double p0, double mu,
+    double sigma)
+{
+
+  return mu + sigma*sample_logistic(s, t, p0);
+}
+
+/**
+ * Deterministically sample from the standard log-logistic
+ * distribution, indexed by a uniform random 32-bit integer s and a
+ * uniform random floating-point number p0 in (0, 1].
+ */
+STATIC double
+sample_log_logistic(uint32_t s, double p0)
+{
+
+  /*
+   * Carve up the interval (0, 1) into (0, 1/2] and [1/2, 1); the
+   * condition numbers of the icdf and the isf coincide at 1/2.
+   */
+  p0 *= 0.5;
+  if ((s & 1) == 0) {
+    /* p = p0 in (0, 1/2] */
+    return p0/(1 - p0);
+  } else {
+    /* p = 1 - p0 in [1/2, 1) */
+    return (1 - p0)/p0;
+  }
+}
+
+/**
+ * Deterministically sample from the log-logistic distribution with
+ * scale alpha and shape beta.
+ */
+static double
+sample_log_logistic_scaleshape(uint32_t s, double p0, double alpha,
+    double beta)
+{
+  double x = sample_log_logistic(s, p0);
+
+  return alpha*pow(x, 1/beta);
+}
+
+/**
+ * Deterministically sample from the standard exponential distribution,
+ * indexed by a uniform random 32-bit integer s and a uniform random
+ * floating-point number p0 in (0, 1].
+ */
+static double
+sample_exponential(uint32_t s, double p0)
+{
+  /*
+   * We would like to evaluate log(p) for p near 0, and log1p(-p)
+   * for p near 1.  Simply carve the interval into (0, 1/2] and
+   * [1/2, 1) by a fair coin toss.
+   */
+  p0 *= 0.5;
+  if ((s & 1) == 0)
+    /* p = p0 in (0, 1/2] */
+    return -log(p0);
+  else
+    /* p = 1 - p0 in [1/2, 1) */
+    return -log1p(-p0);
+}
+
+/**
+ * Deterministically sample from a Weibull distribution with scale
+ * lambda and shape k -- just an exponential with a shape parameter in
+ * addition to a scale parameter.  (Yes, lambda really is the scale,
+ * _not_ the rate.)
+ */
+STATIC double
+sample_weibull(uint32_t s, double p0, double lambda, double k)
+{
+
+  return lambda*pow(sample_exponential(s, p0), 1/k);
+}
+
+/**
+ * Deterministically sample from the generalized Pareto distribution
+ * with shape xi, indexed by a uniform random 32-bit integer s and a
+ * uniform random floating-point number p0 in (0, 1].
+ */
+STATIC double
+sample_genpareto(uint32_t s, double p0, double xi)
+{
+  double x = sample_exponential(s, p0);
+
+  /*
+   * Write f(xi) = (e^{xi x} - 1)/xi for xi near zero as the
+   * absolutely convergent Taylor series
+   *
+   * f(x) = (1/xi) (xi x + \sum_{n=2}^\infty (xi x)^n/n!)
+   *      = x + (1/xi) \sum_{n=2}^\inty (xi x)^n/n!
+   *      = x + \sum_{n=2}^\infty xi^{n-1} x^n/n!
+   *      = x + x \sum_{n=2}^\infty (xi x)^{n-1}/n!
+   *      = x (1 + \sum_{n=2}^\infty (xi x)^{n-1}/n!).
+   *
+   * d = \sum_{n=2}^\infty (xi x)^{n-1}/n! is the relative error
+   * of f(x) from x.  If |xi| < eps/4x, then
+   *
+   * |d| <= \sum_{n=2}^\infty |xi x|^{n-1}/n!
+   *     <= \sum_{n=2}^\infty (eps/4)^{n-1}/n!
+   *     <= \sum_{n=1}^\infty (eps/4)
+   *      = (eps/4) \sum_{n=0}^\infty (eps/4)^n
+   *      = (eps/4)/(1 - eps/4)
+   *      < eps/2,
+   *
+   * for any 0 < eps < 2.  Hence, as long as |xi| < eps/2x, f(xi)
+   * = x (1 + d) for |d| <= eps/2, so x = f(xi) (1 + d') for |d'|
+   * <= eps.  What bound should we use for x?
+   *
+   * - If x is exponentially distributed, x > 200 with
+   *   probability below e^{-200} << 2^{-256}, i.e. never.
+   *
+   * - If x is computed by -log(U) for U in (0, 1], x is
+   *   guaranteed to be below 1000 in IEEE 754 binary64
+   *   floating-point.
+   *
+   * We can safely cut xi off at 1e-20 < eps/4000 and attain an
+   * error bounded by 0.5 ulp for this expression.
+   */
+  return (fabs(xi) < 1e-20 ? x : expm1(xi*x)/xi);
+}
+
+/**
+ * Deterministically sample from a generalized Pareto distribution with
+ * shape xi, scaled by sigma and translated by mu.
+ */
+static double
+sample_genpareto_locscale(uint32_t s, double p0, double mu, double sigma,
+    double xi)
+{
+
+  return mu + sigma*sample_genpareto(s, p0, xi);
+}
+
+/**
+ * Deterministically sample from the geometric distribution with
+ * per-trial success probability p.
+ *
+ * XXX Quantify the error (KL divergence?) of this
+ * ceiling-of-exponential sampler from a true geometric distribution,
+ * which we could get by rejection sampling.  Relevant papers:
+ *
+ *      John F. Monahan, `Accuracy in Random Number Generation',
+ *      Mathematics of Computation 45(172), October 1984, pp. 559--568.
+*https://pdfs.semanticscholar.org/aca6/74b96da1df77b2224e8cfc5dd6d61a471632.pdf
+ *
+ *      Karl Bringmann and Tobias Friedrich, `Exact and Efficient
+ *      Generation of Geometric Random Variates and Random Graphs', in
+ *      Proceedings of the 40th International Colloaquium on Automata,
+ *      Languages, and Programming -- ICALP 2013, Springer LNCS 7965,
+ *      pp.267--278.
+ *      https://doi.org/10.1007/978-3-642-39206-1_23
+ *      https://people.mpi-inf.mpg.de/~kbringma/paper/2013ICALP-1.pdf
+ */
+static double
+sample_geometric(uint32_t s, double p0, double p)
+{
+  double x = sample_exponential(s, p0);
+
+  /* This is actually a check against 1, but we do >= so that the compiler
+     does not raise a -Wfloat-equal */
+  if (p >= 1)
+    return 1;
+
+  return (-x/log1p(-p));
+}
+
+/*******************************************************************/
+
+/** Public API for probability distributions:
+ *
+ *  For each probability distribution we define each public functions
+ *  (sample/cdf/sf/icdf/isf) as part of its dist_ops structure.
+ */
+
+/** Functions for uniform distribution */
+const struct dist_ops uniform_ops = {
+  .name = "uniform",
+  .sample = uniform_sample,
+  .cdf = uniform_cdf,
+  .sf = uniform_sf,
+  .icdf = uniform_icdf,
+  .isf = uniform_isf,
+};
+
+double
+uniform_sample(const struct dist *dist)
+{
+  const struct uniform *U = const_container_of(dist, struct uniform,
+    base);
+  double p0 = random_uniform_01();
+
+  return sample_uniform_interval(p0, U->a, U->b);
+}
+
+double
+uniform_cdf(const struct dist *dist, double x)
+{
+  const struct uniform *U = const_container_of(dist, struct uniform,
+    base);
+
+  if (x < U->a)
+    return 0;
+  else if (x < U->b)
+    return (x - U->a)/(U->b - U->a);
+  else
+    return 1;
+}
+
+double
+uniform_sf(const struct dist *dist, double x)
+{
+  const struct uniform *U = const_container_of(dist, struct uniform,
+    base);
+
+  if (x > U->b)
+    return 0;
+  else if (x > U->a)
+    return (U->b - x)/(U->b - U->a);
+  else
+    return 1;
+}
+
+double
+uniform_icdf(const struct dist *dist, double p)
+{
+  const struct uniform *U = const_container_of(dist, struct uniform,
+    base);
+  double w = U->b - U->a;
+
+  return (p < 0.5 ? (U->a + w*p) : (U->b - w*(1 - p)));
+}
+
+double
+uniform_isf(const struct dist *dist, double p)
+{
+  const struct uniform *U = const_container_of(dist, struct uniform,
+    base);
+  double w = U->b - U->a;
+
+  return (p < 0.5 ? (U->b - w*p) : (U->a + w*(1 - p)));
+}
+
+/** Functions for logistic distribution: */
+const struct dist_ops logistic_ops = {
+  .name = "logistic",
+  .sample = logistic_sample,
+  .cdf = logistic_cdf,
+  .sf = logistic_sf,
+  .icdf = logistic_icdf,
+  .isf = logistic_isf,
+};
+
+double
+logistic_sample(const struct dist *dist)
+{
+  const struct logistic *L = const_container_of(dist, struct logistic,
+    base);
+  uint32_t s = crypto_rand_uint32();
+  double t = random_uniform_01();
+  double p0 = random_uniform_01();
+
+  return sample_logistic_locscale(s, t, p0, L->mu, L->sigma);
+}
+
+double
+logistic_cdf(const struct dist *dist, double x)
+{
+  const struct logistic *L = const_container_of(dist, struct logistic,
+    base);
+
+  return cdf_logistic(x, L->mu, L->sigma);
+}
+
+double
+logistic_sf(const struct dist *dist, double x)
+{
+  const struct logistic *L = const_container_of(dist, struct logistic,
+    base);
+
+  return sf_logistic(x, L->mu, L->sigma);
+}
+
+double
+logistic_icdf(const struct dist *dist, double p)
+{
+  const struct logistic *L = const_container_of(dist, struct logistic,
+    base);
+
+  return icdf_logistic(p, L->mu, L->sigma);
+}
+
+double
+logistic_isf(const struct dist *dist, double p)
+{
+  const struct logistic *L = const_container_of(dist, struct logistic,
+    base);
+
+  return isf_logistic(p, L->mu, L->sigma);
+}
+
+/** Functions for log-logistic distribution: */
+const struct dist_ops log_logistic_ops = {
+  .name = "log logistic",
+  .sample = log_logistic_sample,
+  .cdf = log_logistic_cdf,
+  .sf = log_logistic_sf,
+  .icdf = log_logistic_icdf,
+  .isf = log_logistic_isf,
+};
+
+double
+log_logistic_sample(const struct dist *dist)
+{
+  const struct log_logistic *LL = const_container_of(dist, struct
+    log_logistic, base);
+  uint32_t s = crypto_rand_uint32();
+  double p0 = random_uniform_01();
+
+  return sample_log_logistic_scaleshape(s, p0, LL->alpha, LL->beta);
+}
+
+double
+log_logistic_cdf(const struct dist *dist, double x)
+{
+  const struct log_logistic *LL = const_container_of(dist,
+    struct log_logistic, base);
+
+  return cdf_log_logistic(x, LL->alpha, LL->beta);
+}
+
+double
+log_logistic_sf(const struct dist *dist, double x)
+{
+  const struct log_logistic *LL = const_container_of(dist,
+    struct log_logistic, base);
+
+  return sf_log_logistic(x, LL->alpha, LL->beta);
+}
+
+double
+log_logistic_icdf(const struct dist *dist, double p)
+{
+  const struct log_logistic *LL = const_container_of(dist,
+    struct log_logistic, base);
+
+  return icdf_log_logistic(p, LL->alpha, LL->beta);
+}
+
+double
+log_logistic_isf(const struct dist *dist, double p)
+{
+  const struct log_logistic *LL = const_container_of(dist,
+    struct log_logistic, base);
+
+  return isf_log_logistic(p, LL->alpha, LL->beta);
+}
+
+/** Functions for Weibull distribution */
+const struct dist_ops weibull_ops = {
+  .name = "Weibull",
+  .sample = weibull_sample,
+  .cdf = weibull_cdf,
+  .sf = weibull_sf,
+  .icdf = weibull_icdf,
+  .isf = weibull_isf,
+};
+
+double
+weibull_sample(const struct dist *dist)
+{
+  const struct weibull *W = const_container_of(dist, struct weibull,
+    base);
+  uint32_t s = crypto_rand_uint32();
+  double p0 = random_uniform_01();
+
+  return sample_weibull(s, p0, W->lambda, W->k);
+}
+
+double
+weibull_cdf(const struct dist *dist, double x)
+{
+  const struct weibull *W = const_container_of(dist, struct weibull,
+    base);
+
+  return cdf_weibull(x, W->lambda, W->k);
+}
+
+double
+weibull_sf(const struct dist *dist, double x)
+{
+  const struct weibull *W = const_container_of(dist, struct weibull,
+    base);
+
+  return sf_weibull(x, W->lambda, W->k);
+}
+
+double
+weibull_icdf(const struct dist *dist, double p)
+{
+  const struct weibull *W = const_container_of(dist, struct weibull,
+    base);
+
+  return icdf_weibull(p, W->lambda, W->k);
+}
+
+double
+weibull_isf(const struct dist *dist, double p)
+{
+  const struct weibull *W = const_container_of(dist, struct weibull,
+    base);
+
+  return isf_weibull(p, W->lambda, W->k);
+}
+
+/** Functions for generalized Pareto distributions */
+const struct dist_ops genpareto_ops = {
+  .name = "generalized Pareto",
+  .sample = genpareto_sample,
+  .cdf = genpareto_cdf,
+  .sf = genpareto_sf,
+  .icdf = genpareto_icdf,
+  .isf = genpareto_isf,
+};
+
+double
+genpareto_sample(const struct dist *dist)
+{
+  const struct genpareto *GP = const_container_of(dist, struct genpareto,
+    base);
+  uint32_t s = crypto_rand_uint32();
+  double p0 = random_uniform_01();
+
+  return sample_genpareto_locscale(s, p0, GP->mu, GP->sigma, GP->xi);
+}
+
+double
+genpareto_cdf(const struct dist *dist, double x)
+{
+  const struct genpareto *GP = const_container_of(dist, struct genpareto,
+    base);
+
+  return cdf_genpareto(x, GP->mu, GP->sigma, GP->xi);
+}
+
+double
+genpareto_sf(const struct dist *dist, double x)
+{
+  const struct genpareto *GP = const_container_of(dist, struct genpareto,
+    base);
+
+  return sf_genpareto(x, GP->mu, GP->sigma, GP->xi);
+}
+
+double
+genpareto_icdf(const struct dist *dist, double p)
+{
+  const struct genpareto *GP = const_container_of(dist, struct genpareto,
+    base);
+
+  return icdf_genpareto(p, GP->mu, GP->sigma, GP->xi);
+}
+
+double
+genpareto_isf(const struct dist *dist, double p)
+{
+  const struct genpareto *GP = const_container_of(dist, struct genpareto,
+    base);
+
+  return isf_genpareto(p, GP->mu, GP->sigma, GP->xi);
+}
+
+/* Deterministically sample from the geometric distribution with
+ * per-trial success probability p. */
+double
+geometric_sample(double p)
+{
+  uint32_t s = crypto_rand_uint32();
+  double p0 = random_uniform_01();
+  return sample_geometric(s, p0, p);
+}
+
diff --git a/src/lib/math/prob_distr.h b/src/lib/math/prob_distr.h
new file mode 100644
index 000000000..c2fd6c74b
--- /dev/null
+++ b/src/lib/math/prob_distr.h
@@ -0,0 +1,156 @@
+
+/**
+ * \file prob_distr.h
+ *
+ * \brief Header for prob_distr.c
+ **/
+
+#ifndef TOR_PROB_DISTR_H
+#define TOR_PROB_DISTR_H
+
+#include "lib/cc/compat_compiler.h"
+#include "lib/cc/torint.h"
+#include "lib/testsupport/testsupport.h"
+
+/**
+ * Container for distribution parameters for sampling, CDF, &c.
+ */
+struct dist {
+  const struct dist_ops *ops;
+};
+
+#define DIST_BASE(OPS)  { .ops = (OPS) }
+
+struct dist_ops {
+  const char *name;
+  double (*sample)(const struct dist *);
+  double (*cdf)(const struct dist *, double x);
+  double (*sf)(const struct dist *, double x);
+  double (*icdf)(const struct dist *, double p);
+  double (*isf)(const struct dist *, double p);
+};
+
+/* Geometric distribution */
+
+double geometric_sample(double p);
+
+/* Pareto distribution */
+
+struct genpareto {
+  struct dist base;
+  double mu;
+  double sigma;
+  double xi;
+};
+
+double genpareto_sample(const struct dist *dist);
+double genpareto_cdf(const struct dist *dist, double x);
+double genpareto_sf(const struct dist *dist, double x);
+double genpareto_icdf(const struct dist *dist, double p);
+double genpareto_isf(const struct dist *dist, double p);
+
+extern const struct dist_ops genpareto_ops;
+
+/* Weibull distribution */
+
+struct weibull {
+  struct dist base;
+  double lambda;
+  double k;
+};
+
+double weibull_sample(const struct dist *dist);
+double weibull_cdf(const struct dist *dist, double x);
+double weibull_sf(const struct dist *dist, double x);
+double weibull_icdf(const struct dist *dist, double p);
+double weibull_isf(const struct dist *dist, double p);
+
+extern const struct dist_ops weibull_ops;
+
+/* Log-logistic distribution */
+
+struct log_logistic {
+  struct dist base;
+  double alpha;
+  double beta;
+};
+
+double log_logistic_sample(const struct dist *dist);
+double log_logistic_cdf(const struct dist *dist, double x);
+double log_logistic_sf(const struct dist *dist, double x);
+double log_logistic_icdf(const struct dist *dist, double p);
+double log_logistic_isf(const struct dist *dist, double p);
+
+extern const struct dist_ops log_logistic_ops;
+
+/* Logistic distribution */
+
+struct logistic {
+  struct dist base;
+  double mu;
+  double sigma;
+};
+
+double logistic_sample(const struct dist *dist);
+double logistic_cdf(const struct dist *dist, double x);
+double logistic_sf(const struct dist *dist, double x);
+double logistic_icdf(const struct dist *dist, double p);
+double logistic_isf(const struct dist *dist, double p);
+
+extern const struct dist_ops logistic_ops;
+
+/* Uniform distribution */
+
+struct uniform {
+  struct dist base;
+  double a;
+  double b;
+};
+
+double uniform_sample(const struct dist *dist);
+double uniform_cdf(const struct dist *dist, double x);
+double uniform_sf(const struct dist *dist, double x);
+double uniform_icdf(const struct dist *dist, double p);
+double uniform_isf(const struct dist *dist, double p);
+
+extern const struct dist_ops uniform_ops;
+
+/** Only by unittests */
+
+#ifdef PROB_DISTR_PRIVATE
+
+STATIC double logithalf(double p0);
+STATIC double logit(double p);
+
+STATIC double random_uniform_01(void);
+
+STATIC double logistic(double x);
+STATIC double cdf_logistic(double x, double mu, double sigma);
+STATIC double sf_logistic(double x, double mu, double sigma);
+STATIC double icdf_logistic(double p, double mu, double sigma);
+STATIC double isf_logistic(double p, double mu, double sigma);
+STATIC double sample_logistic(uint32_t s, double t, double p0);
+
+STATIC double cdf_log_logistic(double x, double alpha, double beta);
+STATIC double sf_log_logistic(double x, double alpha, double beta);
+STATIC double icdf_log_logistic(double p, double alpha, double beta);
+STATIC double isf_log_logistic(double p, double alpha, double beta);
+STATIC double sample_log_logistic(uint32_t s, double p0);
+
+STATIC double cdf_weibull(double x, double lambda, double k);
+STATIC double sf_weibull(double x, double lambda, double k);
+STATIC double icdf_weibull(double p, double lambda, double k);
+STATIC double isf_weibull(double p, double lambda, double k);
+STATIC double sample_weibull(uint32_t s, double p0, double lambda, double k);
+
+STATIC double sample_uniform_interval(double p0, double a, double b);
+
+STATIC double cdf_genpareto(double x, double mu, double sigma, double xi);
+STATIC double sf_genpareto(double x, double mu, double sigma, double xi);
+STATIC double icdf_genpareto(double p, double mu, double sigma, double xi);
+STATIC double isf_genpareto(double p, double mu, double sigma, double xi);
+STATIC double sample_genpareto(uint32_t s, double p0, double xi);
+
+#endif
+
+#endif
diff --git a/src/test/include.am b/src/test/include.am
index 4da0b8439..b276500fd 100644
--- a/src/test/include.am
+++ b/src/test/include.am
@@ -157,6 +157,7 @@ src_test_test_SOURCES += \
 	src/test/test_periodic_event.c \
 	src/test/test_policy.c \
 	src/test/test_process.c \
+	src/test/test_prob_distr.c \
 	src/test/test_procmon.c \
 	src/test/test_proto_http.c \
 	src/test/test_proto_misc.c \
@@ -207,6 +208,7 @@ src_test_test_slow_SOURCES += \
 	src/test/test_slow.c \
 	src/test/test_crypto_slow.c \
 	src/test/test_process_slow.c \
+	src/test/test_prob_distr.c \
 	src/test/testing_common.c \
 	src/test/testing_rsakeys.c \
 	src/ext/tinytest.c
diff --git a/src/test/prob_distr_mpfr_ref.c b/src/test/prob_distr_mpfr_ref.c
new file mode 100644
index 000000000..4e64d731c
--- /dev/null
+++ b/src/test/prob_distr_mpfr_ref.c
@@ -0,0 +1,64 @@
+/* Copyright 2012-2018, The Tor Project, Inc
+ * See LICENSE for licensing information */
+
+/** prob_distr_mpfr_ref.c
+ *
+ * Example reference file for GNU MPFR vectors tested in test_prob_distr.c .
+ * Code by Riastradh.
+ */
+
+#include <complex.h>
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+
+/* Must come after <stdio.h> so we get mpfr_printf.  */
+#include <mpfr.h>
+
+/*  gcc -o mpfr prob_distr_mpfr_ref.c -lmpfr -lm */
+
+/* Computes logit(p) for p = .49999 */
+int
+main(void)
+{
+  mpfr_t p, q, r;
+  mpfr_init(p);
+  mpfr_set_prec(p, 200);
+  mpfr_init(q);
+  mpfr_set_prec(q, 200);
+  mpfr_init(r);
+  mpfr_set_prec(r, 200);
+  mpfr_set_d(p, .49999, MPFR_RNDN);
+  mpfr_set_d(q, 1, MPFR_RNDN);
+  /* r := q - p = 1 - p */
+  mpfr_sub(r, q, p, MPFR_RNDN);
+  /* q := p/r = p/(1 - p) */
+  mpfr_div(q, p, r, MPFR_RNDN);
+  /* r := log(q) = log(p/(1 - p)) */
+  mpfr_log(r, q, MPFR_RNDN);
+  mpfr_printf("mpfr 200-bit\t%.128Rg\n", r);
+
+  /*
+   * Print a double approximation to logit three different ways.  All
+   * three agree bit for bit on the libms I tried, with the nextafter
+   * adjustment (which is well within the 10 eps relative error bound
+   * advertised).  Apparently I must have used the Goldberg expression
+   * for what I wrote down in the test case.
+   */
+  printf("mpfr 53-bit\t%.17g\n", nextafter(mpfr_get_d(r, MPFR_RNDN), 0), 0);
+  volatile double p0 = .49999;
+  printf("log1p\t\t%.17g\n", nextafter(-log1p((1 - 2*p0)/p0), 0));
+  volatile double x = (1 - 2*p0)/p0;
+  volatile double xp1 = x + 1;
+  printf("Goldberg\t%.17g\n", -x*log(xp1)/(xp1 - 1));
+
+  /*
+   * Print a bad approximation, using the naive expression, to see a
+   * lot of wrong digits, far beyond the 10 eps relative error attained
+   * by -log1p((1 - 2*p)/p).
+   */
+  printf("naive\t\t%.17g\n", log(p0/(1 - p0)));
+
+  fflush(stdout);
+  return ferror(stdout);
+}
diff --git a/src/test/test.c b/src/test/test.c
index a0a138b03..902565dfb 100644
--- a/src/test/test.c
+++ b/src/test/test.c
@@ -901,6 +901,7 @@ struct testgroup_t testgroups[] = {
   { "parsecommon/", parsecommon_tests },
   { "periodic-event/" , periodic_event_tests },
   { "policy/" , policy_tests },
+  { "prob_distr/", prob_distr_tests },
   { "procmon/", procmon_tests },
   { "process/", process_tests },
   { "proto/http/", proto_http_tests },
diff --git a/src/test/test.h b/src/test/test.h
index 9f6eb0a7e..39953e9f7 100644
--- a/src/test/test.h
+++ b/src/test/test.h
@@ -243,6 +243,8 @@ extern struct testcase_t parsecommon_tests[];
 extern struct testcase_t pem_tests[];
 extern struct testcase_t periodic_event_tests[];
 extern struct testcase_t policy_tests[];
+extern struct testcase_t prob_distr_tests[];
+extern struct testcase_t slow_stochastic_prob_distr_tests[];
 extern struct testcase_t procmon_tests[];
 extern struct testcase_t process_tests[];
 extern struct testcase_t proto_http_tests[];
diff --git a/src/test/test_prob_distr.c b/src/test/test_prob_distr.c
new file mode 100644
index 000000000..bf0f9e059
--- /dev/null
+++ b/src/test/test_prob_distr.c
@@ -0,0 +1,1414 @@
+/* Copyright (c) 2018, The Tor Project, Inc. */
+/* See LICENSE for licensing information */
+
+/**
+ * \file test_prob_distr.c
+ * \brief Test probability distributions.
+ * \detail
+ *
+ * For each probability distribution we do two kinds of tests:
+ *
+ * a) We do numerical deterministic testing of their cdf/icdf/sf/isf functions
+ *    and the various relationships between them for each distribution. We also
+ *    do deterministic tests on their sampling functions. Test vectors for
+ *    these tests were computed from alternative implementations and were
+ *    eyeballed to make sure they make sense
+ *    (e.g. src/test/prob_distr_mpfr_ref.c computes logit(p) using GNU mpfr
+ *    with 200-bit precision and is then tested in test_logit_logistic()).
+ *
+ * b) We do stochastic hypothesis testing (G-test) to ensure that sampling from
+ *    the given distributions is distributed properly. The stochastic tests are
+ *    slow and their false positive rate is not well suited for CI, so they are
+ *    currently disabled-by-default and put into 'tests-slow'.
+ */
+
+#define PROB_DISTR_PRIVATE
+
+#include "orconfig.h"
+
+#include "test/test.h"
+
+#include "core/or/or.h"
+
+#include "lib/math/prob_distr.h"
+#include "lib/math/fp.h"
+#include "lib/crypt_ops/crypto_rand.h"
+
+#include <float.h>
+#include <math.h>
+#include <stdbool.h>
+#include <stddef.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/**
+ * Return floor(d) converted to size_t, as a workaround for complaints
+ * under -Wbad-function-cast for (size_t)floor(d).
+ */
+static size_t
+floor_to_size_t(double d)
+{
+  double integral_d = floor(d);
+  return (size_t)integral_d;
+}
+
+/**
+ * Return ceil(d) converted to size_t, as a workaround for complaints
+ * under -Wbad-function-cast for (size_t)ceil(d).
+ */
+static size_t
+ceil_to_size_t(double d)
+{
+  double integral_d = ceil(d);
+  return (size_t)integral_d;
+}
+
+/*
+ * Geometric(p) distribution, supported on {1, 2, 3, ...}.
+ *
+ * Compute the probability mass function Geom(n; p) of the number of
+ * trials before the first success when success has probability p.
+ */
+static double
+logpmf_geometric(unsigned n, double p)
+{
+  /* This is actually a check against 1, but we do >= so that the compiler
+     does not raise a -Wfloat-equal */
+  if (p >= 1) {
+    if (n == 1)
+      return 0;
+    else
+      return -HUGE_VAL;
+  }
+  return (n - 1)*log1p(-p) + log(p);
+}
+
+/**
+ * Compute the logistic function, translated in output by 1/2:
+ * logistichalf(x) = logistic(x) - 1/2.  Well-conditioned on the entire
+ * real plane, with maximum condition number 1 at 0.
+ *
+ * This implementation gives relative error bounded by 5 eps.
+ */
+static double
+logistichalf(double x)
+{
+  /*
+   * Rewrite this with the identity
+   *
+   *  1/(1 + e^{-x}) - 1/2
+   *  = (1 - 1/2 - e^{-x}/2)/(1 + e^{-x})
+   *  = (1/2 - e^{-x}/2)/(1 + e^{-x})
+   *  = (1 - e^{-x})/[2 (1 + e^{-x})]
+   *  = -(e^{-x} - 1)/[2 (1 + e^{-x})],
+   *
+   * which we can evaluate by -expm1(-x)/[2 (1 + exp(-x))].
+   *
+   * Suppose exp has error d0, + has error d1, expm1 has error
+   * d2, and / has error d3, so we evaluate
+   *
+   *  -(1 + d2) (1 + d3) (e^{-x} - 1)
+   *    / [2 (1 + d1) (1 + (1 + d0) e^{-x})].
+   *
+   * In the denominator,
+   *
+   *  1 + (1 + d0) e^{-x}
+   *  = 1 + e^{-x} + d0 e^{-x}
+   *  = (1 + e^{-x}) (1 + d0 e^{-x}/(1 + e^{-x})),
+   *
+   * so the relative error of the numerator is
+   *
+   *  d' = d2 + d3 + d2 d3,
+   * and of the denominator,
+   *  d'' = d1 + d0 e^{-x}/(1 + e^{-x}) + d0 d1 e^{-x}/(1 + e^{-x})
+   *      = d1 + d0 L(-x) + d0 d1 L(-x),
+   *
+   * where L(-x) is logistic(-x).  By Lemma 1 the relative error
+   * of the quotient is bounded by
+   *
+   *  2|d2 + d3 + d2 d3 - d1 - d0 L(x) + d0 d1 L(x)|,
+   *
+   * Since 0 < L(x) < 1, this is bounded by
+   *
+   *  2|d2| + 2|d3| + 2|d2 d3| + 2|d1| + 2|d0| + 2|d0 d1|
+   *  <= 4 eps + 2 eps^2.
+   */
+  if (x < log(DBL_EPSILON/8)) {
+    /*
+     * Avoid overflow in e^{-x}.  When x < log(eps/4), we
+     * we further have x < logit(eps/4), so that
+     * logistic(x) < eps/4.  Hence the relative error of
+     * logistic(x) - 1/2 from -1/2 is bounded by eps/2, and
+     * so the relative error of -1/2 from logistic(x) - 1/2
+     * is bounded by eps.
+     */
+    return -0.5;
+  } else {
+    return -expm1(-x)/(2*(1 + exp(-x)));
+  }
+}
+
+/**
+ * Compute the log of the sum of the exps.  Caller should arrange the
+ * array in descending order to minimize error because I don't want to
+ * deal with using temporary space and the one caller in this file
+ * arranges that anyway.
+ *
+ * Warning: This implementation does not handle infinite or NaN inputs
+ * sensibly, because I don't need that here at the moment.  (NaN, or
+ * -inf and +inf together, should yield NaN; +inf and finite should
+ * yield +inf; otherwise all -inf should be ignored because exp(-inf) =
+ * 0.)
+ */
+static double
+logsumexp(double *A, size_t n)
+{
+  double maximum, sum;
+  size_t i;
+
+  if (n == 0)
+    return log(0);
+
+  maximum = A[0];
+  for (i = 1; i < n; i++) {
+    if (A[i] > maximum)
+      maximum = A[i];
+  }
+
+  sum = 0;
+  for (i = n; i --> 0;)
+    sum += exp(A[i] - maximum);
+
+  return log(sum) + maximum;
+}
+
+/**
+ * Compute log(1 - e^x).  Defined only for negative x so that e^x < 1.
+ * This is the complement of a probability in log space.
+ */
+static double
+log1mexp(double x)
+{
+
+  /*
+   * We want to compute log on [0, 1/2) but log1p on [1/2, +inf),
+   * so partition x at -log(2) = log(1/2).
+   */
+  if (-log(2) < x)
+    return log(-expm1(x));
+  else
+    return log1p(-exp(x));
+}
+
+/*
+ * Tests of numerical errors in computing logit, logistic, and the
+ * various cdfs, sfs, icdfs, and isfs.
+ */
+
+#define arraycount(A) (sizeof(A)/sizeof(A[0]))
+
+/** Return relative error between <b>actual</b> and <b>expected</b>.
+ *  Special cases: If <b>expected</b> is zero or infinite, return 1 if
+ *  <b>actual</b> is equal to <b>expected</b> and 0 if not, since the
+ *  usual notion of relative error is undefined but we only use this
+ *  for testing relerr(e, a) <= bound.  If either is NaN, return NaN,
+ *  which has the property that NaN <= bound is false no matter what
+ *  bound is.
+ *
+ *  Beware: if you test !(relerr(e, a) > bound), then then the result
+ *  is true when a is NaN because NaN > bound is false too.  See
+ *  CHECK_RELERR for correct use to decide when to report failure.
+ */
+static double
+relerr(double expected, double actual)
+{
+  /*
+   * To silence -Wfloat-equal, we have to test for equality using
+   * inequalities: we have (fabs(expected) <= 0) iff (expected == 0),
+   * and (actual <= expected && actual >= expected) iff actual ==
+   * expected whether expected is zero or infinite.
+   */
+  if (fabs(expected) <= 0 || tor_isinf(expected)) {
+    if (actual <= expected && actual >= expected)
+      return 0;
+    else
+      return 1;
+  } else {
+    return fabs((expected - actual)/expected);
+  }
+}
+
+/** Check that relative error of <b>expected</b> and <b>actual</b> is within
+ *  <b>relerr_bound</b>.  Caller must arrange to have i and relerr_bound in
+ *  scope.  */
+#define CHECK_RELERR(expected, actual) do {                                   \
+  double check_expected = (expected);                                         \
+  double check_actual = (actual);                                             \
+  const char *str_expected = #expected;                                       \
+  const char *str_actual = #actual;                                           \
+  double check_relerr = relerr(expected, actual);                             \
+  if (!(relerr(check_expected, check_actual) <= relerr_bound)) {              \
+    log_warn(LD_GENERAL, "%s:%d: case %u: relerr(%s=%.17e, %s=%.17e)"        \
+             " = %.17e > %.17e\n",                                            \
+             __func__, __LINE__, (unsigned) i,                                \
+             str_expected, check_expected,                                    \
+             str_actual, check_actual,                                        \
+             check_relerr, relerr_bound);                                     \
+    ok = false;                                                               \
+  }                                                                           \
+} while (0)
+
+/* Check that a <= b.
+ * Caller must arrange to have i in scope.  */
+#define CHECK_LE(a, b) do {                                                   \
+  double check_a = (a);                                                       \
+  double check_b = (b);                                                       \
+  const char *str_a = #a;                                                     \
+  const char *str_b = #b;                                                     \
+  if (!(check_a <= check_b)) {                                                \
+    log_warn(LD_GENERAL, "%s:%d: case %u: %s=%.17e > %s=%.17e\n",             \
+             __func__, __LINE__, (unsigned) i,                                \
+             str_a, check_a, str_b, check_b);                                 \
+    ok = false;                                                               \
+  }                                                                           \
+} while (0)
+
+/**
+ * Test the logit and logistic functions.  Confirm that they agree with
+ * the cdf, sf, icdf, and isf of the standard Logistic distribution.
+ * Confirm that the sampler for the standard logistic distribution maps
+ * [0, 1] into the right subinterval for the inverse transform, for
+ * this implementation.
+ */
+static void
+test_logit_logistic(void *arg)
+{
+  (void) arg;
+
+  static const struct {
+    double x;                   /* x = logit(p) */
+    double p;                   /* p = logistic(x) */
+    double phalf;               /* p - 1/2 = logistic(x) - 1/2 */
+  } cases[] = {
+    { -HUGE_VAL, 0, -0.5 },
+    { -1000, 0, -0.5 },
+    { -710, 4.47628622567513e-309, -0.5 },
+    { -708, 3.307553003638408e-308, -0.5 },
+    { -2, .11920292202211755, -.3807970779778824 },
+    { -1.0000001, .2689414017088022, -.23105859829119776 },
+    { -1, .2689414213699951, -.23105857863000487 },
+    { -0.9999999, .26894144103118883, -.2310585589688111 },
+    /* see src/test/prob_distr_mpfr_ref.c for computation */
+    { -4.000000000537333e-5, .49999, -1.0000000000010001e-5 },
+    { -4.000000000533334e-5, .49999, -.00001 },
+    { -4.000000108916878e-9, .499999999, -1.0000000272292198e-9 },
+    { -4e-9, .499999999, -1e-9 },
+    { -4e-16, .5, -1e-16 },
+    { -4e-300, .5, -1e-300 },
+    { 0, .5, 0 },
+    { 4e-300, .5, 1e-300 },
+    { 4e-16, .5, 1e-16 },
+    { 3.999999886872274e-9, .500000001, 9.999999717180685e-10 },
+    { 4e-9, .500000001, 1e-9 },
+    { 4.0000000005333336e-5, .50001, .00001 },
+    { 8.000042667076272e-3, .502, .002 },
+    { 0.9999999, .7310585589688111, .2310585589688111 },
+    { 1, .7310585786300049, .23105857863000487 },
+    { 1.0000001, .7310585982911977, .23105859829119774 },
+    { 2, .8807970779778823, .3807970779778824 },
+    { 708, 1, .5 },
+    { 710, 1, .5 },
+    { 1000, 1, .5 },
+    { HUGE_VAL, 1, .5 },
+  };
+  double relerr_bound = 3e-15; /* >10eps */
+  size_t i;
+  bool ok = true;
+
+  for (i = 0; i < arraycount(cases); i++) {
+    double x = cases[i].x;
+    double p = cases[i].p;
+    double phalf = cases[i].phalf;
+
+    /*
+     * cdf is logistic, icdf is logit, and symmetry for
+     * sf/isf.
+     */
+    CHECK_RELERR(logistic(x), cdf_logistic(x, 0, 1));
+    CHECK_RELERR(logistic(-x), sf_logistic(x, 0, 1));
+    CHECK_RELERR(logit(p), icdf_logistic(p, 0, 1));
+    CHECK_RELERR(-logit(p), isf_logistic(p, 0, 1));
+
+    CHECK_RELERR(cdf_logistic(x, 0, 1), cdf_logistic(x*2, 0, 2));
+    CHECK_RELERR(sf_logistic(x, 0, 1), sf_logistic(x*2, 0, 2));
+    CHECK_RELERR(icdf_logistic(p, 0, 1), icdf_logistic(p, 0, 2)/2);
+    CHECK_RELERR(isf_logistic(p, 0, 1), isf_logistic(p, 0, 2)/2);
+
+    CHECK_RELERR(cdf_logistic(x, 0, 1), cdf_logistic(x/2, 0, .5));
+    CHECK_RELERR(sf_logistic(x, 0, 1), sf_logistic(x/2, 0, .5));
+    CHECK_RELERR(icdf_logistic(p, 0, 1), icdf_logistic(p, 0,.5)*2);
+    CHECK_RELERR(isf_logistic(p, 0, 1), isf_logistic(p, 0, .5)*2);
+
+    CHECK_RELERR(cdf_logistic(x, 0, 1), cdf_logistic(x*2 + 1, 1, 2));
+    CHECK_RELERR(sf_logistic(x, 0, 1), sf_logistic(x*2 + 1, 1, 2));
+
+    /*
+     * For p near 0 and p near 1/2, the arithmetic of
+     * translating by 1 loses precision.
+     */
+    if (fabs(p) > DBL_EPSILON && fabs(p) < 0.4) {
+      CHECK_RELERR(icdf_logistic(p, 0, 1),
+          (icdf_logistic(p, 1, 2) - 1)/2);
+      CHECK_RELERR(isf_logistic(p, 0, 1),
+          (isf_logistic(p, 1, 2) - 1)/2);
+    }
+
+    CHECK_RELERR(p, logistic(x));
+    CHECK_RELERR(phalf, logistichalf(x));
+
+    /*
+     * On the interior floating-point numbers, either logit or
+     * logithalf had better give the correct answer.
+     *
+     * For probabilities near 0, we can get much finer resolution with
+     * logit, and for probabilities near 1/2, we can get much finer
+     * resolution with logithalf by representing them using p - 1/2.
+     *
+     * E.g., we can write -.00001 for phalf, and .49999 for p, but the
+     * difference 1/2 - .00001 gives 1.0000000000010001e-5 in binary64
+     * arithmetic.  So test logit(.49999) which should give the same
+     * answer as logithalf(-1.0000000000010001e-5), namely
+     * -4.000000000537333e-5, and also test logithalf(-.00001) which
+     * gives -4.000000000533334e-5 instead -- but don't expect
+     * logit(.49999) to give -4.000000000533334e-5 even though it looks
+     * like 1/2 - .00001.
+     *
+     * A naive implementation of logit will just use log(p/(1 - p)) and
+     * give the answer -4.000000000551673e-05 for .49999, which is
+     * wrong in a lot of digits, which happens because log is
+     * ill-conditioned near 1 and thus amplifies whatever relative
+     * error we made in computing p/(1 - p).
+     */
+    if ((0 < p && p < 1) || tor_isinf(x)) {
+      if (phalf >= p - 0.5 && phalf <= p - 0.5)
+        CHECK_RELERR(x, logit(p));
+      if (p >= 0.5 + phalf && p <= 0.5 + phalf)
+        CHECK_RELERR(x, logithalf(phalf));
+    }
+
+    CHECK_RELERR(-phalf, logistichalf(-x));
+    if (fabs(phalf) < 0.5 || tor_isinf(x))
+      CHECK_RELERR(-x, logithalf(-phalf));
+    if (p < 1 || tor_isinf(x)) {
+      CHECK_RELERR(1 - p, logistic(-x));
+      if (p > .75 || tor_isinf(x))
+        CHECK_RELERR(-x, logit(1 - p));
+    } else {
+      CHECK_LE(logistic(-x), 1e-300);
+    }
+  }
+
+  for (i = 0; i <= 100; i++) {
+    double p0 = (double)i/100;
+
+    CHECK_RELERR(logit(p0/(1 + M_E)), sample_logistic(0, 0, p0));
+    CHECK_RELERR(-logit(p0/(1 + M_E)), sample_logistic(1, 0, p0));
+    CHECK_RELERR(logithalf(p0*(0.5 - 1/(1 + M_E))),
+        sample_logistic(0, 1, p0));
+    CHECK_RELERR(-logithalf(p0*(0.5 - 1/(1 + M_E))),
+        sample_logistic(1, 1, p0));
+  }
+
+  if (!ok)
+    printf("fail logit/logistic / logistic cdf/sf\n");
+
+  tt_assert(ok);
+
+ done:
+  ;
+}
+
+/**
+ * Test the cdf, sf, icdf, and isf of the LogLogistic distribution.
+ */
+static void
+test_log_logistic(void *arg)
+{
+  (void) arg;
+
+  static const struct {
+    /* x is a point in the support of the LogLogistic distribution */
+    double x;
+    /* 'p' is the probability that a random variable X for a given LogLogistic
+     * probability ditribution will take value less-or-equal to x */
+    double p;
+    /* 'np' is the probability that a random variable X for a given LogLogistic
+     * probability distribution will take value greater-or-equal to x. */
+    double np;
+  } cases[] = {
+    { 0, 0, 1 },
+    { 1e-300, 1e-300, 1 },
+    { 1e-17, 1e-17, 1 },
+    { 1e-15, 1e-15, .999999999999999 },
+    { .1, .09090909090909091, .90909090909090909 },
+    { .25, .2, .8 },
+    { .5, .33333333333333333, .66666666666666667 },
+    { .75, .42857142857142855, .5714285714285714 },
+    { .9999, .49997499874993756, .5000250012500626 },
+    { .99999999, .49999999749999996, .5000000025 },
+    { .999999999999999, .49999999999999994, .5000000000000002 },
+    { 1, .5, .5 },
+  };
+  double relerr_bound = 3e-15;
+  size_t i;
+  bool ok = true;
+
+  for (i = 0; i < arraycount(cases); i++) {
+    double x = cases[i].x;
+    double p = cases[i].p;
+    double np = cases[i].np;
+
+    CHECK_RELERR(p, cdf_log_logistic(x, 1, 1));
+    CHECK_RELERR(p, cdf_log_logistic(x/2, .5, 1));
+    CHECK_RELERR(p, cdf_log_logistic(x*2, 2, 1));
+    CHECK_RELERR(p, cdf_log_logistic(sqrt(x), 1, 2));
+    CHECK_RELERR(p, cdf_log_logistic(sqrt(x)/2, .5, 2));
+    CHECK_RELERR(p, cdf_log_logistic(sqrt(x)*2, 2, 2));
+    if (2*sqrt(DBL_MIN) < x) {
+      CHECK_RELERR(p, cdf_log_logistic(x*x, 1, .5));
+      CHECK_RELERR(p, cdf_log_logistic(x*x/2, .5, .5));
+      CHECK_RELERR(p, cdf_log_logistic(x*x*2, 2, .5));
+    }
+
+    CHECK_RELERR(np, sf_log_logistic(x, 1, 1));
+    CHECK_RELERR(np, sf_log_logistic(x/2, .5, 1));
+    CHECK_RELERR(np, sf_log_logistic(x*2, 2, 1));
+    CHECK_RELERR(np, sf_log_logistic(sqrt(x), 1, 2));
+    CHECK_RELERR(np, sf_log_logistic(sqrt(x)/2, .5, 2));
+    CHECK_RELERR(np, sf_log_logistic(sqrt(x)*2, 2, 2));
+    if (2*sqrt(DBL_MIN) < x) {
+      CHECK_RELERR(np, sf_log_logistic(x*x, 1, .5));
+      CHECK_RELERR(np, sf_log_logistic(x*x/2, .5, .5));
+      CHECK_RELERR(np, sf_log_logistic(x*x*2, 2, .5));
+    }
+
+    CHECK_RELERR(np, cdf_log_logistic(1/x, 1, 1));
+    CHECK_RELERR(np, cdf_log_logistic(1/(2*x), .5, 1));
+    CHECK_RELERR(np, cdf_log_logistic(2/x, 2, 1));
+    CHECK_RELERR(np, cdf_log_logistic(1/sqrt(x), 1, 2));
+    CHECK_RELERR(np, cdf_log_logistic(1/(2*sqrt(x)), .5, 2));
+    CHECK_RELERR(np, cdf_log_logistic(2/sqrt(x), 2, 2));
+    if (2*sqrt(DBL_MIN) < x && x < 1/(2*sqrt(DBL_MIN))) {
+      CHECK_RELERR(np, cdf_log_logistic(1/(x*x), 1, .5));
+      CHECK_RELERR(np, cdf_log_logistic(1/(2*x*x), .5, .5));
+      CHECK_RELERR(np, cdf_log_logistic(2/(x*x), 2, .5));
+    }
+
+    CHECK_RELERR(p, sf_log_logistic(1/x, 1, 1));
+    CHECK_RELERR(p, sf_log_logistic(1/(2*x), .5, 1));
+    CHECK_RELERR(p, sf_log_logistic(2/x, 2, 1));
+    CHECK_RELERR(p, sf_log_logistic(1/sqrt(x), 1, 2));
+    CHECK_RELERR(p, sf_log_logistic(1/(2*sqrt(x)), .5, 2));
+    CHECK_RELERR(p, sf_log_logistic(2/sqrt(x), 2, 2));
+    if (2*sqrt(DBL_MIN) < x && x < 1/(2*sqrt(DBL_MIN))) {
+      CHECK_RELERR(p, sf_log_logistic(1/(x*x), 1, .5));
+      CHECK_RELERR(p, sf_log_logistic(1/(2*x*x), .5, .5));
+      CHECK_RELERR(p, sf_log_logistic(2/(x*x), 2, .5));
+    }
+
+    CHECK_RELERR(x, icdf_log_logistic(p, 1, 1));
+    CHECK_RELERR(x/2, icdf_log_logistic(p, .5, 1));
+    CHECK_RELERR(x*2, icdf_log_logistic(p, 2, 1));
+    CHECK_RELERR(x, icdf_log_logistic(p, 1, 1));
+    CHECK_RELERR(sqrt(x)/2, icdf_log_logistic(p, .5, 2));
+    CHECK_RELERR(sqrt(x)*2, icdf_log_logistic(p, 2, 2));
+    CHECK_RELERR(sqrt(x), icdf_log_logistic(p, 1, 2));
+    CHECK_RELERR(x*x/2, icdf_log_logistic(p, .5, .5));
+    CHECK_RELERR(x*x*2, icdf_log_logistic(p, 2, .5));
+
+    if (np < .9) {
+      CHECK_RELERR(x, isf_log_logistic(np, 1, 1));
+      CHECK_RELERR(x/2, isf_log_logistic(np, .5, 1));
+      CHECK_RELERR(x*2, isf_log_logistic(np, 2, 1));
+      CHECK_RELERR(sqrt(x), isf_log_logistic(np, 1, 2));
+      CHECK_RELERR(sqrt(x)/2, isf_log_logistic(np, .5, 2));
+      CHECK_RELERR(sqrt(x)*2, isf_log_logistic(np, 2, 2));
+      CHECK_RELERR(x*x, isf_log_logistic(np, 1, .5));
+      CHECK_RELERR(x*x/2, isf_log_logistic(np, .5, .5));
+      CHECK_RELERR(x*x*2, isf_log_logistic(np, 2, .5));
+
+      CHECK_RELERR(1/x, icdf_log_logistic(np, 1, 1));
+      CHECK_RELERR(1/(2*x), icdf_log_logistic(np, .5, 1));
+      CHECK_RELERR(2/x, icdf_log_logistic(np, 2, 1));
+      CHECK_RELERR(1/sqrt(x), icdf_log_logistic(np, 1, 2));
+      CHECK_RELERR(1/(2*sqrt(x)),
+          icdf_log_logistic(np, .5, 2));
+      CHECK_RELERR(2/sqrt(x), icdf_log_logistic(np, 2, 2));
+      CHECK_RELERR(1/(x*x), icdf_log_logistic(np, 1, .5));
+      CHECK_RELERR(1/(2*x*x), icdf_log_logistic(np, .5, .5));
+      CHECK_RELERR(2/(x*x), icdf_log_logistic(np, 2, .5));
+    }
+
+    CHECK_RELERR(1/x, isf_log_logistic(p, 1, 1));
+    CHECK_RELERR(1/(2*x), isf_log_logistic(p, .5, 1));
+    CHECK_RELERR(2/x, isf_log_logistic(p, 2, 1));
+    CHECK_RELERR(1/sqrt(x), isf_log_logistic(p, 1, 2));
+    CHECK_RELERR(1/(2*sqrt(x)), isf_log_logistic(p, .5, 2));
+    CHECK_RELERR(2/sqrt(x), isf_log_logistic(p, 2, 2));
+    CHECK_RELERR(1/(x*x), isf_log_logistic(p, 1, .5));
+    CHECK_RELERR(1/(2*x*x), isf_log_logistic(p, .5, .5));
+    CHECK_RELERR(2/(x*x), isf_log_logistic(p, 2, .5));
+  }
+
+  for (i = 0; i <= 100; i++) {
+    double p0 = (double)i/100;
+
+    CHECK_RELERR(0.5*p0/(1 - 0.5*p0), sample_log_logistic(0, p0));
+    CHECK_RELERR((1 - 0.5*p0)/(0.5*p0),
+        sample_log_logistic(1, p0));
+  }
+
+  if (!ok)
+    printf("fail log logistic cdf/sf\n");
+
+  tt_assert(ok);
+
+ done:
+  ;
+}
+
+/**
+ * Test the cdf, sf, icdf, isf of the Weibull distribution.
+ */
+static void
+test_weibull(void *arg)
+{
+  (void) arg;
+
+  static const struct {
+    /* x is a point in the support of the Weibull distribution */
+    double x;
+    /* 'p' is the probability that a random variable X for a given Weibull
+     * probability ditribution will take value less-or-equal to x */
+    double p;
+    /* 'np' is the probability that a random variable X for a given Weibull
+     * probability distribution will take value greater-or-equal to x. */
+    double np;
+  } cases[] = {
+    { 0, 0, 1 },
+    { 1e-300, 1e-300, 1 },
+    { 1e-17, 1e-17, 1 },
+    { .1, .09516258196404043, .9048374180359595 },
+    { .5, .3934693402873666, .6065306597126334 },
+    { .6931471805599453, .5, .5 },
+    { 1, .6321205588285577, .36787944117144233 },
+    { 10, .9999546000702375, 4.5399929762484854e-5 },
+    { 36, .9999999999999998, 2.319522830243569e-16 },
+    { 37, .9999999999999999, 8.533047625744066e-17 },
+    { 38, 1, 3.1391327920480296e-17 },
+    { 100, 1, 3.720075976020836e-44 },
+    { 708, 1, 3.307553003638408e-308 },
+    { 710, 1, 4.47628622567513e-309 },
+    { 1000, 1, 0 },
+    { HUGE_VAL, 1, 0 },
+  };
+  double relerr_bound = 3e-15;
+  size_t i;
+  bool ok = true;
+
+  for (i = 0; i < arraycount(cases); i++) {
+    double x = cases[i].x;
+    double p = cases[i].p;
+    double np = cases[i].np;
+
+    CHECK_RELERR(p, cdf_weibull(x, 1, 1));
+    CHECK_RELERR(p, cdf_weibull(x/2, .5, 1));
+    CHECK_RELERR(p, cdf_weibull(x*2, 2, 1));
+    /* For 0 < x < sqrt(DBL_MIN), x^2 loses lots of bits.  */
+    if (x <= 0 ||
+        sqrt(DBL_MIN) <= x) {
+      CHECK_RELERR(p, cdf_weibull(x*x, 1, .5));
+      CHECK_RELERR(p, cdf_weibull(x*x/2, .5, .5));
+      CHECK_RELERR(p, cdf_weibull(x*x*2, 2, .5));
+    }
+    CHECK_RELERR(p, cdf_weibull(sqrt(x), 1, 2));
+    CHECK_RELERR(p, cdf_weibull(sqrt(x)/2, .5, 2));
+    CHECK_RELERR(p, cdf_weibull(sqrt(x)*2, 2, 2));
+    CHECK_RELERR(np, sf_weibull(x, 1, 1));
+    CHECK_RELERR(np, sf_weibull(x/2, .5, 1));
+    CHECK_RELERR(np, sf_weibull(x*2, 2, 1));
+    CHECK_RELERR(np, sf_weibull(x*x, 1, .5));
+    CHECK_RELERR(np, sf_weibull(x*x/2, .5, .5));
+    CHECK_RELERR(np, sf_weibull(x*x*2, 2, .5));
+    if (x >= 10) {
+      /*
+       * exp amplifies the error of sqrt(x)^2
+       * proportionally to exp(x); for large inputs
+       * this is significant.
+       */
+      double t = -expm1(-x*(2*DBL_EPSILON + DBL_EPSILON));
+      relerr_bound = t + DBL_EPSILON + t*DBL_EPSILON;
+      if (relerr_bound < 3e-15)
+        /*
+         * The tests are written only to 16
+         * decimal places anyway even if your
+         * `double' is, say, i387 binary80, for
+         * whatever reason.
+         */
+        relerr_bound = 3e-15;
+      CHECK_RELERR(np, sf_weibull(sqrt(x), 1, 2));
+      CHECK_RELERR(np, sf_weibull(sqrt(x)/2, .5, 2));
+      CHECK_RELERR(np, sf_weibull(sqrt(x)*2, 2, 2));
+    }
+
+    if (p <= 0.75) {
+      /*
+       * For p near 1, not enough precision near 1 to
+       * recover x.
+       */
+      CHECK_RELERR(x, icdf_weibull(p, 1, 1));
+      CHECK_RELERR(x/2, icdf_weibull(p, .5, 1));
+      CHECK_RELERR(x*2, icdf_weibull(p, 2, 1));
+    }
+    if (p >= 0.25 && !tor_isinf(x) && np > 0) {
+      /*
+       * For p near 0, not enough precision in np
+       * near 1 to recover x.  For 0, isf gives inf,
+       * even if p is precise enough for the icdf to
+       * work.
+       */
+      CHECK_RELERR(x, isf_weibull(np, 1, 1));
+      CHECK_RELERR(x/2, isf_weibull(np, .5, 1));
+      CHECK_RELERR(x*2, isf_weibull(np, 2, 1));
+    }
+  }
+
+  for (i = 0; i <= 100; i++) {
+    double p0 = (double)i/100;
+
+    CHECK_RELERR(3*sqrt(-log(p0/2)), sample_weibull(0, p0, 3, 2));
+    CHECK_RELERR(3*sqrt(-log1p(-p0/2)),
+        sample_weibull(1, p0, 3, 2));
+  }
+
+  if (!ok)
+    printf("fail Weibull cdf/sf\n");
+
+  tt_assert(ok);
+
+ done:
+  ;
+}
+
+/**
+ * Test the cdf, sf, icdf, and isf of the generalized Pareto
+ * distribution.
+ */
+static void
+test_genpareto(void *arg)
+{
+  (void) arg;
+
+  struct {
+    /* xi is the 'xi' parameter of the generalized Pareto distribution, and the
+     * rest are the same as in the above tests */
+    double xi, x, p, np;
+  } cases[] = {
+    { 0, 0, 0, 1 },
+    { 1e-300, .004, 3.992010656008528e-3, .9960079893439915 },
+    { 1e-300, .1, .09516258196404043, .9048374180359595 },
+    { 1e-300, 1, .6321205588285577, .36787944117144233 },
+    { 1e-300, 10, .9999546000702375, 4.5399929762484854e-5 },
+    { 1e-200, 1e-16, 9.999999999999999e-17, .9999999999999999 },
+    { 1e-16, 1e-200, 9.999999999999998e-201, 1 },
+    { 1e-16, 1e-16, 1e-16, 1 },
+    { 1e-16, .004, 3.992010656008528e-3, .9960079893439915 },
+    { 1e-16, .1, .09516258196404043, .9048374180359595 },
+    { 1e-16, 1, .6321205588285577, .36787944117144233 },
+    { 1e-16, 10, .9999546000702375, 4.539992976248509e-5 },
+    { 1e-10, 1e-6, 9.999995000001667e-7, .9999990000005 },
+    { 1e-8, 1e-8, 9.999999950000001e-9, .9999999900000001 },
+    { 1, 1e-300, 1e-300, 1 },
+    { 1, 1e-16, 1e-16, .9999999999999999 },
+    { 1, .1, .09090909090909091, .9090909090909091 },
+    { 1, 1, .5, .5 },
+    { 1, 10, .9090909090909091, .0909090909090909 },
+    { 1, 100, .9900990099009901, .0099009900990099 },
+    { 1, 1000, .999000999000999, 9.990009990009992e-4 },
+    { 10, 1e-300, 1e-300, 1 },
+    { 10, 1e-16, 9.999999999999995e-17, .9999999999999999 },
+    { 10, .1, .06696700846319258, .9330329915368074 },
+    { 10, 1, .21320655780322778, .7867934421967723 },
+    { 10, 10, .3696701667040189, .6303298332959811 },
+    { 10, 100, .49886285755007337, .5011371424499267 },
+    { 10, 1000, .6018968102992647, .3981031897007353 },
+  };
+  double xi_array[] = { -1.5, -1, -1e-30, 0, 1e-30, 1, 1.5 };
+  size_t i, j;
+  double relerr_bound = 3e-15;
+  bool ok = true;
+
+  for (i = 0; i < arraycount(cases); i++) {
+    double xi = cases[i].xi;
+    double x = cases[i].x;
+    double p = cases[i].p;
+    double np = cases[i].np;
+
+    CHECK_RELERR(p, cdf_genpareto(x, 0, 1, xi));
+    CHECK_RELERR(p, cdf_genpareto(x*2, 0, 2, xi));
+    CHECK_RELERR(p, cdf_genpareto(x/2, 0, .5, xi));
+    CHECK_RELERR(np, sf_genpareto(x, 0, 1, xi));
+    CHECK_RELERR(np, sf_genpareto(x*2, 0, 2, xi));
+    CHECK_RELERR(np, sf_genpareto(x/2, 0, .5, xi));
+
+    if (p < .5) {
+      CHECK_RELERR(x, icdf_genpareto(p, 0, 1, xi));
+      CHECK_RELERR(x*2, icdf_genpareto(p, 0, 2, xi));
+      CHECK_RELERR(x/2, icdf_genpareto(p, 0, .5, xi));
+    }
+    if (np < .5) {
+      CHECK_RELERR(x, isf_genpareto(np, 0, 1, xi));
+      CHECK_RELERR(x*2, isf_genpareto(np, 0, 2, xi));
+      CHECK_RELERR(x/2, isf_genpareto(np, 0, .5, xi));
+    }
+  }
+
+  for (i = 0; i < arraycount(xi_array); i++) {
+    for (j = 0; j <= 100; j++) {
+      double p0 = (j == 0 ? 2*DBL_MIN : (double)j/100);
+
+      /* This is actually a check against 0, but we do <= so that the compiler
+         does not raise a -Wfloat-equal */
+      if (fabs(xi_array[i]) <= 0) {
+        /*
+         * When xi == 0, the generalized Pareto
+         * distribution reduces to an
+         * exponential distribution.
+         */
+        CHECK_RELERR(-log(p0/2),
+            sample_genpareto(0, p0, 0));
+        CHECK_RELERR(-log1p(-p0/2),
+            sample_genpareto(1, p0, 0));
+      } else {
+        CHECK_RELERR(expm1(-xi_array[i]*log(p0/2))/xi_array[i],
+            sample_genpareto(0, p0, xi_array[i]));
+        CHECK_RELERR((j == 0 ? DBL_MIN :
+                expm1(-xi_array[i]*log1p(-p0/2))/xi_array[i]),
+            sample_genpareto(1, p0, xi_array[i]));
+      }
+
+      CHECK_RELERR(isf_genpareto(p0/2, 0, 1, xi_array[i]),
+          sample_genpareto(0, p0, xi_array[i]));
+      CHECK_RELERR(icdf_genpareto(p0/2, 0, 1, xi_array[i]),
+          sample_genpareto(1, p0, xi_array[i]));
+    }
+  }
+
+  tt_assert(ok);
+
+ done:
+  ;
+}
+
+/**
+ * Test the deterministic sampler for uniform distribution on [a, b].
+ *
+ * This currently only tests whether the outcome lies within [a, b].
+ */
+static void
+test_uniform_interval(void *arg)
+{
+  (void) arg;
+  struct {
+    /* Sample from a uniform distribution with parameters 'a' and 'b', using
+     * 't' as the sampling index. */
+    double t, a, b;
+  } cases[] = {
+    { 0, 0, 0 },
+    { 0, 0, 1 },
+    { 0, 1.0000000000000007, 3.999999999999995 },
+    { 0, 4000, 4000 },
+    { 0.42475836677491291, 4000, 4000 },
+    { 0, -DBL_MAX, DBL_MAX },
+    { 0.25, -DBL_MAX, DBL_MAX },
+    { 0.5, -DBL_MAX, DBL_MAX },
+  };
+  size_t i = 0;
+  bool ok = true;
+
+  for (i = 0; i < arraycount(cases); i++) {
+    double t = cases[i].t;
+    double a = cases[i].a;
+    double b = cases[i].b;
+
+    CHECK_LE(a, sample_uniform_interval(t, a, b));
+    CHECK_LE(sample_uniform_interval(t, a, b), b);
+
+    CHECK_LE(a, sample_uniform_interval(1 - t, a, b));
+    CHECK_LE(sample_uniform_interval(1 - t, a, b), b);
+
+    CHECK_LE(sample_uniform_interval(t, -b, -a), -a);
+    CHECK_LE(-b, sample_uniform_interval(t, -b, -a));
+
+    CHECK_LE(sample_uniform_interval(1 - t, -b, -a), -a);
+    CHECK_LE(-b, sample_uniform_interval(1 - t, -b, -a));
+  }
+
+  tt_assert(ok);
+
+ done:
+  ;
+}
+
+/********************** Stochastic tests ****************************/
+
+/*
+ * Psi test, sometimes also called G-test.  The psi test statistic,
+ * suitably scaled, has chi^2 distribution, but the psi test tends to
+ * have better statistical power in practice to detect deviations than
+ * the chi^2 test does.  (The chi^2 test statistic is the first term of
+ * the Taylor expansion of the psi test statistic.)  The psi test is
+ * generic, for any CDF; particular distributions might have higher-
+ * power tests to distinguish them from predictable deviations or bugs.
+ *
+ * We choose the psi critical value so that a single psi test has
+ * probability below alpha = 1% of spuriously failing even if all the
+ * code is correct.  But the false positive rate for a suite of n tests
+ * is higher: 1 - Binom(0; n, alpha) = 1 - (1 - alpha)^n.  For n = 10,
+ * this is about 10%, and for n = 100 it is well over 50%.
+ *
+ * We can drive it down by running each test twice, and accepting it if
+ * it passes at least once; in that case, it is as if we used Binom(2;
+ * 2, alpha) = alpha^2 as the false positive rate for each test, and
+ * for n = 10 tests, it would be 0.1%, and for n = 100 tests, still
+ * only 1%.
+ *
+ * The critical value for a chi^2 distribution with 100 degrees of
+ * freedom and false positive rate alpha = 1% was taken from:
+ *
+ *  NIST/SEMATECH e-Handbook of Statistical Methods, Section
+ *  1.3.6.7.4 `Critical Values of the Chi-Square Distribution',
+ *  <http://www.itl.nist.gov/div898/handbook/eda/section3/eda3674.htm>,
+ *  retrieved 2018-10-28.
+ */
+
+static const size_t NSAMPLES = 100000;
+/* Number of chances we give to the test to succeed. */
+static const unsigned NTRIALS = 2;
+/* Number of times we want the test to pass per NTRIALS. */
+static const unsigned NPASSES_MIN = 1;
+
+#define PSI_DF 100                          /* degrees of freedom */
+static const double PSI_CRITICAL = 135.807; /* critical value, alpha = .01 */
+
+/**
+ * Perform a psi test on an array of sample counts, C, adding up to N
+ * samples, and an array of log expected probabilities, logP,
+ * representing the null hypothesis for the distribution of samples
+ * counted.  Return false if the psi test rejects the null hypothesis,
+ * true if otherwise.
+ */
+static bool
+psi_test(const size_t C[PSI_DF], const double logP[PSI_DF], size_t N)
+{
+  double psi = 0;
+  double c = 0;                 /* Kahan compensation */
+  double t, u;
+  size_t i;
+
+  for (i = 0; i < PSI_DF; i++) {
+    /*
+     * c*log(c/(n*p)) = (1/n) * f*log(f/p) where f = c/n is
+     * the frequency, and f*log(f/p) ---> 0 as f ---> 0, so
+     * this is a reasonable choice.  Further, any mass that
+     * _fails_ to turn up in this bin will inflate another
+     * bin instead, so we don't really lose anything by
+     * ignoring empty bins even if they have high
+     * probability.
+     */
+    if (C[i] == 0)
+      continue;
+    t = C[i]*(log((double)C[i]/N) - logP[i]) - c;
+    u = psi + t;
+    c = (u - psi) - t;
+    psi = u;
+  }
+  psi *= 2;
+
+  return psi <= PSI_CRITICAL;
+}
+
+static bool
+test_stochastic_geometric_impl(double p)
+{
+  double logP[PSI_DF] = {0};
+  unsigned ntry = NTRIALS, npass = 0;
+  unsigned i;
+  size_t j;
+
+  /* Compute logP[i] = Geom(i + 1; p).  */
+  for (i = 0; i < PSI_DF - 1; i++)
+    logP[i] = logpmf_geometric(i + 1, p);
+
+  /* Compute logP[n-1] = log (1 - (P[0] + P[1] + ... + P[n-2])).  */
+  logP[PSI_DF - 1] = log1mexp(logsumexp(logP, PSI_DF - 1));
+
+  while (ntry --> 0) {
+    size_t C[PSI_DF] = {0};
+
+    for (j = 0; j < NSAMPLES; j++) {
+      double n_tmp = ceil(geometric_sample(p));
+      unsigned n = (unsigned) n_tmp;
+
+      if (n > PSI_DF)
+        n = PSI_DF;
+      C[n - 1]++;
+    }
+
+    if (psi_test(C, logP, NSAMPLES)) {
+      if (++npass >= NPASSES_MIN)
+        break;
+    }
+  }
+
+  if (npass >= NPASSES_MIN) {
+    /* printf("pass %s sampler\n", "geometric"); */
+    return true;
+  } else {
+    printf("fail %s sampler\n", "geometric");
+    return false;
+  }
+}
+
+/**
+ * Divide the support of <b>dist</b> into histogram bins in <b>logP</b>. Start
+ * at the 1st percentile and ending at the 99th percentile. Pick the bin
+ * boundaries using linear interpolation so that they are uniformly spaced.
+ *
+ * In each bin logP[i] we insert the expected log-probability that a sampled
+ * value will fall into that bin. We will use this as the null hypothesis of
+ * the psi test.
+ *
+ * Set logP[i] = log(CDF(x_i) - CDF(x_{i-1})), where x_-1 = -inf, x_n =
+ * +inf, and x_i = i*(hi - lo)/(n - 2).
+ */
+static void
+bin_cdfs(const struct dist *dist, double lo, double hi, double *logP, size_t n)
+{
+#define CDF(x)  dist->ops->cdf(dist, x)
+#define SF(x)   dist->ops->sf(dist, x)
+  const double w = (hi - lo)/(n - 2);
+  double halfway = dist->ops->icdf(dist, 0.5);
+  double x_0, x_1;
+  size_t i;
+  size_t n2 = ceil_to_size_t((halfway - lo)/w);
+
+  tor_assert(lo <= halfway);
+  tor_assert(halfway <= hi);
+  tor_assert(n2 <= n);
+
+  x_1 = lo;
+  logP[0] = log(CDF(x_1) - 0); /* 0 = CDF(-inf) */
+  for (i = 1; i < n2; i++) {
+    x_0 = x_1;
+    /* do the linear interpolation */
+    x_1 = (i <= n/2 ? lo + i*w : hi - (n - 2 - i)*w);
+    /* set the expected log-probability */
+    logP[i] = log(CDF(x_1) - CDF(x_0));
+  }
+  x_0 = hi;
+  logP[n - 1] = log(SF(x_0) - 0); /* 0 = SF(+inf) = 1 - CDF(+inf) */
+
+  /* In this loop we are filling out the high part of the array. We are using
+   * SF because in these cases the CDF is near 1 where precision is lower. So
+   * instead we are using SF near 0 where the precision is higher. We have
+   * SF(t) = 1 - CDF(t).  */
+  for (i = 1; i < n - n2; i++) {
+    x_1 = x_0;
+    /* do the linear interpolation */
+    x_0 = (i <= n/2 ? hi - i*w : lo + (n - 2 - i)*w);
+    /* set the expected log-probability */
+    logP[n - i - 1] = log(SF(x_0) - SF(x_1));
+  }
+#undef SF
+#undef CDF
+}
+
+/**
+ * Draw NSAMPLES samples from dist, counting the number of samples x in
+ * the ith bin C[i] if x_{i-1} <= x < x_i, where x_-1 = -inf, x_n =
+ * +inf, and x_i = i*(hi - lo)/(n - 2).
+ */
+static void
+bin_samples(const struct dist *dist, double lo, double hi, size_t *C, size_t n)
+{
+  const double w = (hi - lo)/(n - 2);
+  size_t i;
+
+  for (i = 0; i < NSAMPLES; i++) {
+    double x = dist->ops->sample(dist);
+    size_t bin;
+
+    if (x < lo)
+      bin = 0;
+    else if (x < hi)
+      bin = 1 + floor_to_size_t((x - lo)/w);
+    else
+      bin = n - 1;
+    tor_assert(bin < n);
+    C[bin]++;
+  }
+}
+
+/**
+ * Carry out a Psi test on <b>dist</b>.
+ *
+ * Sample NSAMPLES from dist, putting them in bins from -inf to lo to
+ * hi to +inf, and apply up to two psi tests.  True if at least one psi
+ * test passes; false if not.  False positive rate should be bounded by
+ * 0.01^2 = 0.0001.
+ */
+static bool
+test_psi_dist_sample(const struct dist *dist)
+{
+  double logP[PSI_DF] = {0};
+  unsigned ntry = NTRIALS, npass = 0;
+  double lo = dist->ops->icdf(dist, 1/(double)(PSI_DF + 2));
+  double hi = dist->ops->isf(dist, 1/(double)(PSI_DF + 2));
+
+  /* Create the null hypothesis in logP */
+  bin_cdfs(dist, lo, hi, logP, PSI_DF);
+
+  /* Now run the test */
+  while (ntry --> 0) {
+    size_t C[PSI_DF] = {0};
+    bin_samples(dist, lo, hi, C, PSI_DF);
+    if (psi_test(C, logP, NSAMPLES)) {
+      if (++npass >= NPASSES_MIN)
+        break;
+    }
+  }
+
+  /* Did we fail or succeed? */
+  if (npass >= NPASSES_MIN) {
+    /* printf("pass %s sampler\n", dist->ops->name);*/
+    return true;
+  } else {
+    printf("fail %s sampler\n", dist->ops->name);
+    return false;
+  }
+}
+
+/* This is the seed of the deterministic randomness */
+static uint32_t deterministic_rand_counter;
+
+/** Initialize the seed of the deterministic randomness. */
+static void
+init_deterministic_rand(void)
+{
+  deterministic_rand_counter = crypto_rand_uint32();
+}
+
+/** Produce deterministic randomness for the stochastic tests using the global
+ *  deterministic_rand_counter seed
+ *
+ *  This function produces deterministic data over multiple calls iff it's
+ *  called in the same call order with the same 'n' parameter (which is the
+ *  case for the psi test). If not, outputs will deviate. */
+static void
+crypto_rand_deterministic(char *out, size_t n)
+{
+  /* Use a XOF to squeeze bytes out of that silly counter */
+  crypto_xof_t *xof = crypto_xof_new();
+  tor_assert(xof);
+  crypto_xof_add_bytes(xof, (uint8_t*)&deterministic_rand_counter,
+                       sizeof(deterministic_rand_counter));
+  crypto_xof_squeeze_bytes(xof, (uint8_t*)out, n);
+  crypto_xof_free(xof);
+
+  /* Increase counter for next run */
+  deterministic_rand_counter++;
+}
+
+static void
+test_stochastic_uniform(void *arg)
+{
+  (void) arg;
+
+  const struct uniform uniform01 = {
+    .base = DIST_BASE(&uniform_ops),
+    .a = 0,
+    .b = 1,
+  };
+  const struct uniform uniform_pos = {
+    .base = DIST_BASE(&uniform_ops),
+    .a = 1.23,
+    .b = 4.56,
+  };
+  const struct uniform uniform_neg = {
+    .base = DIST_BASE(&uniform_ops),
+    .a = -10,
+    .b = -1,
+  };
+  const struct uniform uniform_cross = {
+    .base = DIST_BASE(&uniform_ops),
+    .a = -1.23,
+    .b = 4.56,
+  };
+  const struct uniform uniform_subnormal = {
+    .base = DIST_BASE(&uniform_ops),
+    .a = 4e-324,
+    .b = 4e-310,
+  };
+  const struct uniform uniform_subnormal_cross = {
+    .base = DIST_BASE(&uniform_ops),
+    .a = -4e-324,
+    .b = 4e-310,
+  };
+  bool ok = true;
+
+  init_deterministic_rand();
+  MOCK(crypto_rand, crypto_rand_deterministic);
+
+  ok &= test_psi_dist_sample(&uniform01.base);
+  ok &= test_psi_dist_sample(&uniform_pos.base);
+  ok &= test_psi_dist_sample(&uniform_neg.base);
+  ok &= test_psi_dist_sample(&uniform_cross.base);
+  ok &= test_psi_dist_sample(&uniform_subnormal.base);
+  ok &= test_psi_dist_sample(&uniform_subnormal_cross.base);
+
+  tt_assert(ok);
+
+ done:
+    ;
+}
+
+static bool
+test_stochastic_logistic_impl(double mu, double sigma)
+{
+  const struct logistic dist = {
+    .base = DIST_BASE(&logistic_ops),
+    .mu = mu,
+    .sigma = sigma,
+  };
+
+  /* XXX Consider some fancier logistic test.  */
+  return test_psi_dist_sample(&dist.base);
+}
+
+static bool
+test_stochastic_log_logistic_impl(double alpha, double beta)
+{
+  const struct log_logistic dist = {
+    .base = DIST_BASE(&log_logistic_ops),
+    .alpha = alpha,
+    .beta = beta,
+  };
+
+  /* XXX Consider some fancier log logistic test.  */
+  return test_psi_dist_sample(&dist.base);
+}
+
+static bool
+test_stochastic_weibull_impl(double lambda, double k)
+{
+  const struct weibull dist = {
+    .base = DIST_BASE(&weibull_ops),
+    .lambda = lambda,
+    .k = k,
+  };
+
+/*
+ * XXX Consider applying a Tiku-Singh test:
+ *
+ *    M.L. Tiku and M. Singh, `Testing the two-parameter
+ *    Weibull distribution', Communications in Statistics --
+ *    Theory and Methods A10(9), 1981, 907--918.
+ *https://www.tandfonline.com/doi/pdf/10.1080/03610928108828082?needAccess=true
+ */
+  return test_psi_dist_sample(&dist.base);
+}
+
+static bool
+test_stochastic_genpareto_impl(double mu, double sigma, double xi)
+{
+  const struct genpareto dist = {
+    .base = DIST_BASE(&genpareto_ops),
+    .mu = mu,
+    .sigma = sigma,
+    .xi = xi,
+  };
+
+  /* XXX Consider some fancier GPD test.  */
+  return test_psi_dist_sample(&dist.base);
+}
+
+static void
+test_stochastic_genpareto(void *arg)
+{
+  bool ok = 0;
+  bool tests_failed = true;
+  (void) arg;
+
+  init_deterministic_rand();
+  MOCK(crypto_rand, crypto_rand_deterministic);
+
+  ok = test_stochastic_genpareto_impl(0, 1, -0.25);
+  tt_assert(ok);
+  ok = test_stochastic_genpareto_impl(0, 1, -1e-30);
+  tt_assert(ok);
+  ok = test_stochastic_genpareto_impl(0, 1, 0);
+  tt_assert(ok);
+  ok = test_stochastic_genpareto_impl(0, 1, 1e-30);
+  tt_assert(ok);
+  ok = test_stochastic_genpareto_impl(0, 1, 0.25);
+  tt_assert(ok);
+  ok = test_stochastic_genpareto_impl(-1, 1, -0.25);
+  tt_assert(ok);
+  ok = test_stochastic_genpareto_impl(1, 2, 0.25);
+  tt_assert(ok);
+
+  tests_failed = false;
+
+ done:
+  if (tests_failed) {
+    printf("seed: %"PRIu32, deterministic_rand_counter);
+  }
+  UNMOCK(crypto_rand);
+}
+
+static void
+test_stochastic_geometric(void *arg)
+{
+  bool ok = 0;
+  bool tests_failed = true;
+
+  (void) arg;
+
+  init_deterministic_rand();
+  MOCK(crypto_rand, crypto_rand_deterministic);
+
+  ok = test_stochastic_geometric_impl(0.1);
+  tt_assert(ok);
+  ok = test_stochastic_geometric_impl(0.5);
+  tt_assert(ok);
+  ok = test_stochastic_geometric_impl(0.9);
+  tt_assert(ok);
+  ok = test_stochastic_geometric_impl(1);
+  tt_assert(ok);
+
+  tests_failed = false;
+
+ done:
+  if (tests_failed) {
+    printf("seed: %"PRIu32, deterministic_rand_counter);
+  }
+  UNMOCK(crypto_rand);
+}
+
+static void
+test_stochastic_logistic(void *arg)
+{
+  bool ok = 0;
+  bool tests_failed = true;
+  (void) arg;
+
+  init_deterministic_rand();
+  MOCK(crypto_rand, crypto_rand_deterministic);
+
+  ok = test_stochastic_logistic_impl(0, 1);
+  tt_assert(ok);
+  ok = test_stochastic_logistic_impl(0, 1e-16);
+  tt_assert(ok);
+  ok = test_stochastic_logistic_impl(1, 10);
+  tt_assert(ok);
+  ok = test_stochastic_logistic_impl(-10, 100);
+  tt_assert(ok);
+
+  tests_failed = false;
+
+ done:
+  if (tests_failed) {
+    printf("seed: %"PRIu32, deterministic_rand_counter);
+  }
+  UNMOCK(crypto_rand);
+}
+
+static void
+test_stochastic_log_logistic(void *arg)
+{
+  bool ok = 0;
+  bool tests_failed = true;
+  (void) arg;
+
+  init_deterministic_rand();
+  MOCK(crypto_rand, crypto_rand_deterministic);
+
+  ok = test_stochastic_log_logistic_impl(1, 1);
+  tt_assert(ok);
+  ok = test_stochastic_log_logistic_impl(1, 10);
+  tt_assert(ok);
+  ok = test_stochastic_log_logistic_impl(M_E, 1e-1);
+  tt_assert(ok);
+  ok = test_stochastic_log_logistic_impl(exp(-10), 1e-2);
+  tt_assert(ok);
+
+  tests_failed = false;
+
+ done:
+  if (tests_failed) {
+    printf("seed: %"PRIu32, deterministic_rand_counter);
+  }
+  UNMOCK(crypto_rand);
+}
+
+static void
+test_stochastic_weibull(void *arg)
+{
+  bool ok = 0;
+  bool tests_failed = true;
+  (void) arg;
+
+  init_deterministic_rand();
+  MOCK(crypto_rand, crypto_rand_deterministic);
+
+  ok = test_stochastic_weibull_impl(1, 0.5);
+  tt_assert(ok);
+  ok = test_stochastic_weibull_impl(1, 1);
+  tt_assert(ok);
+  ok = test_stochastic_weibull_impl(1, 1.5);
+  tt_assert(ok);
+  ok = test_stochastic_weibull_impl(1, 2);
+  tt_assert(ok);
+  ok = test_stochastic_weibull_impl(10, 1);
+  tt_assert(ok);
+
+  tests_failed = false;
+
+ done:
+  if (tests_failed) {
+    printf("seed: %"PRIu32, deterministic_rand_counter);
+  }
+  UNMOCK(crypto_rand);
+}
+
+struct testcase_t prob_distr_tests[] = {
+  { "logit_logistics", test_logit_logistic, TT_FORK, NULL, NULL },
+  { "log_logistic", test_log_logistic, TT_FORK, NULL, NULL },
+  { "weibull", test_weibull, TT_FORK, NULL, NULL },
+  { "genpareto", test_genpareto, TT_FORK, NULL, NULL },
+  { "uniform_interval", test_uniform_interval, TT_FORK, NULL, NULL },
+  END_OF_TESTCASES
+};
+
+struct testcase_t slow_stochastic_prob_distr_tests[] = {
+  { "stochastic_genpareto", test_stochastic_genpareto, TT_FORK, NULL, NULL },
+  { "stochastic_geometric", test_stochastic_geometric, TT_FORK, NULL, NULL },
+  { "stochastic_uniform", test_stochastic_uniform, TT_FORK, NULL, NULL },
+  { "stochastic_logistic", test_stochastic_logistic, TT_FORK, NULL, NULL },
+  { "stochastic_log_logistic", test_stochastic_log_logistic, TT_FORK, NULL,
+    NULL },
+  { "stochastic_weibull", test_stochastic_weibull, TT_FORK, NULL, NULL },
+  END_OF_TESTCASES
+};
diff --git a/src/test/test_slow.c b/src/test/test_slow.c
index 97c2912af..39a203c72 100644
--- a/src/test/test_slow.c
+++ b/src/test/test_slow.c
@@ -21,6 +21,7 @@
 struct testgroup_t testgroups[] = {
   { "slow/crypto/", slow_crypto_tests },
   { "slow/process/", slow_process_tests },
+  { "slow/prob_distr/", slow_stochastic_prob_distr_tests },
   END_OF_GROUPS
 };
 



_______________________________________________
tor-commits mailing list
tor-commits@xxxxxxxxxxxxxxxxxxxx
https://lists.torproject.org/cgi-bin/mailman/listinfo/tor-commits