diff --git a/changes/laplace-edge-cases b/changes/laplace-edge-cases
new file mode 100644
index 0000000000..f8a36e1de3
--- /dev/null
+++ b/changes/laplace-edge-cases
@@ -0,0 +1,7 @@
+ o Minor bugfixes:
+ - Handle edge cases in the laplace functions:
+ * avoid division by zero
+ * avoid taking the log of zero
+ * silence clang type conversion warnings using round and trunc
+ - Add tests for laplace edge cases.
+ Related to HS stats in #13192.
diff --git a/src/common/util.c b/src/common/util.c
index f8d1b7be4b..52b3e04946 100644
--- a/src/common/util.c
+++ b/src/common/util.c
@@ -536,13 +536,18 @@ int64_t
sample_laplace_distribution(double mu, double b, double p)
{
double result;
-
tor_assert(p >= 0.0 && p < 1.0);
+
/* This is the "inverse cumulative distribution function" from:
* http://en.wikipedia.org/wiki/Laplace_distribution */
+ if (p == 0.0) {
+ /* Avoid taking log(0.0) == -INFINITY, as some processors or compiler
+ * options can cause the program to trap. */
+ return INT64_MIN;
+ }
+
result = mu - b * (p > 0.5 ? 1.0 : -1.0)
* tor_mathlog(1.0 - 2.0 * fabs(p - 0.5));
-
if (result >= INT64_MAX)
return INT64_MAX;
else if (result <= INT64_MIN)
@@ -551,18 +556,28 @@ sample_laplace_distribution(double mu, double b, double p)
return (int64_t) result;
}
-/** Add random noise between INT64_MIN and INT64_MAX coming from a
- * Laplace distribution with mu = 0 and b = delta_f/epsilon
- * to signal based on the provided random value in
- * [0.0, 1.0[. */
+/** Add random noise between INT64_MIN and INT64_MAX coming from a Laplace
+ * distribution with mu = 0 and b = delta_f/epsilon to
+ * signal based on the provided random value in [0.0, 1.0[.
+ * The epislon value must be between ]0.0, 1.0]. delta_f must be greater
+ * than 0. */
int64_t
add_laplace_noise(int64_t signal, double random, double delta_f,
double epsilon)
{
- int64_t noise = sample_laplace_distribution(
- 0.0, /* just add noise, no further signal */
- delta_f / epsilon, random);
+ int64_t noise;
+ /* epsilon MUST be between ]0.0, 1.0] */
+ tor_assert(epsilon > 0.0 && epsilon <= 1.0);
+ /* delta_f MUST be greater than 0. */
+ tor_assert(delta_f > 0.0);
+
+ /* Just add noise, no further signal */
+ noise = sample_laplace_distribution(0.0,
+ delta_f / epsilon,
+ random);
+
+ /* Clip (signal + noise) to [INT64_MIN, INT64_MAX] */
if (noise > 0 && INT64_MAX - noise < signal)
return INT64_MAX;
else if (noise < 0 && INT64_MIN - noise > signal)
diff --git a/src/test/test_util.c b/src/test/test_util.c
index 6a4c3ec072..940f9bde52 100644
--- a/src/test/test_util.c
+++ b/src/test/test_util.c
@@ -4100,12 +4100,201 @@ test_util_laplace(void *arg)
*/
tt_i64_op(INT64_MIN + 20, ==,
add_laplace_noise(20, 0.0, delta_f, epsilon));
- tt_i64_op(-60, ==, add_laplace_noise(20, 0.1, delta_f, epsilon));
- tt_i64_op(-14, ==, add_laplace_noise(20, 0.25, delta_f, epsilon));
- tt_i64_op(20, ==, add_laplace_noise(20, 0.5, delta_f, epsilon));
- tt_i64_op(54, ==, add_laplace_noise(20, 0.75, delta_f, epsilon));
- tt_i64_op(100, ==, add_laplace_noise(20, 0.9, delta_f, epsilon));
- tt_i64_op(215, ==, add_laplace_noise(20, 0.99, delta_f, epsilon));
+
+ tt_assert(-60 == add_laplace_noise(20, 0.1, delta_f, epsilon));
+ tt_assert(-14 == add_laplace_noise(20, 0.25, delta_f, epsilon));
+ tt_assert(20 == add_laplace_noise(20, 0.5, delta_f, epsilon));
+ tt_assert(54 == add_laplace_noise(20, 0.75, delta_f, epsilon));
+ tt_assert(100 == add_laplace_noise(20, 0.9, delta_f, epsilon));
+ tt_assert(215 == add_laplace_noise(20, 0.99, delta_f, epsilon));
+
+ /* Test extreme values of signal with maximally negative values of noise
+ * 1.0000000000000002 is the smallest number > 1
+ * 0.0000000000000002 is the double epsilon (error when calculating near 1)
+ * this is approximately 1/(2^52)
+ * per https://en.wikipedia.org/wiki/Double_precision
+ * (let's not descend into the world of subnormals)
+ * >>> laplace.ppf([0, 0.0000000000000002], loc = 0, scale = 1)
+ * array([ -inf, -35.45506713])
+ */
+ const double noscale_df = 1.0, noscale_eps = 1.0;
+
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0, 0.0, noscale_df, noscale_eps));
+
+ /* is it clipped to INT64_MIN? */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(-1, 0.0, noscale_df, noscale_eps));
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(INT64_MIN, 0.0,
+ noscale_df, noscale_eps));
+ /* ... even when scaled? */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0, 0.0, delta_f, epsilon));
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0, 0.0,
+ INT64_MAX, 1));
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(INT64_MIN, 0.0,
+ INT64_MAX, 1));
+
+ /* does it play nice with INT64_MAX? */
+ tt_assert((INT64_MIN + INT64_MAX) ==
+ add_laplace_noise(INT64_MAX, 0.0,
+ noscale_df, noscale_eps));
+
+ /* do near-zero fractional values work? */
+ const double min_dbl_error = 0.0000000000000002;
+
+ tt_assert(-35 ==
+ add_laplace_noise(0, min_dbl_error,
+ noscale_df, noscale_eps));
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(INT64_MIN, min_dbl_error,
+ noscale_df, noscale_eps));
+ tt_assert((-35 + INT64_MAX) ==
+ add_laplace_noise(INT64_MAX, min_dbl_error,
+ noscale_df, noscale_eps));
+ /* ... even when scaled? */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0, min_dbl_error,
+ INT64_MIN, -35));
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0, min_dbl_error,
+ INT64_MIN, -34));
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(0, min_dbl_error,
+ INT64_MIN, 1));
+ tt_assert((INT64_MIN + INT64_MAX) ==
+ add_laplace_noise(INT64_MIN, min_dbl_error,
+ INT64_MIN, 1));
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(INT64_MAX, min_dbl_error,
+ INT64_MIN, 1));
+
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(0, min_dbl_error,
+ INT64_MAX, -35));
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(0, min_dbl_error,
+ INT64_MAX, -34));
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0, min_dbl_error,
+ INT64_MAX, 1));
+ tt_assert((INT64_MAX + INT64_MIN) ==
+ add_laplace_noise(INT64_MAX, min_dbl_error,
+ INT64_MAX, 1));
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(INT64_MIN, min_dbl_error,
+ INT64_MAX, 1));
+
+ /* does it play nice with INT64_MAX? */
+ tt_assert((INT64_MAX - 35) ==
+ add_laplace_noise(INT64_MAX, min_dbl_error,
+ noscale_df, noscale_eps));
+
+ /* Test extreme values of signal with maximally positive values of noise
+ * 1.0000000000000002 is the smallest number > 1
+ * 0.9999999999999998 is the greatest number < 1 by calculation
+ * per https://en.wikipedia.org/wiki/Double_precision
+ * >>> laplace.ppf([1.0, 0.9999999999999998], loc = 0, scale = 1)
+ * array([inf, 35.35050621])
+ * but the function rejects p == 1.0, so we just use max_dbl_lt_one
+ */
+ const double max_dbl_lt_one = 0.9999999999999998;
+
+ /* do near-one fractional values work? */
+ tt_assert(35 ==
+ add_laplace_noise(0, max_dbl_lt_one, noscale_df, noscale_eps));
+
+ /* is it clipped to INT64_MAX? */
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(INT64_MAX - 35, max_dbl_lt_one,
+ noscale_df, noscale_eps));
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(INT64_MAX - 34, max_dbl_lt_one,
+ noscale_df, noscale_eps));
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(INT64_MAX, max_dbl_lt_one,
+ noscale_df, noscale_eps));
+ /* ... even when scaled? */
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(INT64_MAX, max_dbl_lt_one,
+ delta_f, epsilon));
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(0, max_dbl_lt_one,
+ INT64_MAX, 35));
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(0, max_dbl_lt_one,
+ INT64_MAX, 34));
+ tt_assert((INT64_MIN + INT64_MAX) ==
+ add_laplace_noise(INT64_MIN, max_dbl_lt_one,
+ INT64_MAX, 1));
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(INT64_MAX, max_dbl_lt_one,
+ INT64_MAX, 1));
+ tt_assert((INT64_MAX + INT64_MIN) ==
+ add_laplace_noise(INT64_MAX, max_dbl_lt_one,
+ INT64_MIN, 1));
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(INT64_MIN, max_dbl_lt_one,
+ INT64_MIN, 1));
+
+ /* does it play nice with INT64_MIN? */
+ tt_assert((INT64_MIN + 35) ==
+ add_laplace_noise(INT64_MIN, max_dbl_lt_one,
+ noscale_df, noscale_eps));
+
+ /* Test extreme values of b = delta_f / epsilon
+ * >>> laplace.ppf([0.5], loc = 0, scale = 1)
+ * array([0. ])
+ */
+
+ /* Make sure edge cases don't depend on architecture,
+ * optimisation level, or other compiler flags.
+ * XXXX Are these edge cases important enough to make consistent? */
+
+ /* b = positive zero, p yields positive zero */
+ tt_assert(0.0 ==
+ add_laplace_noise(0.0, 0.5, 0.0, 1.0))
+ /* b = negative zero, p yields positive zero */
+ tt_assert(0.0 ==
+ add_laplace_noise(0.0, 0.5, 0.0, -1.0))
+ /* b = positive infinity, p yields positive zero, result is -NaN -> -Inf */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0.0, 0.5, 1.0, 0.0))
+ /* b = negative infinity, p yields positive zero, result is -NaN -> -Inf */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0.0, 0.5, -1.0, 0.0))
+ /* b = positive NaN (rounded to -Inf), p yields positive zero,
+ * result is -NaN -> -Inf */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0.0, 0.5, -0.0, -0.0))
+ /* b = negative NaN (rounded to -Inf), p yields positive zero,
+ * result is -NaN -> -Inf*/
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0.0, 0.5, -0.0, 0.0))
+
+ /* b = positive zero, p yields negative infinity, result is -NaN -> -Inf */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0.0, 0.0, 0.0, 1.0))
+ /* b = negative zero, p yields negative infinity, result is -NaN -> -Inf */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0.0, 0.0, 0.0, -1.0))
+ /* b = positive infinity, p yields negative infinity */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0.0, 0.0, 1.0, 0.0))
+ /* b = negative infinity, p yields negative infinity */
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(0.0, 0.0, -1.0, 0.0))
+ /* b = positive NaN (rounded to -Inf), p yields negative infinity,
+ * result is -NaN -> -Inf */
+ tt_assert(INT64_MIN ==
+ add_laplace_noise(0.0, 0.0, -0.0, -0.0))
+ /* b = negative NaN (rounded to -Inf), p yields negative infinity,
+ * result is NaN -> Inf */
+ tt_assert(INT64_MAX ==
+ add_laplace_noise(0.0, 0.0, -0.0, 0.0))
done:
;