commit bdea94a6653240e0e477eda186ceccbe6be00535 Author: Nick Mathewson <nickm@torproject.org> Date: Thu Jun 28 12:24:45 2018 -0400 Move floating-point math functions into a new lib/math --- .gitignore | 2 + Makefile.am | 2 + src/common/util.c | 158 ---------------------------------------------- src/common/util.h | 7 -- src/include.am | 1 + src/lib/math/.may_include | 5 ++ src/lib/math/fp.c | 113 +++++++++++++++++++++++++++++++++ src/lib/math/fp.h | 17 +++++ src/lib/math/include.am | 20 ++++++ src/lib/math/laplace.c | 67 ++++++++++++++++++++ src/lib/math/laplace.h | 16 +++++ src/or/circpathbias.c | 3 +- src/or/circuitstats.c | 2 +- src/or/circuituse.c | 2 +- src/or/entrynodes.c | 1 + src/or/rephist.c | 2 + src/or/routerlist.c | 1 + src/or/scheduler_kist.c | 2 +- src/test/test_util.c | 2 + 19 files changed, 254 insertions(+), 169 deletions(-) diff --git a/.gitignore b/.gitignore index 8f2d66b0d..c940166aa 100644 --- a/.gitignore +++ b/.gitignore @@ -187,6 +187,8 @@ uptime-*.json /src/lib/libtor-log-testing.a /src/lib/libtor-malloc.a /src/lib/libtor-malloc-testing.a +/src/lib/libtor-math.a +/src/lib/libtor-math-testing.a /src/lib/libtor-memarea.a /src/lib/libtor-memarea-testing.a /src/lib/libtor-net.a diff --git a/Makefile.am b/Makefile.am index f28eae3bb..5a7271e30 100644 --- a/Makefile.am +++ b/Makefile.am @@ -48,6 +48,7 @@ TOR_UTIL_LIBS = \ src/lib/libtor-net.a \ src/lib/libtor-thread.a \ src/lib/libtor-memarea.a \ + src/lib/libtor-math.a \ src/lib/libtor-log.a \ src/lib/libtor-lock.a \ src/lib/libtor-fdio.a \ @@ -71,6 +72,7 @@ TOR_UTIL_TESTING_LIBS = \ src/lib/libtor-net-testing.a \ src/lib/libtor-thread-testing.a \ src/lib/libtor-memarea-testing.a \ + src/lib/libtor-math-testing.a \ src/lib/libtor-log-testing.a \ src/lib/libtor-lock-testing.a \ src/lib/libtor-fdio-testing.a \ diff --git a/src/common/util.c b/src/common/util.c index 6a557c956..e1fae0dea 100644 --- a/src/common/util.c +++ b/src/common/util.c @@ -127,100 +127,6 @@ ENABLE_GCC_WARNING(aggregate-return) * Math * ===== */ -/** - * Returns the natural logarithm of d base e. We defined this wrapper here so - * to avoid conflicts with old versions of tor_log(), which were named log(). - */ -double -tor_mathlog(double d) -{ - return log(d); -} - -/** Return the long integer closest to <b>d</b>. We define this wrapper - * here so that not all users of math.h need to use the right incantations - * to get the c99 functions. */ -long -tor_lround(double d) -{ -#if defined(HAVE_LROUND) - return lround(d); -#elif defined(HAVE_RINT) - return (long)rint(d); -#else - return (long)(d > 0 ? d + 0.5 : ceil(d - 0.5)); -#endif /* defined(HAVE_LROUND) || ... */ -} - -/** Return the 64-bit integer closest to d. We define this wrapper here so - * that not all users of math.h need to use the right incantations to get the - * c99 functions. */ -int64_t -tor_llround(double d) -{ -#if defined(HAVE_LLROUND) - return (int64_t)llround(d); -#elif defined(HAVE_RINT) - return (int64_t)rint(d); -#else - return (int64_t)(d > 0 ? d + 0.5 : ceil(d - 0.5)); -#endif /* defined(HAVE_LLROUND) || ... */ -} - -/** 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>. Truncate the final result - * to be an integer in [INT64_MIN, INT64_MAX]. */ -int64_t -sample_laplace_distribution(double mu, double b, double p) -{ - double result; - tor_assert(p >= 0.0 && p < 1.0); - - /* This is the "inverse cumulative distribution function" from: - * http://en.wikipedia.org/wiki/Laplace_distribution */ - if (p <= 0.0) { - /* Avoid taking log(0.0) == -INFINITY, as some processors or compiler - * options can cause the program to trap. */ - return INT64_MIN; - } - - result = mu - b * (p > 0.5 ? 1.0 : -1.0) - * tor_mathlog(1.0 - 2.0 * fabs(p - 0.5)); - - return clamp_double_to_int64(result); -} - -/** 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[. - * The epsilon value must be between ]0.0, 1.0]. delta_f must be greater - * than 0. */ -int64_t -add_laplace_noise(int64_t signal_, double random_, double delta_f, - double epsilon) -{ - int64_t noise; - - /* epsilon MUST be between ]0.0, 1.0] */ - tor_assert(epsilon > 0.0 && epsilon <= 1.0); - /* delta_f MUST be greater than 0. */ - tor_assert(delta_f > 0.0); - - /* Just add noise, no further signal */ - noise = sample_laplace_distribution(0.0, - delta_f / epsilon, - random_); - - /* Clip (signal + noise) to [INT64_MIN, INT64_MAX] */ - 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; -} - /* ===== * String manipulation * ===== */ @@ -389,67 +295,3 @@ load_windows_system_library(const TCHAR *library_name) return LoadLibrary(path); } #endif /* defined(_WIN32) */ - -/** Cast a given double value to a int64_t. Return 0 if number is NaN. - * Returns either INT64_MIN or INT64_MAX if number is outside of the int64_t - * range. */ -int64_t -clamp_double_to_int64(double number) -{ - int exponent; - -#if defined(MINGW_ANY) && GCC_VERSION >= 409 -/* - Mingw's math.h uses gcc's __builtin_choose_expr() facility to declare - isnan, isfinite, and signbit. But as implemented in at least some - versions of gcc, __builtin_choose_expr() can generate type warnings - even from branches that are not taken. So, suppress those warnings. -*/ -#define PROBLEMATIC_FLOAT_CONVERSION_WARNING -DISABLE_GCC_WARNING(float-conversion) -#endif /* defined(MINGW_ANY) && GCC_VERSION >= 409 */ - -/* - With clang 4.0 we apparently run into "double promotion" warnings here, - since clang thinks we're promoting a double to a long double. - */ -#if defined(__clang__) -#if __has_warning("-Wdouble-promotion") -#define PROBLEMATIC_DOUBLE_PROMOTION_WARNING -DISABLE_GCC_WARNING(double-promotion) -#endif -#endif /* defined(__clang__) */ - - /* NaN is a special case that can't be used with the logic below. */ - if (isnan(number)) { - return 0; - } - - /* Time to validate if result can overflows a int64_t value. Fun with - * float! Find that exponent exp such that - * number == x * 2^exp - * for some x with abs(x) in [0.5, 1.0). Note that this implies that the - * magnitude of number is strictly less than 2^exp. - * - * If number is infinite, the call to frexp is legal but the contents of - * are exponent unspecified. */ - frexp(number, &exponent); - - /* If the magnitude of number is strictly less than 2^63, the truncated - * version of number is guaranteed to be representable. The only - * representable integer for which this is not the case is INT64_MIN, but - * it is covered by the logic below. */ - if (isfinite(number) && exponent <= 63) { - return (int64_t)number; - } - - /* Handle infinities and finite numbers with magnitude >= 2^63. */ - return signbit(number) ? INT64_MIN : INT64_MAX; - -#ifdef PROBLEMATIC_DOUBLE_PROMOTION_WARNING -ENABLE_GCC_WARNING(double-promotion) -#endif -#ifdef PROBLEMATIC_FLOAT_CONVERSION_WARNING -ENABLE_GCC_WARNING(float-conversion) -#endif -} diff --git a/src/common/util.h b/src/common/util.h index 0f991f8c4..3068f023a 100644 --- a/src/common/util.h +++ b/src/common/util.h @@ -62,13 +62,6 @@ void tor_log_mallinfo(int severity); #define bool_neq(a,b) (!(a)!=!(b)) /* Math functions */ -double tor_mathlog(double d) ATTR_CONST; -long tor_lround(double d) ATTR_CONST; -int64_t tor_llround(double d) ATTR_CONST; -int64_t sample_laplace_distribution(double mu, double b, double p); -int64_t add_laplace_noise(int64_t signal, double random, double delta_f, - double epsilon); -int64_t clamp_double_to_int64(double number); /* String manipulation */ diff --git a/src/include.am b/src/include.am index 621c906d9..8345813b7 100644 --- a/src/include.am +++ b/src/include.am @@ -14,6 +14,7 @@ include src/lib/include.libdonna.am include src/lib/intmath/include.am include src/lib/lock/include.am include src/lib/log/include.am +include src/lib/math/include.am include src/lib/memarea/include.am include src/lib/malloc/include.am include src/lib/net/include.am diff --git a/src/lib/math/.may_include b/src/lib/math/.may_include new file mode 100644 index 000000000..1fd26864d --- /dev/null +++ b/src/lib/math/.may_include @@ -0,0 +1,5 @@ +orconfig.h + +lib/cc/*.h +lib/log/*.h +lib/math/*.h diff --git a/src/lib/math/fp.c b/src/lib/math/fp.c new file mode 100644 index 000000000..d1c442825 --- /dev/null +++ b/src/lib/math/fp.c @@ -0,0 +1,113 @@ +/* Copyright (c) 2003, Roger Dingledine + * Copyright (c) 2004-2006, Roger Dingledine, Nick Mathewson. + * Copyright (c) 2007-2018, The Tor Project, Inc. */ +/* See LICENSE for licensing information */ + +#include "orconfig.h" +#include "lib/math/fp.h" + +#include <math.h> + +/** + * Returns the natural logarithm of d base e. We defined this wrapper here so + * to avoid conflicts with old versions of tor_log(), which were named log(). + */ +double +tor_mathlog(double d) +{ + return log(d); +} + +/** Return the long integer closest to <b>d</b>. We define this wrapper + * here so that not all users of math.h need to use the right incantations + * to get the c99 functions. */ +long +tor_lround(double d) +{ +#if defined(HAVE_LROUND) + return lround(d); +#elif defined(HAVE_RINT) + return (long)rint(d); +#else + return (long)(d > 0 ? d + 0.5 : ceil(d - 0.5)); +#endif /* defined(HAVE_LROUND) || ... */ +} + +/** Return the 64-bit integer closest to d. We define this wrapper here so + * that not all users of math.h need to use the right incantations to get the + * c99 functions. */ +int64_t +tor_llround(double d) +{ +#if defined(HAVE_LLROUND) + return (int64_t)llround(d); +#elif defined(HAVE_RINT) + return (int64_t)rint(d); +#else + return (int64_t)(d > 0 ? d + 0.5 : ceil(d - 0.5)); +#endif /* defined(HAVE_LLROUND) || ... */ +} + +/** Cast a given double value to a int64_t. Return 0 if number is NaN. + * Returns either INT64_MIN or INT64_MAX if number is outside of the int64_t + * range. */ +int64_t +clamp_double_to_int64(double number) +{ + int exponent; + +#if defined(MINGW_ANY) && GCC_VERSION >= 409 +/* + Mingw's math.h uses gcc's __builtin_choose_expr() facility to declare + isnan, isfinite, and signbit. But as implemented in at least some + versions of gcc, __builtin_choose_expr() can generate type warnings + even from branches that are not taken. So, suppress those warnings. +*/ +#define PROBLEMATIC_FLOAT_CONVERSION_WARNING +DISABLE_GCC_WARNING(float-conversion) +#endif /* defined(MINGW_ANY) && GCC_VERSION >= 409 */ + +/* + With clang 4.0 we apparently run into "double promotion" warnings here, + since clang thinks we're promoting a double to a long double. + */ +#if defined(__clang__) +#if __has_warning("-Wdouble-promotion") +#define PROBLEMATIC_DOUBLE_PROMOTION_WARNING +DISABLE_GCC_WARNING(double-promotion) +#endif +#endif /* defined(__clang__) */ + + /* NaN is a special case that can't be used with the logic below. */ + if (isnan(number)) { + return 0; + } + + /* Time to validate if result can overflows a int64_t value. Fun with + * float! Find that exponent exp such that + * number == x * 2^exp + * for some x with abs(x) in [0.5, 1.0). Note that this implies that the + * magnitude of number is strictly less than 2^exp. + * + * If number is infinite, the call to frexp is legal but the contents of + * are exponent unspecified. */ + frexp(number, &exponent); + + /* If the magnitude of number is strictly less than 2^63, the truncated + * version of number is guaranteed to be representable. The only + * representable integer for which this is not the case is INT64_MIN, but + * it is covered by the logic below. */ + if (isfinite(number) && exponent <= 63) { + return (int64_t)number; + } + + /* Handle infinities and finite numbers with magnitude >= 2^63. */ + return signbit(number) ? INT64_MIN : INT64_MAX; + +#ifdef PROBLEMATIC_DOUBLE_PROMOTION_WARNING +ENABLE_GCC_WARNING(double-promotion) +#endif +#ifdef PROBLEMATIC_FLOAT_CONVERSION_WARNING +ENABLE_GCC_WARNING(float-conversion) +#endif +} diff --git a/src/lib/math/fp.h b/src/lib/math/fp.h new file mode 100644 index 000000000..b35c18a1c --- /dev/null +++ b/src/lib/math/fp.h @@ -0,0 +1,17 @@ +/* Copyright (c) 2003, Roger Dingledine + * Copyright (c) 2004-2006, Roger Dingledine, Nick Mathewson. + * Copyright (c) 2007-2018, The Tor Project, Inc. */ +/* See LICENSE for licensing information */ + +#ifndef TOR_FP_H +#define TOR_FP_H + +#include "lib/cc/compat_compiler.h" +#include "lib/cc/torint.h" + +double tor_mathlog(double d) ATTR_CONST; +long tor_lround(double d) ATTR_CONST; +int64_t tor_llround(double d) ATTR_CONST; +int64_t clamp_double_to_int64(double number); + +#endif diff --git a/src/lib/math/include.am b/src/lib/math/include.am new file mode 100644 index 000000000..b088b3f3c --- /dev/null +++ b/src/lib/math/include.am @@ -0,0 +1,20 @@ + +noinst_LIBRARIES += src/lib/libtor-math.a + +if UNITTESTS_ENABLED +noinst_LIBRARIES += src/lib/libtor-math-testing.a +endif + +src_lib_libtor_math_a_SOURCES = \ + src/lib/math/fp.c \ + src/lib/math/laplace.c + + +src_lib_libtor_math_testing_a_SOURCES = \ + $(src_lib_libtor_math_a_SOURCES) +src_lib_libtor_math_testing_a_CPPFLAGS = $(AM_CPPFLAGS) $(TEST_CPPFLAGS) +src_lib_libtor_math_testing_a_CFLAGS = $(AM_CFLAGS) $(TEST_CFLAGS) + +noinst_HEADERS += \ + src/lib/math/fp.h \ + src/lib/math/laplace.h diff --git a/src/lib/math/laplace.c b/src/lib/math/laplace.c new file mode 100644 index 000000000..8e45a1fb3 --- /dev/null +++ b/src/lib/math/laplace.c @@ -0,0 +1,67 @@ +/* Copyright (c) 2003, Roger Dingledine + * Copyright (c) 2004-2006, Roger Dingledine, Nick Mathewson. + * Copyright (c) 2007-2018, The Tor Project, Inc. */ +/* See LICENSE for licensing information */ + +#include "orconfig.h" +#include "lib/math/laplace.h" +#include "lib/math/fp.h" + +#include "lib/log/util_bug.h" + +#include <math.h> +#include <stdlib.h> + +/** 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>. Truncate the final result + * to be an integer in [INT64_MIN, INT64_MAX]. */ +int64_t +sample_laplace_distribution(double mu, double b, double p) +{ + double result; + tor_assert(p >= 0.0 && p < 1.0); + + /* This is the "inverse cumulative distribution function" from: + * http://en.wikipedia.org/wiki/Laplace_distribution */ + if (p <= 0.0) { + /* Avoid taking log(0.0) == -INFINITY, as some processors or compiler + * options can cause the program to trap. */ + return INT64_MIN; + } + + result = mu - b * (p > 0.5 ? 1.0 : -1.0) + * tor_mathlog(1.0 - 2.0 * fabs(p - 0.5)); + + return clamp_double_to_int64(result); +} + +/** 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[. + * The epsilon value must be between ]0.0, 1.0]. delta_f must be greater + * than 0. */ +int64_t +add_laplace_noise(int64_t signal_, double random_, double delta_f, + double epsilon) +{ + int64_t noise; + + /* epsilon MUST be between ]0.0, 1.0] */ + tor_assert(epsilon > 0.0 && epsilon <= 1.0); + /* delta_f MUST be greater than 0. */ + tor_assert(delta_f > 0.0); + + /* Just add noise, no further signal */ + noise = sample_laplace_distribution(0.0, + delta_f / epsilon, + random_); + + /* Clip (signal + noise) to [INT64_MIN, INT64_MAX] */ + 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; +} diff --git a/src/lib/math/laplace.h b/src/lib/math/laplace.h new file mode 100644 index 000000000..b22862e64 --- /dev/null +++ b/src/lib/math/laplace.h @@ -0,0 +1,16 @@ +/* Copyright (c) 2003, Roger Dingledine + * Copyright (c) 2004-2006, Roger Dingledine, Nick Mathewson. + * Copyright (c) 2007-2018, The Tor Project, Inc. */ +/* See LICENSE for licensing information */ + +#ifndef TOR_LAPLACE_H +#define TOR_LAPLACE_H + +#include "lib/cc/compat_compiler.h" +#include "lib/cc/torint.h" + +int64_t sample_laplace_distribution(double mu, double b, double p); +int64_t add_laplace_noise(int64_t signal, double random, double delta_f, + double epsilon); + +#endif diff --git a/src/or/circpathbias.c b/src/or/circpathbias.c index 908b76b48..32b3212d3 100644 --- a/src/or/circpathbias.c +++ b/src/or/circpathbias.c @@ -34,6 +34,8 @@ #include "or/entrynodes.h" #include "or/networkstatus.h" #include "or/relay.h" +#include "lib/math/fp.h" +#include "lib/math/laplace.h" #include "or/cell_st.h" #include "or/cpath_build_state_st.h" @@ -1574,4 +1576,3 @@ pathbias_scale_use_rates(entry_guard_t *guard) entry_guards_changed(); } } - diff --git a/src/or/circuitstats.c b/src/or/circuitstats.c index 34605457d..08186ca9a 100644 --- a/src/or/circuitstats.c +++ b/src/or/circuitstats.c @@ -40,6 +40,7 @@ #include "or/statefile.h" #include "or/circuitlist.h" #include "or/circuituse.h" +#include "lib/math/fp.h" #include "or/crypt_path_st.h" #include "or/origin_circuit_st.h" @@ -1946,4 +1947,3 @@ cbt_control_event_buildtimeout_set(const circuit_build_times_t *cbt, tor_free(args); } - diff --git a/src/or/circuituse.c b/src/or/circuituse.c index 0ba3f34b4..3669eb7f1 100644 --- a/src/or/circuituse.c +++ b/src/or/circuituse.c @@ -56,6 +56,7 @@ #include "or/rephist.h" #include "or/router.h" #include "or/routerlist.h" +#include "lib/math/fp.h" #include "or/cpath_build_state_st.h" #include "or/dir_connection_st.h" @@ -3152,4 +3153,3 @@ circuit_read_valid_data(origin_circuit_t *circ, uint16_t relay_body_len) tor_add_u32_nowrap(circ->n_overhead_read_circ_bw, RELAY_PAYLOAD_SIZE-relay_body_len); } - diff --git a/src/or/entrynodes.c b/src/or/entrynodes.c index 47e689b5c..191933622 100644 --- a/src/or/entrynodes.c +++ b/src/or/entrynodes.c @@ -138,6 +138,7 @@ #include "or/routerset.h" #include "or/transports.h" #include "or/statefile.h" +#include "lib/math/fp.h" #include "or/node_st.h" #include "or/origin_circuit_st.h" diff --git a/src/or/rephist.c b/src/or/rephist.c index 2103eecdf..907b01d68 100644 --- a/src/or/rephist.c +++ b/src/or/rephist.c @@ -94,6 +94,8 @@ #include "lib/container/bloomfilt.h" #include "lib/container/order.h" +#include "lib/math/fp.h" +#include "lib/math/laplace.h" static void bw_arrays_init(void); static void predicted_ports_alloc(void); diff --git a/src/or/routerlist.c b/src/or/routerlist.c index a86e29adb..68fd76344 100644 --- a/src/or/routerlist.c +++ b/src/or/routerlist.c @@ -121,6 +121,7 @@ #include "or/routerset.h" #include "lib/sandbox/sandbox.h" #include "or/torcert.h" +#include "lib/math/fp.h" #include "or/dirauth/dirvote.h" #include "or/dirauth/mode.h" diff --git a/src/or/scheduler_kist.c b/src/or/scheduler_kist.c index 6f07458d4..b62281d58 100644 --- a/src/or/scheduler_kist.c +++ b/src/or/scheduler_kist.c @@ -13,6 +13,7 @@ #include "or/channeltls.h" #define SCHEDULER_PRIVATE_ #include "or/scheduler.h" +#include "lib/math/fp.h" #include "or/or_connection_st.h" @@ -835,4 +836,3 @@ scheduler_can_use_kist(void) } #endif /* defined(HAVE_KIST_SUPPORT) */ - diff --git a/src/test/test_util.c b/src/test/test_util.c index bdc6fca7d..8fe308826 100644 --- a/src/test/test_util.c +++ b/src/test/test_util.c @@ -29,6 +29,8 @@ #include "lib/process/subprocess.h" #include "lib/intmath/weakrng.h" #include "lib/thread/numcpus.h" +#include "lib/math/fp.h" +#include "lib/math/laplace.h" #ifdef HAVE_PWD_H #include <pwd.h>