Tor  0.4.4.0-alpha-dev
prob_distr.c
Go to the documentation of this file.
1 /* Copyright (c) 2018-2020, The Tor Project, Inc. */
2 /* See LICENSE for licensing information */
3 
4 /**
5  * \file prob_distr.c
6  *
7  * \brief
8  * Implements various probability distributions.
9  * Almost all code is courtesy of Riastradh.
10  *
11  * \details
12  * Here are some details that might help you understand this file:
13  *
14  * - Throughout this file, `eps' means the largest relative error of a
15  * correctly rounded floating-point operation, which in binary64
16  * floating-point arithmetic is 2^-53. Here the relative error of a
17  * true value x from a computed value y is |x - y|/|x|. This
18  * definition of epsilon is conventional for numerical analysts when
19  * writing error analyses. (If your libm doesn't provide correctly
20  * rounded exp and log, their relative error is usually below 2*2^-53
21  * and probably closer to 1.1*2^-53 instead.)
22  *
23  * The C constant DBL_EPSILON is actually twice this, and should
24  * perhaps rather be named ulp(1) -- that is, it is the distance from
25  * 1 to the next greater floating-point number, which is usually of
26  * more interest to programmers and hardware engineers.
27  *
28  * Since this file is concerned mainly with error bounds rather than
29  * with low-level bit-hacking of floating-point numbers, we adopt the
30  * numerical analysts' definition in the comments, though we do use
31  * DBL_EPSILON in a handful of places where it is convenient to use
32  * some function of eps = DBL_EPSILON/2 in a case analysis.
33  *
34  * - In various functions (e.g. sample_log_logistic()) we jump through hoops so
35  * that we can use reals closer to 0 than closer to 1, since we achieve much
36  * greater accuracy for floating point numbers near 0. In particular, we can
37  * represent differences as small as 10^-300 for numbers near 0, but of no
38  * less than 10^-16 for numbers near 1.
39  **/
40 
41 #define PROB_DISTR_PRIVATE
42 
43 #include "orconfig.h"
44 
45 #include "lib/math/prob_distr.h"
46 
48 #include "lib/cc/ctassert.h"
49 #include "lib/log/util_bug.h"
50 
51 #include <float.h>
52 #include <math.h>
53 #include <stddef.h>
54 
55 #ifndef COCCI
56 /** Declare a function that downcasts from a generic dist struct to the actual
57  * subtype probablity distribution it represents. */
58 #define DECLARE_PROB_DISTR_DOWNCAST_FN(name) \
59  static inline \
60  const struct name##_t * \
61  dist_to_const_##name(const struct dist_t *obj) { \
62  tor_assert(obj->ops == &name##_ops); \
63  return SUBTYPE_P(obj, struct name ## _t, base); \
64  }
71 #endif /* !defined(COCCI) */
72 
73 /**
74  * Count number of one bits in 32-bit word.
75  */
76 static unsigned
77 bitcount32(uint32_t x)
78 {
79 
80  /* Count two-bit groups. */
81  x -= (x >> 1) & UINT32_C(0x55555555);
82 
83  /* Count four-bit groups. */
84  x = ((x >> 2) & UINT32_C(0x33333333)) + (x & UINT32_C(0x33333333));
85 
86  /* Count eight-bit groups. */
87  x = (x + (x >> 4)) & UINT32_C(0x0f0f0f0f);
88 
89  /* Sum all eight-bit groups, and extract the sum. */
90  return (x * UINT32_C(0x01010101)) >> 24;
91 }
92 
93 /**
94  * Count leading zeros in 32-bit word.
95  */
96 static unsigned
97 clz32(uint32_t x)
98 {
99 
100  /* Round up to a power of two. */
101  x |= x >> 1;
102  x |= x >> 2;
103  x |= x >> 4;
104  x |= x >> 8;
105  x |= x >> 16;
106 
107  /* Subtract count of one bits from 32. */
108  return (32 - bitcount32(x));
109 }
110 
111 /*
112  * Some lemmas that will be used throughout this file to prove various error
113  * bounds:
114  *
115  * Lemma 1. If |d| <= 1/2, then 1/(1 + d) <= 2.
116  *
117  * Proof. If 0 <= d <= 1/2, then 1 + d >= 1, so that 1/(1 + d) <= 1.
118  * If -1/2 <= d <= 0, then 1 + d >= 1/2, so that 1/(1 + d) <= 2. QED.
119  *
120  * Lemma 2. If b = a*(1 + d)/(1 + d') for |d'| < 1/2 and nonzero a, b,
121  * then b = a*(1 + e) for |e| <= 2|d' - d|.
122  *
123  * Proof. |a - b|/|a|
124  * = |a - a*(1 + d)/(1 + d')|/|a|
125  * = |1 - (1 + d)/(1 + d')|
126  * = |(1 + d' - 1 - d)/(1 + d')|
127  * = |(d' - d)/(1 + d')|
128  * <= 2|d' - d|, by Lemma 1,
129  *
130  * QED.
131  *
132  * Lemma 3. For |d|, |d'| < 1/4,
133  *
134  * |log((1 + d)/(1 + d'))| <= 4|d - d'|.
135  *
136  * Proof. Write
137  *
138  * log((1 + d)/(1 + d'))
139  * = log(1 + (1 + d)/(1 + d') - 1)
140  * = log(1 + (1 + d - 1 - d')/(1 + d')
141  * = log(1 + (d - d')/(1 + d')).
142  *
143  * By Lemma 1, |(d - d')/(1 + d')| < 2|d' - d| < 1, so the Taylor
144  * series of log(1 + x) converges absolutely for (d - d')/(1 + d'),
145  * and thus we have
146  *
147  * |log(1 + (d - d')/(1 + d'))|
148  * = |\sum_{n=1}^\infty ((d - d')/(1 + d'))^n/n|
149  * <= \sum_{n=1}^\infty |(d - d')/(1 + d')|^n/n
150  * <= \sum_{n=1}^\infty |2(d' - d)|^n/n
151  * <= \sum_{n=1}^\infty |2(d' - d)|^n
152  * = 1/(1 - |2(d' - d)|)
153  * <= 4|d' - d|,
154  *
155  * QED.
156  *
157  * Lemma 4. If 1/e <= 1 + x <= e, then
158  *
159  * log(1 + (1 + d) x) = (1 + d') log(1 + x)
160  *
161  * for |d'| < 8|d|.
162  *
163  * Proof. Write
164  *
165  * log(1 + (1 + d) x)
166  * = log(1 + x + x*d)
167  * = log((1 + x) (1 + x + x*d)/(1 + x))
168  * = log(1 + x) + log((1 + x + x*d)/(1 + x))
169  * = log(1 + x) (1 + log((1 + x + x*d)/(1 + x))/log(1 + x)).
170  *
171  * The relative error is bounded by
172  *
173  * |log((1 + x + x*d)/(1 + x))/log(1 + x)|
174  * <= 4|x + x*d - x|/|log(1 + x)|, by Lemma 3,
175  * = 4|x*d|/|log(1 + x)|
176  * < 8|d|,
177  *
178  * since in this range 0 < 1 - 1/e < x/log(1 + x) <= e - 1 < 2. QED.
179  */
180 
181 /**
182  * Compute the logistic function: f(x) = 1/(1 + e^{-x}) = e^x/(1 + e^x).
183  * Maps a log-odds-space probability in [-infinity, +infinity] into a
184  * direct-space probability in [0,1]. Inverse of logit.
185  *
186  * Ill-conditioned for large x; the identity logistic(-x) = 1 -
187  * logistic(x) and the function logistichalf(x) = logistic(x) - 1/2 may
188  * help to rearrange a computation.
189  *
190  * This implementation gives relative error bounded by 7 eps.
191  */
192 STATIC double
193 logistic(double x)
194 {
195  if (x <= log(DBL_EPSILON/2)) {
196  /*
197  * If x <= log(DBL_EPSILON/2) = log(eps), then e^x <= eps. In this case
198  * we will approximate the logistic() function with e^x because the
199  * relative error is less than eps. Here is a calculation of the
200  * relative error between the logistic() function and e^x and a proof
201  * that it's less than eps:
202  *
203  * |e^x - e^x/(1 + e^x)|/|e^x/(1 + e^x)|
204  * <= |1 - 1/(1 + e^x)|*|1 + e^x|
205  * = |e^x/(1 + e^x)|*|1 + e^x|
206  * = |e^x|
207  * <= eps.
208  */
209  return exp(x); /* return e^x */
210  } else if (x <= -log(DBL_EPSILON/2)) {
211  /*
212  * e^{-x} > 0, so 1 + e^{-x} > 1, and 0 < 1/(1 +
213  * e^{-x}) < 1; further, since e^{-x} < 1 + e^{-x}, we
214  * also have 0 < 1/(1 + e^{-x}) < 1. Thus, if exp has
215  * relative error d0, + has relative error d1, and /
216  * has relative error d2, then we get
217  *
218  * (1 + d2)/[(1 + (1 + d0) e^{-x})(1 + d1)]
219  * = (1 + d0)/[1 + e^{-x} + d0 e^{-x}
220  * + d1 + d1 e^{-x} + d0 d1 e^{-x}]
221  * = (1 + d0)/[(1 + e^{-x})
222  * * (1 + d0 e^{-x}/(1 + e^{-x})
223  * + d1/(1 + e^{-x})
224  * + d0 d1 e^{-x}/(1 + e^{-x}))].
225  * = (1 + d0)/[(1 + e^{-x})(1 + d')]
226  * = [1/(1 + e^{-x})] (1 + d0)/(1 + d')
227  *
228  * where
229  *
230  * d' = d0 e^{-x}/(1 + e^{-x})
231  * + d1/(1 + e^{-x})
232  * + d0 d1 e^{-x}/(1 + e^{-x}).
233  *
234  * By Lemma 2 this relative error is bounded by
235  *
236  * 2|d0 - d'|
237  * = 2|d0 - d0 e^{-x}/(1 + e^{-x})
238  * - d1/(1 + e^{-x})
239  * - d0 d1 e^{-x}/(1 + e^{-x})|
240  * <= 2|d0| + 2|d0 e^{-x}/(1 + e^{-x})|
241  * + 2|d1/(1 + e^{-x})|
242  * + 2|d0 d1 e^{-x}/(1 + e^{-x})|
243  * <= 2|d0| + 2|d0| + 2|d1| + 2|d0 d1|
244  * <= 4|d0| + 2|d1| + 2|d0 d1|
245  * <= 6 eps + 2 eps^2.
246  */
247  return 1/(1 + exp(-x));
248  } else {
249  /*
250  * e^{-x} <= eps, so the relative error of 1 from 1/(1
251  * + e^{-x}) is
252  *
253  * |1/(1 + e^{-x}) - 1|/|1/(1 + e^{-x})|
254  * = |e^{-x}/(1 + e^{-x})|/|1/(1 + e^{-x})|
255  * = |e^{-x}|
256  * <= eps.
257  *
258  * This computation avoids an intermediate overflow
259  * exception, although the effect on the result is
260  * harmless.
261  *
262  * XXX Should maybe raise inexact here.
263  */
264  return 1;
265  }
266 }
267 
268 /**
269  * Compute the logit function: log p/(1 - p). Defined on [0,1]. Maps
270  * a direct-space probability in [0,1] to a log-odds-space probability
271  * in [-infinity, +infinity]. Inverse of logistic.
272  *
273  * Ill-conditioned near 1/2 and 1; the identity logit(1 - p) =
274  * -logit(p) and the function logithalf(p0) = logit(1/2 + p0) may help
275  * to rearrange a computation for p in [1/(1 + e), 1 - 1/(1 + e)].
276  *
277  * This implementation gives relative error bounded by 10 eps.
278  */
279 STATIC double
280 logit(double p)
281 {
282 
283  /* logistic(-1) <= p <= logistic(+1) */
284  if (1/(1 + exp(1)) <= p && p <= 1/(1 + exp(-1))) {
285  /*
286  * For inputs near 1/2, we want to compute log1p(near
287  * 0) rather than log(near 1), so write this as:
288  *
289  * log(p/(1 - p)) = -log((1 - p)/p)
290  * = -log(1 + (1 - p)/p - 1)
291  * = -log(1 + (1 - p - p)/p)
292  * = -log(1 + (1 - 2p)/p).
293  *
294  * Since p = 2p/2 <= 1 <= 2*2p = 4p, the floating-point
295  * evaluation of 1 - 2p is exact; the only error arises
296  * from division and log1p. First, note that if
297  * logistic(-1) <= p <= logistic(+1), (1 - 2p)/p lies
298  * in the bounds of Lemma 4.
299  *
300  * If division has relative error d0 and log1p has
301  * relative error d1, the outcome is
302  *
303  * -(1 + d1) log(1 + (1 - 2p) (1 + d0)/p)
304  * = -(1 + d1) (1 + d') log(1 + (1 - 2p)/p)
305  * = -(1 + d1 + d' + d1 d') log(1 + (1 - 2p)/p).
306  *
307  * where |d'| < 8|d0| by Lemma 4. The relative error
308  * is then bounded by
309  *
310  * |d1 + d' + d1 d'|
311  * <= |d1| + 8|d0| + 8|d1 d0|
312  * <= 9 eps + 8 eps^2.
313  */
314  return -log1p((1 - 2*p)/p);
315  } else {
316  /*
317  * For inputs near 0, although 1 - p may be rounded to
318  * 1, it doesn't matter much because the magnitude of
319  * the result is so much larger. For inputs near 1, we
320  * can compute 1 - p exactly, although the precision on
321  * the input is limited so we won't ever get more than
322  * about 700 for the output.
323  *
324  * If - has relative error d0, / has relative error d1,
325  * and log has relative error d2, then
326  *
327  * (1 + d2) log((1 + d0) p/[(1 - p)(1 + d1)])
328  * = (1 + d2) [log(p/(1 - p)) + log((1 + d0)/(1 + d1))]
329  * = log(p/(1 - p)) + d2 log(p/(1 - p))
330  * + (1 + d2) log((1 + d0)/(1 + d1))
331  * = log(p/(1 - p))*[1 + d2 +
332  * + (1 + d2) log((1 + d0)/(1 + d1))/log(p/(1 - p))]
333  *
334  * Since 0 <= p < logistic(-1) or logistic(+1) < p <=
335  * 1, we have |log(p/(1 - p))| > 1. Hence this error
336  * is bounded by
337  *
338  * |d2 + (1 + d2) log((1 + d0)/(1 + d1))/log(p/(1 - p))|
339  * <= |d2| + |(1 + d2) log((1 + d0)/(1 + d1))
340  * / log(p/(1 - p))|
341  * <= |d2| + |(1 + d2) log((1 + d0)/(1 + d1))|
342  * <= |d2| + 4|(1 + d2) (d0 - d1)|, by Lemma 3,
343  * <= |d2| + 4|d0 - d1 + d2 d0 - d1 d0|
344  * <= |d2| + 4|d0| + 4|d1| + 4|d2 d0| + 4|d1 d0|
345  * <= 9 eps + 8 eps^2.
346  */
347  return log(p/(1 - p));
348  }
349 }
350 
351 /**
352  * Compute the logit function, translated in input by 1/2: logithalf(p)
353  * = logit(1/2 + p). Defined on [-1/2, 1/2]. Inverse of logistichalf.
354  *
355  * Ill-conditioned near +/-1/2. If |p0| > 1/2 - 1/(1 + e), it may be
356  * better to compute 1/2 + p0 or -1/2 - p0 and to use logit instead.
357  * This implementation gives relative error bounded by 34 eps.
358  */
359 STATIC double
360 logithalf(double p0)
361 {
362 
363  if (fabs(p0) <= 0.5 - 1/(1 + exp(1))) {
364  /*
365  * logit(1/2 + p0)
366  * = log((1/2 + p0)/(1 - (1/2 + p0)))
367  * = log((1/2 + p0)/(1/2 - p0))
368  * = log(1 + (1/2 + p0)/(1/2 - p0) - 1)
369  * = log(1 + (1/2 + p0 - (1/2 - p0))/(1/2 - p0))
370  * = log(1 + (1/2 + p0 - 1/2 + p0)/(1/2 - p0))
371  * = log(1 + 2 p0/(1/2 - p0))
372  *
373  * If the error of subtraction is d0, the error of
374  * division is d1, and the error of log1p is d2, then
375  * what we compute is
376  *
377  * (1 + d2) log(1 + (1 + d1) 2 p0/[(1 + d0) (1/2 - p0)])
378  * = (1 + d2) log(1 + (1 + d') 2 p0/(1/2 - p0))
379  * = (1 + d2) (1 + d'') log(1 + 2 p0/(1/2 - p0))
380  * = (1 + d2 + d'' + d2 d'') log(1 + 2 p0/(1/2 - p0)),
381  *
382  * where |d'| < 2|d0 - d1| <= 4 eps by Lemma 2, and
383  * |d''| < 8|d'| < 32 eps by Lemma 4 since
384  *
385  * 1/e <= 1 + 2*p0/(1/2 - p0) <= e
386  *
387  * when |p0| <= 1/2 - 1/(1 + e). Hence the relative
388  * error is bounded by
389  *
390  * |d2 + d'' + d2 d''|
391  * <= |d2| + |d''| + |d2 d''|
392  * <= |d1| + 32 |d0| + 32 |d1 d0|
393  * <= 33 eps + 32 eps^2.
394  */
395  return log1p(2*p0/(0.5 - p0));
396  } else {
397  /*
398  * We have a choice of computing logit(1/2 + p0) or
399  * -logit(1 - (1/2 + p0)) = -logit(1/2 - p0). It
400  * doesn't matter which way we do this: either way,
401  * since 1/2 p0 <= 1/2 <= 2 p0, the sum and difference
402  * are computed exactly. So let's do the one that
403  * skips the final negation.
404  *
405  * The result is
406  *
407  * (1 + d1) log((1 + d0) (1/2 + p0)/[(1 + d2) (1/2 - p0)])
408  * = (1 + d1) (1 + log((1 + d0)/(1 + d2))
409  * / log((1/2 + p0)/(1/2 - p0)))
410  * * log((1/2 + p0)/(1/2 - p0))
411  * = (1 + d') log((1/2 + p0)/(1/2 - p0))
412  * = (1 + d') logit(1/2 + p0)
413  *
414  * where
415  *
416  * d' = d1 + log((1 + d0)/(1 + d2))/logit(1/2 + p0)
417  * + d1 log((1 + d0)/(1 + d2))/logit(1/2 + p0).
418  *
419  * For |p| > 1/2 - 1/(1 + e), logit(1/2 + p0) > 1.
420  * Provided |d0|, |d2| < 1/4, by Lemma 3 we have
421  *
422  * |log((1 + d0)/(1 + d2))| <= 4|d0 - d2|.
423  *
424  * Hence the relative error is bounded by
425  *
426  * |d'| <= |d1| + 4|d0 - d2| + 4|d1| |d0 - d2|
427  * <= |d1| + 4|d0| + 4|d2| + 4|d1 d0| + 4|d1 d2|
428  * <= 9 eps + 8 eps^2.
429  */
430  return log((0.5 + p0)/(0.5 - p0));
431  }
432 }
433 
434 /*
435  * The following random_uniform_01 is tailored for IEEE 754 binary64
436  * floating-point or smaller. It can be adapted to larger
437  * floating-point formats like i387 80-bit or IEEE 754 binary128, but
438  * it may require sampling more bits.
439  */
440 CTASSERT(FLT_RADIX == 2);
441 CTASSERT(-DBL_MIN_EXP <= 1021);
442 CTASSERT(DBL_MANT_DIG <= 53);
443 
444 /**
445  * Draw a floating-point number in [0, 1] with uniform distribution.
446  *
447  * Note that the probability of returning 0 is less than 2^-1074, so
448  * callers need not check for it. However, callers that cannot handle
449  * rounding to 1 must deal with that, because it occurs with
450  * probability 2^-54, which is small but nonnegligible.
451  */
452 STATIC double
454 {
455  uint32_t z, x, hi, lo;
456  double s;
457 
458  /*
459  * Draw an exponent, geometrically distributed, but give up if
460  * we get a run of more than 1088 zeros, which really means the
461  * system is broken.
462  */
463  z = 0;
464  while ((x = crypto_fast_rng_get_u32(get_thread_fast_rng())) == 0) {
465  if (z >= 1088)
466  /* Your bit sampler is broken. Go home. */
467  return 0;
468  z += 32;
469  }
470  z += clz32(x);
471 
472  /*
473  * Pick 32-bit halves of an odd normalized significand.
474  * Picking it odd breaks ties in the subsequent rounding, which
475  * occur only with measure zero in the uniform distribution on
476  * [0, 1].
477  */
478  hi = crypto_fast_rng_get_u32(get_thread_fast_rng()) | UINT32_C(0x80000000);
479  lo = crypto_fast_rng_get_u32(get_thread_fast_rng()) | UINT32_C(0x00000001);
480 
481  /* Round to nearest scaled significand in [2^63, 2^64]. */
482  s = hi*(double)4294967296 + lo;
483 
484  /* Rescale into [1/2, 1] and apply exponent in one swell foop. */
485  return s * ldexp(1, -(64 + z));
486 }
487 
488 /*******************************************************************/
489 
490 /* Functions for specific probability distributions start here: */
491 
492 /*
493  * Logistic(mu, sigma) distribution, supported on (-infinity,+infinity)
494  *
495  * This is the uniform distribution on [0,1] mapped into log-odds
496  * space, scaled by sigma and translated by mu.
497  *
498  * pdf(x) = e^{-(x - mu)/sigma} sigma (1 + e^{-(x - mu)/sigma})^2
499  * cdf(x) = 1/(1 + e^{-(x - mu)/sigma}) = logistic((x - mu)/sigma)
500  * sf(x) = 1 - cdf(x) = 1 - logistic((x - mu)/sigma = logistic(-(x - mu)/sigma)
501  * icdf(p) = mu + sigma log p/(1 - p) = mu + sigma logit(p)
502  * isf(p) = mu + sigma log (1 - p)/p = mu - sigma logit(p)
503  */
504 
505 /**
506  * Compute the CDF of the Logistic(mu, sigma) distribution: the
507  * logistic function. Well-conditioned for negative inputs and small
508  * positive inputs; ill-conditioned for large positive inputs.
509  */
510 STATIC double
511 cdf_logistic(double x, double mu, double sigma)
512 {
513  return logistic((x - mu)/sigma);
514 }
515 
516 /**
517  * Compute the SF of the Logistic(mu, sigma) distribution: the logistic
518  * function reflected over the y axis. Well-conditioned for positive
519  * inputs and small negative inputs; ill-conditioned for large negative
520  * inputs.
521  */
522 STATIC double
523 sf_logistic(double x, double mu, double sigma)
524 {
525  return logistic(-(x - mu)/sigma);
526 }
527 
528 /**
529  * Compute the inverse of the CDF of the Logistic(mu, sigma)
530  * distribution: the logit function. Well-conditioned near 0;
531  * ill-conditioned near 1/2 and 1.
532  */
533 STATIC double
534 icdf_logistic(double p, double mu, double sigma)
535 {
536  return mu + sigma*logit(p);
537 }
538 
539 /**
540  * Compute the inverse of the SF of the Logistic(mu, sigma)
541  * distribution: the -logit function. Well-conditioned near 0;
542  * ill-conditioned near 1/2 and 1.
543  */
544 STATIC double
545 isf_logistic(double p, double mu, double sigma)
546 {
547  return mu - sigma*logit(p);
548 }
549 
550 /*
551  * LogLogistic(alpha, beta) distribution, supported on (0, +infinity).
552  *
553  * This is the uniform distribution on [0,1] mapped into odds space,
554  * scaled by positive alpha and shaped by positive beta.
555  *
556  * Equivalent to computing exp of a Logistic(log alpha, 1/beta) sample.
557  * (Name arises because the pdf has LogLogistic(x; alpha, beta) =
558  * Logistic(log x; log alpha, 1/beta) and mathematicians got their
559  * covariance contravariant.)
560  *
561  * pdf(x) = (beta/alpha) (x/alpha)^{beta - 1}/(1 + (x/alpha)^beta)^2
562  * = (1/e^mu sigma) (x/e^mu)^{1/sigma - 1} /
563  * (1 + (x/e^mu)^{1/sigma})^2
564  * cdf(x) = 1/(1 + (x/alpha)^-beta) = 1/(1 + (x/e^mu)^{-1/sigma})
565  * = 1/(1 + (e^{log x}/e^mu)^{-1/sigma})
566  * = 1/(1 + (e^{log x - mu})^{-1/sigma})
567  * = 1/(1 + e^{-(log x - mu)/sigma})
568  * = logistic((log x - mu)/sigma)
569  * = logistic((log x - log alpha)/(1/beta))
570  * sf(x) = 1 - 1/(1 + (x/alpha)^-beta)
571  * = (x/alpha)^-beta/(1 + (x/alpha)^-beta)
572  * = 1/((x/alpha)^beta + 1)
573  * = 1/(1 + (x/alpha)^beta)
574  * icdf(p) = alpha (p/(1 - p))^{1/beta}
575  * = alpha e^{logit(p)/beta}
576  * = e^{mu + sigma logit(p)}
577  * isf(p) = alpha ((1 - p)/p)^{1/beta}
578  * = alpha e^{-logit(p)/beta}
579  * = e^{mu - sigma logit(p)}
580  */
581 
582 /**
583  * Compute the CDF of the LogLogistic(alpha, beta) distribution.
584  * Well-conditioned for all x and alpha, and the condition number
585  *
586  * -beta/[1 + (x/alpha)^{-beta}]
587  *
588  * grows linearly with beta.
589  *
590  * Loosely, the relative error of this implementation is bounded by
591  *
592  * 4 eps + 2 eps^2 + O(beta eps),
593  *
594  * so don't bother trying this for beta anywhere near as large as
595  * 1/eps, around which point it levels off at 1.
596  */
597 STATIC double
598 cdf_log_logistic(double x, double alpha, double beta)
599 {
600  /*
601  * Let d0 be the error of x/alpha; d1, of pow; d2, of +; and
602  * d3, of the final quotient. The exponentiation gives
603  *
604  * ((1 + d0) x/alpha)^{-beta}
605  * = (x/alpha)^{-beta} (1 + d0)^{-beta}
606  * = (x/alpha)^{-beta} (1 + (1 + d0)^{-beta} - 1)
607  * = (x/alpha)^{-beta} (1 + d')
608  *
609  * where d' = (1 + d0)^{-beta} - 1. If y = (x/alpha)^{-beta},
610  * the denominator is
611  *
612  * (1 + d2) (1 + (1 + d1) (1 + d') y)
613  * = (1 + d2) (1 + y + (d1 + d' + d1 d') y)
614  * = 1 + y + (1 + d2) (d1 + d' + d1 d') y
615  * = (1 + y) (1 + (1 + d2) (d1 + d' + d1 d') y/(1 + y))
616  * = (1 + y) (1 + d''),
617  *
618  * where d'' = (1 + d2) (d1 + d' + d1 d') y/(1 + y). The
619  * final result is
620  *
621  * (1 + d3) / [(1 + d2) (1 + d'') (1 + y)]
622  * = (1 + d''') / (1 + y)
623  *
624  * for |d'''| <= 2|d3 - d''| by Lemma 2 as long as |d''| < 1/2
625  * (which may not be the case for very large beta). This
626  * relative error is therefore bounded by
627  *
628  * |d'''|
629  * <= 2|d3 - d''|
630  * <= 2|d3| + 2|(1 + d2) (d1 + d' + d1 d') y/(1 + y)|
631  * <= 2|d3| + 2|(1 + d2) (d1 + d' + d1 d')|
632  * = 2|d3| + 2|d1 + d' + d1 d' + d2 d1 + d2 d' + d2 d1 d'|
633  * <= 2|d3| + 2|d1| + 2|d'| + 2|d1 d'| + 2|d2 d1| + 2|d2 d'|
634  * + 2|d2 d1 d'|
635  * <= 4 eps + 2 eps^2 + (2 + 2 eps + 2 eps^2) |d'|.
636  *
637  * Roughly, |d'| = |(1 + d0)^{-beta} - 1| grows like beta eps,
638  * until it levels off at 1.
639  */
640  return 1/(1 + pow(x/alpha, -beta));
641 }
642 
643 /**
644  * Compute the SF of the LogLogistic(alpha, beta) distribution.
645  * Well-conditioned for all x and alpha, and the condition number
646  *
647  * beta/[1 + (x/alpha)^beta]
648  *
649  * grows linearly with beta.
650  *
651  * Loosely, the relative error of this implementation is bounded by
652  *
653  * 4 eps + 2 eps^2 + O(beta eps)
654  *
655  * so don't bother trying this for beta anywhere near as large as
656  * 1/eps, beyond which point it grows unbounded.
657  */
658 STATIC double
659 sf_log_logistic(double x, double alpha, double beta)
660 {
661  /*
662  * The error analysis here is essentially the same as in
663  * cdf_log_logistic, except that rather than levelling off at
664  * 1, |(1 + d0)^beta - 1| grows unbounded.
665  */
666  return 1/(1 + pow(x/alpha, beta));
667 }
668 
669 /**
670  * Compute the inverse of the CDF of the LogLogistic(alpha, beta)
671  * distribution. Ill-conditioned for p near 1 and beta near 0 with
672  * condition number 1/[beta (1 - p)].
673  */
674 STATIC double
675 icdf_log_logistic(double p, double alpha, double beta)
676 {
677  return alpha*pow(p/(1 - p), 1/beta);
678 }
679 
680 /**
681  * Compute the inverse of the SF of the LogLogistic(alpha, beta)
682  * distribution. Ill-conditioned for p near 1 and for large beta, with
683  * condition number -1/[beta (1 - p)].
684  */
685 STATIC double
686 isf_log_logistic(double p, double alpha, double beta)
687 {
688  return alpha*pow((1 - p)/p, 1/beta);
689 }
690 
691 /*
692  * Weibull(lambda, k) distribution, supported on (0, +infinity).
693  *
694  * pdf(x) = (k/lambda) (x/lambda)^{k - 1} e^{-(x/lambda)^k}
695  * cdf(x) = 1 - e^{-(x/lambda)^k}
696  * icdf(p) = lambda * (-log (1 - p))^{1/k}
697  * sf(x) = e^{-(x/lambda)^k}
698  * isf(p) = lambda * (-log p)^{1/k}
699  */
700 
701 /**
702  * Compute the CDF of the Weibull(lambda, k) distribution.
703  * Well-conditioned for small x and k, and for large lambda --
704  * condition number
705  *
706  * -k (x/lambda)^k exp(-(x/lambda)^k)/[exp(-(x/lambda)^k) - 1]
707  *
708  * grows linearly with k, x^k, and lambda^{-k}.
709  */
710 STATIC double
711 cdf_weibull(double x, double lambda, double k)
712 {
713  return -expm1(-pow(x/lambda, k));
714 }
715 
716 /**
717  * Compute the SF of the Weibull(lambda, k) distribution.
718  * Well-conditioned for small x and k, and for large lambda --
719  * condition number
720  *
721  * -k (x/lambda)^k
722  *
723  * grows linearly with k, x^k, and lambda^{-k}.
724  */
725 STATIC double
726 sf_weibull(double x, double lambda, double k)
727 {
728  return exp(-pow(x/lambda, k));
729 }
730 
731 /**
732  * Compute the inverse of the CDF of the Weibull(lambda, k)
733  * distribution. Ill-conditioned for p near 1, and for k near 0;
734  * condition number is
735  *
736  * (p/(1 - p))/(k log(1 - p)).
737  */
738 STATIC double
739 icdf_weibull(double p, double lambda, double k)
740 {
741  return lambda*pow(-log1p(-p), 1/k);
742 }
743 
744 /**
745  * Compute the inverse of the SF of the Weibull(lambda, k)
746  * distribution. Ill-conditioned for p near 0, and for k near 0;
747  * condition number is
748  *
749  * 1/(k log(p)).
750  */
751 STATIC double
752 isf_weibull(double p, double lambda, double k)
753 {
754  return lambda*pow(-log(p), 1/k);
755 }
756 
757 /*
758  * GeneralizedPareto(mu, sigma, xi), supported on (mu, +infinity) for
759  * nonnegative xi, or (mu, mu - sigma/xi) for negative xi.
760  *
761  * Samples:
762  * = mu - sigma log U, if xi = 0;
763  * = mu + sigma (U^{-xi} - 1)/xi = mu + sigma*expm1(-xi log U)/xi, if xi =/= 0,
764  * where U is uniform on (0,1].
765  * = mu + sigma (e^{xi X} - 1)/xi,
766  * where X has standard exponential distribution.
767  *
768  * pdf(x) = sigma^{-1} (1 + xi (x - mu)/sigma)^{-(1 + 1/xi)}
769  * cdf(x) = 1 - (1 + xi (x - mu)/sigma)^{-1/xi}
770  * = 1 - e^{-log(1 + xi (x - mu)/sigma)/xi}
771  * --> 1 - e^{-(x - mu)/sigma} as xi --> 0
772  * sf(x) = (1 + xi (x - mu)/sigma)^{-1/xi}
773  * --> e^{-(x - mu)/sigma} as xi --> 0
774  * icdf(p) = mu + sigma*(p^{-xi} - 1)/xi
775  * = mu + sigma*expm1(-xi log p)/xi
776  * --> mu + sigma*log p as xi --> 0
777  * isf(p) = mu + sigma*((1 - p)^{xi} - 1)/xi
778  * = mu + sigma*expm1(-xi log1p(-p))/xi
779  * --> mu + sigma*log1p(-p) as xi --> 0
780  */
781 
782 /**
783  * Compute the CDF of the GeneralizedPareto(mu, sigma, xi)
784  * distribution. Well-conditioned everywhere. For standard
785  * distribution (mu=0, sigma=1), condition number
786  *
787  * (x/(1 + x xi)) / ((1 + x xi)^{1/xi} - 1)
788  *
789  * is bounded by 1, attained only at x = 0.
790  */
791 STATIC double
792 cdf_genpareto(double x, double mu, double sigma, double xi)
793 {
794  double x_0 = (x - mu)/sigma;
795 
796  /*
797  * log(1 + xi x_0)/xi
798  * = (-1/xi) \sum_{n=1}^infinity (-xi x_0)^n/n
799  * = (-1/xi) (-xi x_0 + \sum_{n=2}^infinity (-xi x_0)^n/n)
800  * = x_0 - (1/xi) \sum_{n=2}^infinity (-xi x_0)^n/n
801  * = x_0 - x_0 \sum_{n=2}^infinity (-xi x_0)^{n-1}/n
802  * = x_0 (1 - d),
803  *
804  * where d = \sum_{n=2}^infinity (-xi x_0)^{n-1}/n. If |xi| <
805  * eps/4|x_0|, then
806  *
807  * |d| <= \sum_{n=2}^infinity (eps/4)^{n-1}/n
808  * <= \sum_{n=2}^infinity (eps/4)^{n-1}
809  * = \sum_{n=1}^infinity (eps/4)^n
810  * = (eps/4) \sum_{n=0}^infinity (eps/4)^n
811  * = (eps/4)/(1 - eps/4)
812  * < eps/2
813  *
814  * for any 0 < eps < 2. Thus, the relative error of x_0 from
815  * log(1 + xi x_0)/xi is bounded by eps.
816  */
817  if (fabs(xi) < 1e-17/x_0)
818  return -expm1(-x_0);
819  else
820  return -expm1(-log1p(xi*x_0)/xi);
821 }
822 
823 /**
824  * Compute the SF of the GeneralizedPareto(mu, sigma, xi) distribution.
825  * For standard distribution (mu=0, sigma=1), ill-conditioned for xi
826  * near 0; condition number
827  *
828  * -x (1 + x xi)^{(-1 - xi)/xi}/(1 + x xi)^{-1/xi}
829  * = -x (1 + x xi)^{-1/xi - 1}/(1 + x xi)^{-1/xi}
830  * = -(x/(1 + x xi)) (1 + x xi)^{-1/xi}/(1 + x xi)^{-1/xi}
831  * = -x/(1 + x xi)
832  *
833  * is bounded by 1/xi.
834  */
835 STATIC double
836 sf_genpareto(double x, double mu, double sigma, double xi)
837 {
838  double x_0 = (x - mu)/sigma;
839 
840  if (fabs(xi) < 1e-17/x_0)
841  return exp(-x_0);
842  else
843  return exp(-log1p(xi*x_0)/xi);
844 }
845 
846 /**
847  * Compute the inverse of the CDF of the GeneralizedPareto(mu, sigma,
848  * xi) distribution. Ill-conditioned for p near 1; condition number is
849  *
850  * xi (p/(1 - p))/(1 - (1 - p)^xi)
851  */
852 STATIC double
853 icdf_genpareto(double p, double mu, double sigma, double xi)
854 {
855  /*
856  * To compute f(xi) = (U^{-xi} - 1)/xi = (e^{-xi log U} - 1)/xi
857  * for xi near zero (note f(xi) --> -log U as xi --> 0), write
858  * the absolutely convergent Taylor expansion
859  *
860  * f(xi) = (1/xi)*(-xi log U + \sum_{n=2}^infinity (-xi log U)^n/n!
861  * = -log U + (1/xi)*\sum_{n=2}^infinity (-xi log U)^n/n!
862  * = -log U + \sum_{n=2}^infinity xi^{n-1} (-log U)^n/n!
863  * = -log U - log U \sum_{n=2}^infinity (-xi log U)^{n-1}/n!
864  * = -log U (1 + \sum_{n=2}^infinity (-xi log U)^{n-1}/n!).
865  *
866  * Let d = \sum_{n=2}^infinity (-xi log U)^{n-1}/n!. What do we
867  * lose if we discard it and use -log U as an approximation to
868  * f(xi)? If |xi| < eps/-4log U, then
869  *
870  * |d| <= \sum_{n=2}^infinity |xi log U|^{n-1}/n!
871  * <= \sum_{n=2}^infinity (eps/4)^{n-1}/n!
872  * <= \sum_{n=1}^infinity (eps/4)^n
873  * = (eps/4) \sum_{n=0}^infinity (eps/4)^n
874  * = (eps/4)/(1 - eps/4)
875  * < eps/2,
876  *
877  * for any 0 < eps < 2. Hence, as long as |xi| < eps/-2log U,
878  * f(xi) = -log U (1 + d) for |d| <= eps/2. |d| is the
879  * relative error of f(xi) from -log U; from this bound, the
880  * relative error of -log U from f(xi) is at most (eps/2)/(1 -
881  * eps/2) = eps/2 + (eps/2)^2 + (eps/2)^3 + ... < eps for 0 <
882  * eps < 1. Since -log U < 1000 for all U in (0, 1] in
883  * binary64 floating-point, we can safely cut xi off at 1e-20 <
884  * eps/4000 and attain <1ulp error from series truncation.
885  */
886  if (fabs(xi) <= 1e-20)
887  return mu - sigma*log1p(-p);
888  else
889  return mu + sigma*expm1(-xi*log1p(-p))/xi;
890 }
891 
892 /**
893  * Compute the inverse of the SF of the GeneralizedPareto(mu, sigma,
894  * xi) distribution. Ill-conditioned for p near 1; conditon number is
895  *
896  * -xi/(1 - p^{-xi})
897  */
898 STATIC double
899 isf_genpareto(double p, double mu, double sigma, double xi)
900 {
901  if (fabs(xi) <= 1e-20)
902  return mu - sigma*log(p);
903  else
904  return mu + sigma*expm1(-xi*log(p))/xi;
905 }
906 
907 /*******************************************************************/
908 
909 /**
910  * Deterministic samplers, parametrized by uniform integer and (0,1]
911  * samples. No guarantees are made about _which_ mapping from the
912  * integer and (0,1] samples these use; all that is guaranteed is the
913  * distribution of the outputs conditioned on a uniform distribution on
914  * the inputs. The automatic tests in test_prob_distr.c double-check
915  * the particular mappings we use.
916  *
917  * Beware: Unlike random_uniform_01(), these are not guaranteed to be
918  * supported on all possible outputs. See Ilya Mironov, `On the
919  * Significance of the Least Significant Bits for Differential
920  * Privacy', for an example of what can go wrong if you try to use
921  * these to conceal information from an adversary but you expose the
922  * specific full-precision floating-point values.
923  *
924  * Note: None of these samplers use rejection sampling; they are all
925  * essentially inverse-CDF transforms with tweaks. If you were to add,
926  * say, a Gamma sampler with the Marsaglia-Tsang method, you would have
927  * to parametrize it by a potentially infinite stream of uniform (and
928  * perhaps normal) samples rather than a fixed number, which doesn't
929  * make for quite as nice automatic testing as for these.
930  */
931 
932 /**
933  * Deterministically sample from the interval [a, b], indexed by a
934  * uniform random floating-point number p0 in (0, 1].
935  *
936  * Note that even if p0 is nonzero, the result may be equal to a, if
937  * ulp(a)/2 is nonnegligible, e.g. if a = 1. For maximum resolution,
938  * arrange |a| <= |b|.
939  */
940 STATIC double
941 sample_uniform_interval(double p0, double a, double b)
942 {
943  /*
944  * XXX Prove that the distribution is, in fact, uniform on
945  * [a,b], particularly around p0 = 1, or at least has very
946  * small deviation from uniform, quantified appropriately
947  * (e.g., like in Monahan 1984, or by KL divergence). It
948  * almost certainly does but it would be nice to quantify the
949  * error.
950  */
951  if ((a <= 0 && 0 <= b) || (b <= 0 && 0 <= a)) {
952  /*
953  * When ab < 0, (1 - t) a + t b is monotonic, since for
954  * a <= b it is a sum of nondecreasing functions of t,
955  * and for b <= a, of nonincreasing functions of t.
956  * Further, clearly at 0 and 1 it attains a and b,
957  * respectively. Hence it is bounded within [a, b].
958  */
959  return (1 - p0)*a + p0*b;
960  } else {
961  /*
962  * a + (b - a) t is monotonic -- it is obviously a
963  * nondecreasing function of t for a <= b. Further, it
964  * attains a at 0, and while it may overshoot b at 1,
965  * we have a
966  *
967  * Theorem. If 0 <= t < 1, then the floating-point
968  * evaluation of a + (b - a) t is bounded in [a, b].
969  *
970  * Lemma 1. If 0 <= t < 1 is a floating-point number,
971  * then for any normal floating-point number x except
972  * the smallest in magnitude, |round(x*t)| < |x|.
973  *
974  * Proof. WLOG, assume x >= 0. Since the rounding
975  * function and t |---> x*t are nondecreasing, their
976  * composition t |---> round(x*t) is also
977  * nondecreasing, so it suffices to consider the
978  * largest floating-point number below 1, in particular
979  * t = 1 - ulp(1)/2.
980  *
981  * Case I: If x is a power of two, then the next
982  * floating-point number below x is x - ulp(x)/2 = x -
983  * x*ulp(1)/2 = x*(1 - ulp(1)/2) = x*t, so, since x*t
984  * is a floating-point number, multiplication is exact,
985  * and thus round(x*t) = x*t < x.
986  *
987  * Case II: If x is not a power of two, then the
988  * greatest lower bound of real numbers rounded to x is
989  * x - ulp(x)/2 = x - ulp(T(x))/2 = x - T(x)*ulp(1)/2,
990  * where T(X) is the largest power of two below x.
991  * Anything below this bound is rounded to a
992  * floating-point number smaller than x, and x*t = x*(1
993  * - ulp(1)/2) = x - x*ulp(1)/2 < x - T(x)*ulp(1)/2
994  * since T(x) < x, so round(x*t) < x*t < x. QED.
995  *
996  * Lemma 2. If x and y are subnormal, then round(x +
997  * y) = x + y.
998  *
999  * Proof. It is a matter of adding the significands,
1000  * since if we treat subnormals as having an implicit
1001  * zero bit before the `binary' point, their exponents
1002  * are all the same. There is at most one carry/borrow
1003  * bit, which can always be acommodated either in a
1004  * subnormal, or, at largest, in the implicit one bit
1005  * of a normal.
1006  *
1007  * Lemma 3. Let x and y be floating-point numbers. If
1008  * round(x - y) is subnormal or zero, then it is equal
1009  * to x - y.
1010  *
1011  * Proof. Case I (equal): round(x - y) = 0 iff x = y;
1012  * hence if round(x - y) = 0, then round(x - y) = 0 = x
1013  * - y.
1014  *
1015  * Case II (subnormal/subnormal): If x and y are both
1016  * subnormal, this follows directly from Lemma 2.
1017  *
1018  * Case IIIa (normal/subnormal): If x is normal and y
1019  * is subnormal, then x and y must share sign, or else
1020  * x - y would be larger than x and thus rounded to
1021  * normal. If s is the smallest normal positive
1022  * floating-point number, |x| < 2s since by
1023  * construction 2s - |y| is normal for all subnormal y.
1024  * This means that x and y must have the same exponent,
1025  * so the difference is the difference of significands,
1026  * which is exact.
1027  *
1028  * Case IIIb (subnormal/normal): Same as case IIIa for
1029  * -(y - x).
1030  *
1031  * Case IV (normal/normal): If x and y are both normal,
1032  * then they must share sign, or else x - y would be
1033  * larger than x and thus rounded to normal. Note that
1034  * |y| < 2|x|, for if |y| >= 2|x|, then |x| - |y| <=
1035  * -|x| but -|x| is normal like x. Also, |x|/2 < |y|:
1036  * if |x|/2 is subnormal, it must hold because y is
1037  * normal; if |x|/2 is normal, then |x|/2 >= s, so
1038  * since |x| - |y| < s,
1039  *
1040  * |x|/2 = |x| - |x|/2 <= |x| - s <= |y|;
1041  *
1042  * that is, |x|/2 < |y| < 2|x|, so by the Sterbenz
1043  * lemma, round(x - y) = x - y. QED.
1044  *
1045  * Proof of theorem. WLOG, assume 0 <= a <= b. Since
1046  * round(a + round(round(b - a)*t) is nondecreasing in
1047  * t and attains a at 0, the lower end of the bound is
1048  * trivial; we must show the upper end of the bound
1049  * strictly. It suffices to show this for the largest
1050  * floating-point number below 1, namely 1 - ulp(1)/2.
1051  *
1052  * Case I: round(b - a) is normal. Then it is at most
1053  * the smallest floating-point number above b - a. By
1054  * Lemma 1, round(round(b - a)*t) < round(b - a).
1055  * Since the inequality is strict, and since
1056  * round(round(b - a)*t) is a floating-point number
1057  * below round(b - a), and since there are no
1058  * floating-point numbers between b - a and round(b -
1059  * a), we must have round(round(b - a)*t) < b - a.
1060  * Then since y |---> round(a + y) is nondecreasing, we
1061  * must have
1062  *
1063  * round(a + round(round(b - a)*t))
1064  * <= round(a + (b - a))
1065  * = round(b) = b.
1066  *
1067  * Case II: round(b - a) is subnormal. In this case,
1068  * Lemma 1 falls apart -- we are not guaranteed the
1069  * strict inequality. However, by Lemma 3, the
1070  * difference is exact: round(b - a) = b - a. Thus,
1071  *
1072  * round(a + round(round(b - a)*t))
1073  * <= round(a + round((b - a)*t))
1074  * <= round(a + (b - a))
1075  * = round(b)
1076  * = b,
1077  *
1078  * QED.
1079  */
1080 
1081  /* p0 is restricted to [0,1], but we use >= to silence -Wfloat-equal. */
1082  if (p0 >= 1)
1083  return b;
1084  return a + (b - a)*p0;
1085  }
1086 }
1087 
1088 /**
1089  * Deterministically sample from the standard logistic distribution,
1090  * indexed by a uniform random 32-bit integer s and uniform random
1091  * floating-point numbers t and p0 in (0, 1].
1092  */
1093 STATIC double
1094 sample_logistic(uint32_t s, double t, double p0)
1095 {
1096  double sign = (s & 1) ? -1 : +1;
1097  double r;
1098 
1099  /*
1100  * We carve up the interval (0, 1) into subregions to compute
1101  * the inverse CDF precisely:
1102  *
1103  * A = (0, 1/(1 + e)] ---> (-infinity, -1]
1104  * B = [1/(1 + e), 1/2] ---> [-1, 0]
1105  * C = [1/2, 1 - 1/(1 + e)] ---> [0, 1]
1106  * D = [1 - 1/(1 + e), 1) ---> [1, +infinity)
1107  *
1108  * Cases D and C are mirror images of cases A and B,
1109  * respectively, so we choose between them by the sign chosen
1110  * by a fair coin toss. We choose between cases A and B by a
1111  * coin toss weighted by
1112  *
1113  * 2/(1 + e) = 1 - [1/2 - 1/(1 + e)]/(1/2):
1114  *
1115  * if it comes up heads, scale p0 into a uniform (0, 1/(1 + e)]
1116  * sample p; if it comes up tails, scale p0 into a uniform (0,
1117  * 1/2 - 1/(1 + e)] sample and compute the inverse CDF of p =
1118  * 1/2 - p0.
1119  */
1120  if (t <= 2/(1 + exp(1))) {
1121  /* p uniform in (0, 1/(1 + e)], represented by p. */
1122  p0 /= 1 + exp(1);
1123  r = logit(p0);
1124  } else {
1125  /*
1126  * p uniform in [1/(1 + e), 1/2), actually represented
1127  * by p0 = 1/2 - p uniform in (0, 1/2 - 1/(1 + e)], so
1128  * that p = 1/2 - p.
1129  */
1130  p0 *= 0.5 - 1/(1 + exp(1));
1131  r = logithalf(p0);
1132  }
1133 
1134  /*
1135  * We have chosen from the negative half of the standard
1136  * logistic distribution, which is symmetric with the positive
1137  * half. Now use the sign to choose uniformly between them.
1138  */
1139  return sign*r;
1140 }
1141 
1142 /**
1143  * Deterministically sample from the logistic distribution scaled by
1144  * sigma and translated by mu.
1145  */
1146 static double
1147 sample_logistic_locscale(uint32_t s, double t, double p0, double mu,
1148  double sigma)
1149 {
1150 
1151  return mu + sigma*sample_logistic(s, t, p0);
1152 }
1153 
1154 /**
1155  * Deterministically sample from the standard log-logistic
1156  * distribution, indexed by a uniform random 32-bit integer s and a
1157  * uniform random floating-point number p0 in (0, 1].
1158  */
1159 STATIC double
1160 sample_log_logistic(uint32_t s, double p0)
1161 {
1162 
1163  /*
1164  * Carve up the interval (0, 1) into (0, 1/2] and [1/2, 1); the
1165  * condition numbers of the icdf and the isf coincide at 1/2.
1166  */
1167  p0 *= 0.5;
1168  if ((s & 1) == 0) {
1169  /* p = p0 in (0, 1/2] */
1170  return p0/(1 - p0);
1171  } else {
1172  /* p = 1 - p0 in [1/2, 1) */
1173  return (1 - p0)/p0;
1174  }
1175 }
1176 
1177 /**
1178  * Deterministically sample from the log-logistic distribution with
1179  * scale alpha and shape beta.
1180  */
1181 static double
1182 sample_log_logistic_scaleshape(uint32_t s, double p0, double alpha,
1183  double beta)
1184 {
1185  double x = sample_log_logistic(s, p0);
1186 
1187  return alpha*pow(x, 1/beta);
1188 }
1189 
1190 /**
1191  * Deterministically sample from the standard exponential distribution,
1192  * indexed by a uniform random 32-bit integer s and a uniform random
1193  * floating-point number p0 in (0, 1].
1194  */
1195 static double
1196 sample_exponential(uint32_t s, double p0)
1197 {
1198  /*
1199  * We would like to evaluate log(p) for p near 0, and log1p(-p)
1200  * for p near 1. Simply carve the interval into (0, 1/2] and
1201  * [1/2, 1) by a fair coin toss.
1202  */
1203  p0 *= 0.5;
1204  if ((s & 1) == 0)
1205  /* p = p0 in (0, 1/2] */
1206  return -log(p0);
1207  else
1208  /* p = 1 - p0 in [1/2, 1) */
1209  return -log1p(-p0);
1210 }
1211 
1212 /**
1213  * Deterministically sample from a Weibull distribution with scale
1214  * lambda and shape k -- just an exponential with a shape parameter in
1215  * addition to a scale parameter. (Yes, lambda really is the scale,
1216  * _not_ the rate.)
1217  */
1218 STATIC double
1219 sample_weibull(uint32_t s, double p0, double lambda, double k)
1220 {
1221 
1222  return lambda*pow(sample_exponential(s, p0), 1/k);
1223 }
1224 
1225 /**
1226  * Deterministically sample from the generalized Pareto distribution
1227  * with shape xi, indexed by a uniform random 32-bit integer s and a
1228  * uniform random floating-point number p0 in (0, 1].
1229  */
1230 STATIC double
1231 sample_genpareto(uint32_t s, double p0, double xi)
1232 {
1233  double x = sample_exponential(s, p0);
1234 
1235  /*
1236  * Write f(xi) = (e^{xi x} - 1)/xi for xi near zero as the
1237  * absolutely convergent Taylor series
1238  *
1239  * f(x) = (1/xi) (xi x + \sum_{n=2}^infinity (xi x)^n/n!)
1240  * = x + (1/xi) \sum_{n=2}^infinity (xi x)^n/n!
1241  * = x + \sum_{n=2}^infinity xi^{n-1} x^n/n!
1242  * = x + x \sum_{n=2}^infinity (xi x)^{n-1}/n!
1243  * = x (1 + \sum_{n=2}^infinity (xi x)^{n-1}/n!).
1244  *
1245  * d = \sum_{n=2}^infinity (xi x)^{n-1}/n! is the relative error
1246  * of f(x) from x. If |xi| < eps/4x, then
1247  *
1248  * |d| <= \sum_{n=2}^infinity |xi x|^{n-1}/n!
1249  * <= \sum_{n=2}^infinity (eps/4)^{n-1}/n!
1250  * <= \sum_{n=1}^infinity (eps/4)
1251  * = (eps/4) \sum_{n=0}^infinity (eps/4)^n
1252  * = (eps/4)/(1 - eps/4)
1253  * < eps/2,
1254  *
1255  * for any 0 < eps < 2. Hence, as long as |xi| < eps/2x, f(xi)
1256  * = x (1 + d) for |d| <= eps/2, so x = f(xi) (1 + d') for |d'|
1257  * <= eps. What bound should we use for x?
1258  *
1259  * - If x is exponentially distributed, x > 200 with
1260  * probability below e^{-200} << 2^{-256}, i.e. never.
1261  *
1262  * - If x is computed by -log(U) for U in (0, 1], x is
1263  * guaranteed to be below 1000 in IEEE 754 binary64
1264  * floating-point.
1265  *
1266  * We can safely cut xi off at 1e-20 < eps/4000 and attain an
1267  * error bounded by 0.5 ulp for this expression.
1268  */
1269  return (fabs(xi) < 1e-20 ? x : expm1(xi*x)/xi);
1270 }
1271 
1272 /**
1273  * Deterministically sample from a generalized Pareto distribution with
1274  * shape xi, scaled by sigma and translated by mu.
1275  */
1276 static double
1277 sample_genpareto_locscale(uint32_t s, double p0, double mu, double sigma,
1278  double xi)
1279 {
1280 
1281  return mu + sigma*sample_genpareto(s, p0, xi);
1282 }
1283 
1284 /**
1285  * Deterministically sample from the geometric distribution with
1286  * per-trial success probability p.
1287  *
1288  * XXX Quantify the error (KL divergence?) of this
1289  * ceiling-of-exponential sampler from a true geometric distribution,
1290  * which we could get by rejection sampling. Relevant papers:
1291  *
1292  * John F. Monahan, `Accuracy in Random Number Generation',
1293  * Mathematics of Computation 45(172), October 1984, pp. 559--568.
1294 *https://pdfs.semanticscholar.org/aca6/74b96da1df77b2224e8cfc5dd6d61a471632.pdf
1295  *
1296  * Karl Bringmann and Tobias Friedrich, `Exact and Efficient
1297  * Generation of Geometric Random Variates and Random Graphs', in
1298  * Proceedings of the 40th International Colloaquium on Automata,
1299  * Languages, and Programming -- ICALP 2013, Springer LNCS 7965,
1300  * pp.267--278.
1301  * https://doi.org/10.1007/978-3-642-39206-1_23
1302  * https://people.mpi-inf.mpg.de/~kbringma/paper/2013ICALP-1.pdf
1303  */
1304 static double
1305 sample_geometric(uint32_t s, double p0, double p)
1306 {
1307  double x = sample_exponential(s, p0);
1308 
1309  /* This is actually a check against 1, but we do >= so that the compiler
1310  does not raise a -Wfloat-equal */
1311  if (p >= 1)
1312  return 1;
1313 
1314  return ceil(-x/log1p(-p));
1315 }
1316 
1317 /*******************************************************************/
1318 
1319 /** Public API for probability distributions:
1320  *
1321  * These are wrapper functions on top of the various probability distribution
1322  * operations using the generic <b>dist</b> structure.
1323 
1324  * These are the functions that should be used by consumers of this API.
1325  */
1326 
1327 /** Returns the name of the distribution in <b>dist</b>. */
1328 const char *
1329 dist_name(const struct dist_t *dist)
1330 {
1331  return dist->ops->name;
1332 }
1333 
1334 /* Sample a value from <b>dist</b> and return it. */
1335 double
1336 dist_sample(const struct dist_t *dist)
1337 {
1338  return dist->ops->sample(dist);
1339 }
1340 
1341 /** Compute the CDF of <b>dist</b> at <b>x</b>. */
1342 double
1343 dist_cdf(const struct dist_t *dist, double x)
1344 {
1345  return dist->ops->cdf(dist, x);
1346 }
1347 
1348 /** Compute the SF (Survival function) of <b>dist</b> at <b>x</b>. */
1349 double
1350 dist_sf(const struct dist_t *dist, double x)
1351 {
1352  return dist->ops->sf(dist, x);
1353 }
1354 
1355 /** Compute the iCDF (Inverse CDF) of <b>dist</b> at <b>x</b>. */
1356 double
1357 dist_icdf(const struct dist_t *dist, double p)
1358 {
1359  return dist->ops->icdf(dist, p);
1360 }
1361 
1362 /** Compute the iSF (Inverse Survival function) of <b>dist</b> at <b>x</b>. */
1363 double
1364 dist_isf(const struct dist_t *dist, double p)
1365 {
1366  return dist->ops->isf(dist, p);
1367 }
1368 
1369 /** Functions for uniform distribution */
1370 
1371 static double
1372 uniform_sample(const struct dist_t *dist)
1373 {
1374  const struct uniform_t *U = dist_to_const_uniform(dist);
1375  double p0 = random_uniform_01();
1376 
1377  return sample_uniform_interval(p0, U->a, U->b);
1378 }
1379 
1380 static double
1381 uniform_cdf(const struct dist_t *dist, double x)
1382 {
1383  const struct uniform_t *U = dist_to_const_uniform(dist);
1384  if (x < U->a)
1385  return 0;
1386  else if (x < U->b)
1387  return (x - U->a)/(U->b - U->a);
1388  else
1389  return 1;
1390 }
1391 
1392 static double
1393 uniform_sf(const struct dist_t *dist, double x)
1394 {
1395  const struct uniform_t *U = dist_to_const_uniform(dist);
1396 
1397  if (x > U->b)
1398  return 0;
1399  else if (x > U->a)
1400  return (U->b - x)/(U->b - U->a);
1401  else
1402  return 1;
1403 }
1404 
1405 static double
1406 uniform_icdf(const struct dist_t *dist, double p)
1407 {
1408  const struct uniform_t *U = dist_to_const_uniform(dist);
1409  double w = U->b - U->a;
1410 
1411  return (p < 0.5 ? (U->a + w*p) : (U->b - w*(1 - p)));
1412 }
1413 
1414 static double
1415 uniform_isf(const struct dist_t *dist, double p)
1416 {
1417  const struct uniform_t *U = dist_to_const_uniform(dist);
1418  double w = U->b - U->a;
1419 
1420  return (p < 0.5 ? (U->b - w*p) : (U->a + w*(1 - p)));
1421 }
1422 
1423 const struct dist_ops_t uniform_ops = {
1424  .name = "uniform",
1425  .sample = uniform_sample,
1426  .cdf = uniform_cdf,
1427  .sf = uniform_sf,
1428  .icdf = uniform_icdf,
1429  .isf = uniform_isf,
1430 };
1431 
1432 /*******************************************************************/
1433 
1434 /** Private functions for each probability distribution. */
1435 
1436 /** Functions for logistic distribution: */
1437 
1438 static double
1439 logistic_sample(const struct dist_t *dist)
1440 {
1441  const struct logistic_t *L = dist_to_const_logistic(dist);
1443  double t = random_uniform_01();
1444  double p0 = random_uniform_01();
1445 
1446  return sample_logistic_locscale(s, t, p0, L->mu, L->sigma);
1447 }
1448 
1449 static double
1450 logistic_cdf(const struct dist_t *dist, double x)
1451 {
1452  const struct logistic_t *L = dist_to_const_logistic(dist);
1453  return cdf_logistic(x, L->mu, L->sigma);
1454 }
1455 
1456 static double
1457 logistic_sf(const struct dist_t *dist, double x)
1458 {
1459  const struct logistic_t *L = dist_to_const_logistic(dist);
1460  return sf_logistic(x, L->mu, L->sigma);
1461 }
1462 
1463 static double
1464 logistic_icdf(const struct dist_t *dist, double p)
1465 {
1466  const struct logistic_t *L = dist_to_const_logistic(dist);
1467  return icdf_logistic(p, L->mu, L->sigma);
1468 }
1469 
1470 static double
1471 logistic_isf(const struct dist_t *dist, double p)
1472 {
1473  const struct logistic_t *L = dist_to_const_logistic(dist);
1474  return isf_logistic(p, L->mu, L->sigma);
1475 }
1476 
1477 const struct dist_ops_t logistic_ops = {
1478  .name = "logistic",
1479  .sample = logistic_sample,
1480  .cdf = logistic_cdf,
1481  .sf = logistic_sf,
1482  .icdf = logistic_icdf,
1483  .isf = logistic_isf,
1484 };
1485 
1486 /** Functions for log-logistic distribution: */
1487 
1488 static double
1489 log_logistic_sample(const struct dist_t *dist)
1490 {
1491  const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1493  double p0 = random_uniform_01();
1494 
1495  return sample_log_logistic_scaleshape(s, p0, LL->alpha, LL->beta);
1496 }
1497 
1498 static double
1499 log_logistic_cdf(const struct dist_t *dist, double x)
1500 {
1501  const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1502  return cdf_log_logistic(x, LL->alpha, LL->beta);
1503 }
1504 
1505 static double
1506 log_logistic_sf(const struct dist_t *dist, double x)
1507 {
1508  const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1509  return sf_log_logistic(x, LL->alpha, LL->beta);
1510 }
1511 
1512 static double
1513 log_logistic_icdf(const struct dist_t *dist, double p)
1514 {
1515  const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1516  return icdf_log_logistic(p, LL->alpha, LL->beta);
1517 }
1518 
1519 static double
1520 log_logistic_isf(const struct dist_t *dist, double p)
1521 {
1522  const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1523  return isf_log_logistic(p, LL->alpha, LL->beta);
1524 }
1525 
1526 const struct dist_ops_t log_logistic_ops = {
1527  .name = "log logistic",
1528  .sample = log_logistic_sample,
1529  .cdf = log_logistic_cdf,
1530  .sf = log_logistic_sf,
1531  .icdf = log_logistic_icdf,
1532  .isf = log_logistic_isf,
1533 };
1534 
1535 /** Functions for Weibull distribution */
1536 
1537 static double
1538 weibull_sample(const struct dist_t *dist)
1539 {
1540  const struct weibull_t *W = dist_to_const_weibull(dist);
1542  double p0 = random_uniform_01();
1543 
1544  return sample_weibull(s, p0, W->lambda, W->k);
1545 }
1546 
1547 static double
1548 weibull_cdf(const struct dist_t *dist, double x)
1549 {
1550  const struct weibull_t *W = dist_to_const_weibull(dist);
1551  return cdf_weibull(x, W->lambda, W->k);
1552 }
1553 
1554 static double
1555 weibull_sf(const struct dist_t *dist, double x)
1556 {
1557  const struct weibull_t *W = dist_to_const_weibull(dist);
1558  return sf_weibull(x, W->lambda, W->k);
1559 }
1560 
1561 static double
1562 weibull_icdf(const struct dist_t *dist, double p)
1563 {
1564  const struct weibull_t *W = dist_to_const_weibull(dist);
1565  return icdf_weibull(p, W->lambda, W->k);
1566 }
1567 
1568 static double
1569 weibull_isf(const struct dist_t *dist, double p)
1570 {
1571  const struct weibull_t *W = dist_to_const_weibull(dist);
1572  return isf_weibull(p, W->lambda, W->k);
1573 }
1574 
1575 const struct dist_ops_t weibull_ops = {
1576  .name = "Weibull",
1577  .sample = weibull_sample,
1578  .cdf = weibull_cdf,
1579  .sf = weibull_sf,
1580  .icdf = weibull_icdf,
1581  .isf = weibull_isf,
1582 };
1583 
1584 /** Functions for generalized Pareto distributions */
1585 
1586 static double
1587 genpareto_sample(const struct dist_t *dist)
1588 {
1589  const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1591  double p0 = random_uniform_01();
1592 
1593  return sample_genpareto_locscale(s, p0, GP->mu, GP->sigma, GP->xi);
1594 }
1595 
1596 static double
1597 genpareto_cdf(const struct dist_t *dist, double x)
1598 {
1599  const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1600  return cdf_genpareto(x, GP->mu, GP->sigma, GP->xi);
1601 }
1602 
1603 static double
1604 genpareto_sf(const struct dist_t *dist, double x)
1605 {
1606  const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1607  return sf_genpareto(x, GP->mu, GP->sigma, GP->xi);
1608 }
1609 
1610 static double
1611 genpareto_icdf(const struct dist_t *dist, double p)
1612 {
1613  const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1614  return icdf_genpareto(p, GP->mu, GP->sigma, GP->xi);
1615 }
1616 
1617 static double
1618 genpareto_isf(const struct dist_t *dist, double p)
1619 {
1620  const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1621  return isf_genpareto(p, GP->mu, GP->sigma, GP->xi);
1622 }
1623 
1624 const struct dist_ops_t genpareto_ops = {
1625  .name = "generalized Pareto",
1626  .sample = genpareto_sample,
1627  .cdf = genpareto_cdf,
1628  .sf = genpareto_sf,
1629  .icdf = genpareto_icdf,
1630  .isf = genpareto_isf,
1631 };
1632 
1633 /** Functions for geometric distribution on number of trials before success */
1634 
1635 static double
1636 geometric_sample(const struct dist_t *dist)
1637 {
1638  const struct geometric_t *G = dist_to_const_geometric(dist);
1640  double p0 = random_uniform_01();
1641 
1642  return sample_geometric(s, p0, G->p);
1643 }
1644 
1645 static double
1646 geometric_cdf(const struct dist_t *dist, double x)
1647 {
1648  const struct geometric_t *G = dist_to_const_geometric(dist);
1649 
1650  if (x < 1)
1651  return 0;
1652  /* 1 - (1 - p)^floor(x) = 1 - e^{floor(x) log(1 - p)} */
1653  return -expm1(floor(x)*log1p(-G->p));
1654 }
1655 
1656 static double
1657 geometric_sf(const struct dist_t *dist, double x)
1658 {
1659  const struct geometric_t *G = dist_to_const_geometric(dist);
1660 
1661  if (x < 1)
1662  return 0;
1663  /* (1 - p)^floor(x) = e^{ceil(x) log(1 - p)} */
1664  return exp(floor(x)*log1p(-G->p));
1665 }
1666 
1667 static double
1668 geometric_icdf(const struct dist_t *dist, double p)
1669 {
1670  const struct geometric_t *G = dist_to_const_geometric(dist);
1671 
1672  return log1p(-p)/log1p(-G->p);
1673 }
1674 
1675 static double
1676 geometric_isf(const struct dist_t *dist, double p)
1677 {
1678  const struct geometric_t *G = dist_to_const_geometric(dist);
1679 
1680  return log(p)/log1p(-G->p);
1681 }
1682 
1683 const struct dist_ops_t geometric_ops = {
1684  .name = "geometric (1-based)",
1685  .sample = geometric_sample,
1686  .cdf = geometric_cdf,
1687  .sf = geometric_sf,
1688  .icdf = geometric_icdf,
1689  .isf = geometric_isf,
1690 };
STATIC double logithalf(double p0)
Definition: prob_distr.c:360
Common functions for using (pseudo-)random number generators.
STATIC double sf_logistic(double x, double mu, double sigma)
Definition: prob_distr.c:523
STATIC double cdf_genpareto(double x, double mu, double sigma, double xi)
Definition: prob_distr.c:792
double dist_icdf(const struct dist_t *dist, double p)
Definition: prob_distr.c:1357
STATIC double isf_genpareto(double p, double mu, double sigma, double xi)
Definition: prob_distr.c:899
#define CTASSERT(x)
Definition: ctassert.h:44
static double sample_exponential(uint32_t s, double p0)
Definition: prob_distr.c:1196
static double sample_logistic_locscale(uint32_t s, double t, double p0, double mu, double sigma)
Definition: prob_distr.c:1147
STATIC double isf_logistic(double p, double mu, double sigma)
Definition: prob_distr.c:545
static unsigned bitcount32(uint32_t x)
Definition: prob_distr.c:77
STATIC double icdf_genpareto(double p, double mu, double sigma, double xi)
Definition: prob_distr.c:853
STATIC double logit(double p)
Definition: prob_distr.c:280
STATIC double logistic(double x)
Definition: prob_distr.c:193
STATIC double sf_weibull(double x, double lambda, double k)
Definition: prob_distr.c:726
STATIC double sample_log_logistic(uint32_t s, double p0)
Definition: prob_distr.c:1160
STATIC double icdf_log_logistic(double p, double alpha, double beta)
Definition: prob_distr.c:675
STATIC double sf_log_logistic(double x, double alpha, double beta)
Definition: prob_distr.c:659
#define STATIC
Definition: testsupport.h:32
STATIC double cdf_weibull(double x, double lambda, double k)
Definition: prob_distr.c:711
static double log_logistic_sample(const struct dist_t *dist)
Definition: prob_distr.c:1489
STATIC double cdf_logistic(double x, double mu, double sigma)
Definition: prob_distr.c:511
STATIC double sf_genpareto(double x, double mu, double sigma, double xi)
Definition: prob_distr.c:836
STATIC double cdf_log_logistic(double x, double alpha, double beta)
Definition: prob_distr.c:598
double dist_sf(const struct dist_t *dist, double x)
Definition: prob_distr.c:1350
static double sample_geometric(uint32_t s, double p0, double p)
Definition: prob_distr.c:1305
STATIC double sample_uniform_interval(double p0, double a, double b)
Definition: prob_distr.c:941
static unsigned clz32(uint32_t x)
Definition: prob_distr.c:97
Header for prob_distr.c.
static double uniform_sample(const struct dist_t *dist)
Definition: prob_distr.c:1372
STATIC double icdf_logistic(double p, double mu, double sigma)
Definition: prob_distr.c:534
#define DECLARE_PROB_DISTR_DOWNCAST_FN(name)
Definition: prob_distr.c:58
STATIC double isf_weibull(double p, double lambda, double k)
Definition: prob_distr.c:752
Compile-time assertions: CTASSERT(expression).
double dist_isf(const struct dist_t *dist, double p)
Definition: prob_distr.c:1364
crypto_fast_rng_t * get_thread_fast_rng(void)
STATIC double sample_genpareto(uint32_t s, double p0, double xi)
Definition: prob_distr.c:1231
static double sample_log_logistic_scaleshape(uint32_t s, double p0, double alpha, double beta)
Definition: prob_distr.c:1182
STATIC double random_uniform_01(void)
Definition: prob_distr.c:453
static double weibull_sample(const struct dist_t *dist)
Definition: prob_distr.c:1538
static double sample_genpareto_locscale(uint32_t s, double p0, double mu, double sigma, double xi)
Definition: prob_distr.c:1277
STATIC double sample_logistic(uint32_t s, double t, double p0)
Definition: prob_distr.c:1094
const char * dist_name(const struct dist_t *dist)
Definition: prob_distr.c:1329
double dist_cdf(const struct dist_t *dist, double x)
Definition: prob_distr.c:1343
STATIC double isf_log_logistic(double p, double alpha, double beta)
Definition: prob_distr.c:686
static double logistic_sample(const struct dist_t *dist)
Definition: prob_distr.c:1439
STATIC double sample_weibull(uint32_t s, double p0, double lambda, double k)
Definition: prob_distr.c:1219
Macros to manage assertions, fatal and non-fatal.
STATIC double icdf_weibull(double p, double lambda, double k)
Definition: prob_distr.c:739
uint32_t crypto_fast_rng_get_u32(crypto_fast_rng_t *rng)
static double genpareto_sample(const struct dist_t *dist)
Definition: prob_distr.c:1587
static double geometric_sample(const struct dist_t *dist)
Definition: prob_distr.c:1636