constexpr_iterated_squares_desc

Summary: A scaling array and associated shrink/growth algorithms for use by floating-point algorithms.

Reviewed By: ispeters

Differential Revision: D48483293

fbshipit-source-id: e9d6e928c4c056d8281ba8326c66c82defae1306

folly/ConstexprMath.h

@@ -17,6 +17,7 @@
 #pragma once
 
 #include <cassert>
+#include <cstddef>
 #include <cstdint>
 #include <functional>
 #include <limits>
@@ -80,6 +81,175 @@ constexpr typename floating_point_integral_constant<T, S, Value>::value_type
     floating_point_integral_constant<T, S, Value>::value;
 #endif
 
+//  ----
+
+namespace detail {
+
+template <typename T>
+constexpr size_t constexpr_iterated_squares_desc_size_(T const base) {
+  using lim = std::numeric_limits<T>;
+  size_t s = 1;
+  auto r = base;
+  while (r <= lim::max() / r) {
+    ++s;
+    r *= r;
+  }
+  return s;
+}
+
+} // namespace detail
+
+/// constexpr_iterated_squares_desc_size_v
+///
+/// Effectively calculates: floor(log(max_exponent)/log(base))
+///
+/// For use with constexpr_iterated_squares_desc below.
+template <typename Base>
+FOLLY_INLINE_VARIABLE constexpr size_t constexpr_iterated_squares_desc_size_v =
+    detail::constexpr_iterated_squares_desc_size_(Base::value);
+
+/// constexpr_iterated_squares_desc
+///
+/// A constexpr scaling array of integer powers-of-powers-of-two, descending,
+/// with the associated powers-of-two.
+///
+/// scaling = [..., {8, b^8}, {4, b^4}, {2, b^2}, {1, b^1}] for b = base
+///
+/// Includes select constexpr scaling algorithms based on the scaling array.
+///
+/// The scaling array and the scaling algorithms are general-purpose, if niche.
+/// They may be used by other constexpr math functions (floating-point) either
+/// to improve runtime performance or to improve numerical approximations.
+///
+/// Some compilers fail to support passing some types as non-type template
+/// params. In particular, long double is not universally supported. Therefore,
+/// this utility takes its base as a type rather than as a value. For floating-
+/// point integral bases, that is, bases of floating-point type but of integral
+/// value, floating_point_integral_constant is the easiest parameterization.
+template <typename T, std::size_t Size>
+struct constexpr_iterated_squares_desc {
+  static_assert(Size > 0, "requires non-zero size");
+
+  using size_type = decltype(Size);
+  using base_type = T;
+
+  struct item_type {
+    size_type power;
+    base_type scale;
+  };
+
+  static constexpr size_type size = Size;
+  base_type base;
+  item_type scaling[size];
+
+ private:
+  using lim = std::numeric_limits<base_type>;
+
+  static_assert(
+      lim::max_exponent < std::numeric_limits<size_type>::max(),
+      "size_type too small for base_type");
+
+ public:
+  explicit constexpr constexpr_iterated_squares_desc(base_type r) noexcept
+      : base{r}, scaling{} {
+    assert(size <= detail::constexpr_iterated_squares_desc_size_(base));
+    size_type i = 0;
+    size_type p = 1;
+    while (true) { // a for-loop might cause multiplication overflow below
+      scaling[size - 1 - i] = {p, r};
+      if (++i == size) {
+        break;
+      }
+      p *= 2;
+      r *= r;
+    }
+  }
+
+  /// shrink
+  ///
+  /// Returns scaling params of the form:
+  ///   item_type{power, scale} with scale = base ^ power
+  /// With power the smallest nonnegative integer such that:
+  ///   abs(num) / scale <= max
+  constexpr item_type shrink(base_type const num, base_type const max) const {
+    assert(max > base_type(0));
+    auto const rmax = max / base;
+    auto const snum = num < base_type(0) ? -num : num;
+    auto power = size_type(0);
+    auto scale = base_type(1);
+    if (!(snum / scale <= max)) {
+      for (auto const& i : scaling) {
+        auto const next = scale * i.scale;
+        auto const div = snum / next;
+        if (div <= rmax) {
+          continue;
+        }
+        power += i.power;
+        scale = next;
+        if (div <= max) {
+          break;
+        }
+      }
+    }
+    assert(snum / scale <= max);
+    return {power, scale};
+  }
+
+  /// growth
+  ///
+  /// Returns scaling params of the form:
+  ///   item_type{power, scale} with scale = base ^ power
+  /// With power the smallest nonnegative integer such that:
+  ///   abs(num) * scale >= min
+  constexpr item_type growth(base_type const num, base_type const min) const {
+    assert(min > base_type(0));
+    auto const rmin = min * base;
+    auto const snum = num < base_type(0) ? -num : num;
+    auto power = size_type(0);
+    auto scale = base_type(1);
+    if (!(snum * scale >= min)) {
+      for (auto const& i : scaling) {
+        auto const next = scale * i.scale;
+        auto const mul = snum * next;
+        if (mul >= rmin) {
+          continue;
+        }
+        power += i.power;
+        scale = next;
+        if (mul >= min) {
+          break;
+        }
+      }
+    }
+    assert(snum * scale >= min);
+    return {power, scale};
+  }
+};
+#if FOLLY_CPLUSPLUS < 201703L
+template <typename T, std::size_t Size>
+constexpr typename constexpr_iterated_squares_desc<T, Size>::size_type
+    constexpr_iterated_squares_desc<T, Size>::size;
+#endif
+
+/// constexpr_iterated_squares_desc_v
+///
+/// An instance of constexpr_iterated_squares_desc of max size with the given
+/// base.
+template <typename Base>
+FOLLY_INLINE_VARIABLE constexpr auto constexpr_iterated_squares_desc_v =
+    constexpr_iterated_squares_desc<
+        typename Base::value_type,
+        constexpr_iterated_squares_desc_size_v<Base>>{Base::value};
+
+/// constexpr_iterated_squares_desc_2_v
+///
+/// An alias for constexpr_iterated_squares_desc_v with base 2, which is the
+/// most common base to use with iterated-squares.
+template <typename T>
+constexpr auto& constexpr_iterated_squares_desc_2_v =
+    constexpr_iterated_squares_desc_v<
+        floating_point_integral_constant<T, int, 2>>;
+
 // TLDR: Prefer using operator< for ordering. And when
 // a and b are equivalent objects, we return b to make
 // sorting stable.

folly/test/ConstexprMathTest.cpp

@@ -21,12 +21,197 @@
 #include <type_traits>
 
 #include <folly/lang/Bits.h>
+#include <folly/portability/GMock.h>
 #include <folly/portability/GTest.h>
 
 namespace {
 class ConstexprMathTest : public testing::Test {};
 } // namespace
 
+TEST_F(ConstexprMathTest, constexpr_iterated_squares_desc_scaling_array) {
+  using lim = std::numeric_limits<float>;
+  constexpr auto& isq = folly::constexpr_iterated_squares_desc_2_v<float>;
+
+  auto get = [](auto const& arr, auto fun) {
+    constexpr auto size = sizeof(arr) / sizeof(arr[0]);
+    using res_t = decltype(fun(arr[0]));
+    std::array<res_t, size> res{};
+    for (size_t i = 0; i < size; ++i) {
+      res[i] = fun(arr[i]);
+    }
+    return res;
+  };
+
+  EXPECT_EQ(7, isq.size);
+  auto apowers = get(isq.scaling, [](auto _) { return _.power; });
+  constexpr size_t epowers[] = {
+      64, 32, 16, 8, 4, 2, 1, //
+  };
+  EXPECT_THAT(apowers, testing::ElementsAreArray(epowers));
+  auto ascales = get(isq.scaling, [](auto _) { return _.scale; });
+  constexpr float escales[] = {
+      0x1p64, 0x1p32, 0x1p16, 0x1p8, 0x1p4, 0x1p2, 0x1p1, //
+  };
+  EXPECT_THAT(ascales, testing::ElementsAreArray(escales));
+  EXPECT_GT(isq.scaling[0].scale, lim::max() / isq.scaling[0].scale);
+}
+
+TEST_F(ConstexprMathTest, constexpr_iterated_squares_desc_shrink) {
+  using lim = std::numeric_limits<double>;
+  constexpr auto& isq = folly::constexpr_iterated_squares_desc_2_v<double>;
+
+  {
+    constexpr auto n = 1.;
+    constexpr auto a = isq.shrink(n, 1.);
+    EXPECT_EQ(0, a.power);
+    EXPECT_EQ(1., a.scale);
+    EXPECT_LE(n / a.scale, 1.);
+    EXPECT_GT(n / a.scale, .5);
+  }
+  {
+    constexpr auto n = 2.;
+    constexpr auto a = isq.shrink(n, 1.);
+    EXPECT_EQ(1, a.power);
+    EXPECT_EQ(2., a.scale);
+    EXPECT_LE(n / a.scale, 1.);
+    EXPECT_GT(n / a.scale, .5);
+  }
+  {
+    constexpr auto n = 4.;
+    constexpr auto a = isq.shrink(n, 1.);
+    EXPECT_EQ(2, a.power);
+    EXPECT_EQ(4., a.scale);
+    EXPECT_LE(n / a.scale, 1.);
+    EXPECT_GT(n / a.scale, .5);
+  }
+  {
+    constexpr auto n = 7.;
+    constexpr auto a = isq.shrink(n, 1.);
+    EXPECT_EQ(3, a.power);
+    EXPECT_EQ(8., a.scale);
+    EXPECT_LE(n / a.scale, 1.);
+    EXPECT_GT(n / a.scale, .5);
+  }
+  {
+    constexpr auto n = 8.;
+    constexpr auto a = isq.shrink(n, 1.);
+    EXPECT_EQ(3, a.power);
+    EXPECT_EQ(8., a.scale);
+    EXPECT_LE(n / a.scale, 1.);
+    EXPECT_GT(n / a.scale, .5);
+  }
+  {
+    constexpr auto n = 9.;
+    constexpr auto a = isq.shrink(n, 1.);
+    EXPECT_EQ(4, a.power);
+    EXPECT_EQ(16., a.scale);
+    EXPECT_LE(n / a.scale, 1.);
+    EXPECT_GT(n / a.scale, .5);
+  }
+  {
+    constexpr auto n = 513.;
+    constexpr auto a = isq.shrink(n, 1.);
+    EXPECT_EQ(10, a.power);
+    EXPECT_EQ(1024., a.scale);
+    EXPECT_LE(n / a.scale, 1.);
+    EXPECT_GT(n / a.scale, .5);
+  }
+  {
+    constexpr auto n = 1023.;
+    constexpr auto a = isq.shrink(n, 1.);
+    EXPECT_EQ(10, a.power);
+    EXPECT_EQ(1024., a.scale);
+    EXPECT_LE(n / a.scale, 1.);
+    EXPECT_GT(n / a.scale, .5);
+  }
+  {
+    constexpr auto n = lim::max();
+    constexpr auto a = isq.shrink(n, 2.);
+    EXPECT_EQ(1023, a.power);
+    EXPECT_EQ(folly::constexpr_pow(2., 1023), a.scale);
+    EXPECT_LE(n / a.scale, 2.);
+    EXPECT_GT(n / a.scale, .5);
+  }
+}
+
+TEST_F(ConstexprMathTest, constexpr_iterated_squares_desc_growth) {
+  using lim = std::numeric_limits<double>;
+  constexpr auto& isq = folly::constexpr_iterated_squares_desc_2_v<double>;
+
+  {
+    constexpr auto n = 1.;
+    constexpr auto a = isq.growth(n, 1.);
+    EXPECT_EQ(0, a.power);
+    EXPECT_EQ(1., a.scale);
+    EXPECT_GE(n * a.scale, 1.);
+    EXPECT_LT(n * a.scale, 2.);
+  }
+  {
+    constexpr auto n = 1. / 2.;
+    constexpr auto a = isq.growth(n, 1.);
+    EXPECT_EQ(1, a.power);
+    EXPECT_EQ(2., a.scale);
+    EXPECT_GE(n * a.scale, 1.);
+    EXPECT_LT(n * a.scale, 2.);
+  }
+  {
+    constexpr auto n = 1. / 4.;
+    constexpr auto a = isq.growth(n, 1.);
+    EXPECT_EQ(2, a.power);
+    EXPECT_EQ(4., a.scale);
+    EXPECT_GE(n * a.scale, 1.);
+    EXPECT_LT(n * a.scale, 2.);
+  }
+  {
+    constexpr auto n = 1. / 7.;
+    constexpr auto a = isq.growth(n, 1.);
+    EXPECT_EQ(3, a.power);
+    EXPECT_EQ(8., a.scale);
+    EXPECT_GE(n * a.scale, 1.);
+    EXPECT_LT(n * a.scale, 2.);
+  }
+  {
+    constexpr auto n = 1. / 8.;
+    constexpr auto a = isq.growth(n, 1.);
+    EXPECT_EQ(3, a.power);
+    EXPECT_EQ(8., a.scale);
+    EXPECT_GE(n * a.scale, 1.);
+    EXPECT_LT(n * a.scale, 2.);
+  }
+  {
+    constexpr auto n = 1. / 9.;
+    constexpr auto a = isq.growth(n, 1.);
+    EXPECT_EQ(4, a.power);
+    EXPECT_EQ(16., a.scale);
+    EXPECT_GE(n * a.scale, 1.);
+    EXPECT_LT(n * a.scale, 2.);
+  }
+  {
+    constexpr auto n = 1. / 513.;
+    constexpr auto a = isq.growth(n, 1.);
+    EXPECT_EQ(10, a.power);
+    EXPECT_EQ(1024., a.scale);
+    EXPECT_GE(n * a.scale, 1.);
+    EXPECT_LT(n * a.scale, 2.);
+  }
+  {
+    constexpr auto n = 1. / 1023.;
+    constexpr auto a = isq.growth(n, 1.);
+    EXPECT_EQ(10, a.power);
+    EXPECT_EQ(1024., a.scale);
+    EXPECT_GE(n * a.scale, 1.);
+    EXPECT_LT(n * a.scale, 2.);
+  }
+  {
+    constexpr auto n = lim::min(); // NOLINT
+    constexpr auto a = isq.growth(n, .5);
+    EXPECT_EQ(1021, a.power);
+    EXPECT_EQ(folly::constexpr_pow(2., 1021), a.scale);
+    EXPECT_GE(n * a.scale, .5);
+    EXPECT_LT(n * a.scale, 1.);
+  }
+}
+
 TEST_F(ConstexprMathTest, constexpr_min) {
   constexpr auto x = uint16_t(3);
   constexpr auto y = uint16_t(7);