src/lib/math/laplace.c

   1 /* Copyright (c) 2003, Roger Dingledine
   2  * Copyright (c) 2004-2006, Roger Dingledine, Nick Mathewson.
   3  * Copyright (c) 2007-2021, The Tor Project, Inc. */
   4 /* See LICENSE for licensing information */
   5
   6 /**
   7  * \file laplace.c
   8  *
   9  * \brief Implements a Laplace distribution, used for adding noise to things.
  10  **/
  11
  12 #include "orconfig.h"
  13 #include "lib/math/laplace.h"
  14 #include "lib/math/fp.h"
  15
  16 #include "lib/log/util_bug.h"
  17
  18 #include <math.h>
  19 #include <stdlib.h>
  20
  21 /** Transform a random value <b>p</b> from the uniform distribution in
  22  * [0.0, 1.0[ into a Laplace distributed value with location parameter
  23  * <b>mu</b> and scale parameter <b>b</b>. Truncate the final result
  24  * to be an integer in [INT64_MIN, INT64_MAX]. */
  25 int64_t
  26 sample_laplace_distribution(double mu, double b, double p)
  27 {
  28   double result;
  29   tor_assert(p >= 0.0 && p < 1.0);
  30
  31   /* This is the "inverse cumulative distribution function" from:
  32    * https://en.wikipedia.org/wiki/Laplace_distribution */
  33   if (p <= 0.0) {
  34     /* Avoid taking log(0.0) == -INFINITY, as some processors or compiler
  35      * options can cause the program to trap. */
  36     return INT64_MIN;
  37   }
  38
  39   result = mu - b * (p > 0.5 ? 1.0 : -1.0)
  40                   * tor_mathlog(1.0 - 2.0 * fabs(p - 0.5));
  41
  42   return clamp_double_to_int64(result);
  43 }
  44
  45 /** Add random noise between INT64_MIN and INT64_MAX coming from a Laplace
  46  * distribution with mu = 0 and b = <b>delta_f</b>/<b>epsilon</b> to
  47  * <b>signal</b> based on the provided <b>random</b> value in [0.0, 1.0[.
  48  * The epsilon value must be between ]0.0, 1.0]. delta_f must be greater
  49  * than 0. */
  50 int64_t
  51 add_laplace_noise(int64_t signal_, double random_, double delta_f,
  52                   double epsilon)
  53 {
  54   int64_t noise;
  55
  56   /* epsilon MUST be between ]0.0, 1.0] */
  57   tor_assert(epsilon > 0.0 && epsilon <= 1.0);
  58   /* delta_f MUST be greater than 0. */
  59   tor_assert(delta_f > 0.0);
  60
  61   /* Just add noise, no further signal */
  62   noise = sample_laplace_distribution(0.0,
  63                                       delta_f / epsilon,
  64                                       random_);
  65
  66   /* Clip (signal + noise) to [INT64_MIN, INT64_MAX] */
  67   if (noise > 0 && INT64_MAX - noise < signal_)
  68     return INT64_MAX;
  69   else if (noise < 0 && INT64_MIN - noise > signal_)
  70     return INT64_MIN;
  71   else
  72     return signal_ + noise;
  73 }