Use std::frexp for log2 approx

Author: geonove

ddsketch/include/linearly_interpolated_mapping_impl.hpp

@@ -32,11 +32,17 @@ inline LinearlyInterpolatedMapping::LinearlyInterpolatedMapping(const double& ga
   LogLikeIndexMapping<LinearlyInterpolatedMapping>(gamma, index_offset) {}
 
 inline double LinearlyInterpolatedMapping::log(const double& value) const {
-  return fast_log2(value);
+  int exponent = 0;
+  const double mantissa = std::frexp(value, &exponent);
+  const double significand = 2 * mantissa - 1;
+  return significand + (exponent - 1);
+
 }
 
 inline double LinearlyInterpolatedMapping::log_inverse(const double& index) const {
-  return fast_log2_inverse(index);
+  int exponent = static_cast<int>(std::floor(index)) + 1;
+  double mantissa = (index - exponent + 2) / 2.0;
+  return std::ldexp(mantissa, exponent);
 }
 
 inline IndexMappingLayout LinearlyInterpolatedMapping::layout() const {

ddsketch/include/quadratically_interpolated_mapping_impl.hpp

@@ -32,32 +32,26 @@ inline QuadraticallyInterpolatedMapping::QuadraticallyInterpolatedMapping(const
 
 
 inline double QuadraticallyInterpolatedMapping::log(const double& value) const {
-  int64_t value_bits;
-  std::memcpy(&value_bits, &value, sizeof(value));
-
-  const int64_t mantissa_plus_one_bits = (value_bits & 0x000FFFFFFFFFFFFFL) | 0x3FF0000000000000L;
-  double mantissa_plus_one;
-  std::memcpy(&mantissa_plus_one, &mantissa_plus_one_bits, sizeof(mantissa_plus_one_bits));
-
-  const double exponent = static_cast<double>(((value_bits & 0x7FF0000000000000L) >> 52) - 1023);
-
-  return exponent - (mantissa_plus_one - 5.0) * (mantissa_plus_one - 1) * ONE_THIRD;
+  // int64_t value_bits;
+  // std::memcpy(&value_bits, &value, sizeof(value));
+  //
+  // const int64_t mantissa_plus_one_bits = (value_bits & 0x000FFFFFFFFFFFFFL) | 0x3FF0000000000000L;
+  // double mantissa_plus_one;
+  // std::memcpy(&mantissa_plus_one, &mantissa_plus_one_bits, sizeof(mantissa_plus_one_bits));
+  //
+  // const double exponent = static_cast<double>(((value_bits & 0x7FF0000000000000L) >> 52) - 1023);
+
+  int exponent = 0;
+  const double mantissa = 2 * std::frexp(value, &exponent);
+
+  return exponent - 1 - (mantissa - 5.0) * (mantissa- 1) * ONE_THIRD;
 }
 
 inline double QuadraticallyInterpolatedMapping::log_inverse(const double& index) const {
-  const int64_t exponent = static_cast<int64_t>(std::floor(index));
+  const int exponent = static_cast<int>(std::floor(index));
   const double mantissa_plus_one = 3.0 - std::sqrt(4.0 - 3.0 * (index - exponent));
 
-  int64_t result_bits = (static_cast<uint64_t>(exponent + 1023) << 52) & 0x7FF0000000000000L;
-
-  uint64_t mantissa_plus_one_bits;
-  std::memcpy(&mantissa_plus_one_bits, &mantissa_plus_one, sizeof(mantissa_plus_one));
-
-  result_bits |= (mantissa_plus_one_bits & 0x000FFFFFFFFFFFFFL);
-
-  double result;
-  std::memcpy(&result, &result_bits, sizeof(result_bits));
-  return result;
+  return std::ldexp(mantissa_plus_one, exponent);
 }
 
 inline IndexMappingLayout QuadraticallyInterpolatedMapping::layout() const {

ddsketch/include/quartically_interpolated_mapping_impl.hpp

@@ -29,21 +29,15 @@ inline QuarticallyInterpolatedMapping::QuarticallyInterpolatedMapping(const doub
   LogLikeIndexMapping<QuarticallyInterpolatedMapping>(gamma, index_offset) {}
 
 inline double QuarticallyInterpolatedMapping::log(const double &value) const {
-  int64_t value_bits;
-  std::memcpy(&value_bits, &value, sizeof(value_bits));
+  int exponent = 0;
+  const double mantissa = 2 * std::frexp(value, &exponent) - 1;
 
-  const int64_t mantissa_plus_one_bits = (value_bits & 0x000FFFFFFFFFFFFFL) | 0x3FF0000000000000L;
-  double mantissa_plus_one;
-  std::memcpy(&mantissa_plus_one, &mantissa_plus_one_bits, sizeof(mantissa_plus_one_bits));
-  double mantissa = mantissa_plus_one - 1.0;
-
-  const double exponent = static_cast<double>(((value_bits & 0x7FF0000000000000L) >> 52) - 1023);
-  return (((A * mantissa + B) * mantissa + C) * mantissa + D) * mantissa + exponent;
+  return (((A * mantissa + B) * mantissa + C) * mantissa + D) * mantissa + exponent - 1;
 }
 
 inline double QuarticallyInterpolatedMapping::log_inverse(const double &index) const {
-  const int64_t exponent = static_cast<int64_t>(std::floor(index));
-  const double e = exponent - index;
+  const auto exponent = static_cast<int>(std::floor(index));
+  const double e = static_cast<double>(exponent) - index;
 
   // Derived from Ferrari's method
   const double alpha = -(3 * B * B) / (8 * A * A) + C / A; // 2.5
@@ -57,16 +51,10 @@ inline double QuarticallyInterpolatedMapping::log_inverse(const double &index) c
   const double w = std::sqrt(alpha + 2 * y);
   const double x = (-B / (4 * A) + (w - std::sqrt(-(3 * alpha + 2 * y + (2 * beta) / w))) / 2) + 1;
 
-  int64_t result_bits = (static_cast<uint64_t>(exponent + 1023) << 52) & 0x7FF0000000000000L;
-
-  uint64_t x_plus_one_bits;
-  std::memcpy(&x_plus_one_bits, &x, sizeof(x));
-
-  result_bits |= (x_plus_one_bits & 0x000FFFFFFFFFFFFFL);
-
-  double result;
-  std::memcpy(&result, &result_bits, sizeof(result_bits));
-  return result;
+  int tmp = 0;
+  double m = std::frexp(x, &tmp);
+  double sig = 2.0 * m;
+  return std::ldexp(sig, exponent);
 }
 
 inline IndexMappingLayout QuarticallyInterpolatedMapping::layout() const {