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@mumble.net>
Co-authored-by: Mike Perry <mikeperry-git@torproject.org>

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 \

src/lib/crypt_ops/crypto_rand.c

@@ -528,6 +528,17 @@ crypto_rand_unmocked(char *to, size_t n)
 #endif
 }
 
+/**
+ * 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

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;

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

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
+}
+

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

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

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

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

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);
+}

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 },

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[];

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
 };