[tor-commits] [tor/master] Implement and test probability distributions used by WTF-PAD.
nickm at torproject.org
nickm at torproject.org
Mon Jan 14 19:50:35 UTC 2019
commit 2ccf3268375cd46e8c948e94ba58e0d2f03fe722
Author: George Kadianakis <desnacked at riseup.net>
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 at mumble.net>
Co-authored-by: Mike Perry <mikeperry-git at torproject.org>
---
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
};
More information about the tor-commits
mailing list