41 #define PROB_DISTR_PRIVATE
58 #define DECLARE_PROB_DISTR_DOWNCAST_FN(name) \
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); \
81 x -= (x >> 1) & UINT32_C(0x55555555);
84 x = ((x >> 2) & UINT32_C(0x33333333)) + (x & UINT32_C(0x33333333));
87 x = (x + (x >> 4)) & UINT32_C(0x0f0f0f0f);
90 return (x * UINT32_C(0x01010101)) >> 24;
195 if (x <= log(DBL_EPSILON/2)) {
210 }
else if (x <= -log(DBL_EPSILON/2)) {
247 return 1/(1 + exp(-x));
284 if (1/(1 + exp(1)) <= p && p <= 1/(1 + exp(-1))) {
314 return -log1p((1 - 2*p)/p);
347 return log(p/(1 - p));
363 if (fabs(p0) <= 0.5 - 1/(1 + exp(1))) {
395 return log1p(2*p0/(0.5 - p0));
430 return log((0.5 + p0)/(0.5 - p0));
455 uint32_t z, x, hi, lo;
482 s = hi*(double)4294967296 + lo;
485 return s * ldexp(1, -(64 + z));
536 return mu + sigma*
logit(p);
547 return mu - sigma*
logit(p);
640 return 1/(1 + pow(x/alpha, -beta));
666 return 1/(1 + pow(x/alpha, beta));
677 return alpha*pow(p/(1 - p), 1/beta);
688 return alpha*pow((1 - p)/p, 1/beta);
713 return -expm1(-pow(x/lambda, k));
728 return exp(-pow(x/lambda, k));
741 return lambda*pow(-log1p(-p), 1/k);
754 return lambda*pow(-log(p), 1/k);
794 double x_0 = (x - mu)/sigma;
817 if (fabs(xi) < 1e-17/x_0)
820 return -expm1(-log1p(xi*x_0)/xi);
838 double x_0 = (x - mu)/sigma;
840 if (fabs(xi) < 1e-17/x_0)
843 return exp(-log1p(xi*x_0)/xi);
886 if (fabs(xi) <= 1e-20)
887 return mu - sigma*log1p(-p);
889 return mu + sigma*expm1(-xi*log1p(-p))/xi;
901 if (fabs(xi) <= 1e-20)
902 return mu - sigma*log(p);
904 return mu + sigma*expm1(-xi*log(p))/xi;
951 if ((a <= 0 && 0 <= b) || (b <= 0 && 0 <= a)) {
959 return (1 - p0)*a + p0*b;
1084 return a + (b - a)*p0;
1096 double sign = (s & 1) ? -1 : +1;
1120 if (t <= 2/(1 + exp(1))) {
1130 p0 *= 0.5 - 1/(1 + exp(1));
1187 return alpha*pow(x, 1/beta);
1269 return (fabs(xi) < 1e-20 ? x : expm1(xi*x)/xi);
1316 return ceil(-x/log1p(-p));
1333 return dist->ops->name;
1338 dist_sample(
const struct dist_t *dist)
1340 return dist->ops->sample(dist);
1347 return dist->ops->cdf(dist, x);
1354 return dist->ops->sf(dist, x);
1361 return dist->ops->icdf(dist, p);
1368 return dist->ops->isf(dist, p);
1376 const struct uniform_t *U = dist_to_const_uniform(dist);
1383 uniform_cdf(
const struct dist_t *dist,
double x)
1385 const struct uniform_t *U = dist_to_const_uniform(dist);
1389 return (x - U->a)/(U->b - U->a);
1395 uniform_sf(
const struct dist_t *dist,
double x)
1397 const struct uniform_t *U = dist_to_const_uniform(dist);
1402 return (U->b - x)/(U->b - U->a);
1408 uniform_icdf(
const struct dist_t *dist,
double p)
1410 const struct uniform_t *U = dist_to_const_uniform(dist);
1411 double w = U->b - U->a;
1413 return (p < 0.5 ? (U->a + w*p) : (U->b - w*(1 - p)));
1417 uniform_isf(
const struct dist_t *dist,
double p)
1419 const struct uniform_t *U = dist_to_const_uniform(dist);
1420 double w = U->b - U->a;
1422 return (p < 0.5 ? (U->b - w*p) : (U->a + w*(1 - p)));
1430 .icdf = uniform_icdf,
1443 const struct logistic_t *L = dist_to_const_logistic(dist);
1452 logistic_cdf(
const struct dist_t *dist,
double x)
1454 const struct logistic_t *L = dist_to_const_logistic(dist);
1459 logistic_sf(
const struct dist_t *dist,
double x)
1461 const struct logistic_t *L = dist_to_const_logistic(dist);
1466 logistic_icdf(
const struct dist_t *dist,
double p)
1468 const struct logistic_t *L = dist_to_const_logistic(dist);
1473 logistic_isf(
const struct dist_t *dist,
double p)
1475 const struct logistic_t *L = dist_to_const_logistic(dist);
1482 .cdf = logistic_cdf,
1484 .icdf = logistic_icdf,
1485 .isf = logistic_isf,
1493 const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1501 log_logistic_cdf(
const struct dist_t *dist,
double x)
1503 const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1508 log_logistic_sf(
const struct dist_t *dist,
double x)
1510 const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1515 log_logistic_icdf(
const struct dist_t *dist,
double p)
1517 const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1522 log_logistic_isf(
const struct dist_t *dist,
double p)
1524 const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1529 .name =
"log logistic",
1531 .cdf = log_logistic_cdf,
1532 .sf = log_logistic_sf,
1533 .icdf = log_logistic_icdf,
1534 .isf = log_logistic_isf,
1542 const struct weibull_t *W = dist_to_const_weibull(dist);
1550 weibull_cdf(
const struct dist_t *dist,
double x)
1552 const struct weibull_t *W = dist_to_const_weibull(dist);
1557 weibull_sf(
const struct dist_t *dist,
double x)
1559 const struct weibull_t *W = dist_to_const_weibull(dist);
1564 weibull_icdf(
const struct dist_t *dist,
double p)
1566 const struct weibull_t *W = dist_to_const_weibull(dist);
1571 weibull_isf(
const struct dist_t *dist,
double p)
1573 const struct weibull_t *W = dist_to_const_weibull(dist);
1582 .icdf = weibull_icdf,
1591 const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1599 genpareto_cdf(
const struct dist_t *dist,
double x)
1601 const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1606 genpareto_sf(
const struct dist_t *dist,
double x)
1608 const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1613 genpareto_icdf(
const struct dist_t *dist,
double p)
1615 const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1620 genpareto_isf(
const struct dist_t *dist,
double p)
1622 const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1627 .name =
"generalized Pareto",
1629 .cdf = genpareto_cdf,
1631 .icdf = genpareto_icdf,
1632 .isf = genpareto_isf,
1640 const struct geometric_t *G = dist_to_const_geometric(dist);
1648 geometric_cdf(
const struct dist_t *dist,
double x)
1650 const struct geometric_t *G = dist_to_const_geometric(dist);
1655 return -expm1(floor(x)*log1p(-G->p));
1659 geometric_sf(
const struct dist_t *dist,
double x)
1661 const struct geometric_t *G = dist_to_const_geometric(dist);
1666 return exp(floor(x)*log1p(-G->p));
1670 geometric_icdf(
const struct dist_t *dist,
double p)
1672 const struct geometric_t *G = dist_to_const_geometric(dist);
1674 return log1p(-p)/log1p(-G->p);
1678 geometric_isf(
const struct dist_t *dist,
double p)
1680 const struct geometric_t *G = dist_to_const_geometric(dist);
1682 return log(p)/log1p(-G->p);
1686 .name =
"geometric (1-based)",
1688 .cdf = geometric_cdf,
1690 .icdf = geometric_icdf,
1691 .isf = geometric_isf,
Common functions for using (pseudo-)random number generators.
crypto_fast_rng_t * get_thread_fast_rng(void)
uint32_t crypto_fast_rng_get_u32(crypto_fast_rng_t *rng)
Compile-time assertions: CTASSERT(expression).
STATIC double icdf_log_logistic(double p, double alpha, double beta)
static unsigned clz32(uint32_t x)
double dist_sf(const struct dist_t *dist, double x)
#define DECLARE_PROB_DISTR_DOWNCAST_FN(name)
STATIC double cdf_weibull(double x, double lambda, double k)
static double sample_logistic_locscale(uint32_t s, double t, double p0, double mu, double sigma)
STATIC double sample_genpareto(uint32_t s, double p0, double xi)
STATIC double random_uniform_01(void)
static double sample_geometric(uint32_t s, double p0, double p)
static double log_logistic_sample(const struct dist_t *dist)
double dist_isf(const struct dist_t *dist, double p)
static double genpareto_sample(const struct dist_t *dist)
STATIC double cdf_log_logistic(double x, double alpha, double beta)
STATIC double sf_weibull(double x, double lambda, double k)
static unsigned bitcount32(uint32_t x)
static double weibull_sample(const struct dist_t *dist)
STATIC double cdf_logistic(double x, double mu, double sigma)
STATIC double sf_log_logistic(double x, double alpha, double beta)
STATIC double isf_weibull(double p, double lambda, double k)
STATIC double isf_logistic(double p, double mu, double sigma)
double dist_icdf(const struct dist_t *dist, double p)
static double sample_genpareto_locscale(uint32_t s, double p0, double mu, double sigma, double xi)
STATIC double sf_logistic(double x, double mu, double sigma)
double dist_cdf(const struct dist_t *dist, double x)
STATIC double logithalf(double p0)
STATIC double isf_genpareto(double p, double mu, double sigma, double xi)
STATIC double icdf_weibull(double p, double lambda, double k)
static double sample_exponential(uint32_t s, double p0)
STATIC double sample_logistic(uint32_t s, double t, double p0)
STATIC double sample_weibull(uint32_t s, double p0, double lambda, double k)
STATIC double logistic(double x)
STATIC double sample_log_logistic(uint32_t s, double p0)
static double sample_log_logistic_scaleshape(uint32_t s, double p0, double alpha, double beta)
static double logistic_sample(const struct dist_t *dist)
STATIC double icdf_logistic(double p, double mu, double sigma)
STATIC double cdf_genpareto(double x, double mu, double sigma, double xi)
STATIC double icdf_genpareto(double p, double mu, double sigma, double xi)
STATIC double logit(double p)
static double uniform_sample(const struct dist_t *dist)
STATIC double sample_uniform_interval(double p0, double a, double b)
STATIC double isf_log_logistic(double p, double alpha, double beta)
static double geometric_sample(const struct dist_t *dist)
const char * dist_name(const struct dist_t *dist)
STATIC double sf_genpareto(double x, double mu, double sigma, double xi)
Macros to manage assertions, fatal and non-fatal.