laplace.c (2304B)
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 }