Program seems to hang in frexp for decimal128_t having large exponent

[ckormanyos]

I haven't quite figured out why, but the program below seems to hang indefinitely in the frexp function. I can't reach the debug line. It seems like there might be some wrong assumptions in the frexp function. I remember myself having quite a bit of trouble with the impl of frexp. I think this is a real issue, but I can't yet figure out what's going on.

Basically, the argument to frexp is

$$ \dfrac{104743{\times}10^{500}}{204745}{\approx}5.11577816308090551661823243546851{\times}10^{499} $$

which should be a valid decimal128_t for frexp, but maybe my assumptions in the test code are faulty?

#include <boost/decimal/decimal128_t.hpp>
#include <boost/decimal/cmath.hpp>
#include <boost/decimal/iostream.hpp>

volatile unsigned debug { };

auto main() -> int
{
  boost::decimal::decimal128_t arg_x { 104743, 500 };

  arg_x /= 204745;

  {
    std::stringstream strm { };

    strm << std::setprecision(std::numeric_limits<boost::decimal::decimal128_t>::digits10)
          << "arg_x: " << arg_x;

    std::cout << strm.str() << std::endl;
  }

  int exp2 { };

  // The program hangs here in frexp.
  const boost::decimal::decimal128_t frexp_result { frexp(arg_x, &exp2) };

  ++debug;
}

[mborland]

In this line result_frexp is multiplied by 0, and then in the loop in line 78 we end up infinitely multiplying zero by a non-zero number hoping it'll eventually be > 1/2

[ckormanyos]

Hi @mborland I can do a preliminary handling of this one based on our offline ideas. We should probably use a pre-pass thorugh frexp10 and then with a reduced number, go into frexp-base2. This should work for all exponents, large and small.

Also the frexp tests are somewhat weka. I'll be adding tests.

I'd like to get a working correction into 1.91, even if we have to optimize it later in the next cycle.

[ckormanyos]

Note to self: Test program two troublesome spot values.

#include <boost/decimal/decimal128_t.hpp>
#include <boost/decimal/cmath.hpp>
#include <boost/decimal/iostream.hpp>

template<typename DecimalType>
auto test(const DecimalType arg_x) -> void
{
  {
    std::stringstream strm { };

    strm << std::setprecision(std::numeric_limits<boost::decimal::decimal128_t>::digits10)
          << "arg_x: " << arg_x;

    std::cout << strm.str() << std::endl;
  }

  int exp2 { };

  // The program hands here in frexp.
  const boost::decimal::decimal128_t frexp_result { frexp(arg_x, &exp2) };

  {
    std::stringstream strm { };

    strm << std::setprecision(std::numeric_limits<boost::decimal::decimal128_t>::digits10)
          << "frexp_result: " << frexp_result;

    std::cout << strm.str() << std::endl;
  }
}

auto main() -> int
{
  {
    boost::decimal::decimal128_t arg_x { 104743, 4000 };

    arg_x /= 204745;

    test(arg_x);
  }

  {
    boost::decimal::decimal128_t arg_x { 104743, -4000 };

    arg_x /= 204745;

    test(arg_x);
  }
}

[ckormanyos]

Note to self: Stress test cbrt-versus-boost::math::cbrt. This program includes skinny Math bindings also.

#include <boost/decimal/decimal128_t.hpp>
#include <boost/decimal/cmath.hpp>
#include <boost/decimal/iostream.hpp>
#include <boost/math/special_functions/cbrt.hpp>

#include <chrono>
#include <iomanip>
#include <iostream>
#include <random>

namespace boost { namespace math { namespace policies {

// Specialization of the precision structure.
template<typename ThisPolicy>
struct precision<boost::decimal::decimal128_t, ThisPolicy>
{
  using local_decimal_type = boost::decimal::decimal128_t;

  using precision_type = typename ThisPolicy::precision_type;

  using local_digits_2 = digits2<((static_cast<std::intmax_t>(std::numeric_limits<local_decimal_type>::digits10) + INTMAX_C(1)) * INTMAX_C(1000)) / INTMAX_C(301)>;

  using type =
    typename std::conditional<((local_digits_2::value <= precision_type::value) || (precision_type::value <= 0)),
                                local_digits_2,
                                precision_type>::type;
};

} // namespace policies

namespace constants { namespace detail {

  template <class T>
  struct constant_pi;

  template<>
  struct constant_pi<boost::decimal::decimal128_t>
  {
  public:
    using result_type = boost::decimal::decimal128_t;

  public:
    template<int N> static auto get(const std::integral_constant<int, N>&) -> const result_type& // NOLINT(hicpp-named-parameter,readability-named-parameter)
    {
      static const result_type result { boost::decimal::numbers::pi_v<result_type> };

      return result;
    }
  };

  #if(__cplusplus >= 201703L)
  } // namespace constants::detail
  #else
  } // namespace detail
  } // namespace constants
  #endif

} // namespace math
} // namespace boost

namespace local
{
  template<typename IntegralTimePointType,
           typename ClockType = std::chrono::high_resolution_clock>
  auto time_point() noexcept -> IntegralTimePointType
  {
    using local_integral_time_point_type = IntegralTimePointType;
    using local_clock_type               = ClockType;

    const auto current_now =
      static_cast<std::uintmax_t>
      (
        std::chrono::duration_cast<std::chrono::nanoseconds>
        (
          local_clock_type::now().time_since_epoch()
        ).count()
      );

    return static_cast<local_integral_time_point_type>(current_now);
  }

  template<typename NumericType>
  auto is_close_fraction(const NumericType& a,
                         const NumericType& b,
                         const NumericType& tol) noexcept -> bool
  {
    using std::fabs;

    const NumericType delta
    {
      (b == static_cast<NumericType>(0)) ? fabs(a - b) : fabs(1 - (a / b))
    };

    const bool result_is_ok { (delta < tol) };

    return result_is_ok;
  }

  auto test_cbrt128_versus_math() -> void
  {
    std::random_device rd;
    std::mt19937_64 gen(rd());

    gen.seed(time_point<typename std::mt19937_64::result_type>());

    auto dis_denom =
      std::uniform_int_distribution<std::uint32_t>
      {
        UINT32_C(104744),
        UINT32_C(999999)
      };

    auto dis_exp10 =
      std::uniform_int_distribution<std::int32_t>
      {
        INT32_C(-4000),
        INT32_C(+4000)
      };

    bool result_is_ok { true };

    std::uint64_t count { };

    while(result_is_ok)
    //for(std::uint32_t index { UINT32_C(0) }; index < unsigned { UINT32_C(10) }; ++index)
    {
      boost::decimal::decimal128_t arg_x { UINT32_C(104743), dis_exp10(gen) };

      arg_x /= dis_denom(gen);

      const boost::decimal::decimal128_t trial { cbrt(arg_x) };
      const boost::decimal::decimal128_t ctrl  { boost::math::cbrt(arg_x) };

      result_is_ok = is_close_fraction(trial, ctrl, 4 * std::numeric_limits<boost::decimal::decimal128_t>::epsilon());

      if((!result_is_ok) || ((count % 512U) == 0U))
      {
        std::stringstream strm { };

        strm << std::setprecision(std::numeric_limits<boost::decimal::decimal128_t>::digits10)
             << "arg_x: " << arg_x
             << "\n  trial: "
             << trial
             << "\n  ctrl:  "
             << ctrl
             << "\n  is_ok: "
             << std::boolalpha
             << result_is_ok
             << "\n  count: "
             << count;

        std::cout << strm.str() << std::endl;
      }

      ++count;
    }
  }
} // namespace local

auto main() -> int
{
  local::test_cbrt128_versus_math();
}

[ckormanyos]

Test frexp no hang

#include <boost/decimal/decimal128_t.hpp>
#include <boost/decimal/cmath.hpp>
#include <boost/decimal/iostream.hpp>

#include <chrono>
#include <iomanip>
#include <iostream>
#include <random>

namespace local
{
  template<typename IntegralTimePointType,
           typename ClockType = std::chrono::high_resolution_clock>
  auto time_point() noexcept -> IntegralTimePointType
  {
    using local_integral_time_point_type = IntegralTimePointType;
    using local_clock_type               = ClockType;

    const auto current_now =
      static_cast<std::uintmax_t>
      (
        std::chrono::duration_cast<std::chrono::nanoseconds>
        (
          local_clock_type::now().time_since_epoch()
        ).count()
      );

    return static_cast<local_integral_time_point_type>(current_now);
  }

  auto test_frexp_nohang() -> void
  {
    std::random_device rd;
    std::mt19937_64 gen(rd());

    gen.seed(time_point<typename std::mt19937_64::result_type>());

    auto dis_denom =
      std::uniform_int_distribution<std::uint32_t>
      {
        UINT32_C(104744),
        UINT32_C(999999)
      };

    auto dis_exp10 =
      std::uniform_int_distribution<std::int32_t>
      {
        INT32_C(-4000),
        INT32_C(+4000)
      };

    bool result_is_ok { true };

    std::uint64_t count { };

    while(result_is_ok)
    {
      boost::decimal::decimal128_t arg_x { UINT32_C(104743), dis_exp10(gen) };

      arg_x /= dis_denom(gen);

      int nexp2 { };
      const boost::decimal::decimal128_t trial { frexp(arg_x, &nexp2) };

      if((count % 512U) == 0U)
      {
        std::stringstream strm { };

        strm << "count: " << count;

        std::cout << strm.str() << std::endl;
      }

      ++count;
    }
  }
} // namespace local

auto main() -> int
{
  local::test_frexp_nohang();
}

[ckormanyos]

See also #1330

[ckormanyos]

Test frexp no-hang for decimal128_t versus cpp_dec_float

#include <boost/decimal/decimal128_t.hpp>
#include <boost/decimal/cmath.hpp>
#include <boost/decimal/iostream.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

#include <chrono>
#include <iomanip>
#include <iostream>
#include <random>

namespace local
{
  template<typename IntegralTimePointType,
           typename ClockType = std::chrono::high_resolution_clock>
  auto time_point() noexcept -> IntegralTimePointType
  {
    using local_integral_time_point_type = IntegralTimePointType;
    using local_clock_type               = ClockType;

    const auto current_now =
      static_cast<std::uintmax_t>
      (
        std::chrono::duration_cast<std::chrono::nanoseconds>
        (
          local_clock_type::now().time_since_epoch()
        ).count()
      );

    return static_cast<local_integral_time_point_type>(current_now);
  }

  template<typename NumericType>
  auto is_close_fraction(const NumericType& a,
                         const NumericType& b,
                         const NumericType& tol) noexcept -> bool
  {
    using std::fabs;

    const NumericType delta
    {
      (b == static_cast<NumericType>(0)) ? fabs(a - b) : fabs(1 - (a / b))
    };

    const bool result_is_ok { (delta < tol) };

    return result_is_ok;
  }

  auto test_frexp_nohang() -> void
  {
    std::random_device rd;
    std::mt19937_64 gen(rd());

    gen.seed(time_point<typename std::mt19937_64::result_type>());

    auto dis_denom =
      std::uniform_int_distribution<std::uint32_t>
      {
        UINT32_C(104744),
        UINT32_C(999999)
      };

    auto dis_exp10 =
      std::uniform_int_distribution<std::int32_t>
      {
        INT32_C(-4000),
        INT32_C(+4000)
      };

    bool result_is_ok { true };

    std::uint64_t count { };

    for(;;)
    {
      using mp_type = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<34>, boost::multiprecision::et_off>;

      const auto exp10 { dis_exp10(gen) };
      const auto denom { dis_denom(gen) };

      boost::decimal::decimal128_t arg_x { UINT32_C(104743), exp10 }; arg_x /= denom;
      mp_type arg_mp { UINT32_C(104743) * pow(mp_type(10), exp10) }; arg_mp /= denom;

      int nexp2 { };
      const boost::decimal::decimal128_t trial { frexp(arg_x, &nexp2) };

      int n_mp { };
      const mp_type ctrl_mp { frexp(arg_mp, &n_mp) };

      std::string str_ctrl { };

      boost::decimal::decimal128_t ctrl_mp_as_dec { };

      {
        std::stringstream strm { };

        strm << std::setprecision(std::numeric_limits<boost::decimal::decimal128_t>::digits10) << ctrl_mp;

        ctrl_mp_as_dec = boost::decimal::decimal128_t(strm.str());
      }

      ++count;

      result_is_ok = { (nexp2 == n_mp) && is_close_fraction(trial, ctrl_mp_as_dec, std::numeric_limits<boost::decimal::decimal128_t>::epsilon() * 64) };

      {
        std::stringstream strm { };

        strm << std::boolalpha << "result_is_ok: " << result_is_ok << ", count: " << count;

        std::cout << strm.str() << std::endl;
      }

      if((count % 1024U) == 0U)
      {
        gen.seed(time_point<typename std::mt19937_64::result_type>());
      }

      if(!result_is_ok)
      {
        break;
      }
    }
  }
} // namespace local

auto main() -> int
{
  local::test_frexp_nohang();
}