Commit bdea94a6 authored by Nick Mathewson's avatar Nick Mathewson 🤹
Browse files

Move floating-point math functions into a new lib/math

parent 6178a9f7

.gitignore

+2 −0
Original line number Diff line number Diff line
@@ -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

Makefile.am

+2 −0
Original line number Diff line number Diff line
@@ -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 \

src/common/util.c

+0 −158
Original line number Diff line number Diff line
@@ -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

}

src/common/util.h

+0 −7
Original line number Diff line number Diff line
@@ -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 */

src/include.am

+1 −0
Original line number Diff line number Diff line
@@ -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