/////////////////////////////////////////////////////////////////////////////// // 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(__clang__) && (__clang_major__ <= 9)) #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; // TBD: constexpr on hash_value? template 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 (say a0 and a1) which satisfy |a1| <= (1 / 2) * ulp(a0). // The type of the floating-point constituents should adhere to IEEE754. 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::digits == 53) || (cpp_df_qf_detail::ccmath::numeric_limits::digits == 64) || (cpp_df_qf_detail::ccmath::numeric_limits::digits == 113) }, "Error: float_type does not fulfil 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(a.my_first()) { // TBD: Maybe specialize this constructor for cases either wider or less wide. operator+=(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) { using local_unsigned_integral_type = UnsignedIntegralType; constexpr local_unsigned_integral_type flt_mask { static_cast ( static_cast ( static_cast(UINT8_C(1)) << static_cast(std::numeric_limits::digits) ) - static_cast(UINT8_C(1)) ) }; data.second = static_cast(0.0F); data.first = static_cast(u & flt_mask); int u_index { 1 }; while (u > static_cast(UINT8_C(0))) { u >>= static_cast(std::numeric_limits::digits); const float_type xhi { static_cast(u & flt_mask) * cpp_df_qf_detail::pow2_maker::value(std::numeric_limits::digits * u_index) }; eval_add(*this, cpp_double_fp_backend(xhi)); ++u_index; } } 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* s) { *this = s; } // 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* v) -> cpp_double_fp_backend& { rd_string(v); return *this; } auto hash() const -> ::std::size_t { // Hash the raw values of the data field with direct-memory access. // Use 16-bit (2 byte) chunks as the data size when hashing. static_assert( ( sizeof(data.first) == sizeof(data.second)) && ( sizeof(float_type) >= sizeof(std::uint16_t)) && ((sizeof(float_type) % sizeof(std::uint16_t)) == std::size_t { UINT8_C(0) }), "Error: float_type size is inappropriate for hashing routine"); auto hash_one { [](std::size_t& res, const float_type& val) { const std::uint16_t* first { reinterpret_cast(&val) }; const std::uint16_t* last { first + std::size_t { sizeof(float_type) / sizeof(std::uint16_t) } }; while (first != last) { boost::multiprecision::detail::hash_combine(res, *first); ++first; } return res; } }; std::size_t result { UINT8_C(0) }; static_cast(hash_one(result, data.first)); static_cast(hash_one(result, data.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 { data.first = -data.first; if (!iszero_unchecked()) { data.second = -data.second; 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); } } const float_type xlo { data.second }; data = arithmetic::two_sum(data.first, v.data.first); const rep_type thi_tlo { arithmetic::two_sum(xlo, v.data.second) }; 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& { 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 float_type xlo { data.second }; data = arithmetic::two_diff(data.first, v.data.first); const rep_type thi_tlo { arithmetic::two_diff(xlo, v.data.second) }; 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 atempting to retain proper rounding corrections. 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 }; float_type hv { }; if (cpp_df_qf_detail::ccmath::isinf(u)) { // Handle overflow by scaling down (and then back up) with the split. hv = cpp_df_qf_detail::ccmath::ldexp ( 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 ); } else { hv = 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; } } auto str(std::streamsize number_of_digits, const std::ios::fmtflags format_flags) const -> std::string { if (number_of_digits == 0) { number_of_digits = cpp_double_fp_backend::my_digits10; } // 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_type f_bin { data.first }; f_bin += 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. constexpr float_type hi_part { (cpp_df_qf_detail::ccmath::numeric_limits::max)() * ( static_cast(1.0F) - static_cast(1.5F) * cpp_df_qf_detail::ccmath::unsafe::sqrt(cpp_df_qf_detail::ccmath::numeric_limits::epsilon()) ) }; constexpr float_type lo_part { cpp_df_qf_detail::ccmath::unsafe::ldexp ( (cpp_df_qf_detail::ccmath::numeric_limits::max)(), -cpp_df_qf_detail::ccmath::numeric_limits::digits ) }; constexpr cpp_double_fp_backend my_value_max_constexpr { arithmetic::two_hilo_sum(hi_part, lo_part) }; static_assert ( eval_gt(my_value_max_constexpr, cpp_double_fp_backend { hi_part }), "Error: maximum value is too small in relation to the maximum of its constituent type" ); return my_value_max_constexpr; } 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. constexpr cpp_double_fp_backend my_value_min_constexpr { cpp_df_qf_detail::ccmath::unsafe::ldexp ( (cpp_df_qf_detail::ccmath::numeric_limits::min)(), cpp_df_qf_detail::ccmath::numeric_limits::digits ) }; static_assert ( eval_gt(my_value_min_constexpr, cpp_double_fp_backend { (cpp_df_qf_detail::ccmath::numeric_limits::min)() }), "Error: minimum value is too small and must exceed the min of its constituent type" ); return my_value_min_constexpr; } static constexpr auto my_value_eps() noexcept -> cpp_double_fp_backend { constexpr cpp_double_fp_backend my_value_eps_constexpr { cpp_df_qf_detail::ccmath::unsafe::ldexp(float_type { 1 }, int { 3 - my_digits }) }; static_assert ( eval_lt(cpp_double_fp_backend { 1 } - my_value_eps_constexpr, cpp_double_fp_backend { 1 }), "Error: epsilon value is too small and must be large enough to differentiate (1 - epsilon) from 1" ); return my_value_eps_constexpr; } static constexpr auto my_value_nan() noexcept -> cpp_double_fp_backend { return cpp_double_fp_backend(static_cast(NAN), float_type { 0.0F }); } static constexpr auto my_value_inf() noexcept -> cpp_double_fp_backend { return cpp_double_fp_backend(static_cast(HUGE_VAL), float_type { 0.0F }); // conversion from double infinity OK } 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>; using cpp_bin_float_read_write_exp_type = typename cpp_bin_float_read_write_backend_type::exponent_type; using cpp_bin_float_read_write_type = boost::multiprecision::number; auto rd_string(const char* pstr) -> bool; 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 atempting to retain proper rounding corrections. 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 }; float_type t1 { cpp_df_qf_detail::ccmath::unsafe::fma(hu, hv, -C) }; t1 = cpp_df_qf_detail::ccmath::unsafe::fma(hu, tv, t1); const float_type tu { data.first - hu }; t1 = cpp_df_qf_detail::ccmath::unsafe::fma(tu, hv, t1); c = cpp_df_qf_detail::ccmath::unsafe::fma(tu, tv, 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 auto cpp_double_fp_backend::rd_string(const char* pstr) -> bool { using local_double_fp_type = cpp_double_fp_backend; cpp_bin_float_read_write_type f_bin { pstr }; const auto fpc = fpclassify(f_bin); const auto is_definitely_nan = (fpc == FP_NAN); if (is_definitely_nan) { static_cast(operator=(local_double_fp_type::my_value_nan())); } else { const bool b_neg { (signbit(f_bin) == 1) }; if (b_neg) { f_bin = -f_bin; } const int expval_from_f_bin { [&f_bin]() { int expval; frexp(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 }; } else { const auto is_definitely_inf = (f_bin > (::std::numeric_limits::max)()); if (is_definitely_inf) { static_cast(operator=(local_double_fp_type::my_value_inf())); if (b_neg) { negate(); } } else { data.first = float_type { 0.0F }; data.second = float_type { 0.0F }; using local_builtin_float_type = typename std::conditional<(sizeof(float_type) <= sizeof(double)), double, float_type>::type; constexpr int pow2_scaling_for_small_input { cpp_df_qf_detail::ccmath::numeric_limits::digits }; const auto 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) { f_bin = ldexp(f_bin, pow2_scaling_for_small_input); } constexpr auto dig_lim = 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 < dig_lim; i = static_cast(i + static_cast(cpp_df_qf_detail::ccmath::numeric_limits::digits))) { const local_builtin_float_type f { static_cast(f_bin) }; f_bin -= cpp_bin_float_read_write_type { f }; eval_add(*this, local_double_fp_type { f }); } 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.rep().first } + cc; } template constexpr auto eval_pow(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, const cpp_double_fp_backend& a) -> void { using double_float_type = cpp_double_fp_backend; constexpr double_float_type zero { 0 }; result = zero; signed long long val_nll { }; eval_convert_to(&val_nll, a); const int na { static_cast(val_nll) }; const int fpc_a { eval_fpclassify(a) }; if((fpc_a == FP_NORMAL) && (a.compare(double_float_type(na)) == 0)) { eval_pow(result, x, na); } else { constexpr double_float_type one { 1 }; const int fpc_x { eval_fpclassify(x) }; const int fpc_a { eval_fpclassify(a) }; if (fpc_a == FP_ZERO) { // pow(base, +/-0) returns 1 for any base, even when base is NaN. result = one; } else if (fpc_a == FP_NAN) { result = double_float_type::my_value_nan(); } else if (fpc_x == FP_ZERO) { if ((fpc_a == FP_NORMAL) || (fpc_a == FP_INFINITE)) { // pow(+/-0, exp), where exp is negative and finite, returns +infinity. // pow(+/-0, exp), where exp is positive non-integer, returns +0. // pow(+/-0, -infinity) returns +infinity. // pow(+/-0, +infinity) returns +0. result = (eval_signbit(a) ? double_float_type::my_value_inf() : zero); } } else if (fpc_x == FP_INFINITE) { if ((fpc_a == FP_NORMAL) || (fpc_a == FP_INFINITE)) { // pow(+infinity, exp) returns +0 for any negative exp. // pow(-infinity, exp) returns +infinity for any positive exp. result = (eval_signbit(a) ? zero : double_float_type::my_value_inf()); } } else if (fpc_x != FP_NORMAL) { result = x; } else { if (fpc_a == FP_INFINITE) { result = ( (x.compare(one) == -1) ? (eval_signbit(a) ? double_float_type::my_value_inf() : zero) : (x.compare(one) == +1) ? (eval_signbit(a) ? zero : double_float_type::my_value_inf()) : one ); } else { double_float_type log_x { }; eval_log(log_x, x); double_float_type a_log_x { }; eval_multiply(a_log_x, a, log_x); eval_exp(result, a_log_x); } } } } template constexpr auto eval_pow(cpp_double_fp_backend& result, const cpp_double_fp_backend