boost/multiprecision/cpp_double_fp.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //  Copyright 2021 - 2025 Fahad Syed.
0003 //  Copyright 2021 - 2025 Christopher Kormanyos.
0004 //  Copyright 2021 - 2025 Janek Kozicki.
0005 //  Copyright 2025 Matt Borland.
0006 //  Distributed under the Boost Software License, Version 1.0.
0007 //  (See accompanying file LICENSE_1_0.txt or copy at
0008 //  http://www.boost.org/LICENSE_1_0.txt)
0009 //
0010
0011 #ifndef BOOST_MP_CPP_DOUBLE_FP_2021_06_05_HPP
0012 #define BOOST_MP_CPP_DOUBLE_FP_2021_06_05_HPP
0013
0014 #include <boost/multiprecision/cpp_bin_float.hpp>
0015 #include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail.hpp>
0016 #include <boost/multiprecision/detail/hash.hpp>
0017 #include <boost/multiprecision/traits/max_digits10.hpp>
0018 #include <boost/multiprecision/traits/std_integer_traits.hpp>
0019
0020 #ifdef BOOST_MP_MATH_AVAILABLE
0021 //
0022 // Headers required for Boost.Math integration:
0023 //
0024 #include <boost/math/policies/policy.hpp>
0025 //
0026 // Some includes we need from Boost.Math, since we rely on that library to provide these functions:
0027 //
0028 #include <boost/math/special_functions/acosh.hpp>
0029 #include <boost/math/special_functions/asinh.hpp>
0030 #include <boost/math/special_functions/atanh.hpp>
0031 #include <boost/math/special_functions/cbrt.hpp>
0032 #include <boost/math/special_functions/expm1.hpp>
0033 #include <boost/math/special_functions/gamma.hpp>
0034 #endif
0035
0036 #include <limits>
0037 #include <string>
0038 #include <type_traits>
0039
0040 #if (defined(BOOST_CLANG) && defined(BOOST_CLANG_VERSION) && (BOOST_CLANG_VERSION <= 90000))
0041 #define BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE struct
0042 #else
0043 #define BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE class
0044 #endif
0045
0046 namespace boost { namespace multiprecision { namespace backends {
0047
0048 template <typename FloatingPointType>
0049 class cpp_double_fp_backend;
0050
0051 template <typename FloatingPointType>
0052 constexpr auto eval_add(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0053 template <typename FloatingPointType>
0054 constexpr auto eval_add(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void;
0055 template <typename FloatingPointType>
0056 constexpr auto eval_subtract(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0057 template <typename FloatingPointType>
0058 constexpr auto eval_subtract(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void;
0059 template <typename FloatingPointType>
0060 constexpr auto eval_multiply(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0061 template <typename FloatingPointType>
0062 constexpr auto eval_multiply(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void;
0063 template <typename FloatingPointType>
0064 constexpr auto eval_divide(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0065 template <typename FloatingPointType>
0066 constexpr auto eval_divide(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void;
0067 template <typename FloatingPointType>
0068 constexpr auto eval_eq(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool;
0069 template <typename FloatingPointType>
0070 constexpr auto eval_lt(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool;
0071 template <typename FloatingPointType>
0072 constexpr auto eval_gt(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool;
0073 template <typename FloatingPointType>
0074 constexpr auto eval_is_zero(const cpp_double_fp_backend<FloatingPointType>& x) -> bool;
0075 template <typename FloatingPointType>
0076 constexpr auto eval_signbit(const cpp_double_fp_backend<FloatingPointType>& x) -> int;
0077
0078 template <typename FloatingPointType>
0079 constexpr auto eval_fabs(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a) -> void;
0080 template <typename FloatingPointType>
0081 constexpr auto eval_frexp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, int* v) -> void;
0082 template <typename FloatingPointType>
0083 constexpr auto eval_ldexp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, int v) -> void;
0084 template <typename FloatingPointType>
0085 constexpr auto eval_floor(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0086 template <typename FloatingPointType>
0087 constexpr auto eval_ceil(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0088 template <typename FloatingPointType>
0089 constexpr auto eval_fpclassify(const cpp_double_fp_backend<FloatingPointType>& o) -> int;
0090
0091 template <typename FloatingPointType>
0092 constexpr auto eval_sqrt(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& o) -> void;
0093
0094 template <typename FloatingPointType>
0095 constexpr auto eval_pow(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x, const cpp_double_fp_backend<FloatingPointType>& a) -> void;
0096
0097 template <typename FloatingPointType,
0098           typename IntegralType>
0099 constexpr auto eval_pow(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x, IntegralType n) -> typename ::std::enable_if<boost::multiprecision::detail::is_integral<IntegralType>::value, void>::type;
0100
0101 template <typename FloatingPointType,
0102           typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) < 16))>::type const* = nullptr>
0103 constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0104
0105 template <typename FloatingPointType,
0106           typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && (((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) <= 36)))>::type const* = nullptr>
0107 constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0108
0109 template <typename FloatingPointType,
0110           typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) > 36))>::type const* = nullptr>
0111 constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0112
0113 template <typename FloatingPointType>
0114 constexpr auto eval_log(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
0115
0116 template <typename FloatingPointType>
0117 constexpr auto eval_convert_to(signed long long* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> void;
0118
0119 template <typename FloatingPointType>
0120 constexpr auto eval_convert_to(unsigned long long* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> void;
0121
0122 #ifdef BOOST_HAS_INT128
0123 template <typename FloatingPointType>
0124 constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type;
0125
0126 template <typename FloatingPointType>
0127 constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<!(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type;
0128
0129 template <typename FloatingPointType>
0130 constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type;
0131
0132 template <typename FloatingPointType>
0133 constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<!(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type;
0134 #endif
0135
0136 template <typename FloatingPointType,
0137           typename OtherFloatingPointType>
0138 constexpr auto eval_convert_to(OtherFloatingPointType* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename ::std::enable_if<cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value>::type;
0139
0140 template <typename FloatingPointType>
0141 constexpr auto hash_value(const cpp_double_fp_backend<FloatingPointType>& a) -> ::std::size_t;
0142
0143 template <typename FloatingPointType>
0144 constexpr auto fabs(const cpp_double_fp_backend<FloatingPointType>& a) -> cpp_double_fp_backend<FloatingPointType>;
0145
0146 } } } // namespace boost::multiprecision::backends
0147
0148 namespace std {
0149
0150 // Foward declarations of various specializations of std::numeric_limits
0151
0152 template <typename FloatingPointType,
0153           const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
0154 BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >;
0155
0156 } // namespace std
0157
0158 namespace boost { namespace multiprecision {
0159
0160 template <typename FloatingPointType>
0161 struct number_category<backends::cpp_double_fp_backend<FloatingPointType>> : public std::integral_constant<int, number_kind_floating_point> { };
0162
0163 namespace backends {
0164
0165 // A cpp_double_fp_backend is represented by an unevaluated sum of two
0166 // floating-point units, a0 and a1, which satisfy |a1| <= (1 / 2) * ulp(a0).
0167 // The type of the floating-point constituents should adhere to IEEE754.
0168 // Although the constituent parts (a0 and a1) satisfy |a1| <= (1 / 2) * ulp(a0),
0169 // the composite type does not adhere to these strict error bounds. Its error
0170 // bounds are larger.
0171
0172 // This class has been tested with floats having single-precision (4 byte),
0173 // double-precision (8 byte) and quad precision (16 byte, such as GCC's __float128).
0174
0175 template <typename FloatingPointType>
0176 class cpp_double_fp_backend
0177 {
0178  public:
0179    using float_type = FloatingPointType;
0180
0181    static_assert
0182    (
0183          cpp_df_qf_detail::is_floating_point<float_type>::value
0184       && bool
0185          {
0186                ((cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits ==  24) && cpp_df_qf_detail::ccmath::numeric_limits<float_type>::is_specialized && cpp_df_qf_detail::ccmath::numeric_limits<float_type>::is_iec559)
0187             || ((cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits ==  53) && cpp_df_qf_detail::ccmath::numeric_limits<float_type>::is_specialized && cpp_df_qf_detail::ccmath::numeric_limits<float_type>::is_iec559)
0188             || ((cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits ==  64) && cpp_df_qf_detail::ccmath::numeric_limits<float_type>::is_specialized && cpp_df_qf_detail::ccmath::numeric_limits<float_type>::is_iec559)
0189             ||  (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits == 113)
0190          }, "Error: float_type does not fulfill the backend requirements of cpp_double_fp_backend"
0191    );
0192
0193    using rep_type   = cpp_df_qf_detail::pair<float_type, float_type>;
0194    using arithmetic = cpp_df_qf_detail::exact_arithmetic<float_type>;
0195
0196    using signed_types   = std::tuple<signed char, signed short, signed int, signed long, signed long long, std::intmax_t>;
0197    using unsigned_types = std::tuple<unsigned char, unsigned short, unsigned int, unsigned long, unsigned long long, std::uintmax_t>;
0198    using float_types    = std::tuple<float, double, long double>;
0199    using exponent_type  = int;
0200
0201    static constexpr int my_digits         = static_cast<int>(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits * static_cast<int>(INT8_C(2)));
0202    static constexpr int my_digits10       = static_cast<int>(boost::multiprecision::detail::calc_digits10    <static_cast<unsigned>(my_digits)>::value);
0203    static constexpr int my_max_digits10   = static_cast<int>(boost::multiprecision::detail::calc_max_digits10<static_cast<unsigned>(my_digits)>::value);
0204    static constexpr int my_max_exponent   = cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max_exponent;
0205    static constexpr int my_max_exponent10 = cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max_exponent10;
0206
0207    static constexpr int my_min_exponent =
0208       static_cast<int>
0209       (
0210            cpp_df_qf_detail::ccmath::numeric_limits<float_type>::min_exponent
0211          + cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits
0212       );
0213
0214    static constexpr int my_min_exponent10 =
0215       static_cast<int>
0216       (
0217          -static_cast<int>
0218           (
0219              boost::multiprecision::detail::calc_digits10<static_cast<unsigned>(-my_min_exponent)>::value
0220           )
0221       );
0222
0223    // Default constructor.
0224    constexpr cpp_double_fp_backend() noexcept { }
0225
0226    // Copy constructor.
0227    constexpr cpp_double_fp_backend(const cpp_double_fp_backend& other) : data(other.data) { }
0228
0229    // Move constructor.
0230    constexpr cpp_double_fp_backend(cpp_double_fp_backend&& other) noexcept : data(static_cast<rep_type&&>(other.data)) { }
0231
0232    // Constructors from other floating-point types.
0233    template <typename OtherFloatType,
0234              typename ::std::enable_if<(    cpp_df_qf_detail::is_floating_point<OtherFloatType>::value
0235                                         && (cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatType>::digits <= cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
0236    constexpr cpp_double_fp_backend(const OtherFloatType& f)
0237       : data(f, static_cast<float_type>(0.0F)) { }
0238
0239    template <typename OtherFloatType,
0240              typename ::std::enable_if<(    cpp_df_qf_detail::is_floating_point<OtherFloatType>::value
0241                                         && (cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatType>::digits > cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
0242    constexpr cpp_double_fp_backend(const OtherFloatType& f)
0243       : data(static_cast<float_type>(f),
0244              static_cast<float_type>(f - static_cast<OtherFloatType>(static_cast<float_type>(f)))) { }
0245
0246    // Construtor from another kind of cpp_double_fp_backend<> object.
0247
0248    template <typename OtherFloatType,
0249              typename ::std::enable_if<(    cpp_df_qf_detail::is_floating_point<OtherFloatType>::value
0250                                         && (!std::is_same<FloatingPointType, OtherFloatType>::value))>::type const* = nullptr>
0251    constexpr cpp_double_fp_backend(const cpp_double_fp_backend<OtherFloatType>& a)
0252       : cpp_double_fp_backend(cpp_double_fp_backend(a.my_first()) += a.my_second()) { }
0253
0254    // Constructors from integers.
0255    template <typename SignedIntegralType,
0256              typename ::std::enable_if<(     boost::multiprecision::detail::is_integral<SignedIntegralType>::value
0257                                         && (!boost::multiprecision::detail::is_unsigned<SignedIntegralType>::value)
0258                                         && ((static_cast<int>(sizeof(SignedIntegralType) * 8) - 1) <= cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
0259    constexpr cpp_double_fp_backend(const SignedIntegralType& n)
0260       : data(static_cast<float_type>(n), static_cast<float_type>(0.0F)) { }
0261
0262    template <typename UnsignedIntegralType,
0263              typename ::std::enable_if<(    boost::multiprecision::detail::is_integral<UnsignedIntegralType>::value
0264                                         &&  boost::multiprecision::detail::is_unsigned<UnsignedIntegralType>::value
0265                                         && (static_cast<int>(sizeof(UnsignedIntegralType) * 8u) <= cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
0266    constexpr cpp_double_fp_backend(const UnsignedIntegralType& u)
0267       : data(static_cast<float_type>(u), static_cast<float_type>(0.0F)) { }
0268
0269    // Constructors from integers which hold more information than *this can contain.
0270    template <typename UnsignedIntegralType,
0271              typename ::std::enable_if<(    boost::multiprecision::detail::is_integral<UnsignedIntegralType>::value
0272                                         &&  boost::multiprecision::detail::is_unsigned<UnsignedIntegralType>::value
0273                                         && (static_cast<int>(sizeof(UnsignedIntegralType) * 8u) > cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
0274    constexpr cpp_double_fp_backend(UnsignedIntegralType u)
0275       : data(static_cast<float_type>(u & cpp_df_qf_detail::float_mask<UnsignedIntegralType, float_type>()),
0276              static_cast<float_type>(0.0F))
0277    {
0278       using local_unsigned_integral_type = UnsignedIntegralType;
0279
0280       if (u > cpp_df_qf_detail::float_mask<local_unsigned_integral_type, float_type>())
0281       {
0282          local_unsigned_integral_type
0283             local_flt_mask
0284             {
0285                cpp_df_qf_detail::float_mask<local_unsigned_integral_type, float_type>()
0286             };
0287
0288          for (int     index_mask_lsb =  cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits;
0289                      (index_mask_lsb <  static_cast<int>(sizeof(local_unsigned_integral_type) * 8u))
0290                   && (local_flt_mask != local_unsigned_integral_type { UINT8_C(0) });
0291                       index_mask_lsb += cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits)
0292          {
0293             local_flt_mask <<= static_cast<unsigned>(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits);
0294
0295             add_unchecked_limb(static_cast<float_type>(u & local_flt_mask));
0296          }
0297       }
0298    }
0299
0300    template <typename SignedIntegralType,
0301              typename ::std::enable_if<(     boost::multiprecision::detail::is_integral<SignedIntegralType>::value
0302                                         && (!boost::multiprecision::detail::is_unsigned<SignedIntegralType>::value)
0303                                         && ((static_cast<int>(sizeof(SignedIntegralType) * 8) - 1) > cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
0304    constexpr cpp_double_fp_backend(SignedIntegralType n)
0305    {
0306       const bool is_neg { (n < SignedIntegralType { INT8_C(0) }) };
0307
0308       using local_unsigned_integral_type = typename boost::multiprecision::detail::make_unsigned<SignedIntegralType>::type;
0309
0310       const local_unsigned_integral_type
0311          u_val
0312          {
0313             (!is_neg)
0314                ? static_cast<local_unsigned_integral_type>(n)
0315                : static_cast<local_unsigned_integral_type>
0316                  (
0317                       static_cast<local_unsigned_integral_type>(~n)
0318                     + static_cast<local_unsigned_integral_type>(UINT8_C(1))
0319                  )
0320          };
0321
0322       data = cpp_double_fp_backend(u_val).data;
0323
0324       if(is_neg) { negate(); }
0325    }
0326
0327    constexpr cpp_double_fp_backend(const float_type& a, const float_type& b) noexcept : data(a, b) { }
0328
0329    constexpr cpp_double_fp_backend(const cpp_df_qf_detail::pair<float_type, float_type>& p) noexcept : data(p) { }
0330
0331    cpp_double_fp_backend(const char* p_str)
0332    {
0333       *this = p_str;
0334    }
0335
0336    // Assignment operator.
0337    constexpr auto operator=(const cpp_double_fp_backend& other) -> cpp_double_fp_backend&
0338    {
0339       if (this != &other)
0340       {
0341          data = other.data;
0342       }
0343
0344       return *this;
0345    }
0346
0347    // Move assignment operator.
0348    constexpr auto operator=(cpp_double_fp_backend&& other) noexcept -> cpp_double_fp_backend&
0349    {
0350       data = static_cast<rep_type&&>(other.data);
0351
0352       return *this;
0353    }
0354
0355    // Assignment operator from another kind of cpp_double_fp_backend<> object.
0356    template <typename OtherFloatType,
0357              typename ::std::enable_if<(   cpp_df_qf_detail::is_floating_point<OtherFloatType>::value
0358                                         && (!std::is_same<FloatingPointType, OtherFloatType>::value))>::type const* = nullptr>
0359    constexpr auto operator=(const cpp_double_fp_backend<OtherFloatType>& other) -> cpp_double_fp_backend&
0360    {
0361      return operator=(cpp_double_fp_backend(other));
0362    }
0363
0364    template <typename OtherFloatType>
0365    constexpr auto operator=(const OtherFloatType f) -> typename ::std::enable_if<cpp_df_qf_detail::is_floating_point<OtherFloatType>::value, cpp_double_fp_backend&>::type
0366    {
0367      return operator=(cpp_double_fp_backend(f));
0368    }
0369
0370    template <typename IntegralType>
0371    constexpr auto operator=(const IntegralType n) -> typename ::std::enable_if<boost::multiprecision::detail::is_integral<IntegralType>::value, cpp_double_fp_backend&>::type
0372    {
0373      return operator=(cpp_double_fp_backend(n));
0374    }
0375
0376    auto operator=(const char* p_str) -> cpp_double_fp_backend&
0377    {
0378       rd_string(p_str);
0379
0380       return *this;
0381    }
0382
0383    constexpr auto hash() const -> ::std::size_t
0384    {
0385       #if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128)
0386       using local_float_type = typename std::conditional<::std::is_same<float_type, ::boost::float128_type>::value,
0387                                                          long double,
0388                                                          float_type>::type;
0389       #else
0390       using local_float_type = float_type;
0391       #endif
0392
0393       std::size_t result { UINT8_C(0) };
0394
0395       int n_first { };
0396       int n_second { };
0397
0398       boost::multiprecision::detail::hash_combine(result, static_cast<local_float_type>(cpp_df_qf_detail::ccmath::frexp(data.first, &n_first)));
0399       boost::multiprecision::detail::hash_combine(result, static_cast<local_float_type>(cpp_df_qf_detail::ccmath::frexp(data.second, &n_second)));
0400       boost::multiprecision::detail::hash_combine(result, n_first);
0401       boost::multiprecision::detail::hash_combine(result, n_second);
0402
0403       return result;
0404    }
0405
0406    // The public methods follow.
0407
0408    constexpr auto isneg_unchecked() const noexcept -> bool { return (data.first < 0); }
0409
0410    constexpr auto iszero_unchecked() const noexcept -> bool { return (data.first  == float_type { 0.0F }); }
0411
0412    constexpr auto is_one() const noexcept -> bool
0413    {
0414       return
0415          (
0416                (data.second == float_type { 0.0F })
0417             && (data.first == float_type { 1.0F })
0418          );
0419    }
0420
0421    constexpr auto negate() -> void
0422    {
0423       const bool isinf_u { (cpp_df_qf_detail::ccmath::isinf)(data.first) };
0424       const bool isnan_u { (cpp_df_qf_detail::ccmath::isnan)(data.first) };
0425
0426       if      (isnan_u) { }
0427       else if (isinf_u)
0428       {
0429          data.first = -data.first;
0430       }
0431       else
0432       {
0433          if (!iszero_unchecked())
0434          {
0435             data = arithmetic::normalize(-data.first, -data.second);
0436          }
0437       }
0438    }
0439
0440    // Getters/Setters
0441    constexpr auto my_first () const noexcept -> const float_type& { return data.first; }
0442    constexpr auto my_second() const noexcept -> const float_type& { return data.second; }
0443
0444    constexpr auto rep() noexcept -> rep_type& { return data; }
0445
0446    constexpr auto rep() const noexcept -> const rep_type& { return data; }
0447
0448    constexpr auto crep() const noexcept -> const rep_type& { return data; }
0449
0450    // Unary add/sub/mul/div follow in the upcoming paragraphs.
0451
0452    constexpr auto operator+=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend&
0453    {
0454       const int fpc_u { eval_fpclassify(*this) };
0455       const int fpc_v { eval_fpclassify(v) };
0456
0457       if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL))
0458       {
0459          // Handle special cases like zero, inf and NaN.
0460
0461          if (fpc_u == FP_NAN)
0462          {
0463             return *this;
0464          }
0465
0466          const bool isinf_v { (fpc_v == FP_INFINITE) };
0467
0468          if (fpc_u == FP_INFINITE)
0469          {
0470             if (isinf_v && (isneg_unchecked() != v.isneg_unchecked()))
0471             {
0472                *this = cpp_double_fp_backend::my_value_nan();
0473             }
0474
0475             return *this;
0476          }
0477
0478          const bool iszero_u { ((fpc_u == FP_ZERO) || (fpc_u == FP_SUBNORMAL)) };
0479          const bool isnan_v  { (fpc_v == FP_NAN) };
0480
0481          if (iszero_u || (isnan_v || isinf_v))
0482          {
0483             if (iszero_u)
0484             {
0485                data.first  = float_type { 0.0F };
0486                data.second = float_type { 0.0F };
0487             }
0488
0489             const bool iszero_v { ((fpc_v == FP_ZERO) || (fpc_v == FP_SUBNORMAL)) };
0490
0491             return ((!iszero_v) ? operator=(v) : *this);
0492          }
0493       }
0494
0495       add_unchecked(v);
0496
0497       return *this;
0498    }
0499
0500    constexpr auto operator-=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend&
0501    {
0502       const int fpc_u { eval_fpclassify(*this) };
0503       const int fpc_v { eval_fpclassify(v) };
0504
0505       if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL))
0506       {
0507          // Handle special cases like zero, inf and NaN.
0508
0509          if (fpc_u == FP_NAN)
0510          {
0511             return *this;
0512          }
0513
0514          const bool isinf_v { (fpc_v == FP_INFINITE) };
0515
0516          if (fpc_u == FP_INFINITE)
0517          {
0518             if (isinf_v && (isneg_unchecked() == v.isneg_unchecked()))
0519             {
0520                *this = cpp_double_fp_backend::my_value_nan();
0521             }
0522
0523             return *this;
0524          }
0525
0526          const bool iszero_u { ((fpc_u == FP_ZERO) || (fpc_u == FP_SUBNORMAL)) };
0527          const bool isnan_v  { (fpc_v == FP_NAN) };
0528
0529          if (iszero_u || (isnan_v || isinf_v))
0530          {
0531             if (iszero_u)
0532             {
0533                data.first  = float_type { 0.0F };
0534                data.second = float_type { 0.0F };
0535             }
0536
0537             const bool iszero_v { ((fpc_v == FP_ZERO) || (fpc_v == FP_SUBNORMAL)) };
0538
0539             return ((!iszero_v) ? operator=(-v) : *this);
0540          }
0541       }
0542
0543       if (this == &v)
0544       {
0545          data.first  = float_type { 0.0F };
0546          data.second = float_type { 0.0F };
0547
0548          return *this;
0549       }
0550
0551       const rep_type thi_tlo { arithmetic::two_diff(data.second, v.data.second) };
0552
0553       data = arithmetic::two_diff(data.first, v.data.first);
0554
0555       if (cpp_df_qf_detail::ccmath::isinf(data.first))
0556       {
0557          // Handle overflow.
0558          const bool b_neg { (data.first < float_type { 0.0F }) };
0559
0560          *this = cpp_double_fp_backend::my_value_inf();
0561
0562          if (b_neg)
0563          {
0564             negate();
0565          }
0566
0567          return *this;
0568       }
0569
0570       data = arithmetic::two_hilo_sum(data.first, data.second + thi_tlo.first);
0571
0572       data = arithmetic::two_hilo_sum(data.first, thi_tlo.second + data.second);
0573
0574       return *this;
0575    }
0576
0577    constexpr auto operator*=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend&
0578    {
0579       // Evaluate the sign of the result.
0580
0581       const int fpc_u { eval_fpclassify(*this) };
0582       const int fpc_v { eval_fpclassify(v) };
0583
0584       if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL))
0585       {
0586          // Handle special cases like zero, inf and NaN.
0587          const bool isinf_u  { (fpc_u == FP_INFINITE) };
0588          const bool isinf_v  { (fpc_v == FP_INFINITE) };
0589          const bool iszero_u { (fpc_u == FP_ZERO) };
0590          const bool iszero_v { (fpc_v == FP_ZERO) };
0591
0592          if (((fpc_u == FP_NAN) || (fpc_v == FP_NAN)) || (isinf_u && iszero_v) || (isinf_v && iszero_u))
0593          {
0594             return operator=( cpp_double_fp_backend::my_value_nan());
0595          }
0596
0597          if (isinf_u || isinf_v)
0598          {
0599             const bool b_neg { (isneg_unchecked() != v.isneg_unchecked()) };
0600
0601             *this = cpp_double_fp_backend::my_value_inf();
0602
0603             if (b_neg)
0604             {
0605                negate();
0606             }
0607
0608             return *this;
0609          }
0610
0611          if (iszero_u || iszero_v)
0612          {
0613             return operator=(cpp_double_fp_backend(0));
0614          }
0615       }
0616
0617       mul_unchecked(v);
0618
0619       return *this;
0620    }
0621
0622    constexpr auto operator/=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend&
0623    {
0624       const int fpc_u { eval_fpclassify(*this) };
0625       const int fpc_v { eval_fpclassify(v) };
0626
0627       if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL))
0628       {
0629          // Handle special cases like zero, inf and NaN.
0630          const bool isnan_u { (fpc_u == FP_NAN) };
0631          const bool isnan_v { (fpc_v == FP_NAN) };
0632
0633          if (isnan_u || isnan_v)
0634          {
0635             return operator=(cpp_double_fp_backend::my_value_nan());
0636          }
0637
0638          const bool iszero_u { (fpc_u == FP_ZERO) };
0639          const bool iszero_v { (fpc_v == FP_ZERO) };
0640
0641          if (iszero_u)
0642          {
0643             if (iszero_v)
0644             {
0645                return operator=(cpp_double_fp_backend::my_value_nan());
0646             }
0647             else
0648             {
0649                return operator=(cpp_double_fp_backend(0));
0650             }
0651          }
0652
0653          // Handle more special cases like zero, inf and NaN.
0654          if (iszero_v)
0655          {
0656             const bool b_neg = isneg_unchecked();
0657
0658             *this = cpp_double_fp_backend::my_value_inf();
0659
0660             if (b_neg)
0661             {
0662                negate();
0663             }
0664
0665             return *this;
0666          }
0667
0668          const bool isinf_v { (fpc_v == FP_INFINITE) };
0669          const bool isinf_u { (fpc_u == FP_INFINITE) };
0670
0671          if (isinf_u)
0672          {
0673             if (isinf_v)
0674             {
0675                return operator=(cpp_double_fp_backend::my_value_nan());
0676             }
0677             else
0678             {
0679                const bool b_neg { isneg_unchecked() };
0680
0681                return operator=((!b_neg) ? cpp_double_fp_backend::my_value_inf() : -cpp_double_fp_backend::my_value_inf());
0682             }
0683          }
0684
0685          if (isinf_v)
0686          {
0687             return operator=(cpp_double_fp_backend(0));
0688          }
0689       }
0690
0691       if (this == &v)
0692       {
0693          data.first  = float_type { 1.0F };
0694          data.second = float_type { 0.0F };
0695
0696          return *this;
0697       }
0698
0699       // The division algorithm has been taken from Victor Shoup,
0700       // package WinNTL-5_3_2. It might originally be related to the
0701       // K. Briggs work. The algorithm has been significantly simplified
0702       // while still attempting to retain proper rounding corrections.
0703       // Checks for overflow and underflow have been added.
0704
0705       const float_type C { data.first / v.data.first };
0706
0707       float_type c { cpp_df_qf_detail::split_maker<float_type>::value * C };
0708
0709       float_type hc { };
0710
0711       if (cpp_df_qf_detail::ccmath::isinf(c))
0712       {
0713          // Handle overflow by scaling down (and then back up) with the split.
0714
0715          hc =
0716             cpp_df_qf_detail::ccmath::ldexp
0717             (
0718                C - float_type { C - cpp_df_qf_detail::ccmath::ldexp(C, -cpp_df_qf_detail::split_maker<float_type>::n_shl) },
0719                cpp_df_qf_detail::split_maker<float_type>::n_shl
0720             );
0721       }
0722       else
0723       {
0724          hc = c - float_type { c - C };
0725       }
0726
0727       float_type u { cpp_df_qf_detail::split_maker<float_type>::value * v.data.first };
0728
0729       const float_type hv =
0730       (
0731          cpp_df_qf_detail::ccmath::isinf(u)
0732             ? cpp_df_qf_detail::ccmath::ldexp
0733               (
0734                  // Handle overflow by scaling down (and then back up) with the split.
0735                  v.data.first - float_type { v.data.first - cpp_df_qf_detail::ccmath::ldexp(v.data.first, -cpp_df_qf_detail::split_maker<float_type>::n_shl) },
0736                  cpp_df_qf_detail::split_maker<float_type>::n_shl
0737               )
0738             : u - float_type { u - v.data.first }
0739       );
0740
0741       const float_type U { C * v.data.first };
0742
0743       u = cpp_df_qf_detail::ccmath::unsafe::fma(hc, hv, -U);
0744
0745       {
0746          const float_type tv { v.data.first - hv };
0747
0748          u = cpp_df_qf_detail::ccmath::unsafe::fma(hc, tv, u);
0749
0750          const float_type tc { C - hc };
0751
0752          u = cpp_df_qf_detail::ccmath::unsafe::fma(tc, hv, u);
0753          u = cpp_df_qf_detail::ccmath::unsafe::fma(tc, tv, u);
0754       }
0755
0756       c = float_type { (data.first - U) - u } + data.second;
0757
0758       c = (c - float_type { C * v.data.second }) / v.data.first;
0759
0760       // Perform even more simplifications compared to Victor Shoup.
0761       data.first  = C + c;
0762       data.second = float_type { C - data.first } + c;
0763
0764       return *this;
0765    }
0766
0767    // Unary minus operator.
0768    constexpr auto operator-() const -> cpp_double_fp_backend
0769    {
0770       cpp_double_fp_backend v { *this };
0771
0772       v.negate();
0773
0774       return v;
0775    }
0776
0777    // Public helper functions.
0778    static constexpr auto pown(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, int p) -> void
0779    {
0780       using local_float_type = cpp_double_fp_backend;
0781
0782       if (p == 2)
0783       {
0784          result = x; result.mul_unchecked(x);
0785       }
0786       else if (p == 3)
0787       {
0788          result = x; result.mul_unchecked(x); result.mul_unchecked(x);
0789       }
0790       else if (p == 4)
0791       {
0792          local_float_type x2 { x };
0793          x2.mul_unchecked(x);
0794
0795          result = x2;
0796          result.mul_unchecked(x2);
0797       }
0798       else if (p == 1)
0799       {
0800          result = x;
0801       }
0802       else if (p < 0)
0803       {
0804          // The case p == 0 is checked in a higher calling layer.
0805
0806          pown(result, local_float_type(1U) / x, -p);
0807       }
0808       else
0809       {
0810          result = local_float_type(1U);
0811
0812          local_float_type y(x);
0813
0814          auto p_local = static_cast<std::uint32_t>(p);
0815
0816          for (;;)
0817          {
0818             if (static_cast<std::uint_fast8_t>(p_local & static_cast<std::uint32_t>(UINT8_C(1))) != static_cast<std::uint_fast8_t>(UINT8_C(0)))
0819             {
0820                result.mul_unchecked(y);
0821             }
0822
0823             p_local >>= 1U;
0824
0825             if (p_local == static_cast<std::uint32_t>(UINT8_C(0)))
0826             {
0827                break;
0828             }
0829             else
0830             {
0831                y.mul_unchecked(y);
0832             }
0833          }
0834       }
0835    }
0836
0837    constexpr auto swap(cpp_double_fp_backend& other) -> void
0838    {
0839       if (this != &other)
0840       {
0841          const rep_type tmp { data };
0842
0843          data = other.data;
0844
0845          other.data = tmp;
0846       }
0847    }
0848
0849    constexpr auto swap(cpp_double_fp_backend&& other) noexcept -> void
0850    {
0851       const rep_type tmp { static_cast<typename cpp_double_fp_backend::rep_type&&>(data) };
0852
0853       data = other.data;
0854
0855       other.data = tmp;
0856    }
0857
0858    constexpr auto compare(const cpp_double_fp_backend& other) const noexcept -> int
0859    {
0860       // Return 1 for *this > other, -1 for *this < other, 0 for *this = other.
0861
0862       if ((cpp_df_qf_detail::ccmath::isnan)(data.first))
0863       {
0864         return -1;
0865       }
0866       else
0867       {
0868         return (my_first() > other.my_first()) ?  1 : (my_first()  < other.my_first())
0869                                                ? -1 : (my_second() > other.my_second())
0870                                                ?  1 : (my_second() < other.my_second())
0871                                                ? -1 : 0;
0872       }
0873    }
0874
0875    // TBD: Exactly what compilers/language-standards are needed to make this constexpr?
0876    // TBD: It odes not really become constexpr until we stop using an intermediate
0877    // cpp_bin_float anyway. But I will leave this comment for future library evolution.
0878
0879    auto str(std::streamsize number_of_digits, const std::ios::fmtflags format_flags) const -> std::string
0880    {
0881       // Use cpp_bin_float when writing to string. This is similar
0882       // to the use of cpp_bin_float when reading from string.
0883
0884       cpp_bin_float_read_write_backend_type f_bin { data.first };
0885
0886       eval_add(f_bin, cpp_bin_float_read_write_backend_type(data.second));
0887
0888       return f_bin.str(number_of_digits, format_flags);
0889    }
0890
0891    constexpr auto order02() const -> int
0892    {
0893       // TBD: Is there another option to get the base-2 log
0894       // that's more unequivocally closer to constexpr?
0895
0896       // TBD: Either way, integrate this (or something like it)
0897       // into any potential implementation of eval_ilogb().
0898
0899       int e2 { };
0900
0901       cpp_double_fp_backend dummy { };
0902
0903       eval_frexp(dummy, *this, &e2);
0904
0905       return e2;
0906    }
0907
0908    static constexpr auto my_value_max() noexcept -> cpp_double_fp_backend
0909    {
0910       // Use the non-normalized sum of two maximum values, where the lower
0911       // value is "shifted" right in the sense of floating-point ldexp.
0912
0913       return
0914          cpp_double_fp_backend
0915          (
0916             arithmetic::two_hilo_sum
0917             (
0918                float_type
0919                (
0920                     (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)()
0921                   * (
0922                          static_cast<float_type>(1.0F)
0923                        - static_cast<float_type>(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::epsilon())
0924                     )
0925                ),
0926                float_type
0927                (
0928                   cpp_df_qf_detail::ccmath::ldexp
0929                   (
0930                      (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)(),
0931                      -cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits
0932                   )
0933                )
0934             )
0935          );
0936    }
0937
0938    static constexpr auto my_value_min() noexcept -> cpp_double_fp_backend
0939    {
0940       // Use the non-normalized minimum value, where the lower value
0941       // is "shifted" left in the sense of floating-point ldexp.
0942
0943       return
0944          cpp_double_fp_backend
0945          (
0946             float_type
0947             (
0948                cpp_df_qf_detail::ccmath::ldexp
0949                (
0950                   (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::min)(),
0951                    cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits
0952                )
0953             )
0954          );
0955    }
0956
0957    static constexpr auto my_value_eps() noexcept -> cpp_double_fp_backend
0958    {
0959       return
0960          cpp_double_fp_backend
0961          (
0962             float_type(cpp_df_qf_detail::ccmath::ldexp(float_type { 1 }, int { 3 - my_digits }))
0963          );
0964    }
0965
0966    static constexpr auto my_value_nan() noexcept -> cpp_double_fp_backend
0967    {
0968       return cpp_double_fp_backend(static_cast<float_type>(NAN), static_cast<float_type>(0.0F));
0969    }
0970
0971    static constexpr auto my_value_inf() noexcept -> cpp_double_fp_backend
0972    {
0973       return cpp_double_fp_backend(static_cast<float_type>(HUGE_VAL), static_cast<float_type>(0.0F)); // conversion from double infinity OK
0974    }
0975
0976    static constexpr auto my_value_logmax() -> cpp_double_fp_backend
0977    {
0978       return
0979          cpp_double_fp_backend
0980          (
0981             cpp_df_qf_detail::ccmath::log
0982             (
0983                float_type
0984                (
0985                     (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)()
0986                   * (
0987                          static_cast<float_type>(1.0F)
0988                        - static_cast<float_type>(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::epsilon())
0989                     )
0990                )
0991             )
0992          );
0993    }
0994
0995    static constexpr auto my_value_logmin() -> cpp_double_fp_backend
0996    {
0997       return
0998          cpp_double_fp_backend
0999          (
1000             cpp_df_qf_detail::ccmath::log
1001             (
1002                float_type
1003                (
1004                   cpp_df_qf_detail::ccmath::ldexp
1005                   (
1006                      (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::min)(),
1007                       cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits
1008                   )
1009                )
1010             )
1011          );
1012    }
1013
1014    constexpr auto add_unchecked_limb(const float_type v_first) -> void
1015    {
1016       const float_type thi { data.second };
1017
1018       data = arithmetic::two_sum(data.first, v_first);
1019
1020       data = arithmetic::two_hilo_sum(data.first, data.second + thi);
1021
1022       data = arithmetic::two_hilo_sum(data.first, data.second);
1023    }
1024
1025 private:
1026    rep_type data;
1027
1028    using cpp_bin_float_read_write_backend_type = boost::multiprecision::backends::cpp_bin_float<static_cast<unsigned>(my_digits), digit_base_2, void, int, cpp_df_qf_detail::ccmath::numeric_limits<float_type>::min_exponent, cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max_exponent>;
1029
1030    constexpr auto rd_string(const char* p_str) -> bool;
1031
1032    constexpr auto add_unchecked(const cpp_double_fp_backend& v) -> void
1033    {
1034       const rep_type thi_tlo { arithmetic::two_sum(data.second, v.data.second) };
1035
1036       data = arithmetic::two_sum(data.first, v.data.first);
1037
1038       if (cpp_df_qf_detail::ccmath::isinf(data.first))
1039       {
1040          // Handle overflow.
1041          const bool b_neg { (data.first < float_type { 0.0F }) };
1042
1043          *this = cpp_double_fp_backend::my_value_inf();
1044
1045          if (b_neg)
1046          {
1047             negate();
1048          }
1049
1050          return;
1051       }
1052
1053       data = arithmetic::two_hilo_sum(data.first, data.second + thi_tlo.first);
1054
1055       data = arithmetic::two_hilo_sum(data.first, thi_tlo.second + data.second);
1056    }
1057
1058    constexpr auto mul_unchecked(const cpp_double_fp_backend& v) -> void
1059    {
1060       // The multiplication algorithm has been taken from Victor Shoup,
1061       // package WinNTL-5_3_2. It might originally be related to the
1062       // K. Briggs work. The algorithm has been significantly simplified
1063       // while still attempting to retain proper rounding corrections.
1064       // Checks for overflow and underflow have been added.
1065
1066       float_type C { cpp_df_qf_detail::split_maker<float_type>::value * data.first };
1067
1068       float_type hu { };
1069
1070       if (cpp_df_qf_detail::ccmath::isinf(C))
1071       {
1072          // Handle overflow by scaling down (and then back up) with the split.
1073
1074          C = data.first - cpp_df_qf_detail::ccmath::ldexp(data.first, -cpp_df_qf_detail::split_maker<float_type>::n_shl);
1075
1076          hu = cpp_df_qf_detail::ccmath::ldexp(data.first - C, cpp_df_qf_detail::split_maker<float_type>::n_shl);
1077       }
1078       else
1079       {
1080          hu = C - float_type { C - data.first };
1081       }
1082
1083       C = data.first * v.data.first;
1084
1085       if (cpp_df_qf_detail::ccmath::isinf(C))
1086       {
1087          // Handle overflow.
1088          const bool b_neg { (isneg_unchecked() != v.isneg_unchecked()) };
1089
1090          *this = cpp_double_fp_backend::my_value_inf();
1091
1092          if (b_neg)
1093          {
1094             negate();
1095          }
1096
1097          return;
1098       }
1099
1100       float_type c { cpp_df_qf_detail::split_maker<float_type>::value * v.data.first };
1101
1102       float_type hv { };
1103
1104       if (cpp_df_qf_detail::ccmath::isinf(c))
1105       {
1106          // Handle overflow by scaling down (and then back up) with the split.
1107
1108          c = v.data.first - cpp_df_qf_detail::ccmath::ldexp(v.data.first, -cpp_df_qf_detail::split_maker<float_type>::n_shl);
1109
1110          hv = cpp_df_qf_detail::ccmath::ldexp(v.data.first - c, cpp_df_qf_detail::split_maker<float_type>::n_shl);
1111       }
1112       else
1113       {
1114          hv = c - float_type { c - v.data.first };
1115       }
1116
1117       {
1118          const float_type tv { v.data.first - hv };
1119
1120          const float_type
1121             t1
1122             {
1123                cpp_df_qf_detail::ccmath::unsafe::fma
1124                (
1125                   hu,
1126                   tv,
1127                   cpp_df_qf_detail::ccmath::unsafe::fma(hu, hv, -C)
1128                )
1129             };
1130
1131          const float_type tu { data.first - hu };
1132
1133          c =    cpp_df_qf_detail::ccmath::unsafe::fma(tu, tv, cpp_df_qf_detail::ccmath::unsafe::fma(tu, hv, t1))
1134              + (data.first * v.data.second)
1135              + (data.second * v.data.first);
1136       }
1137
1138       // Perform even more simplifications compared to Victor Shoup.
1139       data.first  = C + c;
1140       data.second = float_type { C - data.first } + c;
1141    }
1142
1143    template <typename OtherFloatingPointType,
1144              typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits10 * 2) < 16))>::type const*>
1145    friend constexpr auto eval_exp(cpp_double_fp_backend<OtherFloatingPointType>& result, const cpp_double_fp_backend<OtherFloatingPointType>& x) -> void;
1146
1147    template <typename OtherFloatingPointType,
1148              typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value && (((cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits10 * 2) <= 36)))>::type const*>
1149    friend constexpr auto eval_exp(cpp_double_fp_backend<OtherFloatingPointType>& result, const cpp_double_fp_backend<OtherFloatingPointType>& x) -> void;
1150
1151    template <typename OtherFloatingPointType,
1152              typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits10 * 2) > 36))>::type const*>
1153    friend constexpr auto eval_exp(cpp_double_fp_backend<OtherFloatingPointType>& result, const cpp_double_fp_backend<OtherFloatingPointType>& x) -> void;
1154 };
1155
1156 template <typename FloatingPointType>
1157 constexpr auto cpp_double_fp_backend<FloatingPointType>::rd_string(const char* p_str) -> bool
1158 {
1159    // Use an intermediate cpp_bin_float backend type for reading string input.
1160
1161    cpp_bin_float_read_write_backend_type f_bin { };
1162
1163    f_bin = p_str;
1164
1165    const int fpc { eval_fpclassify(f_bin) };
1166
1167    const bool is_definitely_nan { (fpc == FP_NAN) };
1168
1169    using local_double_fp_type = cpp_double_fp_backend<FloatingPointType>;
1170
1171    if (is_definitely_nan)
1172    {
1173       static_cast<void>(operator=(local_double_fp_type::my_value_nan()));
1174
1175       return true;
1176    }
1177
1178    const bool b_neg { (eval_signbit(f_bin) == 1) };
1179
1180    if (b_neg) { f_bin.negate(); }
1181
1182    const int
1183       expval_from_f_bin
1184       {
1185          [&f_bin]()
1186          {
1187             int expval { };
1188
1189             cpp_bin_float_read_write_backend_type dummy { };
1190
1191             eval_frexp(dummy, f_bin, &expval);
1192
1193             return expval;
1194          }()
1195       };
1196
1197    const auto is_zero_or_subnormal =
1198    (
1199          (fpc == FP_ZERO)
1200       || (expval_from_f_bin < static_cast<typename cpp_bin_float_read_write_backend_type::exponent_type>(local_double_fp_type::my_min_exponent))
1201    );
1202
1203    if (is_zero_or_subnormal)
1204    {
1205       data.first  = float_type { 0.0F };
1206       data.second = float_type { 0.0F };
1207
1208       return true;
1209    }
1210
1211    float_type flt_inf_check_first { };
1212
1213    eval_convert_to(&flt_inf_check_first, f_bin);
1214
1215    bool is_definitely_inf { ((fpc == FP_INFINITE) || cpp_df_qf_detail::ccmath::isinf(flt_inf_check_first)) };
1216
1217    if ((!is_definitely_inf) && (flt_inf_check_first > my_value_max().my_first()))
1218    {
1219       cpp_bin_float_read_write_backend_type f_bin_inf_check(f_bin);
1220
1221       eval_subtract(f_bin_inf_check, cpp_bin_float_read_write_backend_type(flt_inf_check_first));
1222
1223       float_type flt_inf_check_second { };
1224
1225       eval_convert_to(&flt_inf_check_second, f_bin_inf_check);
1226
1227       is_definitely_inf = eval_gt(local_double_fp_type(flt_inf_check_first, flt_inf_check_second), my_value_max());
1228    };
1229
1230    if (is_definitely_inf)
1231    {
1232       static_cast<void>(operator=(local_double_fp_type::my_value_inf()));
1233
1234       if (b_neg)
1235       {
1236         negate();
1237       }
1238
1239       return true;
1240    }
1241
1242    // The input string is normal. We will now extract its value.
1243
1244    data.first  = float_type { 0.0F };
1245    data.second = float_type { 0.0F };
1246
1247    // Handle small input values. Scale (and re-scale them below) if possible.
1248
1249    constexpr int pow2_scaling_for_small_input { cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits };
1250
1251    const bool
1252       has_pow2_scaling_for_small_input
1253       {
1254          (expval_from_f_bin < static_cast<int>(local_double_fp_type::my_min_exponent + pow2_scaling_for_small_input))
1255       };
1256
1257    if (has_pow2_scaling_for_small_input)
1258    {
1259       eval_ldexp(f_bin, f_bin, pow2_scaling_for_small_input);
1260    }
1261
1262    using local_builtin_float_type = typename std::conditional<(sizeof(float_type) <= sizeof(double)), double, float_type>::type;
1263
1264    constexpr unsigned
1265       digit_limit
1266       {
1267          static_cast<unsigned>
1268          (
1269             static_cast<int>
1270             (
1271                    (local_double_fp_type::my_digits / cpp_df_qf_detail::ccmath::numeric_limits<local_builtin_float_type>::digits)
1272                + (((local_double_fp_type::my_digits % cpp_df_qf_detail::ccmath::numeric_limits<local_builtin_float_type>::digits) != 0) ? 1 : 0)
1273             )
1274             * cpp_df_qf_detail::ccmath::numeric_limits<local_builtin_float_type>::digits
1275          )
1276       };
1277
1278    for(auto i = static_cast<unsigned>(UINT8_C(0));
1279             i < digit_limit;
1280             i = static_cast<unsigned>(i + static_cast<unsigned>(cpp_df_qf_detail::ccmath::numeric_limits<local_builtin_float_type>::digits)))
1281    {
1282       local_builtin_float_type flt_part { };
1283
1284       eval_convert_to(&flt_part, f_bin);
1285
1286       eval_subtract(f_bin, cpp_bin_float_read_write_backend_type(flt_part));
1287
1288       eval_add(*this, local_double_fp_type { flt_part });
1289    }
1290
1291    if (has_pow2_scaling_for_small_input)
1292    {
1293       eval_ldexp(*this, *this, -pow2_scaling_for_small_input);
1294    }
1295
1296    if (b_neg) { negate(); }
1297
1298    return true;
1299 }
1300
1301 } } } // namespace boost::multiprecision::backends
1302
1303 #include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_constants.hpp>
1304
1305 namespace boost { namespace multiprecision { namespace backends {
1306
1307 template <typename FloatingPointType>
1308 constexpr int cpp_double_fp_backend<FloatingPointType>::my_digits;
1309 template <typename FloatingPointType>
1310 constexpr int cpp_double_fp_backend<FloatingPointType>::my_digits10;
1311 template <typename FloatingPointType>
1312 constexpr int cpp_double_fp_backend<FloatingPointType>::my_max_digits10;
1313 template <typename FloatingPointType>
1314 constexpr int cpp_double_fp_backend<FloatingPointType>::my_max_exponent;
1315 template <typename FloatingPointType>
1316 constexpr int cpp_double_fp_backend<FloatingPointType>::my_min_exponent;
1317 template <typename FloatingPointType>
1318 constexpr int cpp_double_fp_backend<FloatingPointType>::my_max_exponent10;
1319 template <typename FloatingPointType>
1320 constexpr int cpp_double_fp_backend<FloatingPointType>::my_min_exponent10;
1321
1322 template <typename FloatingPointType>
1323 constexpr cpp_double_fp_backend<FloatingPointType> operator+(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) { return cpp_double_fp_backend<FloatingPointType>(a) += b; }
1324 template <typename FloatingPointType>
1325 constexpr cpp_double_fp_backend<FloatingPointType> operator-(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) { return cpp_double_fp_backend<FloatingPointType>(a) -= b; }
1326 template <typename FloatingPointType>
1327 constexpr cpp_double_fp_backend<FloatingPointType> operator*(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) { return cpp_double_fp_backend<FloatingPointType>(a) *= b; }
1328 template <typename FloatingPointType>
1329 constexpr cpp_double_fp_backend<FloatingPointType> operator/(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) { return cpp_double_fp_backend<FloatingPointType>(a) /= b; }
1330
1331 template <typename FloatingPointType>
1332 constexpr auto eval_add(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void { result += x; }
1333 template <typename FloatingPointType>
1334 constexpr auto eval_add(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void { result = cpp_double_fp_backend<FloatingPointType>(a) += b; }
1335 template <typename FloatingPointType>
1336 constexpr auto eval_subtract(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void { result -= x; }
1337 template <typename FloatingPointType>
1338 constexpr auto eval_subtract(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void { result = cpp_double_fp_backend<FloatingPointType>(a) -= b; }
1339 template <typename FloatingPointType>
1340 constexpr auto eval_multiply(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void { result *= x; }
1341 template <typename FloatingPointType>
1342 constexpr auto eval_multiply(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void { result = cpp_double_fp_backend<FloatingPointType>(a) *= b; }
1343 template <typename FloatingPointType>
1344 constexpr auto eval_divide(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void { result /= x; }
1345 template <typename FloatingPointType>
1346 constexpr auto eval_divide(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void { result = cpp_double_fp_backend<FloatingPointType>(a) /= b; }
1347 template <typename FloatingPointType>
1348 constexpr auto eval_eq(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool { return (a.compare(b) == 0); }
1349 template <typename FloatingPointType>
1350 constexpr auto eval_lt(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool { return (a.compare(b) == -1); }
1351 template <typename FloatingPointType>
1352 constexpr auto eval_gt(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool { return (a.compare(b) == 1); }
1353 template <typename FloatingPointType>
1354
1355 constexpr auto eval_is_zero(const cpp_double_fp_backend<FloatingPointType>& x) -> bool
1356 {
1357    return x.iszero_unchecked();
1358 }
1359
1360 template <typename FloatingPointType>
1361 constexpr auto eval_signbit(const cpp_double_fp_backend<FloatingPointType>& x) -> int
1362 {
1363    return (x.isneg_unchecked() ? 1 : 0);
1364 }
1365
1366 template <typename FloatingPointType>
1367 constexpr auto eval_fabs(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a) -> void
1368 {
1369    result = a;
1370
1371    if (a.isneg_unchecked())
1372    {
1373       result.negate();
1374    }
1375 }
1376
1377 template <typename FloatingPointType>
1378 constexpr auto eval_frexp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, int* v) -> void
1379 {
1380    using local_backend_type = cpp_double_fp_backend<FloatingPointType>;
1381    using local_float_type = typename local_backend_type::float_type;
1382
1383    int expptr { };
1384
1385    const local_float_type fhi { cpp_df_qf_detail::ccmath::frexp(a.rep().first, &expptr) };
1386    const local_float_type flo { cpp_df_qf_detail::ccmath::ldexp(a.rep().second, -expptr) };
1387
1388    if (v != nullptr)
1389    {
1390       *v = expptr;
1391    }
1392
1393    result.rep() = local_backend_type::arithmetic::normalize(fhi, flo);
1394 }
1395
1396 template <typename FloatingPointType>
1397 constexpr auto eval_ldexp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, int v) -> void
1398 {
1399    using local_backend_type = cpp_double_fp_backend<FloatingPointType>;
1400    using local_float_type = typename local_backend_type::float_type;
1401
1402    const local_float_type fhi { cpp_df_qf_detail::ccmath::ldexp(a.crep().first,  v) };
1403    const local_float_type flo { cpp_df_qf_detail::ccmath::ldexp(a.crep().second, v) };
1404
1405    result.rep() = local_backend_type::arithmetic::normalize(fhi, flo);
1406 }
1407
1408 template <typename FloatingPointType>
1409 constexpr auto eval_floor(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
1410 {
1411    using local_backend_type = cpp_double_fp_backend<FloatingPointType>;
1412    using local_float_type = typename local_backend_type::float_type;
1413
1414    const local_float_type fhi { cpp_df_qf_detail::ccmath::floor(x.my_first()) };
1415
1416    if (fhi != x.my_first())
1417    {
1418       result.rep().first  = fhi;
1419       result.rep().second = local_float_type { 0 };
1420    }
1421    else
1422    {
1423       const local_float_type flo = { cpp_df_qf_detail::ccmath::floor(x.my_second()) };
1424
1425       result.rep() = local_backend_type::arithmetic::normalize(fhi, flo);
1426    }
1427 }
1428
1429 template <typename FloatingPointType>
1430 constexpr auto eval_ceil(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
1431 {
1432    // Compute -floor(-x);
1433    eval_floor(result, -x);
1434
1435    result.negate();
1436 }
1437
1438 template <typename FloatingPointType>
1439 constexpr auto eval_fpclassify(const cpp_double_fp_backend<FloatingPointType>& o) -> int
1440 {
1441    if      (cpp_df_qf_detail::ccmath::isnan(o.crep().first)) { return FP_NAN; }
1442    else if (cpp_df_qf_detail::ccmath::isinf(o.crep().first)) { return FP_INFINITE; }
1443    else if (eval_is_zero(o)) { return FP_ZERO; }
1444    else
1445    {
1446       using local_backend_type = cpp_double_fp_backend<FloatingPointType>;
1447       using local_float_type = typename local_backend_type::float_type;
1448
1449       const local_float_type fabs_x { cpp_df_qf_detail::ccmath::fabs(o.crep().first) };
1450
1451       if ((fabs_x > 0) && (fabs_x < (cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::min)()))
1452       {
1453          return FP_SUBNORMAL;
1454       }
1455       else
1456       {
1457          return FP_NORMAL;
1458       }
1459    }
1460 }
1461
1462 template <typename FloatingPointType>
1463 constexpr auto eval_sqrt(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& o) -> void
1464 {
1465    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
1466    using local_float_type = typename double_float_type::float_type;
1467
1468    const int fpc { eval_fpclassify(o) };
1469
1470    const bool isneg_o { o.isneg_unchecked() };
1471
1472    if ((fpc != FP_NORMAL) || isneg_o)
1473    {
1474       if (fpc == FP_ZERO)
1475       {
1476          result = double_float_type(0);
1477
1478          return;
1479       }
1480       else if ((fpc == FP_NAN) || isneg_o)
1481       {
1482          result = double_float_type::my_value_nan();
1483
1484          return;
1485       }
1486       else if (fpc == FP_INFINITE)
1487       {
1488          result = double_float_type::my_value_inf();
1489
1490          return;
1491       }
1492    }
1493
1494    // TBD: Do we need any overflow/underflow guards when multiplying
1495    // by the split or when multiplying (hx * tx) and/or (hx * hx)?
1496
1497    const local_float_type c { cpp_df_qf_detail::ccmath::sqrt(o.crep().first) };
1498
1499    local_float_type p { cpp_df_qf_detail::split_maker<local_float_type>::value * c };
1500
1501    const local_float_type hx { local_float_type { c - p } + p };
1502    const local_float_type tx { c  - hx };
1503
1504    local_float_type q  = hx * tx;
1505
1506    q = q + q;
1507
1508    p = hx * hx;
1509
1510    const local_float_type u { p + q };
1511
1512    const local_float_type uu { cpp_df_qf_detail::ccmath::unsafe::fma(tx, tx, local_float_type { p - u } + q) };
1513
1514    const local_float_type
1515       cc
1516       {
1517         local_float_type { local_float_type { o.crep().first - u } - uu + o.crep().second } / local_float_type { c + c }
1518       };
1519
1520    result.rep().first  = c + cc;
1521    result.rep().second = local_float_type { c - result.my_first() } + cc;
1522 }
1523
1524 template <typename FloatingPointType>
1525 constexpr auto eval_pow(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x, const cpp_double_fp_backend<FloatingPointType>& a) -> void
1526 {
1527    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
1528
1529    constexpr double_float_type zero { 0 };
1530
1531    result = zero;
1532
1533    signed long long val_nll { };
1534
1535    eval_convert_to(&val_nll, a);
1536
1537    const int na { static_cast<int>(val_nll) };
1538
1539    const int fpc_a { eval_fpclassify(a) };
1540
1541    if((fpc_a == FP_NORMAL) && (a.compare(double_float_type(na)) == 0))
1542    {
1543       eval_pow(result, x, na);
1544    }
1545    else
1546    {
1547       constexpr double_float_type one { 1 };
1548
1549       const int fpc_x { eval_fpclassify(x) };
1550
1551       if (fpc_a == FP_ZERO)
1552       {
1553          // pow(base, +/-0) returns 1 for any base, even when base is NaN.
1554
1555          result = one;
1556       }
1557       else if (fpc_x == FP_ZERO)
1558       {
1559          if ((fpc_a == FP_NORMAL) || (fpc_a == FP_INFINITE))
1560          {
1561             // pow(+/-0, exp), where exp is negative and finite, returns +infinity.
1562             // pow(+/-0, exp), where exp is positive non-integer, returns +0.
1563
1564             // pow(+/-0, -infinity) returns +infinity.
1565             // pow(+/-0, +infinity) returns +0.
1566
1567             result = (eval_signbit(a) ? double_float_type::my_value_inf() : zero);
1568          }
1569          else if (fpc_a == FP_NAN)
1570          {
1571             result = double_float_type::my_value_nan();
1572          }
1573       }
1574       else if (fpc_x == FP_INFINITE)
1575       {
1576          if ((fpc_a == FP_NORMAL) || (fpc_a == FP_INFINITE))
1577          {
1578             // pow(+infinity, exp) returns +0 for any negative exp.
1579             // pow(-infinity, exp) returns +infinity for any positive exp.
1580
1581             result = (eval_signbit(a) ? zero : double_float_type::my_value_inf());
1582          }
1583          else if (fpc_a == FP_NAN)
1584          {
1585             result = double_float_type::my_value_nan();
1586          }
1587       }
1588       else if (fpc_x != FP_NORMAL)
1589       {
1590          result = x;
1591       }
1592       else
1593       {
1594          if (fpc_a == FP_INFINITE)
1595          {
1596             constexpr double_float_type one_minus { -1 };
1597
1598             if (x.compare(one_minus) == 0)
1599             {
1600                result = one;
1601             }
1602             else
1603             {
1604                double_float_type xabs { };
1605
1606                eval_fabs(xabs, x);
1607
1608                const int compare_one_result { xabs.compare(one) };
1609
1610                result =
1611                   (
1612                        (compare_one_result < 0) ? (eval_signbit(a) ? double_float_type::my_value_inf() : zero)
1613                      : (compare_one_result > 0) ? (eval_signbit(a) ? zero : double_float_type::my_value_inf())
1614                      : one
1615                   );
1616             }
1617          }
1618          else if (fpc_a == FP_NAN)
1619          {
1620             result = (x.compare(one) == 0) ? one : double_float_type::my_value_nan();
1621          }
1622          else
1623          {
1624             double_float_type log_x { };
1625
1626             eval_log(log_x, x);
1627
1628             double_float_type a_log_x { };
1629
1630             eval_multiply(a_log_x, a, log_x);
1631
1632             eval_exp(result, a_log_x);
1633          }
1634       }
1635    }
1636 }
1637
1638 template <typename FloatingPointType,
1639           typename IntegralType>
1640 constexpr auto eval_pow(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x, IntegralType p) -> typename ::std::enable_if<::boost::multiprecision::detail::is_integral<IntegralType>::value, void>::type
1641 {
1642    const int fpc { eval_fpclassify(x) };
1643
1644    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
1645
1646    using local_integral_type = IntegralType;
1647
1648    const bool p_is_odd { (static_cast<local_integral_type>(p & 1) != static_cast<local_integral_type>(0)) };
1649
1650    if (p == static_cast<local_integral_type>(0))
1651    {
1652       // pow(base, +/-0) returns 1 for any base, even when base is NaN.
1653
1654       result = double_float_type { unsigned { UINT8_C(1) } };
1655    }
1656    else if (fpc != FP_NORMAL)
1657    {
1658       if (fpc == FP_ZERO)
1659       {
1660          // pow(  +0, exp), where exp is a negative odd  integer, returns +infinity.
1661          // pow(  -0, exp), where exp is a negative odd  integer, returns +infinity.
1662          // pow(+/-0, exp), where exp is a negative even integer, returns +infinity.
1663
1664          // pow(  +0, exp), where exp is a positive odd  integer, returns +0.
1665          // pow(  -0, exp), where exp is a positive odd  integer, returns -0.
1666          // pow(+/-0, exp), where exp is a positive even integer, returns +0.
1667
1668          result = ((p < static_cast<local_integral_type>(0)) ? double_float_type::my_value_inf() : double_float_type { unsigned { UINT8_C(0) } });
1669       }
1670       else if (fpc == FP_INFINITE)
1671       {
1672          if (eval_signbit(x))
1673          {
1674             if (p < static_cast<local_integral_type>(0))
1675             {
1676                // pow(-infinity, exp) returns -0 if exp is a negative odd integer.
1677                // pow(-infinity, exp) returns +0 if exp is a negative even integer.
1678
1679                result = double_float_type { unsigned { UINT8_C(0) } };
1680             }
1681             else
1682             {
1683                // pow(-infinity, exp) returns -infinity if exp is a positive odd integer.
1684                // pow(-infinity, exp) returns +infinity if exp is a positive even integer.
1685
1686                result = (p_is_odd ? -double_float_type::my_value_inf() : double_float_type::my_value_inf());
1687             }
1688          }
1689          else
1690          {
1691             // pow(+infinity, exp) returns +0 for any negative exp.
1692             // pow(+infinity, exp) returns +infinity for any positive exp.
1693
1694             result = ((p < static_cast<local_integral_type>(0)) ? double_float_type { unsigned { UINT8_C(0) } } : double_float_type::my_value_inf());
1695          }
1696       }
1697       else
1698       {
1699          result = double_float_type::my_value_nan();
1700       }
1701    }
1702    else
1703    {
1704       double_float_type::pown(result, x, p);
1705    }
1706 }
1707
1708 template <typename FloatingPointType,
1709           typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) < 16))>::type const*>
1710 constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
1711 {
1712    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
1713
1714    constexpr double_float_type one { 1 };
1715
1716    const int fpc { eval_fpclassify(x) };
1717
1718    const bool b_neg { x.isneg_unchecked() };
1719
1720    if (fpc == FP_ZERO)
1721    {
1722       result = one;
1723    }
1724    else if (fpc != FP_NORMAL)
1725    {
1726       if (fpc == FP_INFINITE)
1727       {
1728          result = (b_neg ? double_float_type { 0.0F } : double_float_type::my_value_inf());
1729       }
1730       else if (fpc == FP_NAN)
1731       {
1732          result = x;
1733       }
1734    }
1735    else
1736    {
1737       using local_float_type = typename double_float_type::float_type;
1738
1739       // Get a local copy of the argument and force it to be positive.
1740       const double_float_type xx { (!b_neg) ? x : -x };
1741
1742       // Check the range of the input.
1743       if (eval_lt(x, double_float_type::my_value_logmin()))
1744       {
1745          result = double_float_type(0U);
1746       }
1747       else if (eval_gt(xx, double_float_type::my_value_logmax()))
1748       {
1749          result = double_float_type::my_value_inf();
1750       }
1751       else if (xx.is_one())
1752       {
1753          if(!b_neg)
1754          {
1755             result = cpp_df_qf_detail::constant_df_exp1<local_float_type>();
1756          }
1757          else
1758          {
1759             eval_divide(result, one, cpp_df_qf_detail::constant_df_exp1<local_float_type>());
1760          }
1761       }
1762       else
1763       {
1764          // Use an argument reduction algorithm for exp() in classic MPFUN-style.
1765
1766          double_float_type nf { };
1767
1768          // Prepare the scaled variables.
1769          const bool b_scale { (xx.order02() > -1) };
1770
1771          double_float_type r { };
1772
1773          if (b_scale)
1774          {
1775             eval_floor(nf, xx / cpp_df_qf_detail::constant_df_ln_two<local_float_type>());
1776
1777             eval_ldexp(r, xx - (nf * cpp_df_qf_detail::constant_df_ln_two<local_float_type>()), -2);
1778          }
1779          else
1780          {
1781             r = xx;
1782          }
1783
1784          // PadeApproximant[Exp[x] - 1, {x, 0, {6, 6}}]
1785          // FullSimplify[%]
1786          //   (84 x (7920 + 240 x^2 + x^4))
1787          // / (665280 + x (-332640 + x (75600 + x (-10080 + x (840 + (-42 + x) x)))))
1788
1789          constexpr double_float_type n84(84);
1790          constexpr double_float_type n240(240);
1791          constexpr double_float_type n7920(7920);
1792
1793          constexpr double_float_type n665280(665280);
1794          constexpr double_float_type n332640(332640);
1795          constexpr double_float_type n75600(75600);
1796          constexpr double_float_type n10080(10080);
1797          constexpr double_float_type n840(840);
1798          constexpr double_float_type n42(42);
1799
1800          const double_float_type r2 { r * r };
1801
1802          // Use the small-argument Pade approximation having coefficients shown above.
1803          result = (n84 * r * (n7920 + (n240 + r2) * r2));
1804          const double_float_type bot = (n665280 + r * (-n332640 + r * (n75600 + r * (-n10080 + r * (n840 + (-n42 + r) * r)))));
1805
1806          eval_divide(result, bot);
1807          result.add_unchecked_limb(local_float_type { 1.0F });
1808
1809          // Rescale the result.
1810          if (b_scale)
1811          {
1812             result *= result;
1813             result *= result;
1814
1815             signed long long lln { };
1816
1817             eval_convert_to(&lln, nf);
1818
1819             const int n { static_cast<int>(lln) };
1820
1821             if (n > 0)
1822             {
1823                eval_ldexp(result, result, n);
1824             }
1825          }
1826
1827          if (b_neg)
1828          {
1829             eval_divide(result, one, result);
1830          }
1831       }
1832    }
1833 }
1834
1835 template <typename FloatingPointType,
1836           typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && (((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) <= 36)))>::type const*>
1837 constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
1838 {
1839    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
1840
1841    constexpr double_float_type one { 1 };
1842
1843    const int fpc { eval_fpclassify(x) };
1844
1845    const bool b_neg { x.isneg_unchecked() };
1846
1847    if (fpc == FP_ZERO)
1848    {
1849       result = one;
1850    }
1851    else if (fpc != FP_NORMAL)
1852    {
1853       if (fpc == FP_INFINITE)
1854       {
1855          result = (b_neg ? double_float_type(0) : double_float_type::my_value_inf());
1856       }
1857       else if (fpc == FP_NAN)
1858       {
1859          result = x;
1860       }
1861    }
1862    else
1863    {
1864       using local_float_type  = typename double_float_type::float_type;
1865
1866       // Get a local copy of the argument and force it to be positive.
1867       const double_float_type xx { (!b_neg) ? x : -x };
1868
1869       // Check the range of the input.
1870       if (eval_lt(x, double_float_type::my_value_logmin()))
1871       {
1872          result = double_float_type(0U);
1873       }
1874       else if (eval_gt(xx, double_float_type::my_value_logmax()))
1875       {
1876          result = double_float_type::my_value_inf();
1877       }
1878       else if (xx.is_one())
1879       {
1880          if(!b_neg)
1881          {
1882             result = cpp_df_qf_detail::constant_df_exp1<local_float_type>();
1883          }
1884          else
1885          {
1886             eval_divide(result, one, cpp_df_qf_detail::constant_df_exp1<local_float_type>());
1887          }
1888       }
1889       else
1890       {
1891          // Use an argument reduction algorithm for exp() in classic MPFUN-style.
1892
1893          double_float_type nf { };
1894
1895          // Prepare the scaled variables.
1896          const bool b_scale { (xx.order02() > -4) };
1897
1898          double_float_type r { };
1899
1900          if (b_scale)
1901          {
1902             eval_floor(nf, xx / cpp_df_qf_detail::constant_df_ln_two<local_float_type>());
1903
1904             eval_ldexp(r, xx - (nf * cpp_df_qf_detail::constant_df_ln_two<local_float_type>()), -4);
1905          }
1906          else
1907          {
1908             r = xx;
1909          }
1910
1911          // PadeApproximant[Exp[r], {r, 0, 8, 8}]
1912          // FullSimplify[%]
1913
1914          constexpr double_float_type n144(144U);
1915          constexpr double_float_type n3603600(3603600UL);
1916          constexpr double_float_type n120120(120120UL);
1917          constexpr double_float_type n770(770U);
1918
1919          constexpr double_float_type n518918400(518918400UL);
1920          constexpr double_float_type n259459200(259459200UL);
1921          constexpr double_float_type n60540480(60540480UL);
1922          constexpr double_float_type n8648640(8648640UL);
1923          constexpr double_float_type n831600(831600UL);
1924          constexpr double_float_type n55440(55440U);
1925          constexpr double_float_type n2520(2520U);
1926          constexpr double_float_type n72(72U);
1927
1928          const double_float_type r2 { r * r };
1929
1930          result = (n144 * r) * (n3603600 + r2 * (n120120 + r2 * (n770 + r2)));
1931          const double_float_type bot = (n518918400 + r * (-n259459200 + r * (n60540480 + r * (-n8648640 + r * (n831600 + r * (-n55440 + r * (n2520 + r * (-n72 + r))))))));
1932
1933          eval_divide(result, bot);
1934          result.add_unchecked_limb(local_float_type { 1.0F });
1935
1936          // Rescale the result.
1937          if (b_scale)
1938          {
1939             result *= result;
1940             result *= result;
1941             result *= result;
1942             result *= result;
1943
1944             signed long long lln { };
1945
1946             eval_convert_to(&lln, nf);
1947
1948             const int n { static_cast<int>(lln) };
1949
1950             if (n > 0)
1951             {
1952                eval_ldexp(result, result, n);
1953             }
1954          }
1955
1956          if (b_neg)
1957          {
1958             eval_divide(result, one, result);
1959          }
1960       }
1961    }
1962 }
1963
1964 template <typename FloatingPointType,
1965           typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) > 36))>::type const*>
1966 constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
1967 {
1968    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
1969
1970    constexpr double_float_type one { 1 };
1971
1972    const int fpc { eval_fpclassify(x) };
1973
1974    const bool b_neg { x.isneg_unchecked() };
1975
1976    if (fpc == FP_ZERO)
1977    {
1978       result = one;
1979    }
1980    else if (fpc != FP_NORMAL)
1981    {
1982       if (fpc == FP_INFINITE)
1983       {
1984          result = (b_neg ? double_float_type(0) : double_float_type::my_value_inf());
1985       }
1986       else if (fpc == FP_NAN)
1987       {
1988          result = x;
1989       }
1990    }
1991    else
1992    {
1993       using local_float_type = typename double_float_type::float_type;
1994
1995       // Get a local copy of the argument and force it to be positive.
1996       const double_float_type xx { (!b_neg) ? x : -x };
1997
1998       // Check the range of the input.
1999       if (eval_lt(x, double_float_type::my_value_logmin()))
2000       {
2001          result = double_float_type(0U);
2002       }
2003       else if (eval_gt(xx, double_float_type::my_value_logmax()))
2004       {
2005          result = double_float_type::my_value_inf();
2006       }
2007       else if (xx.is_one())
2008       {
2009          if(!b_neg)
2010          {
2011             result = cpp_df_qf_detail::constant_df_exp1<local_float_type>();
2012          }
2013          else
2014          {
2015             eval_divide(result, one, cpp_df_qf_detail::constant_df_exp1<local_float_type>());
2016          }
2017       }
2018       else
2019       {
2020          // Use an argument reduction algorithm for exp() in classic MPFUN-style.
2021
2022          double_float_type nf { };
2023
2024          // Prepare the scaled variables.
2025          const bool b_scale { (xx.order02() > -4) };
2026
2027          double_float_type xh { };
2028
2029          if (b_scale)
2030          {
2031             eval_floor(nf, xx / cpp_df_qf_detail::constant_df_ln_two<local_float_type>());
2032
2033             eval_ldexp(xh, xx - (nf * cpp_df_qf_detail::constant_df_ln_two<local_float_type>()), -4);
2034          }
2035          else
2036          {
2037             xh = xx;
2038          }
2039
2040          double_float_type x_pow_n_div_n_fact(xh);
2041
2042          result = x_pow_n_div_n_fact;
2043          result.add_unchecked_limb(local_float_type { 1.0F });
2044
2045          double_float_type dummy { };
2046
2047          // Use the Taylor series expansion of hypergeometric_0f0(; ; x).
2048          // For this high(er) digit count, a scaled argument with subsequent
2049          // Taylor series expansion is actually more precise than Pade approximation.
2050          for (unsigned n { UINT8_C(2) }; n < unsigned { UINT8_C(64) }; ++n)
2051          {
2052             eval_multiply(x_pow_n_div_n_fact, xh);
2053
2054             eval_divide(x_pow_n_div_n_fact, double_float_type(n));
2055
2056             int n_tol { };
2057
2058             eval_frexp(dummy, x_pow_n_div_n_fact, &n_tol);
2059
2060             if ((n > 4U) && (n_tol < -(double_float_type::my_digits - 1)))
2061             {
2062                break;
2063             }
2064
2065             eval_add(result, x_pow_n_div_n_fact);
2066          }
2067
2068          // Rescale the result.
2069          if (b_scale)
2070          {
2071             result *= result;
2072             result *= result;
2073             result *= result;
2074             result *= result;
2075
2076             signed long long lln { };
2077
2078             eval_convert_to(&lln, nf);
2079
2080             const int n { static_cast<int>(lln) };
2081
2082             if (n > 0)
2083             {
2084                eval_ldexp(result, result, n);
2085             }
2086          }
2087
2088          if (b_neg)
2089          {
2090             eval_divide(result, one, result);
2091          }
2092       }
2093    }
2094 }
2095
2096 template <typename FloatingPointType>
2097 constexpr auto eval_log(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
2098 {
2099    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
2100
2101    constexpr double_float_type one { 1 };
2102
2103    const int result_compare_with_one { x.compare(one) };
2104
2105    const int fpc { eval_fpclassify(x) };
2106
2107    if (fpc == FP_ZERO)
2108    {
2109       result = -double_float_type::my_value_inf();
2110    }
2111    else if (x.isneg_unchecked() || (fpc == FP_NAN))
2112    {
2113       result = double_float_type::my_value_nan();
2114    }
2115    else if (fpc == FP_INFINITE)
2116    {
2117       result = double_float_type::my_value_inf();
2118    }
2119    else if (result_compare_with_one == -1)
2120    {
2121       // Use argument inversion and negation of the result.
2122
2123       double_float_type x_inv { };
2124
2125       eval_divide(x_inv, one, x);
2126       eval_log(result, x_inv);
2127
2128       result.negate();
2129    }
2130    else if (result_compare_with_one == 1)
2131    {
2132       // Optimize to only use eval_frexp if (and only if) the exponent is
2133       // actually large/small in the sense of above/below a defined cutoff.
2134
2135       double_float_type x2 { };
2136
2137       int n2 { };
2138
2139       eval_frexp(x2, x, &n2);
2140
2141       // Get initial estimate using the self-written, detail math function log.
2142       using local_float_type = typename double_float_type::float_type;
2143
2144       const local_float_type s { cpp_df_qf_detail::ccmath::log(x2.my_first()) };
2145
2146       double_float_type E { };
2147
2148       eval_exp(E, double_float_type(s));
2149
2150       // Perform one single step of Newton-Raphson iteration.
2151       // result = s + (x2 - E) / E;
2152
2153       eval_subtract(result, x2, E);
2154       eval_divide(result, E);
2155       result.add_unchecked_limb(s);
2156
2157       double_float_type xn2 { n2 };
2158
2159       eval_multiply(xn2, cpp_df_qf_detail::constant_df_ln_two<local_float_type>());
2160
2161       eval_add(result, xn2);
2162    }
2163    else
2164    {
2165       result = double_float_type { 0 };
2166    }
2167 }
2168
2169 namespace detail {
2170
2171 template<typename DestType, typename FloatingPointType>
2172 constexpr auto extract(cpp_double_fp_backend<FloatingPointType>& source) -> DestType
2173 {
2174    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
2175    using local_float_type =  typename double_float_type::float_type;
2176
2177    using destination_type = DestType;
2178
2179    destination_type result { };
2180
2181    constexpr double_float_type zero { 0 };
2182
2183    result = static_cast<destination_type>(0);
2184
2185    unsigned fail_safe { UINT32_C(32) };
2186
2187    constexpr bool
2188       destination_type_is_longer
2189       {
2190          (std::numeric_limits<signed long long>::digits > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
2191       };
2192
2193    using float_extract_type = typename std::conditional<destination_type_is_longer, local_float_type, float>::type;
2194
2195    while((source.compare(zero) != 0) && (fail_safe > unsigned { UINT8_C(0) }))
2196    {
2197       const float_extract_type next_flt_val { static_cast<float_extract_type>(source.my_first()) };
2198
2199       result += static_cast<destination_type>(next_flt_val);
2200
2201       eval_subtract(source, double_float_type(next_flt_val));
2202
2203       --fail_safe;
2204    }
2205
2206    return result;
2207 }
2208
2209 } // namespace detail
2210
2211 template <typename FloatingPointType>
2212 constexpr auto eval_convert_to(signed long long* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> void
2213 {
2214    const auto fpc = eval_fpclassify(backend);
2215
2216    if (fpc != FP_NORMAL)
2217    {
2218       *result = static_cast<signed long long>(backend.crep().first);
2219    }
2220    else
2221    {
2222       using double_float_type = cpp_double_fp_backend<FloatingPointType>;
2223       using local_float_type =  typename double_float_type::float_type;
2224
2225       static_assert(std::is_same<local_float_type, FloatingPointType>::value, "Error something went wrong with the limb type");
2226
2227       constexpr bool
2228          long_long_is_longer
2229          {
2230             (std::numeric_limits<signed long long>::digits > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
2231          };
2232
2233       using longer_type = typename ::std::conditional<long_long_is_longer, signed long long, double_float_type>::type;
2234
2235       constexpr longer_type my_max_val { static_cast<longer_type>((std::numeric_limits<signed long long>::max)()) };
2236       constexpr longer_type my_min_val { static_cast<longer_type>((std::numeric_limits<signed long long>::min)()) };
2237
2238       constexpr double_float_type my_max_val_dd { static_cast<double_float_type>(my_max_val) };
2239       constexpr double_float_type my_min_val_dd { static_cast<double_float_type>(my_min_val) };
2240
2241       if (backend.compare(my_max_val_dd) >= 0)
2242       {
2243          *result = (std::numeric_limits<signed long long>::max)();
2244       }
2245       else if (backend.compare(my_min_val_dd) <= 0)
2246       {
2247          *result = (std::numeric_limits<signed long long>::min)();
2248       }
2249       else
2250       {
2251          double_float_type source { backend };
2252
2253          *result = detail::extract<signed long long>(source);
2254
2255          #if !defined(__x86_64__) && !defined(_M_X64)
2256
2257          // It has been "empirically found" that non-X64 needs this workaround.
2258          // Even though the same conditions are met for x86_64 on GCC and MSVC,
2259          // this workaround will actually break the long long conversion tests
2260          // on those platforms.
2261          //
2262          // Our assumption is that on x64 there is x87 math (double -> long double)
2263          // being performed in the background. Seemingly these might "aid" the conversion
2264          // of double value to long long. Somehow I get the feeling this issue will arise
2265          // in future evolution of the cpp_double_fp_backend.
2266          //
2267          // This workaround has been tested on: ARM64 (linux and mac), s390x and PPC64LE.
2268
2269          constexpr bool
2270             needs_workaround
2271             {
2272                   (sizeof(signed long long) == 8U)
2273                && (std::is_same<local_float_type, double>::value || (std::is_same<local_float_type, long double>::value))
2274             };
2275
2276          BOOST_IF_CONSTEXPR (needs_workaround)
2277          {
2278             // This is the last value stored in a double as 9223372036854775808
2279             constexpr signed long long upper_bound = 9223372036854775296LL;
2280
2281             if (!eval_signbit(backend) && *result >= upper_bound)
2282             {
2283                // LONG_MAX is stored with .second = -1, so we compensate for the offset
2284                // We also only need this at the upper end where the values aren't exactly
2285                // representable in double. Below a certain point we are fine
2286                *result += static_cast<signed long long>(1);
2287             }
2288          }
2289
2290          #endif // non-x64
2291       }
2292    }
2293 }
2294
2295 template <typename FloatingPointType>
2296 constexpr auto eval_convert_to(unsigned long long* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> void
2297 {
2298    const auto fpc = eval_fpclassify(backend);
2299
2300    if (fpc != FP_NORMAL)
2301    {
2302       *result = static_cast<unsigned long long>(backend.crep().first);
2303    }
2304    else
2305    {
2306       using double_float_type = cpp_double_fp_backend<FloatingPointType>;
2307       using local_float_type =  typename double_float_type::float_type;
2308
2309       static_assert(std::is_same<local_float_type, FloatingPointType>::value, "Error something went wrong with the limb type");
2310
2311       constexpr bool
2312          ulong_long_is_longer
2313          {
2314             (std::numeric_limits<unsigned long long>::digits > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
2315          };
2316
2317       using longer_type = typename ::std::conditional<ulong_long_is_longer, unsigned long long, double_float_type>::type;
2318
2319       constexpr longer_type my_max_val { static_cast<longer_type>((std::numeric_limits<unsigned long long>::max)()) };
2320
2321       constexpr double_float_type my_max_val_dd { static_cast<double_float_type>(my_max_val) };
2322
2323       if (backend.compare(my_max_val_dd) >= 0)
2324       {
2325          *result = (std::numeric_limits<unsigned long long>::max)();
2326       }
2327       else
2328       {
2329          double_float_type source { backend };
2330
2331          *result = detail::extract<unsigned long long>(source);
2332       }
2333    }
2334 }
2335
2336 #ifdef BOOST_HAS_INT128
2337 template <typename FloatingPointType>
2338 constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type
2339 {
2340    const auto fpc = eval_fpclassify(backend);
2341
2342    if (fpc != FP_NORMAL)
2343    {
2344       *result = static_cast<boost::int128_type>(backend.crep().first);
2345    }
2346    else
2347    {
2348       using double_float_type = cpp_double_fp_backend<FloatingPointType>;
2349       using local_float_type =  typename double_float_type::float_type;
2350
2351       static_assert(std::is_same<local_float_type, FloatingPointType>::value, "Error something went wrong with the limb type");
2352
2353       constexpr bool
2354          n128_is_longer
2355          {
2356             ((static_cast<int>(sizeof(boost::int128_type) * static_cast<std::size_t>(CHAR_BIT)) - 1) > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
2357          };
2358
2359       using longer_type = typename ::std::conditional<n128_is_longer, boost::int128_type, double_float_type>::type;
2360
2361       constexpr boost::int128_type my_max_val_n128 = (((static_cast<boost::int128_type>(1) << (sizeof(boost::int128_type) * CHAR_BIT - 2)) - 1) << 1) + 1;
2362       constexpr boost::int128_type my_min_val_n128 = static_cast<boost::int128_type>(-my_max_val_n128 - 1);
2363
2364       constexpr longer_type my_max_val(static_cast<longer_type>(my_max_val_n128));
2365       constexpr longer_type my_min_val(static_cast<longer_type>(my_min_val_n128));
2366
2367       constexpr double_float_type my_max_val_dd(static_cast<double_float_type>(my_max_val));
2368       constexpr double_float_type my_min_val_dd(static_cast<double_float_type>(my_min_val));
2369
2370       if (backend.compare(my_max_val_dd) >= 0)
2371       {
2372          *result = my_max_val;
2373       }
2374       else if (backend.compare(my_min_val_dd) <= 0)
2375       {
2376          *result = my_min_val;
2377       }
2378       else
2379       {
2380          double_float_type source { backend };
2381
2382          *result = detail::extract<boost::int128_type>(source);
2383       }
2384    }
2385 }
2386
2387 template <typename FloatingPointType>
2388 constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<!(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type
2389 {
2390    const auto fpc = eval_fpclassify(backend);
2391
2392    if (fpc != FP_NORMAL)
2393    {
2394       *result = static_cast<boost::int128_type>(backend.crep().first);
2395    }
2396    else
2397    {
2398       using double_float_type = cpp_double_fp_backend<FloatingPointType>;
2399
2400       double_float_type source { backend };
2401
2402       *result = detail::extract<boost::int128_type>(source);
2403    }
2404 }
2405
2406 template <typename FloatingPointType>
2407 constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type
2408 {
2409    const auto fpc = eval_fpclassify(backend);
2410
2411    if (fpc != FP_NORMAL)
2412    {
2413       *result = static_cast<boost::int128_type>(backend.crep().first);
2414    }
2415    else
2416    {
2417       using double_float_type = cpp_double_fp_backend<FloatingPointType>;
2418       using local_float_type =  typename double_float_type::float_type;
2419
2420       static_assert(std::is_same<local_float_type, FloatingPointType>::value, "Error something went wrong with the limb type");
2421
2422       constexpr bool
2423          u128_is_longer
2424          {
2425             (static_cast<int>(sizeof(boost::uint128_type) * static_cast<std::size_t>(CHAR_BIT)) > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
2426          };
2427
2428       using longer_type = typename ::std::conditional<u128_is_longer, boost::uint128_type, double_float_type>::type;
2429
2430       constexpr boost::uint128_type my_max_val_u128 = static_cast<boost::uint128_type>(~static_cast<boost::uint128_type>(0));
2431
2432       constexpr longer_type my_max_val(static_cast<longer_type>(my_max_val_u128));
2433
2434       constexpr double_float_type my_max_val_dd(static_cast<double_float_type>(my_max_val));
2435
2436       if (backend.compare(my_max_val_dd) >= 0)
2437       {
2438          *result = my_max_val;
2439       }
2440       else
2441       {
2442          double_float_type source { backend };
2443
2444          *result = detail::extract<boost::uint128_type>(source);
2445       }
2446    }
2447 }
2448
2449 template <typename FloatingPointType>
2450 constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<!(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type
2451 {
2452    const auto fpc = eval_fpclassify(backend);
2453
2454    if (fpc != FP_NORMAL)
2455    {
2456       *result = static_cast<boost::uint128_type>(backend.crep().first);
2457    }
2458    else
2459    {
2460       using double_float_type = cpp_double_fp_backend<FloatingPointType>;
2461
2462       double_float_type source { backend };
2463
2464       *result = detail::extract<boost::uint128_type>(source);
2465    }
2466 }
2467 #endif
2468
2469 template <typename FloatingPointType,
2470           typename OtherFloatingPointType>
2471 constexpr auto eval_convert_to(OtherFloatingPointType* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename ::std::enable_if<cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value>::type
2472 {
2473    const auto fpc = eval_fpclassify(backend);
2474
2475    // TBD: Implement min/max check for the destination floating-point type result.
2476
2477    if (fpc != FP_NORMAL)
2478    {
2479       *result = static_cast<OtherFloatingPointType>(backend.my_first());
2480    }
2481    else
2482    {
2483       BOOST_IF_CONSTEXPR(cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits > cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits)
2484       {
2485          *result  = static_cast<OtherFloatingPointType>(backend.my_first());
2486          *result += static_cast<OtherFloatingPointType>(backend.my_second());
2487       }
2488       else
2489       {
2490          cpp_double_fp_backend<FloatingPointType> source = backend;
2491
2492          *result = 0;
2493
2494          for(auto digit_count  = static_cast<int>(0);
2495                   digit_count  < cpp_double_fp_backend<FloatingPointType>::my_digits;
2496                   digit_count += cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits)
2497          {
2498             const auto next = static_cast<OtherFloatingPointType>(source.my_first());
2499
2500             *result += next;
2501
2502             eval_subtract(source, cpp_double_fp_backend<FloatingPointType>(next));
2503          }
2504       }
2505    }
2506 }
2507
2508 template <typename FloatingPointType>
2509 constexpr auto hash_value(const cpp_double_fp_backend<FloatingPointType>& a) -> ::std::size_t
2510 {
2511    return a.hash();
2512 }
2513
2514 } // namespace backends
2515
2516 using backends::cpp_double_fp_backend;
2517
2518 using cpp_double_float       = number<cpp_double_fp_backend<float>,       ::boost::multiprecision::et_off>;
2519 using cpp_double_double      = number<cpp_double_fp_backend<double>,      ::boost::multiprecision::et_off>;
2520 using cpp_double_long_double = number<cpp_double_fp_backend<long double>, ::boost::multiprecision::et_off>;
2521 #ifdef BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128
2522 using cpp_double_float128    = number<cpp_double_fp_backend<::boost::float128_type>, ::boost::multiprecision::et_off>;
2523 #endif
2524
2525 } } // namespace boost::multiprecision
2526
2527 namespace std {
2528
2529 // Specialization of numeric_limits for boost::multiprecision::number<cpp_double_fp_backend<>>
2530 template <typename FloatingPointType,
2531           const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2532 BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >
2533 {
2534  private:
2535    using local_float_type = FloatingPointType;
2536    using inner_self_type  = boost::multiprecision::cpp_double_fp_backend<local_float_type>;
2537
2538    using self_type = boost::multiprecision::number<inner_self_type, ExpressionTemplatesOption>;
2539
2540  public:
2541    static constexpr bool                    is_specialized                = true;
2542    static constexpr bool                    is_signed                     = true;
2543    static constexpr bool                    is_integer                    = false;
2544    static constexpr bool                    is_exact                      = false;
2545    static constexpr bool                    is_bounded                    = true;
2546    static constexpr bool                    is_modulo                     = false;
2547    static constexpr bool                    is_iec559                     = false;
2548    static constexpr bool                    has_infinity                  = true;
2549    static constexpr bool                    has_quiet_NaN                 = true;
2550    static constexpr bool                    has_signaling_NaN             = false;
2551
2552    // These members were deprecated in C++23, but only MSVC warns (as of June 25)
2553    #ifdef BOOST_MSVC
2554    #pragma warning(push)
2555    #pragma warning(disable:4996)
2556    #endif
2557
2558    static constexpr float_denorm_style      has_denorm                    = denorm_absent;
2559    static constexpr bool                    has_denorm_loss               = true;
2560
2561    #ifdef BOOST_MSVC
2562    #pragma warning(pop)
2563    #endif // deprecated members
2564
2565    static constexpr bool                    traps                         = false;
2566    static constexpr bool                    tinyness_before               = false;
2567    static constexpr float_round_style       round_style                   = round_to_nearest;
2568
2569    static constexpr int radix                          = 2;
2570    static constexpr int digits                         = inner_self_type::my_digits;
2571    static constexpr int digits10                       = inner_self_type::my_digits10;
2572    static constexpr int max_digits10                   = inner_self_type::my_max_digits10;
2573
2574    static constexpr int max_exponent                   = inner_self_type::my_max_exponent;
2575    static constexpr int min_exponent                   = inner_self_type::my_min_exponent;
2576    static constexpr int max_exponent10                 = inner_self_type::my_max_exponent10;
2577    static constexpr int min_exponent10                 = inner_self_type::my_min_exponent10;
2578
2579    static constexpr auto(min)         () noexcept -> self_type { return static_cast<self_type>(inner_self_type::my_value_min()); }
2580    static constexpr auto(max)         () noexcept -> self_type { return static_cast<self_type>(inner_self_type::my_value_max()); }
2581    static constexpr auto lowest       () noexcept -> self_type { return static_cast<self_type>(-(max)()); }
2582    static constexpr auto epsilon      () noexcept -> self_type { return static_cast<self_type>(inner_self_type::my_value_eps()); }
2583    static constexpr auto round_error  () noexcept -> self_type { return static_cast<self_type>(static_cast<local_float_type>(0.5)); }
2584    static constexpr auto denorm_min   () noexcept -> self_type { return static_cast<self_type>((min)()); }
2585    static constexpr auto infinity     () noexcept -> self_type { return static_cast<self_type>(inner_self_type::my_value_inf()); }
2586    static constexpr auto quiet_NaN    () noexcept -> self_type { return static_cast<self_type>(inner_self_type::my_value_nan()); }
2587    static constexpr auto signaling_NaN() noexcept -> self_type { return static_cast<self_type>(static_cast<local_float_type>(0.0)); }
2588 };
2589
2590 } // namespace std
2591
2592 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2593 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_specialized;
2594 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2595 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_signed;
2596 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2597 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_integer;
2598 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2599 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_exact;
2600 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2601 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_bounded;
2602 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2603 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_modulo;
2604 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2605 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_iec559;
2606
2607 #ifdef BOOST_MSVC
2608 #pragma warning(push)
2609 #pragma warning(disable:4996)
2610 #endif
2611
2612 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2613 constexpr std::float_denorm_style std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_denorm;
2614 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2615 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_denorm_loss;
2616
2617 #ifdef BOOST_MSVC
2618 #pragma warning(pop)
2619 #endif // deprecated members
2620
2621 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2622 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_infinity;
2623 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2624 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_quiet_NaN;
2625 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2626 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_signaling_NaN;
2627 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2628 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::traps;
2629 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2630 constexpr bool                    std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::tinyness_before;
2631 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2632 constexpr std::float_round_style  std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::round_style;
2633 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2634 constexpr int                     std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::radix;
2635 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2636 constexpr int                     std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::digits;
2637 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2638 constexpr int                     std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::digits10;
2639 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2640 constexpr int                     std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::max_digits10;
2641 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2642 constexpr int                     std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::max_exponent;
2643 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2644 constexpr int                     std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::min_exponent;
2645 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2646 constexpr int                     std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::max_exponent10;
2647 template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
2648 constexpr int                     std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::min_exponent10;
2649
2650 #if defined(BOOST_MP_MATH_AVAILABLE)
2651 namespace boost { namespace math { namespace policies {
2652
2653 template <class FloatingPointType, class Policy, boost::multiprecision::expression_template_option ExpressionTemplates>
2654 struct precision<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplates>, Policy>
2655 {
2656 private:
2657    using my_multiprecision_backend_type = boost::multiprecision::cpp_double_fp_backend<FloatingPointType>;
2658
2659    using digits2_type = digits2<my_multiprecision_backend_type::my_digits>;
2660
2661    static constexpr auto use_full_precision() noexcept -> bool
2662    {
2663       return ((digits2_type::value <= precision_type::value) || (precision_type::value <= 0));
2664    }
2665
2666 public:
2667    using precision_type = typename Policy::precision_type;
2668
2669    using type =
2670       typename std::conditional<use_full_precision(),
2671                                 digits2_type,           // This is the default case: Use full precision for FloatingPointType.
2672                                 precision_type>::type;  // Here we find (and use) user-customized precision.
2673 };
2674
2675 } } } // namespace boost::math::policies
2676 #endif
2677
2678 #endif // BOOST_MP_CPP_DOUBLE_FP_2021_06_05_HPP