Merge branch 'asn-karsten-task-13192-5-squashed'

changes/feature13192

@@ -0,0 +1,13 @@
+m  o Major features (hidden services):
+    - Add a HiddenServiceStatistics option that allows Tor relays to
+      gather and publish statistics about hidden service usage, to
+      better understand the size and volume of the hidden service
+      network. Specifically, if a Tor relay is an HSDir it will
+      publish the approximate number of hidden services that have
+      published descriptors to it the past 24 hours. Also, if a relay
+      has acted as a hidden service rendezvous point, it will publish
+      the approximate amount of rendezvous cells it has relayed the
+      past 24 hours. The statistics themselves are obfuscated so that
+      the exact values cannot be derived. For more details see
+      proposal 238 "Better hidden service stats from Tor relays". This
+      feature is currently disabled by default. Implements feature 13192.

doc/tor.1.txt

@@ -1764,6 +1764,13 @@ is non-zero):
     When this option is enabled, Tor writes statistics on the bidirectional use
     of connections to disk every 24 hours. (Default: 0)
 
+[[HiddenServiceStatistics]] **HiddenServiceStatistics** **0**|**1**::
+    When this option is enabled, a Tor relay writes obfuscated
+    statistics on its role as hidden-service directory, introduction
+    point, or rendezvous point to disk every 24 hours. If
+    ExtraInfoStatistics is also enabled, these statistics are further
+    published to the directory authorities. (Default: 0)
+
 [[ExtraInfoStatistics]] **ExtraInfoStatistics** **0**|**1**::
     When this option is enabled, Tor includes previously gathered statistics in
     its extra-info documents that it uploads to the directory authorities.

src/common/crypto.c

@@ -1293,7 +1293,7 @@ crypto_pk_asn1_decode(const char *str, size_t len)
  * Return 0 on success, -1 on failure.
  */
 int
-crypto_pk_get_digest(crypto_pk_t *pk, char *digest_out)
+crypto_pk_get_digest(const crypto_pk_t *pk, char *digest_out)
 {
   unsigned char *buf = NULL;
   int len;

src/common/crypto.h

@@ -180,7 +180,7 @@ int crypto_pk_private_hybrid_decrypt(crypto_pk_t *env, char *to,
 
 int crypto_pk_asn1_encode(crypto_pk_t *pk, char *dest, size_t dest_len);
 crypto_pk_t *crypto_pk_asn1_decode(const char *str, size_t len);
-int crypto_pk_get_digest(crypto_pk_t *pk, char *digest_out);
+int crypto_pk_get_digest(const crypto_pk_t *pk, char *digest_out);
 int crypto_pk_get_all_digests(crypto_pk_t *pk, digests_t *digests_out);
 int crypto_pk_get_fingerprint(crypto_pk_t *pk, char *fp_out,int add_space);
 int crypto_pk_get_hashed_fingerprint(crypto_pk_t *pk, char *fp_out);

src/common/util.c

@@ -513,6 +513,51 @@ round_uint64_to_next_multiple_of(uint64_t number, uint64_t divisor)
   return number;
 }
 
+/** Return the lowest x in [INT64_MIN, INT64_MAX] such that x is at least
+ * <b>number</b>, and x modulo <b>divisor</b> == 0. */
+int64_t
+round_int64_to_next_multiple_of(int64_t number, int64_t divisor)
+{
+  tor_assert(divisor > 0);
+  if (number >= 0 && INT64_MAX - divisor + 1 >= number)
+    number += divisor - 1;
+  number -= number % divisor;
+  return number;
+}
+
+/** Transform a random value <b>p</b> from the uniform distribution in
+ * [0.0, 1.0[ into a Laplace distributed value with location parameter
+ * <b>mu</b> and scale parameter <b>b</b> in [-Inf, Inf[. */
+double
+sample_laplace_distribution(double mu, double b, double p)
+{
+  tor_assert(p >= 0.0 && p < 1.0);
+  /* This is the "inverse cumulative distribution function" from:
+   * http://en.wikipedia.org/wiki/Laplace_distribution */
+  return mu - b * (p > 0.5 ? 1.0 : -1.0)
+                * tor_mathlog(1.0 - 2.0 * fabs(p - 0.5));
+}
+
+/** Add random noise between INT64_MIN and INT64_MAX coming from a
+ * Laplace distribution with mu = 0 and b = <b>delta_f</b>/<b>epsilon</b>
+ * to <b>signal</b> based on the provided <b>random</b> value in
+ * [0.0, 1.0[. */
+int64_t
+add_laplace_noise(int64_t signal, double random, double delta_f,
+                  double epsilon)
+{
+  /* cast to int64_t intended */
+  int64_t noise = sample_laplace_distribution(
+               0.0, /* just add noise, no further signal */
+               delta_f / epsilon, random);
+  if (noise > 0 && INT64_MAX - noise < signal)
+    return INT64_MAX;
+  else if (noise < 0 && INT64_MIN - noise > signal)
+    return INT64_MIN;
+  else
+    return signal + noise;
+}
+
 /** Return the number of bits set in <b>v</b>. */
 int
 n_bits_set_u8(uint8_t v)

src/common/util.h

@@ -172,6 +172,10 @@ uint64_t round_to_power_of_2(uint64_t u64);
 unsigned round_to_next_multiple_of(unsigned number, unsigned divisor);
 uint32_t round_uint32_to_next_multiple_of(uint32_t number, uint32_t divisor);
 uint64_t round_uint64_to_next_multiple_of(uint64_t number, uint64_t divisor);
+int64_t round_int64_to_next_multiple_of(int64_t number, int64_t divisor);
+double sample_laplace_distribution(double mu, double b, double p);
+int64_t add_laplace_noise(int64_t signal, double random, double delta_f,
+                          double epsilon);
 int n_bits_set_u8(uint8_t v);
 
 /* Compute the CEIL of <b>a</b> divided by <b>b</b>, for nonnegative <b>a</b>

src/or/command.c

@@ -438,6 +438,7 @@ command_process_created_cell(cell_t *cell, channel_t *chan)
 static void
 command_process_relay_cell(cell_t *cell, channel_t *chan)
 {
+  const or_options_t *options = get_options();
   circuit_t *circ;
   int reason, direction;
 
@@ -511,6 +512,14 @@ command_process_relay_cell(cell_t *cell, channel_t *chan)
            direction==CELL_DIRECTION_OUT?"forward":"backward");
     circuit_mark_for_close(circ, -reason);
   }
+
+  /* If this is a cell in an RP circuit, count it as part of the
+     hidden service stats */
+  if (options->HiddenServiceStatistics &&
+      !CIRCUIT_IS_ORIGIN(circ) &&
+      TO_OR_CIRCUIT(circ)->circuit_carries_hs_traffic_stats) {
+    rep_hist_seen_new_rp_cell();
+  }
 }
 
 /** Process a 'destroy' <b>cell</b> that just arrived from

src/or/config.c

@@ -268,6 +268,7 @@ static config_var_t option_vars_[] = {
   VAR("HiddenServicePort",   LINELIST_S, RendConfigLines,    NULL),
   VAR("HiddenServiceVersion",LINELIST_S, RendConfigLines,    NULL),
   VAR("HiddenServiceAuthorizeClient",LINELIST_S,RendConfigLines, NULL),
+  V(HiddenServiceStatistics,     BOOL,     "0"),
   V(HidServAuth,                 LINELIST, NULL),
   V(CloseHSClientCircuitsImmediatelyOnTimeout, BOOL, "0"),
   V(CloseHSServiceRendCircuitsImmediatelyOnTimeout, BOOL, "0"),
@@ -1714,6 +1715,7 @@ options_act(const or_options_t *old_options)
   if (options->CellStatistics || options->DirReqStatistics ||
       options->EntryStatistics || options->ExitPortStatistics ||
       options->ConnDirectionStatistics ||
+      options->HiddenServiceStatistics ||
       options->BridgeAuthoritativeDir) {
     time_t now = time(NULL);
     int print_notice = 0;
@@ -1722,6 +1724,7 @@ options_act(const or_options_t *old_options)
     if (!public_server_mode(options)) {
       options->CellStatistics = 0;
       options->EntryStatistics = 0;
+      options->HiddenServiceStatistics = 0;
       options->ExitPortStatistics = 0;
     }
 
@@ -1767,6 +1770,11 @@ options_act(const or_options_t *old_options)
         options->ConnDirectionStatistics) {
       rep_hist_conn_stats_init(now);
     }
+    if ((!old_options || !old_options->HiddenServiceStatistics) &&
+        options->HiddenServiceStatistics) {
+      log_info(LD_CONFIG, "Configured to measure hidden service statistics.");
+      rep_hist_hs_stats_init(now);
+    }
     if ((!old_options || !old_options->BridgeAuthoritativeDir) &&
         options->BridgeAuthoritativeDir) {
       rep_hist_desc_stats_init(now);
@@ -1778,6 +1786,8 @@ options_act(const or_options_t *old_options)
                  "data directory in 24 hours from now.");
   }
 
+  /* If we used to have statistics enabled but we just disabled them,
+     stop gathering them.  */
   if (old_options && old_options->CellStatistics &&
       !options->CellStatistics)
     rep_hist_buffer_stats_term();
@@ -1787,6 +1797,9 @@ options_act(const or_options_t *old_options)
   if (old_options && old_options->EntryStatistics &&
       !options->EntryStatistics)
     geoip_entry_stats_term();
+  if (old_options && old_options->HiddenServiceStatistics &&
+      !options->HiddenServiceStatistics)
+    rep_hist_hs_stats_term();
   if (old_options && old_options->ExitPortStatistics &&
       !options->ExitPortStatistics)
     rep_hist_exit_stats_term();

src/or/main.c

@@ -1384,6 +1384,11 @@ run_scheduled_events(time_t now)
       if (next_write && next_write < next_time_to_write_stats_files)
         next_time_to_write_stats_files = next_write;
     }
+    if (options->HiddenServiceStatistics) {
+      time_t next_write = rep_hist_hs_stats_write(time_to_write_stats_files);
+      if (next_write && next_write < next_time_to_write_stats_files)
+        next_time_to_write_stats_files = next_write;
+    }
     if (options->ExitPortStatistics) {
       time_t next_write = rep_hist_exit_stats_write(time_to_write_stats_files);
       if (next_write && next_write < next_time_to_write_stats_files)

src/or/or.h

@@ -3204,6 +3204,10 @@ typedef struct or_circuit_t {
   /** True iff this circuit was made with a CREATE_FAST cell. */
   unsigned int is_first_hop : 1;
 
+  /** If set, this circuit carries HS traffic. Consider it in any HS
+   *  statistics. */
+  unsigned int circuit_carries_hs_traffic_stats : 1;
+
   /** Number of cells that were removed from circuit queue; reset every
    * time when writing buffer stats to disk. */
   uint32_t processed_cells;
@@ -3961,6 +3965,10 @@ typedef struct {
   /** If true, the user wants us to collect statistics as entry node. */
   int EntryStatistics;
 
+  /** If true, the user wants us to collect statistics as hidden service
+   * directory, introduction point, or rendezvous point. */
+  int HiddenServiceStatistics;
+
   /** If true, include statistics file contents in extra-info documents. */
   int ExtraInfoStatistics;
 

src/or/rendcommon.c

@@ -924,6 +924,7 @@ rend_cache_lookup_v2_desc_as_dir(const char *desc_id, const char **desc)
 rend_cache_store_status_t
 rend_cache_store_v2_desc_as_dir(const char *desc)
 {
+  const or_options_t *options = get_options();
   rend_service_descriptor_t *parsed;
   char desc_id[DIGEST_LEN];
   char *intro_content;
@@ -1003,6 +1004,12 @@ rend_cache_store_v2_desc_as_dir(const char *desc)
     log_info(LD_REND, "Successfully stored service descriptor with desc ID "
                       "'%s' and len %d.",
              safe_str(desc_id_base32), (int)encoded_size);
+
+    /* Statistics: Note down this potentially new HS. */
+    if (options->HiddenServiceStatistics) {
+      rep_hist_stored_maybe_new_hs(e->parsed->pk);
+    }
+
     number_stored++;
     goto advance;
   skip:

src/or/rendmid.c

@@ -281,6 +281,7 @@ int
 rend_mid_rendezvous(or_circuit_t *circ, const uint8_t *request,
                     size_t request_len)
 {
+  const or_options_t *options = get_options();
   or_circuit_t *rend_circ;
   char hexid[9];
   int reason = END_CIRC_REASON_INTERNAL;
@@ -316,6 +317,12 @@ rend_mid_rendezvous(or_circuit_t *circ, const uint8_t *request,
     goto err;
   }
 
+  /* Statistics: Mark this circuit as an RP circuit so that we collect
+     stats from it. */
+  if (options->HiddenServiceStatistics) {
+    circ->circuit_carries_hs_traffic_stats = 1;
+  }
+
   /* Send the RENDEZVOUS2 cell to Alice. */
   if (relay_send_command_from_edge(0, TO_CIRCUIT(rend_circ),
                                    RELAY_COMMAND_RENDEZVOUS2,

src/or/rephist.c

@@ -2908,11 +2908,227 @@ rep_hist_log_circuit_handshake_stats(time_t now)
   memset(onion_handshakes_requested, 0, sizeof(onion_handshakes_requested));
 }
 
+/* Hidden service statistics section */
+
+/** Start of the current hidden service stats interval or 0 if we're
+ * not collecting hidden service statistics. */
+static time_t start_of_hs_stats_interval;
+
+/** Carries the various hidden service statistics, and any other
+ *  information needed. */
+typedef struct hs_stats_t {
+  /** How many relay cells have we seen as rendezvous points? */
+  int64_t rp_relay_cells_seen;
+
+  /** Set of unique public key digests we've seen this stat period
+   * (could also be implemented as sorted smartlist). */
+  digestmap_t *onions_seen_this_period;
+} hs_stats_t;
+
+/** Our statistics structure singleton. */
+static hs_stats_t *hs_stats = NULL;
+
+/** Allocate, initialize and return an hs_stats_t structure. */
+static hs_stats_t *
+hs_stats_new(void)
+{
+  hs_stats_t * hs_stats = tor_malloc_zero(sizeof(hs_stats_t));
+  hs_stats->onions_seen_this_period = digestmap_new();
+
+  return hs_stats;
+}
+
+/** Free an hs_stats_t structure. */
+static void
+hs_stats_free(hs_stats_t *hs_stats)
+{
+  if (!hs_stats) {
+    return;
+  }
+
+  digestmap_free(hs_stats->onions_seen_this_period, NULL);
+  tor_free(hs_stats);
+}
+
+/** Initialize hidden service statistics. */
+void
+rep_hist_hs_stats_init(time_t now)
+{
+  if (!hs_stats) {
+    hs_stats = hs_stats_new();
+  }
+
+  start_of_hs_stats_interval = now;
+}
+
+/** Clear history of hidden service statistics and set the measurement
+ * interval start to <b>now</b>. */
+static void
+rep_hist_reset_hs_stats(time_t now)
+{
+  if (!hs_stats) {
+    hs_stats = hs_stats_new();
+  }
+
+  hs_stats->rp_relay_cells_seen = 0;
+
+  digestmap_free(hs_stats->onions_seen_this_period, NULL);
+  hs_stats->onions_seen_this_period = digestmap_new();
+
+  start_of_hs_stats_interval = now;
+}
+
+/** Stop collecting hidden service stats in a way that we can re-start
+ * doing so in rep_hist_buffer_stats_init(). */
+void
+rep_hist_hs_stats_term(void)
+{
+  rep_hist_reset_hs_stats(0);
+}
+
+/** We saw a new HS relay cell, Count it! */
+void
+rep_hist_seen_new_rp_cell(void)
+{
+  if (!hs_stats) {
+    return; // We're not collecting stats
+  }
+
+  hs_stats->rp_relay_cells_seen++;
+}
+
+/** As HSDirs, we saw another hidden service with public key
+ *  <b>pubkey</b>. Check whether we have counted it before, if not
+ *  count it now! */
+void
+rep_hist_stored_maybe_new_hs(const crypto_pk_t *pubkey)
+{
+  char pubkey_hash[DIGEST_LEN];
+
+  if (!hs_stats) {
+    return; // We're not collecting stats
+  }
+
+  /* Get the digest of the pubkey which will be used to detect whether
+     we've seen this hidden service before or not.  */
+  if (crypto_pk_get_digest(pubkey, pubkey_hash) < 0) {
+    /*  This fail should not happen; key has been validated by
+        descriptor parsing code first. */
+    return;
+  }
+
+  /* Check if this is the first time we've seen this hidden
+     service. If it is, count it as new. */
+  if (!digestmap_get(hs_stats->onions_seen_this_period,
+                     pubkey_hash)) {
+    digestmap_set(hs_stats->onions_seen_this_period,
+                  pubkey_hash, (void*)(uintptr_t)1);
+  }
+}
+
+/* The number of cells that are supposed to be hidden from the adversary
+ * by adding noise from the Laplace distribution.  This value, divided by
+ * EPSILON, is Laplace parameter b. */
+#define REND_CELLS_DELTA_F 2048
+/* Security parameter for obfuscating number of cells with a value between
+ * 0 and 1.  Smaller values obfuscate observations more, but at the same
+ * time make statistics less usable. */
+#define REND_CELLS_EPSILON 0.3
+/* The number of cells that are supposed to be hidden from the adversary
+ * by rounding up to the next multiple of this number. */
+#define REND_CELLS_BIN_SIZE 1024
+/* The number of service identities that are supposed to be hidden from
+ * the adversary by adding noise from the Laplace distribution.  This
+ * value, divided by EPSILON, is Laplace parameter b. */
+#define ONIONS_SEEN_DELTA_F 8
+/* Security parameter for obfuscating number of service identities with a
+ * value between 0 and 1.  Smaller values obfuscate observations more, but
+ * at the same time make statistics less usable. */
+#define ONIONS_SEEN_EPSILON 0.3
+/* The number of service identities that are supposed to be hidden from
+ * the adversary by rounding up to the next multiple of this number. */
+#define ONIONS_SEEN_BIN_SIZE 8
+
+/** Allocate and return a string containing hidden service stats that
+ *  are meant to be placed in the extra-info descriptor. */
+static char *
+rep_hist_format_hs_stats(time_t now)
+{
+  char t[ISO_TIME_LEN+1];
+  char *hs_stats_string;
+  int64_t obfuscated_cells_seen;
+  int64_t obfuscated_onions_seen;
+
+  obfuscated_cells_seen = round_int64_to_next_multiple_of(
+                          hs_stats->rp_relay_cells_seen,
+                          REND_CELLS_BIN_SIZE);
+  obfuscated_cells_seen = add_laplace_noise(obfuscated_cells_seen,
+                          crypto_rand_double(),
+                          REND_CELLS_DELTA_F, REND_CELLS_EPSILON);
+  obfuscated_onions_seen = round_int64_to_next_multiple_of(digestmap_size(
+                           hs_stats->onions_seen_this_period),
+                           ONIONS_SEEN_BIN_SIZE);
+  obfuscated_onions_seen = add_laplace_noise(obfuscated_onions_seen,
+                           crypto_rand_double(), ONIONS_SEEN_DELTA_F,
+                           ONIONS_SEEN_EPSILON);
+
+  format_iso_time(t, now);
+  tor_asprintf(&hs_stats_string, "hidserv-stats-end %s (%d s)\n"
+               "hidserv-rend-relayed-cells "I64_FORMAT" delta_f=%d "
+                                           "epsilon=%.2f bin_size=%d\n"
+               "hidserv-dir-onions-seen "I64_FORMAT" delta_f=%d "
+                                        "epsilon=%.2f bin_size=%d\n",
+               t, (unsigned) (now - start_of_hs_stats_interval),
+               I64_PRINTF_ARG(obfuscated_cells_seen), REND_CELLS_DELTA_F,
+               REND_CELLS_EPSILON, REND_CELLS_BIN_SIZE,
+               I64_PRINTF_ARG(obfuscated_onions_seen),
+               ONIONS_SEEN_DELTA_F,
+               ONIONS_SEEN_EPSILON, ONIONS_SEEN_BIN_SIZE);
+
+  return hs_stats_string;
+}
+
+/** If 24 hours have passed since the beginning of the current HS
+ * stats period, write buffer stats to $DATADIR/stats/hidserv-stats
+ * (possibly overwriting an existing file) and reset counters.  Return
+ * when we would next want to write buffer stats or 0 if we never want to
+ * write. */
+time_t
+rep_hist_hs_stats_write(time_t now)
+{
+  char *str = NULL;
+
+  if (!start_of_hs_stats_interval) {
+    return 0; /* Not initialized. */
+  }
+
+  if (start_of_hs_stats_interval + WRITE_STATS_INTERVAL > now) {
+    goto done; /* Not ready to write */
+  }
+
+  /* Generate history string. */
+  str = rep_hist_format_hs_stats(now);
+
+  /* Reset HS history. */
+  rep_hist_reset_hs_stats(now);
+
+  /* Try to write to disk. */
+  if (!check_or_create_data_subdir("stats")) {
+    write_to_data_subdir("stats", "hidserv-stats", str,
+                         "hidden service stats");
+  }
+
+ done:
+  tor_free(str);
+  return start_of_hs_stats_interval + WRITE_STATS_INTERVAL;
+}
+
 /** Free all storage held by the OR/link history caches, by the
  * bandwidth history arrays, by the port history, or by statistics . */
 void
 rep_hist_free_all(void)
 {
+  hs_stats_free(hs_stats);
   digestmap_free(history_map, free_or_history);
   tor_free(read_array);
   tor_free(write_array);

src/or/rephist.h

@@ -99,6 +99,13 @@ void rep_hist_note_circuit_handshake_requested(uint16_t type);
 void rep_hist_note_circuit_handshake_assigned(uint16_t type);
 void rep_hist_log_circuit_handshake_stats(time_t now);
 
+void rep_hist_hs_stats_init(time_t now);
+void rep_hist_hs_stats_term(void);
+time_t rep_hist_hs_stats_write(time_t now);
+char *rep_hist_get_hs_stats_string(void);
+void rep_hist_seen_new_rp_cell(void);
+void rep_hist_stored_maybe_new_hs(const crypto_pk_t *pubkey);
+
 void rep_hist_free_all(void);
 
 #endif

src/or/router.c

@@ -2658,6 +2658,11 @@ extrainfo_dump_to_string(char **s_out, extrainfo_t *extrainfo,
                         "dirreq-stats-end", now, &contents) > 0) {
       smartlist_add(chunks, contents);
     }
+    if (options->HiddenServiceStatistics &&
+        load_stats_file("stats"PATH_SEPARATOR"hidserv-stats",
+                        "hidserv-stats-end", now, &contents) > 0) {
+      smartlist_add(chunks, contents);
+    }
     if (options->EntryStatistics &&
         load_stats_file("stats"PATH_SEPARATOR"entry-stats",
                         "entry-stats-end", now, &contents) > 0) {

src/test/test_util.c

@@ -4619,6 +4619,58 @@ test_util_round_to_next_multiple_of(void *arg)
   tt_assert(round_uint64_to_next_multiple_of(99,7) == 105);
   tt_assert(round_uint64_to_next_multiple_of(99,9) == 99);
 
+  tt_assert(round_int64_to_next_multiple_of(0,1) == 0);
+  tt_assert(round_int64_to_next_multiple_of(0,7) == 0);
+
+  tt_assert(round_int64_to_next_multiple_of(99,1) == 99);
+  tt_assert(round_int64_to_next_multiple_of(99,7) == 105);
+  tt_assert(round_int64_to_next_multiple_of(99,9) == 99);
+
+  tt_assert(round_int64_to_next_multiple_of(-99,1) == -99);
+  tt_assert(round_int64_to_next_multiple_of(-99,7) == -98);
+  tt_assert(round_int64_to_next_multiple_of(-99,9) == -99);
+
+  tt_assert(round_int64_to_next_multiple_of(INT64_MIN,2) == INT64_MIN);
+  tt_assert(round_int64_to_next_multiple_of(INT64_MAX,2) ==
+                                            INT64_MAX-INT64_MAX%2);
+ done:
+  ;
+}
+
+static void
+test_util_laplace(void *arg)
+{
+  /* Sample values produced using Python's SciPy:
+   *
+   * >>> from scipy.stats import laplace
+   * >>> laplace.ppf([-0.01, 0.0, 0.01, 0.5, 0.51, 0.99, 1.0, 1.01],
+     ...             loc = 24, scale = 24)
+   * array([          nan,          -inf,  -69.88855213,   24.        ,
+   *          24.48486498,  117.88855213,           inf,           nan])
+   */
+  const double mu = 24.0, b = 24.0;
+  const double delta_f = 15.0, epsilon = 0.3; /* b = 15.0 / 0.3 = 50.0 */
+  (void)arg;
+
+  tt_assert(isinf(sample_laplace_distribution(mu, b, 0.0)));
+  test_feq(-69.88855213, sample_laplace_distribution(mu, b, 0.01));
+  test_feq(24.0, sample_laplace_distribution(mu, b, 0.5));
+  test_feq(24.48486498, sample_laplace_distribution(mu, b, 0.51));
+  test_feq(117.88855213, sample_laplace_distribution(mu, b, 0.99));
+
+  /* >>> laplace.ppf([0.0, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99],
+   * ...             loc = 0, scale = 50)
+   * array([         -inf,  -80.47189562,  -34.65735903,    0.        ,
+   *          34.65735903,   80.47189562,  195.60115027])
+   */
+  tt_assert(LONG_MIN + 20 ==
+            add_laplace_noise(20, 0.0, 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));
  done:
   ;
 }
@@ -4880,6 +4932,7 @@ struct testcase_t util_tests[] = {
   UTIL_LEGACY(strtok),
   UTIL_LEGACY(di_ops),
   UTIL_TEST(round_to_next_multiple_of, 0),
+  UTIL_TEST(laplace, 0),
   UTIL_TEST(strclear, 0),
   UTIL_TEST(find_str_at_start_of_line, 0),
   UTIL_TEST(string_is_C_identifier, 0),