diff options
Diffstat (limited to 'src/lib/math/prob_distr.c')
| -rw-r--r-- | src/lib/math/prob_distr.c | 165 |
1 files changed, 68 insertions, 97 deletions
diff --git a/src/lib/math/prob_distr.c b/src/lib/math/prob_distr.c index c952dadc06..d44dc28265 100644 --- a/src/lib/math/prob_distr.c +++ b/src/lib/math/prob_distr.c @@ -46,26 +46,27 @@ #include "lib/crypt_ops/crypto_rand.h" #include "lib/cc/ctassert.h" +#include "lib/log/util_bug.h" #include <float.h> #include <math.h> #include <stddef.h> -/** Validators for downcasting macros below */ -#define validate_container_of(PTR, TYPE, FIELD) \ - (0 * sizeof((PTR) - &((TYPE *)(((char *)(PTR)) - \ - offsetof(TYPE, FIELD)))->FIELD)) -#define validate_const_container_of(PTR, TYPE, FIELD) \ - (0 * sizeof((PTR) - &((const TYPE *)(((const char *)(PTR)) - \ - offsetof(TYPE, FIELD)))->FIELD)) -/** Downcasting macro */ -#define container_of(PTR, TYPE, FIELD) \ - ((TYPE *)(((char *)(PTR)) - offsetof(TYPE, FIELD)) \ - + validate_container_of(PTR, TYPE, FIELD)) -/** Constified downcasting macro */ -#define const_container_of(PTR, TYPE, FIELD) \ - ((const TYPE *)(((const char *)(PTR)) - offsetof(TYPE, FIELD)) \ - + validate_const_container_of(PTR, TYPE, FIELD)) +/** Declare a function that downcasts from a generic dist struct to the actual + * subtype probablity distribution it represents. */ +#define DECLARE_PROB_DISTR_DOWNCAST_FN(name) \ + static inline \ + const struct name * \ + dist_to_const_##name(const struct dist *obj) { \ + tor_assert(obj->ops == &name##_ops); \ + return SUBTYPE_P(obj, struct name, base); \ + } +DECLARE_PROB_DISTR_DOWNCAST_FN(uniform) +DECLARE_PROB_DISTR_DOWNCAST_FN(geometric) +DECLARE_PROB_DISTR_DOWNCAST_FN(logistic) +DECLARE_PROB_DISTR_DOWNCAST_FN(log_logistic) +DECLARE_PROB_DISTR_DOWNCAST_FN(genpareto) +DECLARE_PROB_DISTR_DOWNCAST_FN(weibull) /** * Count number of one bits in 32-bit word. @@ -458,7 +459,7 @@ random_uniform_01(void) * system is broken. */ z = 0; - while ((x = crypto_rand_u32()) == 0) { + while ((x = crypto_fast_rng_get_u32(get_thread_fast_rng())) == 0) { if (z >= 1088) /* Your bit sampler is broken. Go home. */ return 0; @@ -472,8 +473,8 @@ random_uniform_01(void) * occur only with measure zero in the uniform distribution on * [0, 1]. */ - hi = crypto_rand_u32() | UINT32_C(0x80000000); - lo = crypto_rand_u32() | UINT32_C(0x00000001); + hi = crypto_fast_rng_get_u32(get_thread_fast_rng()) | UINT32_C(0x80000000); + lo = crypto_fast_rng_get_u32(get_thread_fast_rng()) | UINT32_C(0x00000001); /* Round to nearest scaled significand in [2^63, 2^64]. */ s = hi*(double)4294967296 + lo; @@ -1315,40 +1316,48 @@ sample_geometric(uint32_t s, double p0, double p) /** Public API for probability distributions: * - * For each probability distribution we define each public functions - * (sample/cdf/sf/icdf/isf) as part of its dist_ops structure. + * These are wrapper functions on top of the various probability distribution + * operations using the generic <b>dist</b> structure. + + * These are the functions that should be used by consumers of this API. */ +/** Returns the name of the distribution in <b>dist</b>. */ const char * dist_name(const struct dist *dist) { return dist->ops->name; } +/* Sample a value from <b>dist</b> and return it. */ double dist_sample(const struct dist *dist) { return dist->ops->sample(dist); } +/** Compute the CDF of <b>dist</b> at <b>x</b>. */ double dist_cdf(const struct dist *dist, double x) { return dist->ops->cdf(dist, x); } +/** Compute the SF (Survival function) of <b>dist</b> at <b>x</b>. */ double dist_sf(const struct dist *dist, double x) { return dist->ops->sf(dist, x); } +/** Compute the iCDF (Inverse CDF) of <b>dist</b> at <b>x</b>. */ double dist_icdf(const struct dist *dist, double p) { return dist->ops->icdf(dist, p); } +/** Compute the iSF (Inverse Survival function) of <b>dist</b> at <b>x</b>. */ double dist_isf(const struct dist *dist, double p) { @@ -1360,8 +1369,7 @@ dist_isf(const struct dist *dist, double p) static double uniform_sample(const struct dist *dist) { - const struct uniform *U = const_container_of(dist, struct uniform, - base); + const struct uniform *U = dist_to_const_uniform(dist); double p0 = random_uniform_01(); return sample_uniform_interval(p0, U->a, U->b); @@ -1370,9 +1378,7 @@ uniform_sample(const struct dist *dist) static double uniform_cdf(const struct dist *dist, double x) { - const struct uniform *U = const_container_of(dist, struct uniform, - base); - + const struct uniform *U = dist_to_const_uniform(dist); if (x < U->a) return 0; else if (x < U->b) @@ -1384,8 +1390,7 @@ uniform_cdf(const struct dist *dist, double x) static double uniform_sf(const struct dist *dist, double x) { - const struct uniform *U = const_container_of(dist, struct uniform, - base); + const struct uniform *U = dist_to_const_uniform(dist); if (x > U->b) return 0; @@ -1398,8 +1403,7 @@ uniform_sf(const struct dist *dist, double x) static double uniform_icdf(const struct dist *dist, double p) { - const struct uniform *U = const_container_of(dist, struct uniform, - base); + const struct uniform *U = dist_to_const_uniform(dist); double w = U->b - U->a; return (p < 0.5 ? (U->a + w*p) : (U->b - w*(1 - p))); @@ -1408,8 +1412,7 @@ uniform_icdf(const struct dist *dist, double p) static double uniform_isf(const struct dist *dist, double p) { - const struct uniform *U = const_container_of(dist, struct uniform, - base); + const struct uniform *U = dist_to_const_uniform(dist); double w = U->b - U->a; return (p < 0.5 ? (U->b - w*p) : (U->a + w*(1 - p))); @@ -1424,14 +1427,17 @@ const struct dist_ops uniform_ops = { .isf = uniform_isf, }; +/*******************************************************************/ + +/** Private functions for each probability distribution. */ + /** Functions for logistic distribution: */ static double logistic_sample(const struct dist *dist) { - const struct logistic *L = const_container_of(dist, struct logistic, - base); - uint32_t s = crypto_rand_u32(); + const struct logistic *L = dist_to_const_logistic(dist); + uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng()); double t = random_uniform_01(); double p0 = random_uniform_01(); @@ -1441,36 +1447,28 @@ logistic_sample(const struct dist *dist) static double logistic_cdf(const struct dist *dist, double x) { - const struct logistic *L = const_container_of(dist, struct logistic, - base); - + const struct logistic *L = dist_to_const_logistic(dist); return cdf_logistic(x, L->mu, L->sigma); } static double logistic_sf(const struct dist *dist, double x) { - const struct logistic *L = const_container_of(dist, struct logistic, - base); - + const struct logistic *L = dist_to_const_logistic(dist); return sf_logistic(x, L->mu, L->sigma); } static double logistic_icdf(const struct dist *dist, double p) { - const struct logistic *L = const_container_of(dist, struct logistic, - base); - + const struct logistic *L = dist_to_const_logistic(dist); return icdf_logistic(p, L->mu, L->sigma); } static double logistic_isf(const struct dist *dist, double p) { - const struct logistic *L = const_container_of(dist, struct logistic, - base); - + const struct logistic *L = dist_to_const_logistic(dist); return isf_logistic(p, L->mu, L->sigma); } @@ -1488,9 +1486,8 @@ const struct dist_ops logistic_ops = { static double log_logistic_sample(const struct dist *dist) { - const struct log_logistic *LL = const_container_of(dist, struct - log_logistic, base); - uint32_t s = crypto_rand_u32(); + const struct log_logistic *LL = dist_to_const_log_logistic(dist); + uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng()); double p0 = random_uniform_01(); return sample_log_logistic_scaleshape(s, p0, LL->alpha, LL->beta); @@ -1499,36 +1496,28 @@ log_logistic_sample(const struct dist *dist) static double log_logistic_cdf(const struct dist *dist, double x) { - const struct log_logistic *LL = const_container_of(dist, - struct log_logistic, base); - + const struct log_logistic *LL = dist_to_const_log_logistic(dist); return cdf_log_logistic(x, LL->alpha, LL->beta); } static double log_logistic_sf(const struct dist *dist, double x) { - const struct log_logistic *LL = const_container_of(dist, - struct log_logistic, base); - + const struct log_logistic *LL = dist_to_const_log_logistic(dist); return sf_log_logistic(x, LL->alpha, LL->beta); } static double log_logistic_icdf(const struct dist *dist, double p) { - const struct log_logistic *LL = const_container_of(dist, - struct log_logistic, base); - + const struct log_logistic *LL = dist_to_const_log_logistic(dist); return icdf_log_logistic(p, LL->alpha, LL->beta); } static double log_logistic_isf(const struct dist *dist, double p) { - const struct log_logistic *LL = const_container_of(dist, - struct log_logistic, base); - + const struct log_logistic *LL = dist_to_const_log_logistic(dist); return isf_log_logistic(p, LL->alpha, LL->beta); } @@ -1546,9 +1535,8 @@ const struct dist_ops log_logistic_ops = { static double weibull_sample(const struct dist *dist) { - const struct weibull *W = const_container_of(dist, struct weibull, - base); - uint32_t s = crypto_rand_u32(); + const struct weibull *W = dist_to_const_weibull(dist); + uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng()); double p0 = random_uniform_01(); return sample_weibull(s, p0, W->lambda, W->k); @@ -1557,36 +1545,28 @@ weibull_sample(const struct dist *dist) static double weibull_cdf(const struct dist *dist, double x) { - const struct weibull *W = const_container_of(dist, struct weibull, - base); - + const struct weibull *W = dist_to_const_weibull(dist); return cdf_weibull(x, W->lambda, W->k); } static double weibull_sf(const struct dist *dist, double x) { - const struct weibull *W = const_container_of(dist, struct weibull, - base); - + const struct weibull *W = dist_to_const_weibull(dist); return sf_weibull(x, W->lambda, W->k); } static double weibull_icdf(const struct dist *dist, double p) { - const struct weibull *W = const_container_of(dist, struct weibull, - base); - + const struct weibull *W = dist_to_const_weibull(dist); return icdf_weibull(p, W->lambda, W->k); } static double weibull_isf(const struct dist *dist, double p) { - const struct weibull *W = const_container_of(dist, struct weibull, - base); - + const struct weibull *W = dist_to_const_weibull(dist); return isf_weibull(p, W->lambda, W->k); } @@ -1604,9 +1584,8 @@ const struct dist_ops weibull_ops = { static double genpareto_sample(const struct dist *dist) { - const struct genpareto *GP = const_container_of(dist, struct genpareto, - base); - uint32_t s = crypto_rand_u32(); + const struct genpareto *GP = dist_to_const_genpareto(dist); + uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng()); double p0 = random_uniform_01(); return sample_genpareto_locscale(s, p0, GP->mu, GP->sigma, GP->xi); @@ -1615,36 +1594,28 @@ genpareto_sample(const struct dist *dist) static double genpareto_cdf(const struct dist *dist, double x) { - const struct genpareto *GP = const_container_of(dist, struct genpareto, - base); - + const struct genpareto *GP = dist_to_const_genpareto(dist); return cdf_genpareto(x, GP->mu, GP->sigma, GP->xi); } static double genpareto_sf(const struct dist *dist, double x) { - const struct genpareto *GP = const_container_of(dist, struct genpareto, - base); - + const struct genpareto *GP = dist_to_const_genpareto(dist); return sf_genpareto(x, GP->mu, GP->sigma, GP->xi); } static double genpareto_icdf(const struct dist *dist, double p) { - const struct genpareto *GP = const_container_of(dist, struct genpareto, - base); - + const struct genpareto *GP = dist_to_const_genpareto(dist); return icdf_genpareto(p, GP->mu, GP->sigma, GP->xi); } static double genpareto_isf(const struct dist *dist, double p) { - const struct genpareto *GP = const_container_of(dist, struct genpareto, - base); - + const struct genpareto *GP = dist_to_const_genpareto(dist); return isf_genpareto(p, GP->mu, GP->sigma, GP->xi); } @@ -1662,8 +1633,8 @@ const struct dist_ops genpareto_ops = { static double geometric_sample(const struct dist *dist) { - const struct geometric *G = const_container_of(dist, struct geometric, base); - uint32_t s = crypto_rand_u32(); + const struct geometric *G = dist_to_const_geometric(dist); + uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng()); double p0 = random_uniform_01(); return sample_geometric(s, p0, G->p); @@ -1672,7 +1643,7 @@ geometric_sample(const struct dist *dist) static double geometric_cdf(const struct dist *dist, double x) { - const struct geometric *G = const_container_of(dist, struct geometric, base); + const struct geometric *G = dist_to_const_geometric(dist); if (x < 1) return 0; @@ -1683,7 +1654,7 @@ geometric_cdf(const struct dist *dist, double x) static double geometric_sf(const struct dist *dist, double x) { - const struct geometric *G = const_container_of(dist, struct geometric, base); + const struct geometric *G = dist_to_const_geometric(dist); if (x < 1) return 0; @@ -1694,7 +1665,7 @@ geometric_sf(const struct dist *dist, double x) static double geometric_icdf(const struct dist *dist, double p) { - const struct geometric *G = const_container_of(dist, struct geometric, base); + const struct geometric *G = dist_to_const_geometric(dist); return log1p(-p)/log1p(-G->p); } @@ -1702,7 +1673,7 @@ geometric_icdf(const struct dist *dist, double p) static double geometric_isf(const struct dist *dist, double p) { - const struct geometric *G = const_container_of(dist, struct geometric, base); + const struct geometric *G = dist_to_const_geometric(dist); return log(p)/log1p(-G->p); } |