/////////////////////////////////////////////////////////////////////////////// // Copyright 2021 - 2025 Fahad Syed. // Copyright 2021 - 2025 Christopher Kormanyos. // Copyright 2021 - 2025 Janek Kozicki. // Copyright 2025 Matt Borland. // Distributed under the Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) // #ifndef BOOST_MP_CPP_DOUBLE_FP_2021_06_05_HPP #define BOOST_MP_CPP_DOUBLE_FP_2021_06_05_HPP #include #include #include #include #include #ifdef BOOST_MP_MATH_AVAILABLE // // Headers required for Boost.Math integration: // #include // // Some includes we need from Boost.Math, since we rely on that library to provide these functions: // #include #include #include #include #include #include #endif #include #include #include #if (defined(BOOST_CLANG) && defined(BOOST_CLANG_VERSION) && (BOOST_CLANG_VERSION <= 90000)) #define BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE struct #else #define BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE class #endif namespace boost { namespace multiprecision { namespace backends { template class cpp_double_fp_backend; template constexpr auto eval_add(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template constexpr auto eval_add(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void; template constexpr auto eval_subtract(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template constexpr auto eval_subtract(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void; template constexpr auto eval_multiply(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template constexpr auto eval_multiply(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void; template constexpr auto eval_divide(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template constexpr auto eval_divide(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void; template constexpr auto eval_eq(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool; template constexpr auto eval_lt(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool; template constexpr auto eval_gt(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool; template constexpr auto eval_is_zero(const cpp_double_fp_backend& x) -> bool; template constexpr auto eval_signbit(const cpp_double_fp_backend& x) -> int; template constexpr auto eval_fabs(cpp_double_fp_backend& result, const cpp_double_fp_backend& a) -> void; template constexpr auto eval_frexp(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, int* v) -> void; template constexpr auto eval_ldexp(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, int v) -> void; template constexpr auto eval_floor(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template constexpr auto eval_ceil(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template constexpr auto eval_fpclassify(const cpp_double_fp_backend& o) -> int; template constexpr auto eval_sqrt(cpp_double_fp_backend& result, const cpp_double_fp_backend& o) -> void; template constexpr auto eval_pow(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, const cpp_double_fp_backend& a) -> void; template constexpr auto eval_pow(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, IntegralType n) -> typename ::std::enable_if::value, void>::type; template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) < 16))>::type const* = nullptr> constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template ::value && (((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) <= 36)))>::type const* = nullptr> constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) > 36))>::type const* = nullptr> constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template constexpr auto eval_log(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template constexpr auto eval_convert_to(signed long long* result, const cpp_double_fp_backend& backend) -> void; template constexpr auto eval_convert_to(unsigned long long* result, const cpp_double_fp_backend& backend) -> void; #ifdef BOOST_HAS_INT128 template constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits::digits > 24), void>::type; template constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if::digits > 24), void>::type; template constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits::digits > 24), void>::type; template constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if::digits > 24), void>::type; #endif template constexpr auto eval_convert_to(OtherFloatingPointType* result, const cpp_double_fp_backend& backend) -> typename ::std::enable_if::value>::type; template constexpr auto hash_value(const cpp_double_fp_backend& a) -> ::std::size_t; template constexpr auto fabs(const cpp_double_fp_backend& a) -> cpp_double_fp_backend; } } } // namespace boost::multiprecision::backends namespace std { // Foward declarations of various specializations of std::numeric_limits template BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE numeric_limits, ExpressionTemplatesOption> >; } // namespace std namespace boost { namespace multiprecision { template struct number_category> : public std::integral_constant { }; namespace backends { // A cpp_double_fp_backend is represented by an unevaluated sum of two // floating-point units, a0 and a1, which satisfy |a1| <= (1 / 2) * ulp(a0). // The type of the floating-point constituents should adhere to IEEE754. // Although the constituent parts (a0 and a1) satisfy |a1| <= (1 / 2) * ulp(a0), // the composite type does not adhere to these strict error bounds. Its error // bounds are larger. // This class has been tested with floats having single-precision (4 byte), // double-precision (8 byte) and quad precision (16 byte, such as GCC's __float128). template class cpp_double_fp_backend { public: using float_type = FloatingPointType; static_assert ( cpp_df_qf_detail::is_floating_point::value && bool { ((cpp_df_qf_detail::ccmath::numeric_limits::digits == 24) && cpp_df_qf_detail::ccmath::numeric_limits::is_specialized && cpp_df_qf_detail::ccmath::numeric_limits::is_iec559) || ((cpp_df_qf_detail::ccmath::numeric_limits::digits == 53) && cpp_df_qf_detail::ccmath::numeric_limits::is_specialized && cpp_df_qf_detail::ccmath::numeric_limits::is_iec559) || ((cpp_df_qf_detail::ccmath::numeric_limits::digits == 64) && cpp_df_qf_detail::ccmath::numeric_limits::is_specialized && cpp_df_qf_detail::ccmath::numeric_limits::is_iec559) || (cpp_df_qf_detail::ccmath::numeric_limits::digits == 113) }, "Error: float_type does not fulfill the backend requirements of cpp_double_fp_backend" ); using rep_type = cpp_df_qf_detail::pair; using arithmetic = cpp_df_qf_detail::exact_arithmetic; using signed_types = std::tuple; using unsigned_types = std::tuple; using float_types = std::tuple; using exponent_type = int; static constexpr int my_digits = static_cast(cpp_df_qf_detail::ccmath::numeric_limits::digits * static_cast(INT8_C(2))); static constexpr int my_digits10 = static_cast(boost::multiprecision::detail::calc_digits10 (my_digits)>::value); static constexpr int my_max_digits10 = static_cast(boost::multiprecision::detail::calc_max_digits10(my_digits)>::value); static constexpr int my_max_exponent = cpp_df_qf_detail::ccmath::numeric_limits::max_exponent; static constexpr int my_max_exponent10 = cpp_df_qf_detail::ccmath::numeric_limits::max_exponent10; static constexpr int my_min_exponent = static_cast ( cpp_df_qf_detail::ccmath::numeric_limits::min_exponent + cpp_df_qf_detail::ccmath::numeric_limits::digits ); static constexpr int my_min_exponent10 = static_cast ( -static_cast ( boost::multiprecision::detail::calc_digits10(-my_min_exponent)>::value ) ); // Default constructor. constexpr cpp_double_fp_backend() noexcept { } // Copy constructor. constexpr cpp_double_fp_backend(const cpp_double_fp_backend& other) : data(other.data) { } // Move constructor. constexpr cpp_double_fp_backend(cpp_double_fp_backend&& other) noexcept : data(static_cast(other.data)) { } // Constructors from other floating-point types. template ::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits <= cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> constexpr cpp_double_fp_backend(const OtherFloatType& f) : data(f, static_cast(0.0F)) { } template ::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits > cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> constexpr cpp_double_fp_backend(const OtherFloatType& f) : data(static_cast(f), static_cast(f - static_cast(static_cast(f)))) { } // Construtor from another kind of cpp_double_fp_backend<> object. template ::value && (!std::is_same::value))>::type const* = nullptr> constexpr cpp_double_fp_backend(const cpp_double_fp_backend& a) : cpp_double_fp_backend(cpp_double_fp_backend(a.my_first()) += a.my_second()) { } // Constructors from integers. template ::value && (!boost::multiprecision::detail::is_unsigned::value) && ((static_cast(sizeof(SignedIntegralType) * 8) - 1) <= cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> constexpr cpp_double_fp_backend(const SignedIntegralType& n) : data(static_cast(n), static_cast(0.0F)) { } template ::value && boost::multiprecision::detail::is_unsigned::value && (static_cast(sizeof(UnsignedIntegralType) * 8u) <= cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> constexpr cpp_double_fp_backend(const UnsignedIntegralType& u) : data(static_cast(u), static_cast(0.0F)) { } // Constructors from integers which hold more information than *this can contain. template ::value && boost::multiprecision::detail::is_unsigned::value && (static_cast(sizeof(UnsignedIntegralType) * 8u) > cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> constexpr cpp_double_fp_backend(UnsignedIntegralType u) : data(static_cast(u & cpp_df_qf_detail::float_mask()), static_cast(0.0F)) { using local_unsigned_integral_type = UnsignedIntegralType; if (u > cpp_df_qf_detail::float_mask()) { local_unsigned_integral_type local_flt_mask { cpp_df_qf_detail::float_mask() }; for (int index_mask_lsb = cpp_df_qf_detail::ccmath::numeric_limits::digits; (index_mask_lsb < static_cast(sizeof(local_unsigned_integral_type) * 8u)) && (local_flt_mask != local_unsigned_integral_type { UINT8_C(0) }); index_mask_lsb += cpp_df_qf_detail::ccmath::numeric_limits::digits) { local_flt_mask <<= static_cast(cpp_df_qf_detail::ccmath::numeric_limits::digits); add_unchecked_limb(static_cast(u & local_flt_mask)); } } } template ::value && (!boost::multiprecision::detail::is_unsigned::value) && ((static_cast(sizeof(SignedIntegralType) * 8) - 1) > cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> constexpr cpp_double_fp_backend(SignedIntegralType n) { const bool is_neg { (n < SignedIntegralType { INT8_C(0) }) }; using local_unsigned_integral_type = typename boost::multiprecision::detail::make_unsigned::type; const local_unsigned_integral_type u_val { (!is_neg) ? static_cast(n) : static_cast ( static_cast(~n) + static_cast(UINT8_C(1)) ) }; data = cpp_double_fp_backend(u_val).data; if(is_neg) { negate(); } } constexpr cpp_double_fp_backend(const float_type& a, const float_type& b) noexcept : data(a, b) { } constexpr cpp_double_fp_backend(const cpp_df_qf_detail::pair& p) noexcept : data(p) { } cpp_double_fp_backend(const char* p_str) { *this = p_str; } // Assignment operator. constexpr auto operator=(const cpp_double_fp_backend& other) -> cpp_double_fp_backend& { if (this != &other) { data = other.data; } return *this; } // Move assignment operator. constexpr auto operator=(cpp_double_fp_backend&& other) noexcept -> cpp_double_fp_backend& { data = static_cast(other.data); return *this; } // Assignment operator from another kind of cpp_double_fp_backend<> object. template ::value && (!std::is_same::value))>::type const* = nullptr> constexpr auto operator=(const cpp_double_fp_backend& other) -> cpp_double_fp_backend& { return operator=(cpp_double_fp_backend(other)); } template constexpr auto operator=(const OtherFloatType f) -> typename ::std::enable_if::value, cpp_double_fp_backend&>::type { return operator=(cpp_double_fp_backend(f)); } template constexpr auto operator=(const IntegralType n) -> typename ::std::enable_if::value, cpp_double_fp_backend&>::type { return operator=(cpp_double_fp_backend(n)); } auto operator=(const char* p_str) -> cpp_double_fp_backend& { rd_string(p_str); return *this; } constexpr auto hash() const -> ::std::size_t { #if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128) using local_float_type = typename std::conditional<::std::is_same::value, long double, float_type>::type; #else using local_float_type = float_type; #endif std::size_t result { UINT8_C(0) }; int n_first { }; int n_second { }; boost::multiprecision::detail::hash_combine(result, static_cast(cpp_df_qf_detail::ccmath::frexp(data.first, &n_first))); boost::multiprecision::detail::hash_combine(result, static_cast(cpp_df_qf_detail::ccmath::frexp(data.second, &n_second))); boost::multiprecision::detail::hash_combine(result, n_first); boost::multiprecision::detail::hash_combine(result, n_second); return result; } // The public methods follow. constexpr auto isneg_unchecked() const noexcept -> bool { return (data.first < 0); } constexpr auto iszero_unchecked() const noexcept -> bool { return (data.first == float_type { 0.0F }); } constexpr auto is_one() const noexcept -> bool { return ( (data.second == float_type { 0.0F }) && (data.first == float_type { 1.0F }) ); } constexpr auto negate() -> void { const bool isinf_u { (cpp_df_qf_detail::ccmath::isinf)(data.first) }; const bool isnan_u { (cpp_df_qf_detail::ccmath::isnan)(data.first) }; if (isnan_u) { } else if (isinf_u) { data.first = -data.first; } else { if (!iszero_unchecked()) { data = arithmetic::normalize(-data.first, -data.second); } } } // Getters/Setters constexpr auto my_first () const noexcept -> const float_type& { return data.first; } constexpr auto my_second() const noexcept -> const float_type& { return data.second; } constexpr auto rep() noexcept -> rep_type& { return data; } constexpr auto rep() const noexcept -> const rep_type& { return data; } constexpr auto crep() const noexcept -> const rep_type& { return data; } // Unary add/sub/mul/div follow in the upcoming paragraphs. constexpr auto operator+=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend& { const int fpc_u { eval_fpclassify(*this) }; const int fpc_v { eval_fpclassify(v) }; if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL)) { // Handle special cases like zero, inf and NaN. if (fpc_u == FP_NAN) { return *this; } const bool isinf_v { (fpc_v == FP_INFINITE) }; if (fpc_u == FP_INFINITE) { if (isinf_v && (isneg_unchecked() != v.isneg_unchecked())) { *this = cpp_double_fp_backend::my_value_nan(); } return *this; } const bool iszero_u { ((fpc_u == FP_ZERO) || (fpc_u == FP_SUBNORMAL)) }; const bool isnan_v { (fpc_v == FP_NAN) }; if (iszero_u || (isnan_v || isinf_v)) { if (iszero_u) { data.first = float_type { 0.0F }; data.second = float_type { 0.0F }; } const bool iszero_v { ((fpc_v == FP_ZERO) || (fpc_v == FP_SUBNORMAL)) }; return ((!iszero_v) ? operator=(v) : *this); } } add_unchecked(v); return *this; } constexpr auto operator-=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend& { const int fpc_u { eval_fpclassify(*this) }; const int fpc_v { eval_fpclassify(v) }; if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL)) { // Handle special cases like zero, inf and NaN. if (fpc_u == FP_NAN) { return *this; } const bool isinf_v { (fpc_v == FP_INFINITE) }; if (fpc_u == FP_INFINITE) { if (isinf_v && (isneg_unchecked() == v.isneg_unchecked())) { *this = cpp_double_fp_backend::my_value_nan(); } return *this; } const bool iszero_u { ((fpc_u == FP_ZERO) || (fpc_u == FP_SUBNORMAL)) }; const bool isnan_v { (fpc_v == FP_NAN) }; if (iszero_u || (isnan_v || isinf_v)) { if (iszero_u) { data.first = float_type { 0.0F }; data.second = float_type { 0.0F }; } const bool iszero_v { ((fpc_v == FP_ZERO) || (fpc_v == FP_SUBNORMAL)) }; return ((!iszero_v) ? operator=(-v) : *this); } } if (this == &v) { data.first = float_type { 0.0F }; data.second = float_type { 0.0F }; return *this; } const rep_type thi_tlo { arithmetic::two_diff(data.second, v.data.second) }; data = arithmetic::two_diff(data.first, v.data.first); if (cpp_df_qf_detail::ccmath::isinf(data.first)) { // Handle overflow. const bool b_neg { (data.first < float_type { 0.0F }) }; *this = cpp_double_fp_backend::my_value_inf(); if (b_neg) { negate(); } return *this; } data = arithmetic::two_hilo_sum(data.first, data.second + thi_tlo.first); data = arithmetic::two_hilo_sum(data.first, thi_tlo.second + data.second); return *this; } constexpr auto operator*=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend& { // Evaluate the sign of the result. const int fpc_u { eval_fpclassify(*this) }; const int fpc_v { eval_fpclassify(v) }; if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL)) { // Handle special cases like zero, inf and NaN. const bool isinf_u { (fpc_u == FP_INFINITE) }; const bool isinf_v { (fpc_v == FP_INFINITE) }; const bool iszero_u { (fpc_u == FP_ZERO) }; const bool iszero_v { (fpc_v == FP_ZERO) }; if (((fpc_u == FP_NAN) || (fpc_v == FP_NAN)) || (isinf_u && iszero_v) || (isinf_v && iszero_u)) { return operator=( cpp_double_fp_backend::my_value_nan()); } if (isinf_u || isinf_v) { const bool b_neg { (isneg_unchecked() != v.isneg_unchecked()) }; *this = cpp_double_fp_backend::my_value_inf(); if (b_neg) { negate(); } return *this; } if (iszero_u || iszero_v) { return operator=(cpp_double_fp_backend(0)); } } mul_unchecked(v); return *this; } constexpr auto operator/=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend& { const int fpc_u { eval_fpclassify(*this) }; const int fpc_v { eval_fpclassify(v) }; if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL)) { // Handle special cases like zero, inf and NaN. const bool isnan_u { (fpc_u == FP_NAN) }; const bool isnan_v { (fpc_v == FP_NAN) }; if (isnan_u || isnan_v) { return operator=(cpp_double_fp_backend::my_value_nan()); } const bool iszero_u { (fpc_u == FP_ZERO) }; const bool iszero_v { (fpc_v == FP_ZERO) }; if (iszero_u) { if (iszero_v) { return operator=(cpp_double_fp_backend::my_value_nan()); } else { return operator=(cpp_double_fp_backend(0)); } } // Handle more special cases like zero, inf and NaN. if (iszero_v) { const bool b_neg = isneg_unchecked(); *this = cpp_double_fp_backend::my_value_inf(); if (b_neg) { negate(); } return *this; } const bool isinf_v { (fpc_v == FP_INFINITE) }; const bool isinf_u { (fpc_u == FP_INFINITE) }; if (isinf_u) { if (isinf_v) { return operator=(cpp_double_fp_backend::my_value_nan()); } else { const bool b_neg { isneg_unchecked() }; return operator=((!b_neg) ? cpp_double_fp_backend::my_value_inf() : -cpp_double_fp_backend::my_value_inf()); } } if (isinf_v) { return operator=(cpp_double_fp_backend(0)); } } if (this == &v) { data.first = float_type { 1.0F }; data.second = float_type { 0.0F }; return *this; } // The division algorithm has been taken from Victor Shoup, // package WinNTL-5_3_2. It might originally be related to the // K. Briggs work. The algorithm has been significantly simplified // while still attempting to retain proper rounding corrections. // Checks for overflow and underflow have been added. const float_type C { data.first / v.data.first }; float_type c { cpp_df_qf_detail::split_maker::value * C }; float_type hc { }; if (cpp_df_qf_detail::ccmath::isinf(c)) { // Handle overflow by scaling down (and then back up) with the split. hc = cpp_df_qf_detail::ccmath::ldexp ( C - float_type { C - cpp_df_qf_detail::ccmath::ldexp(C, -cpp_df_qf_detail::split_maker::n_shl) }, cpp_df_qf_detail::split_maker::n_shl ); } else { hc = c - float_type { c - C }; } float_type u { cpp_df_qf_detail::split_maker::value * v.data.first }; const float_type hv = ( cpp_df_qf_detail::ccmath::isinf(u) ? cpp_df_qf_detail::ccmath::ldexp ( // Handle overflow by scaling down (and then back up) with the split. v.data.first - float_type { v.data.first - cpp_df_qf_detail::ccmath::ldexp(v.data.first, -cpp_df_qf_detail::split_maker::n_shl) }, cpp_df_qf_detail::split_maker::n_shl ) : u - float_type { u - v.data.first } ); const float_type U { C * v.data.first }; u = cpp_df_qf_detail::ccmath::unsafe::fma(hc, hv, -U); { const float_type tv { v.data.first - hv }; u = cpp_df_qf_detail::ccmath::unsafe::fma(hc, tv, u); const float_type tc { C - hc }; u = cpp_df_qf_detail::ccmath::unsafe::fma(tc, hv, u); u = cpp_df_qf_detail::ccmath::unsafe::fma(tc, tv, u); } c = float_type { (data.first - U) - u } + data.second; c = (c - float_type { C * v.data.second }) / v.data.first; // Perform even more simplifications compared to Victor Shoup. data.first = C + c; data.second = float_type { C - data.first } + c; return *this; } // Unary minus operator. constexpr auto operator-() const -> cpp_double_fp_backend { cpp_double_fp_backend v { *this }; v.negate(); return v; } // Public helper functions. static constexpr auto pown(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, int p) -> void { using local_float_type = cpp_double_fp_backend; if (p == 2) { result = x; result.mul_unchecked(x); } else if (p == 3) { result = x; result.mul_unchecked(x); result.mul_unchecked(x); } else if (p == 4) { local_float_type x2 { x }; x2.mul_unchecked(x); result = x2; result.mul_unchecked(x2); } else if (p == 1) { result = x; } else if (p < 0) { // The case p == 0 is checked in a higher calling layer. pown(result, local_float_type(1U) / x, -p); } else { result = local_float_type(1U); local_float_type y(x); auto p_local = static_cast(p); for (;;) { if (static_cast(p_local & static_cast(UINT8_C(1))) != static_cast(UINT8_C(0))) { result.mul_unchecked(y); } p_local >>= 1U; if (p_local == static_cast(UINT8_C(0))) { break; } else { y.mul_unchecked(y); } } } } constexpr auto swap(cpp_double_fp_backend& other) -> void { if (this != &other) { const rep_type tmp { data }; data = other.data; other.data = tmp; } } constexpr auto swap(cpp_double_fp_backend&& other) noexcept -> void { const rep_type tmp { static_cast(data) }; data = other.data; other.data = tmp; } constexpr auto compare(const cpp_double_fp_backend& other) const noexcept -> int { // Return 1 for *this > other, -1 for *this < other, 0 for *this = other. if ((cpp_df_qf_detail::ccmath::isnan)(data.first)) { return -1; } else { return (my_first() > other.my_first()) ? 1 : (my_first() < other.my_first()) ? -1 : (my_second() > other.my_second()) ? 1 : (my_second() < other.my_second()) ? -1 : 0; } } // TBD: Exactly what compilers/language-standards are needed to make this constexpr? // TBD: It odes not really become constexpr until we stop using an intermediate // cpp_bin_float anyway. But I will leave this comment for future library evolution. auto str(std::streamsize number_of_digits, const std::ios::fmtflags format_flags) const -> std::string { // Use cpp_bin_float when writing to string. This is similar // to the use of cpp_bin_float when reading from string. cpp_bin_float_read_write_backend_type f_bin { data.first }; eval_add(f_bin, cpp_bin_float_read_write_backend_type(data.second)); return f_bin.str(number_of_digits, format_flags); } constexpr auto order02() const -> int { // TBD: Is there another option to get the base-2 log // that's more unequivocally closer to constexpr? // TBD: Either way, integrate this (or something like it) // into any potential implementation of eval_ilogb(). int e2 { }; cpp_double_fp_backend dummy { }; eval_frexp(dummy, *this, &e2); return e2; } static constexpr auto my_value_max() noexcept -> cpp_double_fp_backend { // Use the non-normalized sum of two maximum values, where the lower // value is "shifted" right in the sense of floating-point ldexp. return cpp_double_fp_backend ( arithmetic::two_hilo_sum ( float_type ( (cpp_df_qf_detail::ccmath::numeric_limits::max)() * ( static_cast(1.0F) - static_cast(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits::epsilon()) ) ), float_type ( cpp_df_qf_detail::ccmath::ldexp ( (cpp_df_qf_detail::ccmath::numeric_limits::max)(), -cpp_df_qf_detail::ccmath::numeric_limits::digits ) ) ) ); } static constexpr auto my_value_min() noexcept -> cpp_double_fp_backend { // Use the non-normalized minimum value, where the lower value // is "shifted" left in the sense of floating-point ldexp. return cpp_double_fp_backend ( float_type ( cpp_df_qf_detail::ccmath::ldexp ( (cpp_df_qf_detail::ccmath::numeric_limits::min)(), cpp_df_qf_detail::ccmath::numeric_limits::digits ) ) ); } static constexpr auto my_value_eps() noexcept -> cpp_double_fp_backend { return cpp_double_fp_backend ( float_type(cpp_df_qf_detail::ccmath::ldexp(float_type { 1 }, int { 3 - my_digits })) ); } static constexpr auto my_value_nan() noexcept -> cpp_double_fp_backend { return cpp_double_fp_backend(static_cast(NAN), static_cast(0.0F)); } static constexpr auto my_value_inf() noexcept -> cpp_double_fp_backend { return cpp_double_fp_backend(static_cast(HUGE_VAL), static_cast(0.0F)); // conversion from double infinity OK } static constexpr auto my_value_logmax() -> cpp_double_fp_backend { return cpp_double_fp_backend ( cpp_df_qf_detail::ccmath::log ( float_type ( (cpp_df_qf_detail::ccmath::numeric_limits::max)() * ( static_cast(1.0F) - static_cast(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits::epsilon()) ) ) ) ); } static constexpr auto my_value_logmin() -> cpp_double_fp_backend { return cpp_double_fp_backend ( cpp_df_qf_detail::ccmath::log ( float_type ( cpp_df_qf_detail::ccmath::ldexp ( (cpp_df_qf_detail::ccmath::numeric_limits::min)(), cpp_df_qf_detail::ccmath::numeric_limits::digits ) ) ) ); } constexpr auto add_unchecked_limb(const float_type v_first) -> void { const float_type thi { data.second }; data = arithmetic::two_sum(data.first, v_first); data = arithmetic::two_hilo_sum(data.first, data.second + thi); data = arithmetic::two_hilo_sum(data.first, data.second); } private: rep_type data; using cpp_bin_float_read_write_backend_type = boost::multiprecision::backends::cpp_bin_float(my_digits), digit_base_2, void, int, cpp_df_qf_detail::ccmath::numeric_limits::min_exponent, cpp_df_qf_detail::ccmath::numeric_limits::max_exponent>; constexpr auto rd_string(const char* p_str) -> bool; constexpr auto add_unchecked(const cpp_double_fp_backend& v) -> void { const rep_type thi_tlo { arithmetic::two_sum(data.second, v.data.second) }; data = arithmetic::two_sum(data.first, v.data.first); if (cpp_df_qf_detail::ccmath::isinf(data.first)) { // Handle overflow. const bool b_neg { (data.first < float_type { 0.0F }) }; *this = cpp_double_fp_backend::my_value_inf(); if (b_neg) { negate(); } return; } data = arithmetic::two_hilo_sum(data.first, data.second + thi_tlo.first); data = arithmetic::two_hilo_sum(data.first, thi_tlo.second + data.second); } constexpr auto mul_unchecked(const cpp_double_fp_backend& v) -> void { // The multiplication algorithm has been taken from Victor Shoup, // package WinNTL-5_3_2. It might originally be related to the // K. Briggs work. The algorithm has been significantly simplified // while still attempting to retain proper rounding corrections. // Checks for overflow and underflow have been added. float_type C { cpp_df_qf_detail::split_maker::value * data.first }; float_type hu { }; if (cpp_df_qf_detail::ccmath::isinf(C)) { // Handle overflow by scaling down (and then back up) with the split. C = data.first - cpp_df_qf_detail::ccmath::ldexp(data.first, -cpp_df_qf_detail::split_maker::n_shl); hu = cpp_df_qf_detail::ccmath::ldexp(data.first - C, cpp_df_qf_detail::split_maker::n_shl); } else { hu = C - float_type { C - data.first }; } C = data.first * v.data.first; if (cpp_df_qf_detail::ccmath::isinf(C)) { // Handle overflow. const bool b_neg { (isneg_unchecked() != v.isneg_unchecked()) }; *this = cpp_double_fp_backend::my_value_inf(); if (b_neg) { negate(); } return; } float_type c { cpp_df_qf_detail::split_maker::value * v.data.first }; float_type hv { }; if (cpp_df_qf_detail::ccmath::isinf(c)) { // Handle overflow by scaling down (and then back up) with the split. c = v.data.first - cpp_df_qf_detail::ccmath::ldexp(v.data.first, -cpp_df_qf_detail::split_maker::n_shl); hv = cpp_df_qf_detail::ccmath::ldexp(v.data.first - c, cpp_df_qf_detail::split_maker::n_shl); } else { hv = c - float_type { c - v.data.first }; } { const float_type tv { v.data.first - hv }; const float_type t1 { cpp_df_qf_detail::ccmath::unsafe::fma ( hu, tv, cpp_df_qf_detail::ccmath::unsafe::fma(hu, hv, -C) ) }; const float_type tu { data.first - hu }; c = cpp_df_qf_detail::ccmath::unsafe::fma(tu, tv, cpp_df_qf_detail::ccmath::unsafe::fma(tu, hv, t1)) + (data.first * v.data.second) + (data.second * v.data.first); } // Perform even more simplifications compared to Victor Shoup. data.first = C + c; data.second = float_type { C - data.first } + c; } template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) < 16))>::type const*> friend constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template ::value && (((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) <= 36)))>::type const*> friend constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) > 36))>::type const*> friend constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; }; template constexpr auto cpp_double_fp_backend::rd_string(const char* p_str) -> bool { // Use an intermediate cpp_bin_float backend type for reading string input. cpp_bin_float_read_write_backend_type f_bin { }; f_bin = p_str; const int fpc { eval_fpclassify(f_bin) }; const bool is_definitely_nan { (fpc == FP_NAN) }; using local_double_fp_type = cpp_double_fp_backend; if (is_definitely_nan) { static_cast(operator=(local_double_fp_type::my_value_nan())); return true; } const bool b_neg { (eval_signbit(f_bin) == 1) }; if (b_neg) { f_bin.negate(); } const int expval_from_f_bin { [&f_bin]() { int expval { }; cpp_bin_float_read_write_backend_type dummy { }; eval_frexp(dummy, f_bin, &expval); return expval; }() }; const auto is_zero_or_subnormal = ( (fpc == FP_ZERO) || (expval_from_f_bin < static_cast(local_double_fp_type::my_min_exponent)) ); if (is_zero_or_subnormal) { data.first = float_type { 0.0F }; data.second = float_type { 0.0F }; return true; } float_type flt_inf_check_first { }; eval_convert_to(&flt_inf_check_first, f_bin); bool is_definitely_inf { ((fpc == FP_INFINITE) || cpp_df_qf_detail::ccmath::isinf(flt_inf_check_first)) }; if ((!is_definitely_inf) && (flt_inf_check_first > my_value_max().my_first())) { cpp_bin_float_read_write_backend_type f_bin_inf_check(f_bin); eval_subtract(f_bin_inf_check, cpp_bin_float_read_write_backend_type(flt_inf_check_first)); float_type flt_inf_check_second { }; eval_convert_to(&flt_inf_check_second, f_bin_inf_check); is_definitely_inf = eval_gt(local_double_fp_type(flt_inf_check_first, flt_inf_check_second), my_value_max()); }; if (is_definitely_inf) { static_cast(operator=(local_double_fp_type::my_value_inf())); if (b_neg) { negate(); } return true; } // The input string is normal. We will now extract its value. data.first = float_type { 0.0F }; data.second = float_type { 0.0F }; // Handle small input values. Scale (and re-scale them below) if possible. constexpr int pow2_scaling_for_small_input { cpp_df_qf_detail::ccmath::numeric_limits::digits }; const bool has_pow2_scaling_for_small_input { (expval_from_f_bin < static_cast(local_double_fp_type::my_min_exponent + pow2_scaling_for_small_input)) }; if (has_pow2_scaling_for_small_input) { eval_ldexp(f_bin, f_bin, pow2_scaling_for_small_input); } using local_builtin_float_type = typename std::conditional<(sizeof(float_type) <= sizeof(double)), double, float_type>::type; constexpr unsigned digit_limit { static_cast ( static_cast ( (local_double_fp_type::my_digits / cpp_df_qf_detail::ccmath::numeric_limits::digits) + (((local_double_fp_type::my_digits % cpp_df_qf_detail::ccmath::numeric_limits::digits) != 0) ? 1 : 0) ) * cpp_df_qf_detail::ccmath::numeric_limits::digits ) }; for(auto i = static_cast(UINT8_C(0)); i < digit_limit; i = static_cast(i + static_cast(cpp_df_qf_detail::ccmath::numeric_limits::digits))) { local_builtin_float_type flt_part { }; eval_convert_to(&flt_part, f_bin); eval_subtract(f_bin, cpp_bin_float_read_write_backend_type(flt_part)); eval_add(*this, local_double_fp_type { flt_part }); } if (has_pow2_scaling_for_small_input) { eval_ldexp(*this, *this, -pow2_scaling_for_small_input); } if (b_neg) { negate(); } return true; } } } } // namespace boost::multiprecision::backends #include namespace boost { namespace multiprecision { namespace backends { template constexpr int cpp_double_fp_backend::my_digits; template constexpr int cpp_double_fp_backend::my_digits10; template constexpr int cpp_double_fp_backend::my_max_digits10; template constexpr int cpp_double_fp_backend::my_max_exponent; template constexpr int cpp_double_fp_backend::my_min_exponent; template constexpr int cpp_double_fp_backend::my_max_exponent10; template constexpr int cpp_double_fp_backend::my_min_exponent10; template constexpr cpp_double_fp_backend operator+(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) { return cpp_double_fp_backend(a) += b; } template constexpr cpp_double_fp_backend operator-(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) { return cpp_double_fp_backend(a) -= b; } template constexpr cpp_double_fp_backend operator*(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) { return cpp_double_fp_backend(a) *= b; } template constexpr cpp_double_fp_backend operator/(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) { return cpp_double_fp_backend(a) /= b; } template constexpr auto eval_add(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { result += x; } template constexpr auto eval_add(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void { result = cpp_double_fp_backend(a) += b; } template constexpr auto eval_subtract(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { result -= x; } template constexpr auto eval_subtract(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void { result = cpp_double_fp_backend(a) -= b; } template constexpr auto eval_multiply(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { result *= x; } template constexpr auto eval_multiply(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void { result = cpp_double_fp_backend(a) *= b; } template constexpr auto eval_divide(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { result /= x; } template constexpr auto eval_divide(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void { result = cpp_double_fp_backend(a) /= b; } template constexpr auto eval_eq(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool { return (a.compare(b) == 0); } template constexpr auto eval_lt(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool { return (a.compare(b) == -1); } template constexpr auto eval_gt(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool { return (a.compare(b) == 1); } template constexpr auto eval_is_zero(const cpp_double_fp_backend& x) -> bool { return x.iszero_unchecked(); } template constexpr auto eval_signbit(const cpp_double_fp_backend& x) -> int { return (x.isneg_unchecked() ? 1 : 0); } template constexpr auto eval_fabs(cpp_double_fp_backend& result, const cpp_double_fp_backend& a) -> void { result = a; if (a.isneg_unchecked()) { result.negate(); } } template constexpr auto eval_frexp(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, int* v) -> void { using local_backend_type = cpp_double_fp_backend; using local_float_type = typename local_backend_type::float_type; int expptr { }; const local_float_type fhi { cpp_df_qf_detail::ccmath::frexp(a.rep().first, &expptr) }; const local_float_type flo { cpp_df_qf_detail::ccmath::ldexp(a.rep().second, -expptr) }; if (v != nullptr) { *v = expptr; } result.rep() = local_backend_type::arithmetic::normalize(fhi, flo); } template constexpr auto eval_ldexp(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, int v) -> void { using local_backend_type = cpp_double_fp_backend; using local_float_type = typename local_backend_type::float_type; const local_float_type fhi { cpp_df_qf_detail::ccmath::ldexp(a.crep().first, v) }; const local_float_type flo { cpp_df_qf_detail::ccmath::ldexp(a.crep().second, v) }; result.rep() = local_backend_type::arithmetic::normalize(fhi, flo); } template constexpr auto eval_floor(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { using local_backend_type = cpp_double_fp_backend; using local_float_type = typename local_backend_type::float_type; const local_float_type fhi { cpp_df_qf_detail::ccmath::floor(x.my_first()) }; if (fhi != x.my_first()) { result.rep().first = fhi; result.rep().second = local_float_type { 0 }; } else { const local_float_type flo = { cpp_df_qf_detail::ccmath::floor(x.my_second()) }; result.rep() = local_backend_type::arithmetic::normalize(fhi, flo); } } template constexpr auto eval_ceil(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { // Compute -floor(-x); eval_floor(result, -x); result.negate(); } template constexpr auto eval_fpclassify(const cpp_double_fp_backend& o) -> int { if (cpp_df_qf_detail::ccmath::isnan(o.crep().first)) { return FP_NAN; } else if (cpp_df_qf_detail::ccmath::isinf(o.crep().first)) { return FP_INFINITE; } else if (eval_is_zero(o)) { return FP_ZERO; } else { using local_backend_type = cpp_double_fp_backend; using local_float_type = typename local_backend_type::float_type; const local_float_type fabs_x { cpp_df_qf_detail::ccmath::fabs(o.crep().first) }; if ((fabs_x > 0) && (fabs_x < (cpp_df_qf_detail::ccmath::numeric_limits::min)())) { return FP_SUBNORMAL; } else { return FP_NORMAL; } } } template constexpr auto eval_sqrt(cpp_double_fp_backend& result, const cpp_double_fp_backend& o) -> void { using double_float_type = cpp_double_fp_backend; using local_float_type = typename double_float_type::float_type; const int fpc { eval_fpclassify(o) }; const bool isneg_o { o.isneg_unchecked() }; if ((fpc != FP_NORMAL) || isneg_o) { if (fpc == FP_ZERO) { result = double_float_type(0); return; } else if ((fpc == FP_NAN) || isneg_o) { result = double_float_type::my_value_nan(); return; } else if (fpc == FP_INFINITE) { result = double_float_type::my_value_inf(); return; } } // TBD: Do we need any overflow/underflow guards when multiplying // by the split or when multiplying (hx * tx) and/or (hx * hx)? const local_float_type c { cpp_df_qf_detail::ccmath::sqrt(o.crep().first) }; local_float_type p { cpp_df_qf_detail::split_maker::value * c }; const local_float_type hx { local_float_type { c - p } + p }; const local_float_type tx { c - hx }; local_float_type q = hx * tx; q = q + q; p = hx * hx; const local_float_type u { p + q }; const local_float_type uu { cpp_df_qf_detail::ccmath::unsafe::fma(tx, tx, local_float_type { p - u } + q) }; const local_float_type cc { local_float_type { local_float_type { o.crep().first - u } - uu + o.crep().second } / local_float_type { c + c } }; result.rep().first = c + cc; result.rep().second = local_float_type { c - result.my_first() } + cc; } template constexpr auto eval_pow(cpp_double_fp_backend