diff --git a/CMakeLists.txt b/CMakeLists.txt index ea70e16..26731a5 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -320,6 +320,13 @@ endif() target_compile_definitions(infinipaint PRIVATE BOOST_MP_STANDALONE EIGEN_MPL2_ONLY) +include(CheckCCompilerFlag) +check_c_compiler_flag(-Wdeprecated-declarations HAS_DEPRECATED_DECLARATIONS) + +if(HAS_DEPRECATED_DECLARATIONS) + target_compile_options(infinipaint PRIVATE -Wno-deprecated-declarations) +endif() + if(APPLE OR WIN32) install(TARGETS infinipaint DESTINATION .) diff --git a/include/boost/config/stdlib/libcpp.hpp b/include/boost/config/stdlib/libcpp.hpp index 0e9f244..777228f 100644 --- a/include/boost/config/stdlib/libcpp.hpp +++ b/include/boost/config/stdlib/libcpp.hpp @@ -177,4 +177,9 @@ #endif #endif +#if _LIBCPP_VERSION <= 170006 +// no std::ranges::join_view +# define BOOST_NO_CXX20_HDR_RANGES +#endif + // --- end --- diff --git a/include/boost/multiprecision/cpp_bin_float.hpp b/include/boost/multiprecision/cpp_bin_float.hpp index 9eb8f71..1472352 100644 --- a/include/boost/multiprecision/cpp_bin_float.hpp +++ b/include/boost/multiprecision/cpp_bin_float.hpp @@ -1,6 +1,6 @@ //////////////////////////////////////////////////////////////// -// Copyright 2013 - 2022 John Maddock. -// Copyright 2022 Christopher Kormanyos. +// Copyright 2013 - 2025 John Maddock. +// Copyright 2022 - 2025 Christopher Kormanyos. // Distributed under the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt // or copy at https://www.boost.org/LICENSE_1_0.txt) @@ -1614,6 +1614,10 @@ inline void convert_to_unsigned_int(I* res, const cpp_bin_float::exponent_nan: BOOST_MP_THROW_EXCEPTION(std::runtime_error("Could not convert NaN to integer.")); + #if defined(BOOST_NO_EXCEPTIONS) + *res = max_val; + return; + #endif case cpp_bin_float::exponent_infinity: *res = max_val; return; @@ -1622,15 +1626,34 @@ inline void convert_to_unsigned_int(I* res, const cpp_bin_float::rep_type man(arg.bits()); using shift_type = typename std::conditional::exponent_type) < sizeof(int), int, typename cpp_bin_float::exponent_type>::type; - shift_type shift = (shift_type)cpp_bin_float::bit_count - 1 - arg.exponent(); - if (shift > (shift_type)cpp_bin_float::bit_count - 1) + + const shift_type shift = (shift_type)cpp_bin_float::bit_count - 1 - arg.exponent(); + + if (arg.sign()) + { + using si_type = typename boost::multiprecision::detail::make_signed::type; + + si_type val; + + convert_to_signed_int(&val, arg); + + *res = static_cast(val); + + return; + } + else if (shift > (shift_type)cpp_bin_float::bit_count - 1) { *res = 0; return; } + else if (arg.compare(max_val) >= 0) + { + *res = max_val; + return; + } else if (shift < 0) { - if (cpp_bin_float::bit_count - shift <= digits) + if ((shift_type)cpp_bin_float::bit_count - shift <= (shift_type)digits) { // We have more bits in ulong_long_type than the float, so it's OK to left shift: eval_convert_to(res, man); @@ -1640,7 +1663,7 @@ inline void convert_to_unsigned_int(I* res, const cpp_bin_float(shift)); eval_convert_to(res, man); } @@ -2031,18 +2054,18 @@ inline void eval_ceil(cpp_bin_float(shift)); if (fractional && !res.sign()) { eval_increment(res.bits()); - if ((std::ptrdiff_t)eval_msb(res.bits()) != cpp_bin_float::bit_count - 1 - shift) + if ((std::ptrdiff_t)eval_msb(res.bits()) != (std::ptrdiff_t)(static_cast(cpp_bin_float::bit_count) - 1 - shift)) { // Must have extended result by one bit in the increment: --shift; ++res.exponent(); } } - eval_left_shift(res.bits(), shift); + eval_left_shift(res.bits(), static_cast(shift)); } template diff --git a/include/boost/multiprecision/cpp_dec_float.hpp b/include/boost/multiprecision/cpp_dec_float.hpp index 6270556..c9b958c 100644 --- a/include/boost/multiprecision/cpp_dec_float.hpp +++ b/include/boost/multiprecision/cpp_dec_float.hpp @@ -1,6 +1,6 @@ /////////////////////////////////////////////////////////////////////////////// -// Copyright Christopher Kormanyos 2002 - 2021. -// Copyright 2011 -2021 John Maddock. Distributed under the Boost +// Copyright Christopher Kormanyos 2002 - 2025. +// Copyright 2011 -2025 John Maddock. Distributed under the Boost // Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) // @@ -163,7 +163,7 @@ class cpp_dec_float a.prec_elem = cpp_dec_float_elem_number; return a; - } + } // LCOV_EXCL_LINE public: // Public Constructors @@ -501,7 +501,9 @@ class cpp_dec_float prec_elem = (std::min)(cpp_dec_float_elem_number, (std::max)(elems, static_cast(2))); } + static cpp_dec_float pow2(long long i); + exponent_type order() const { const bool bo_order_is_zero = ((!(isfinite)()) || (data[0] == static_cast(0u))); @@ -515,12 +517,7 @@ class cpp_dec_float if (data[0] >= 10000000UL) { if (data[0] >= 100000000UL) - { - if (data[0] >= 1000000000UL) - prefix = 9; - else - prefix = 8; - } + prefix = 8; else prefix = 7; } @@ -1164,7 +1161,7 @@ cpp_dec_float& cpp_dec_float& cpp_dec_float(n); @@ -1645,7 +1642,7 @@ long double cpp_dec_float::extract_long_doubl // Check if *this cpp_dec_float exceeds the maximum of double. if (xx.compare(long_double_max()) > 0) { - return ((!neg) ? std::numeric_limits::infinity() + return ((!neg) ? std::numeric_limits::infinity() : -std::numeric_limits::infinity()); } @@ -1713,7 +1710,7 @@ long long cpp_dec_float::extract_signed_long_ // See https://svn.boost.org/trac/boost/ticket/9740. // long long sval = static_cast(val - 1); - sval = -sval; + sval = -sval; --sval; return sval; } @@ -1749,14 +1746,14 @@ unsigned long long cpp_dec_float::extract_uns else { // Extract the data into an unsigned long long value. - val = static_cast(xn.data[0]); + val = static_cast(xn.data[std::size_t { 0U }]); const std::int32_t imax = (std::min)(static_cast(static_cast(xn.exp) / cpp_dec_float_elem_digits10), static_cast(cpp_dec_float_elem_number - static_cast(1))); for (std::int32_t i = static_cast(1); i <= imax; i++) { val *= static_cast(cpp_dec_float_elem_mask); - val += static_cast(xn.data[i]); + val += static_cast(xn.data[static_cast(i)]); } } @@ -1795,7 +1792,7 @@ int128_type cpp_dec_float::extract_signed_int } else { - // Extract the data into an unsigned long long value. + // Extract the data into an (unsigned) boost::uint128_type value. cpp_dec_float xn(extract_integer_part()); if (xn.isneg()) xn.negate(); @@ -1817,9 +1814,9 @@ int128_type cpp_dec_float::extract_signed_int } else { - // This strange expression avoids a hardware trap in the corner case - // that val is the most negative value permitted in long long. - // See https://svn.boost.org/trac/boost/ticket/9740. + // This strange expression avoids a hardware trap in the corner case that + // val is the most negative value permitted in (signed) boost::int128_type. + // See also https://svn.boost.org/trac/boost/ticket/9740. // int128_type sval = static_cast(val - 1); sval = -sval; @@ -3157,7 +3154,7 @@ cpp_dec_float cpp_dec_float::two(), static_cast(p), std::integral_constant()); return t; -} +} // LCOV_EXCL_LINE template inline void eval_add(cpp_dec_float& result, const cpp_dec_float& o) @@ -3592,7 +3589,7 @@ class numeric_limits, ExpressionTemplates>(min)() { return (boost::multiprecision::cpp_dec_float::min)(); } static constexpr boost::multiprecision::number, ExpressionTemplates>(max)() { return (boost::multiprecision::cpp_dec_float::max)(); } - static constexpr boost::multiprecision::number, ExpressionTemplates> lowest() { return boost::multiprecision::cpp_dec_float::zero(); } + static constexpr boost::multiprecision::number, ExpressionTemplates> lowest() { return -(max)(); } static constexpr boost::multiprecision::number, ExpressionTemplates> epsilon() { return boost::multiprecision::cpp_dec_float::eps(); } static constexpr boost::multiprecision::number, ExpressionTemplates> round_error() { return 0.5L; } static constexpr boost::multiprecision::number, ExpressionTemplates> infinity() { return boost::multiprecision::cpp_dec_float::inf(); } diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail.hpp new file mode 100644 index 0000000..a95f247 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail.hpp @@ -0,0 +1,169 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2021 Fahad Syed. +// Copyright 2021 - 2025 Christopher Kormanyos. +// Copyright 2021 Janek Kozicki. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_2023_01_02_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_2023_01_02_HPP + +#include + +#if defined(BOOST_HAS_FLOAT128) +#if defined(__has_include) +#if __has_include() + +#include + +#if !defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128) +#define BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128 +#endif + +#endif // __has_include() +#endif // defined(__has_include) +#endif // defined(BOOST_HAS_FLOAT128) + +#include +#include +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { + +template +struct pair +{ + static_assert(std::is_same::value, "Error: floating point types A and B must be identical"); + + using float_type = FloatingPointTypeA; + + float_type first; + float_type second; + + constexpr pair() : first { }, second { } { }; + constexpr pair(float_type a, float_type b) : first { a }, second { b } { } + constexpr pair(const pair& other) : first { other.first }, second { other.second } { } + constexpr pair(pair&& other) noexcept : first { other.first }, second { other.second } { } + + constexpr auto operator=(const pair& other) -> pair& + { + if (this != &other) + { + first = other.first; + second = other.second; + } + + return *this; + } + + constexpr auto operator=(pair&& other) noexcept -> pair& + { + first = other.first; + second = other.second; + + return *this; + } +}; + +template +struct is_floating_point +{ + static constexpr auto value = ::std::is_same::value + || ::std::is_same::value + || ::std::is_same::value +#if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128) + || ::std::is_same::value +#endif + ; +}; + +template +struct split_maker +{ +private: + using float_type = FloatType; + +public: + static constexpr int + n_shl + { + static_cast((ccmath::numeric_limits::digits + 1) / 2) + }; + + static_assert(n_shl < std::numeric_limits::digits, + "Error: Left-shift amount for split does not fin in std::uint64_t"); + + static constexpr float_type + value + { + static_cast + ( + std::uint64_t + { + UINT64_C(1) + + std::uint64_t { UINT64_C(1) << static_cast(n_shl) } + } + ) + }; +}; + +template +struct exact_arithmetic +{ + // The exact_arithmetic<> struct implements a few extended + // precision algorithms that are used in cpp_double_fp_backend. + + static_assert(is_floating_point::value, "Error: exact_arithmetic<> invoked with unknown floating-point type"); + + using float_type = FloatingPointType; + using float_pair = pair; + + static constexpr auto two_sum(const float_type a, const float_type b) -> float_pair + { + const float_type hi { a + b }; + const float_type a1 { hi - b }; + + return { hi, float_type { (a - a1) + (b - float_type { hi - a1 }) } }; + } + + static constexpr auto two_diff(const float_type a, const float_type b) -> float_pair + { + const float_type hi { a - b }; + const float_type a1 { hi + b }; + + return { hi, float_type { (a - a1) - (b + float_type { hi - a1 }) } }; + } + + static constexpr auto two_hilo_sum(const float_type a, const float_type b) -> float_pair + { + const float_type hi { a + b }; + + return { hi, float_type { b - (hi - a) } }; + } + + static constexpr auto normalize(float_type a, float_type b) -> float_pair + { + const float_type u { a + b }; + + return + { + u, + float_type { a - u } + b + }; + } +}; + +template +struct pow2_maker +{ + static constexpr auto value(const int power_value) noexcept -> ArithmeticType + { + return ((power_value == 0) ? ArithmeticType { 1 } : ArithmeticType { 2 } * pow2_maker::value(power_value - 1)); + } +}; + +} } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_2023_01_02_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath.hpp new file mode 100644 index 0000000..e5e05e5 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath.hpp @@ -0,0 +1,25 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2023 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_2023_01_07_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_2023_01_07_HPP + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_2023_01_07_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fabs.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fabs.hpp new file mode 100644 index 0000000..cfd3b52 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fabs.hpp @@ -0,0 +1,27 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2023 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FABS_2023_01_07_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FABS_2023_01_07_HPP + +#include +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +template +constexpr auto fabs(Real x) -> Real +{ + return (cpp_df_qf_detail::ccmath::isnan(x)) ? cpp_df_qf_detail::ccmath::numeric_limits::quiet_NaN() + : (x == static_cast(-0)) ? static_cast(0) + : (x >= 0) ? x + : -x; +}; + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FABS_2023_01_07_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_floor.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_floor.hpp new file mode 100644 index 0000000..4df2f46 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_floor.hpp @@ -0,0 +1,37 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2024 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FLOOR_2024_12_30_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FLOOR_2024_12_30_HPP + +#include +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +namespace detail { + +template +constexpr auto floor_impl(T x) -> T +{ + // Default to the regular floor function. + using std::floor; + + return floor(x); +} + +} // namespace detail + +template +constexpr auto floor(Real x) -> Real +{ + return cpp_df_qf_detail::ccmath::detail::floor_impl(x); +} + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FLOOR_2024_12_30_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fma.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fma.hpp new file mode 100644 index 0000000..f7117c8 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fma.hpp @@ -0,0 +1,61 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2024 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FMA_2024_12_17_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FMA_2024_12_17_HPP + +#include +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +namespace unsafe { + +namespace detail { + +template +constexpr auto fma_impl(T x, T y, T z) noexcept -> T +{ + // Default to the written-out operations ... + // ... and hope the compiler chooses FMA if available). + + return (x * y) + z; +} + +#if defined(__GNUC__) && !defined(__clang__) && !defined(__INTEL_COMPILER) && !defined(__INTEL_LLVM_COMPILER) +template <> +constexpr auto fma_impl(float x, float y, float z) noexcept -> float +{ + return __builtin_fmaf(x, y, z); +} + +template <> +constexpr auto fma_impl(double x, double y, double z) noexcept -> double +{ + return __builtin_fma(x, y, z); +} + +template <> +constexpr auto fma_impl(long double x, long double y, long double z) noexcept -> long double +{ + return __builtin_fmal(x, y, z); +} +#endif + +} // namespace detail + +template +constexpr auto fma(Real x, Real y, Real z) noexcept -> Real +{ + return detail::fma_impl(x, y, z); +} + +} // namespace unsafe + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FMA_2024_12_17_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fpclassify.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fpclassify.hpp new file mode 100644 index 0000000..7824e99 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fpclassify.hpp @@ -0,0 +1,51 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2024 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FPCLASSIFY_2024_12_26_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FPCLASSIFY_2024_12_26_HPP + +#include +#include +#include + +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +template +constexpr auto fpclassify(T x) -> int +{ + if ((::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::isnan)(x)) + { + return FP_NAN; + } + else if ((::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::isinf)(x)) + { + return FP_INFINITE; + } + else + { + const T fabs_x { ::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::fabs(x) }; + + if (fabs_x == T(0)) + { + return FP_ZERO; + } + else if ((fabs_x > 0) && (fabs_x < (::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits::min)())) + { + return FP_SUBNORMAL; + } + else + { + return FP_NORMAL; + } + } +} + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FPCLASSIFY_2024_12_26_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_frexp.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_frexp.hpp new file mode 100644 index 0000000..2ca53f3 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_frexp.hpp @@ -0,0 +1,38 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2023 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FREXP_2023_01_07_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FREXP_2023_01_07_HPP + +#include +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +namespace detail +{ + +template +constexpr auto frexp_impl(T arg, int* expptr) -> T +{ + // Default to the regular frexp function. + using std::frexp; + + return frexp(arg, expptr); +} + +} // namespace detail + +template +constexpr auto frexp(Real arg, int* expptr) -> Real +{ + return detail::frexp_impl(arg, expptr); +} + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FREXP_2023_01_07_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_isinf.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_isinf.hpp new file mode 100644 index 0000000..fdd5ce2 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_isinf.hpp @@ -0,0 +1,24 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2023 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_ISINF_2023_01_07_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_ISINF_2023_01_07_HPP + +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +template +constexpr auto isinf(FloatingPointType x) noexcept -> bool +{ + return ( ( x == cpp_df_qf_detail::ccmath::numeric_limits::infinity()) + || (-x == cpp_df_qf_detail::ccmath::numeric_limits::infinity())); +} + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_ISINF_2023_01_07_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_isnan.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_isnan.hpp new file mode 100644 index 0000000..468a97b --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_isnan.hpp @@ -0,0 +1,21 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2023 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_ISNAN_2023_01_07_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_ISNAN_2023_01_07_HPP + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +template +constexpr auto isnan(FloatingPointType x) -> bool +{ + return (x != x); +} + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_ISNAN_2023_01_07_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_ldexp.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_ldexp.hpp new file mode 100644 index 0000000..da25109 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_ldexp.hpp @@ -0,0 +1,58 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2023 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LDEXP_2023_01_07_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LDEXP_2023_01_07_HPP + +#include +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +namespace detail { + +template +constexpr auto ldexp_impl(T arg, int expval) -> T +{ + // Default to the regular ldexp function. + using std::ldexp; + + return ldexp(arg, expval); +} + +} // Namespace detail + +template +constexpr auto ldexp(Real arg, int expval) -> Real +{ + return detail::ldexp_impl(arg, expval); +} + +namespace unsafe { + +template +constexpr auto ldexp(Real arg, int exp) noexcept -> Real +{ + while(exp > 0) + { + arg *= 2; + --exp; + } + while(exp < 0) + { + arg /= 2; + ++exp; + } + + return arg; +} + +} // namespace unsafe + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LDEXP_2023_01_07_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_limits.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_limits.hpp new file mode 100644 index 0000000..92593a1 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_limits.hpp @@ -0,0 +1,194 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2023 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LIMITS_2023_01_07_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LIMITS_2023_01_07_HPP + +#include +#include +#include +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +template +struct numeric_limits +{ + static constexpr bool is_specialized = false; +}; + +template +struct numeric_limits::value + || std::is_same::value + || std::is_same::value)>::type> +{ +private: + using float_type = FloatingPointType; + +public: + static constexpr bool is_specialized = std::numeric_limits::is_specialized; + static constexpr bool is_signed = std::numeric_limits::is_signed; + static constexpr bool is_integer = std::numeric_limits::is_integer; + static constexpr bool is_exact = std::numeric_limits::is_exact; + static constexpr bool is_bounded = std::numeric_limits::is_bounded; + static constexpr bool is_modulo = std::numeric_limits::is_modulo; + static constexpr bool is_iec559 = std::numeric_limits::is_iec559; + static constexpr std::float_denorm_style has_denorm = std::numeric_limits::has_denorm; + static constexpr bool has_infinity = std::numeric_limits::has_infinity; + static constexpr bool has_quiet_NaN = std::numeric_limits::has_quiet_NaN; + static constexpr bool has_signaling_NaN = std::numeric_limits::has_signaling_NaN; + static constexpr bool has_denorm_loss = std::numeric_limits::has_denorm_loss; + static constexpr bool traps = std::numeric_limits::traps; + static constexpr bool tinyness_before = std::numeric_limits::tinyness_before; + static constexpr std::float_round_style round_style = std::numeric_limits::round_style; + + static constexpr int radix = std::numeric_limits::radix; + static constexpr int digits = std::numeric_limits::digits; + static constexpr int digits10 = std::numeric_limits::digits10; + static constexpr int max_digits10 = std::numeric_limits::max_digits10; + + static constexpr int max_exponent = std::numeric_limits::max_exponent; + static constexpr int min_exponent = std::numeric_limits::min_exponent; + static constexpr int max_exponent10 = std::numeric_limits::max_exponent10; + static constexpr int min_exponent10 = std::numeric_limits::min_exponent10; + + // LCOV_EXCL_START + static constexpr auto (min) () noexcept -> float_type { return (std::numeric_limits::min) (); } + static constexpr auto (max) () noexcept -> float_type { return (std::numeric_limits::max) (); } + static constexpr auto lowest () noexcept -> float_type { return std::numeric_limits::lowest (); } + static constexpr auto epsilon () noexcept -> float_type { return std::numeric_limits::epsilon (); } + static constexpr auto round_error () noexcept -> float_type { return std::numeric_limits::round_error (); } + static constexpr auto denorm_min () noexcept -> float_type { return std::numeric_limits::denorm_min (); } + static constexpr auto infinity () noexcept -> float_type { return std::numeric_limits::infinity (); } + static constexpr auto quiet_NaN () noexcept -> float_type { return std::numeric_limits::quiet_NaN (); } + static constexpr auto signaling_NaN() noexcept -> float_type { return std::numeric_limits::signaling_NaN(); } + // LCOV_EXCL_STOP +}; + +#if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128) +template +struct numeric_limits::value>::type> +{ +private: + using self_type = ::boost::float128_type; + +public: + static constexpr bool is_specialized = true; + static constexpr bool is_signed = true; + static constexpr bool is_integer = false; + static constexpr bool is_exact = false; + static constexpr bool is_bounded = true; + static constexpr bool is_modulo = false; + static constexpr bool is_iec559 = true; + static constexpr std::float_denorm_style has_denorm = std::denorm_present; + static constexpr bool has_infinity = true; + static constexpr bool has_quiet_NaN = true; + static constexpr bool has_signaling_NaN = false; + static constexpr bool has_denorm_loss = true; + static constexpr bool traps = false; + static constexpr bool tinyness_before = false; + static constexpr std::float_round_style round_style = std::round_to_nearest; + + static constexpr int radix = 2; + static constexpr int digits = 113; + static constexpr int digits10 = 33; + static constexpr int max_digits10 = 36; + + static constexpr int max_exponent = 16384; + static constexpr int min_exponent = -16381; + static constexpr int max_exponent10 = static_cast(static_cast(static_cast(max_exponent) * 301LL) / 1000LL); + static constexpr int min_exponent10 = static_cast(static_cast(static_cast(min_exponent) * 301LL) / 1000LL); + + static constexpr auto (min)() noexcept -> self_type + { + return static_cast(1) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) / 1073741824; + } + + static constexpr auto (max)() noexcept -> self_type + { + // This has one bit set only. + constexpr double dbl_mult = 8.9884656743115795386e+307; + + return (static_cast(1) - 9.62964972193617926527988971292463659e-35) // This now has all bits sets to 1 + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) + * static_cast(dbl_mult) * 65536; + } + + static constexpr auto lowest() noexcept -> self_type { return -(max)(); } + + static constexpr auto epsilon() -> self_type + { + // This double value has only one bit set and so is exact. + return 1.92592994438723585305597794258492732e-34; + } + + static constexpr auto round_error() noexcept -> self_type { return static_cast(0.5F); } + + static constexpr auto denorm_min() noexcept -> self_type + { + return static_cast(1) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + * static_cast(DBL_MIN) + / 5.5751862996326557854e+42; + } + + static constexpr auto infinity () noexcept -> self_type { return HUGE_VAL; } + static constexpr auto quiet_NaN () noexcept -> self_type { return NAN; } + static constexpr auto signaling_NaN() noexcept -> self_type { return 0; } +}; +#endif + +} } } } } // namespace boost::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LIMITS_2023_01_07_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_log.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_log.hpp new file mode 100644 index 0000000..64cd301 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_log.hpp @@ -0,0 +1,37 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2024 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LOG_2024_12_30_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LOG_2024_12_30_HPP + +#include +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +namespace detail { + +template +constexpr auto log_impl(T x) -> T +{ + // Default to the regular log function. + using std::log; + + return log(x); +} + +} // namespace detail + +template +constexpr auto log(Real x) -> Real +{ + return cpp_df_qf_detail::ccmath::detail::log_impl(x); +} + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LOG_2024_12_30_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_sqrt.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_sqrt.hpp new file mode 100644 index 0000000..381570b --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_sqrt.hpp @@ -0,0 +1,69 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2023 - 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_SQRT_2023_01_07_HPP +#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_SQRT_2023_01_07_HPP + +#include +#include + +namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath { + +namespace detail { + +template +constexpr auto sqrt_impl(T x) -> T +{ + // Default to the regular sqrt function. + using std::sqrt; + + return sqrt(x); +} + +} // namespace detail + +template +constexpr auto sqrt(Real x) -> Real +{ + return cpp_df_qf_detail::ccmath::detail::sqrt_impl(x); +} + +namespace unsafe { + +namespace detail { + +template +constexpr auto sqrt_impl_2(Real x, Real s, Real s2) noexcept -> Real +{ + return !(s < s2) ? s2 : sqrt_impl_2(x, (x / s + s) / 2, s); +} + +template +constexpr auto sqrt_impl_1(Real x, Real s) noexcept -> Real +{ + return sqrt_impl_2(x, (x / s + s) / 2, s); +} + +template +constexpr auto sqrt_impl(Real x) noexcept -> Real +{ + return sqrt_impl_1(x, x > 1 ? x : Real(1)); +} + +} // namespace detail + +template +constexpr auto sqrt(Real x) noexcept -> Real +{ + return detail::sqrt_impl(x); +} + +} // namespace unsafe + +} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_SQRT_2023_01_07_HPP diff --git a/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_constants.hpp b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_constants.hpp new file mode 100644 index 0000000..fed4682 --- /dev/null +++ b/include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_constants.hpp @@ -0,0 +1,80 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright Christopher Kormanyos 2025. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CONSTANTS_2025_06_23_HPP + #define BOOST_MP_CPP_DF_QF_DETAIL_CONSTANTS_2025_06_23_HPP + + #include + + #include + + #if (defined(__GNUC__) && defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128)) + // + // This is the only way we can avoid + // warning: non-standard suffix on floating constant [-Wpedantic] + // when building with -Wall -pedantic. Neither __extension__ + // nor #pragma diagnostic ignored work :( + // + #pragma GCC system_header + #endif + + namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { + + // N[Pi, 101] + // 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680 + + // 3.14159250259, 1.50995788317e-07 + // 3.141592653589793116, 1.2246467991473529607e-16 + // 3.14159265358979323851281, -5.01655761266833202345176e-20 + // 3.14159265358979323846264338327950279748, 8.67181013012378102479704402604335225411e-35 + + template constexpr auto constant_df_pi() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 24)), cpp_double_fp_backend>::type { return { static_cast(3.14159250259L), static_cast(1.50995788317e-07L) }; } + template constexpr auto constant_df_pi() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 53)), cpp_double_fp_backend>::type { return { static_cast(3.141592653589793116L), static_cast(1.2246467991473529607e-16L) }; } + template constexpr auto constant_df_pi() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 64)), cpp_double_fp_backend>::type { return { static_cast(3.14159265358979323851281L), static_cast(-5.01655761266833202345176e-20L) }; } + #if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128) + template constexpr auto constant_df_pi() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 113)), cpp_double_fp_backend>::type { return { static_cast(3.14159265358979323846264338327950279748Q), static_cast(8.67181013012378102479704402604335225411e-35Q) }; } + #else + template constexpr auto constant_df_pi() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 113)), cpp_double_fp_backend>::type { return { static_cast(3.14159265358979323846264338327950279748L), static_cast(8.67181013012378102479704402604335225411e-35L) }; } + #endif + + // N[Log[2], 101] + // 0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754 + + // 0.69314712286, 5.76999887869e-08 + // 0.6931471805599451752, 1.3421277060097865271e-16 + // 0.69314718055994530942869, -1.14583527267987328094768e-20 + // 0.693147180559945309417232121458176575084, -7.00813947454958516341266200877162272784e-36 + + template constexpr auto constant_df_ln_two() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 24)), cpp_double_fp_backend>::type { return { static_cast(0.69314712286L), static_cast(5.76999887869e-08L) }; } + template constexpr auto constant_df_ln_two() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 53)), cpp_double_fp_backend>::type { return { static_cast(0.6931471805599451752L), static_cast(1.3421277060097865271e-16L) }; } + template constexpr auto constant_df_ln_two() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 64)), cpp_double_fp_backend>::type { return { static_cast(0.69314718055994530942869L), static_cast(-1.14583527267987328094768e-20L) }; } + #if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128) + template constexpr auto constant_df_ln_two() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 113)), cpp_double_fp_backend>::type { return { static_cast(0.693147180559945309417232121458176575084Q), static_cast(-7.00813947454958516341266200877162272784e-36Q) }; } + #else + template constexpr auto constant_df_ln_two() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 113)), cpp_double_fp_backend>::type { return { static_cast(0.693147180559945309417232121458176575084L), static_cast(-7.00813947454958516341266200877162272784e-36L) }; } + #endif + + // N[Exp[1], 101] + // 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274 + + // 2.71828174591, 8.25483965627e-08 + // 2.7182818284590450908, 1.4456468917292501578e-16 + // 2.71828182845904523521133, 1.4895979785582304563159e-19 + // 2.71828182845904523536028747135266231436, 1.83398825226506410712297736767396397644e-34 + + template constexpr auto constant_df_exp1() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 24)), cpp_double_fp_backend>::type { return { static_cast(2.71828174591L), static_cast(8.25483965627e-08L) }; } + template constexpr auto constant_df_exp1() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 53)), cpp_double_fp_backend>::type { return { static_cast(2.7182818284590450908L), static_cast(1.4456468917292501578e-16L) }; } + template constexpr auto constant_df_exp1() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 64)), cpp_double_fp_backend>::type { return { static_cast(2.71828182845904523521133L), static_cast(1.4895979785582304563159e-19L) }; } + #if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128) + template constexpr auto constant_df_exp1() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 113)), cpp_double_fp_backend>::type { return { static_cast(2.71828182845904523536028747135266231436Q), static_cast(1.83398825226506410712297736767396397644e-34Q) }; } + #else + template constexpr auto constant_df_exp1() -> typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point::value && (cpp_df_qf_detail::ccmath::numeric_limits::digits == 113)), cpp_double_fp_backend>::type { return { static_cast(2.71828182845904523536028747135266231436L), static_cast(1.83398825226506410712297736767396397644e-34L) }; } + #endif + + } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail + +#endif // BOOST_MP_CPP_DF_QF_DETAIL_CONSTANTS_2025_06_23_HPP diff --git a/include/boost/multiprecision/cpp_double_fp.hpp b/include/boost/multiprecision/cpp_double_fp.hpp new file mode 100644 index 0000000..4e8245f --- /dev/null +++ b/include/boost/multiprecision/cpp_double_fp.hpp @@ -0,0 +1,2634 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2021 - 2025 Fahad Syed. +// Copyright 2021 - 2025 Christopher Kormanyos. +// Copyright 2021 - 2025 Janek Kozicki. +// Copyright 2025 Matt Borland. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// + +#ifndef BOOST_MP_CPP_DOUBLE_FP_2021_06_05_HPP +#define BOOST_MP_CPP_DOUBLE_FP_2021_06_05_HPP + +#include +#include +#include +#include +#include + +#ifdef BOOST_MP_MATH_AVAILABLE +// +// Headers required for Boost.Math integration: +// +#include +// +// Some includes we need from Boost.Math, since we rely on that library to provide these functions: +// +#include +#include +#include +#include +#include +#include +#endif + +#include +#include +#include + +#if (defined(__clang__) && (__clang_major__ <= 9)) +#define BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE struct +#else +#define BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE class +#endif + +namespace boost { namespace multiprecision { namespace backends { + +template +class cpp_double_fp_backend; + +template +constexpr auto eval_add(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; +template +constexpr auto eval_add(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void; +template +constexpr auto eval_subtract(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; +template +constexpr auto eval_subtract(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void; +template +constexpr auto eval_multiply(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; +template +constexpr auto eval_multiply(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void; +template +constexpr auto eval_divide(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; +template +constexpr auto eval_divide(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void; +template +constexpr auto eval_eq(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool; +template +constexpr auto eval_lt(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool; +template +constexpr auto eval_gt(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool; +template +constexpr auto eval_is_zero(const cpp_double_fp_backend& x) -> bool; +template +constexpr auto eval_signbit(const cpp_double_fp_backend& x) -> int; + +template +constexpr auto eval_fabs(cpp_double_fp_backend& result, const cpp_double_fp_backend& a) -> void; +template +constexpr auto eval_frexp(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, int* v) -> void; +template +constexpr auto eval_ldexp(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, int v) -> void; +template +constexpr auto eval_floor(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; +template +constexpr auto eval_ceil(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; +template +constexpr auto eval_fpclassify(const cpp_double_fp_backend& o) -> int; + +template +constexpr auto eval_sqrt(cpp_double_fp_backend& result, const cpp_double_fp_backend& o) -> void; + +template +constexpr auto eval_pow(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, const cpp_double_fp_backend& a) -> void; + +template +constexpr auto eval_pow(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, IntegralType n) -> typename ::std::enable_if::value, void>::type; + +template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) < 16))>::type const* = nullptr> +constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; + +template ::value && (((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) <= 36)))>::type const* = nullptr> +constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; + +template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) > 36))>::type const* = nullptr> +constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; + +template +constexpr auto eval_log(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; + +template +constexpr auto eval_convert_to(signed long long* result, const cpp_double_fp_backend& backend) -> void; + +template +constexpr auto eval_convert_to(unsigned long long* result, const cpp_double_fp_backend& backend) -> void; + +#ifdef BOOST_HAS_INT128 +template +constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits::digits > 24), void>::type; + +template +constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if::digits > 24), void>::type; + +template +constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits::digits > 24), void>::type; + +template +constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if::digits > 24), void>::type; +#endif + +template +constexpr auto eval_convert_to(OtherFloatingPointType* result, const cpp_double_fp_backend& backend) -> typename ::std::enable_if::value>::type; + +// TBD: constexpr on hash_value? +template +auto hash_value(const cpp_double_fp_backend& a) -> ::std::size_t; + +template +constexpr auto fabs(const cpp_double_fp_backend& a) -> cpp_double_fp_backend; + +} } } // namespace boost::multiprecision::backends + +namespace std { + +// Foward declarations of various specializations of std::numeric_limits + +template +BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE numeric_limits, ExpressionTemplatesOption> >; + +} // namespace std + +namespace boost { namespace multiprecision { + +template +struct number_category> : public std::integral_constant { }; + +namespace backends { + +// A cpp_double_fp_backend is represented by an unevaluated sum of two floating-point +// units (say a0 and a1) which satisfy |a1| <= (1 / 2) * ulp(a0). +// The type of the floating-point constituents should adhere to IEEE754. + +template +class cpp_double_fp_backend +{ + public: + using float_type = FloatingPointType; + + static_assert + ( + cpp_df_qf_detail::is_floating_point::value + && bool + { + (cpp_df_qf_detail::ccmath::numeric_limits::digits == 24) + || (cpp_df_qf_detail::ccmath::numeric_limits::digits == 53) + || (cpp_df_qf_detail::ccmath::numeric_limits::digits == 64) + || (cpp_df_qf_detail::ccmath::numeric_limits::digits == 113) + }, "Error: float_type does not fulfil the backend requirements of cpp_double_fp_backend" + ); + + using rep_type = cpp_df_qf_detail::pair; + using arithmetic = cpp_df_qf_detail::exact_arithmetic; + + using signed_types = std::tuple; + using unsigned_types = std::tuple; + using float_types = std::tuple; + using exponent_type = int; + + static constexpr int my_digits = static_cast(cpp_df_qf_detail::ccmath::numeric_limits::digits * static_cast(INT8_C(2))); + static constexpr int my_digits10 = static_cast(boost::multiprecision::detail::calc_digits10 (my_digits)>::value); + static constexpr int my_max_digits10 = static_cast(boost::multiprecision::detail::calc_max_digits10(my_digits)>::value); + static constexpr int my_max_exponent = cpp_df_qf_detail::ccmath::numeric_limits::max_exponent; + static constexpr int my_max_exponent10 = cpp_df_qf_detail::ccmath::numeric_limits::max_exponent10; + + static constexpr int my_min_exponent = + static_cast + ( + cpp_df_qf_detail::ccmath::numeric_limits::min_exponent + + cpp_df_qf_detail::ccmath::numeric_limits::digits + ); + + static constexpr int my_min_exponent10 = + static_cast + ( + -static_cast + ( + boost::multiprecision::detail::calc_digits10(-my_min_exponent)>::value + ) + ); + + // Default constructor. + constexpr cpp_double_fp_backend() noexcept { } + + // Copy constructor. + constexpr cpp_double_fp_backend(const cpp_double_fp_backend& other) : data(other.data) { } + + // Move constructor. + constexpr cpp_double_fp_backend(cpp_double_fp_backend&& other) noexcept : data(static_cast(other.data)) { } + + // Constructors from other floating-point types. + template ::value + && (cpp_df_qf_detail::ccmath::numeric_limits::digits <= cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> + constexpr cpp_double_fp_backend(const OtherFloatType& f) + : data(f, static_cast(0.0F)) { } + + template ::value + && (cpp_df_qf_detail::ccmath::numeric_limits::digits > cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> + constexpr cpp_double_fp_backend(const OtherFloatType& f) + : data(static_cast(f), + static_cast(f - static_cast(static_cast(f)))) { } + + // Construtor from another kind of cpp_double_fp_backend<> object. + + template ::value + && (!std::is_same::value))>::type const* = nullptr> + constexpr cpp_double_fp_backend(const cpp_double_fp_backend& a) + : cpp_double_fp_backend(a.my_first()) + { + // TBD: Maybe specialize this constructor for cases either wider or less wide. + operator+=(a.my_second()); + } + + // Constructors from integers. + template ::value + && (!boost::multiprecision::detail::is_unsigned::value) + && ((static_cast(sizeof(SignedIntegralType) * 8) - 1) <= cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> + constexpr cpp_double_fp_backend(const SignedIntegralType& n) + : data(static_cast(n), static_cast(0.0F)) { } + + template ::value + && boost::multiprecision::detail::is_unsigned::value + && (static_cast(sizeof(UnsignedIntegralType) * 8u) <= cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> + constexpr cpp_double_fp_backend(const UnsignedIntegralType& u) + : data(static_cast(u), static_cast(0.0F)) { } + + // Constructors from integers which hold more information than *this can contain. + template ::value + && boost::multiprecision::detail::is_unsigned::value + && (static_cast(sizeof(UnsignedIntegralType) * 8u) > cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> + constexpr cpp_double_fp_backend(UnsignedIntegralType u) + { + using local_unsigned_integral_type = UnsignedIntegralType; + + constexpr local_unsigned_integral_type + flt_mask + { + static_cast + ( + static_cast + ( + static_cast(UINT8_C(1)) << static_cast(std::numeric_limits::digits) + ) + - static_cast(UINT8_C(1)) + ) + }; + + data.second = static_cast(0.0F); + data.first = static_cast(u & flt_mask); + + int u_index { 1 }; + + while (u > static_cast(UINT8_C(0))) + { + u >>= static_cast(std::numeric_limits::digits); + + const float_type + xhi + { + static_cast(u & flt_mask) + * cpp_df_qf_detail::pow2_maker::value(std::numeric_limits::digits * u_index) + }; + + eval_add(*this, cpp_double_fp_backend(xhi)); + + ++u_index; + } + } + + template ::value + && (!boost::multiprecision::detail::is_unsigned::value) + && ((static_cast(sizeof(SignedIntegralType) * 8) - 1) > cpp_df_qf_detail::ccmath::numeric_limits::digits))>::type const* = nullptr> + constexpr cpp_double_fp_backend(SignedIntegralType n) + { + const bool is_neg { n < SignedIntegralType { INT8_C(0) } }; + + using local_unsigned_integral_type = typename boost::multiprecision::detail::make_unsigned::type; + + const local_unsigned_integral_type + u_val + { + (!is_neg) + ? static_cast(n) + : static_cast + ( + static_cast(~n) + + static_cast(UINT8_C(1)) + ) + }; + + data = cpp_double_fp_backend(u_val).data; + + if(is_neg) { negate(); } + } + + constexpr cpp_double_fp_backend(const float_type& a, const float_type& b) noexcept : data(a, b) { } + + constexpr cpp_double_fp_backend(const cpp_df_qf_detail::pair& p) noexcept : data(p) { } + + cpp_double_fp_backend(const char* s) + { + *this = s; + } + + // Assignment operator. + constexpr auto operator=(const cpp_double_fp_backend& other) -> cpp_double_fp_backend& + { + if (this != &other) + { + data = other.data; + } + + return *this; + } + + // Move assignment operator. + constexpr auto operator=(cpp_double_fp_backend&& other) noexcept -> cpp_double_fp_backend& + { + data = static_cast(other.data); + + return *this; + } + + // Assignment operator from another kind of cpp_double_fp_backend<> object. + template ::value + && (!std::is_same::value))>::type const* = nullptr> + constexpr auto operator=(const cpp_double_fp_backend& other) -> cpp_double_fp_backend& + { + return operator=(cpp_double_fp_backend(other)); + } + + template + constexpr auto operator=(const OtherFloatType f) -> typename ::std::enable_if::value, cpp_double_fp_backend&>::type + { + return operator=(cpp_double_fp_backend(f)); + } + + template + constexpr auto operator=(const IntegralType n) -> typename ::std::enable_if::value, cpp_double_fp_backend&>::type + { + return operator=(cpp_double_fp_backend(n)); + } + + auto operator=(const char* v) -> cpp_double_fp_backend& + { + rd_string(v); + + return *this; + } + + auto hash() const -> ::std::size_t + { + // Hash the raw values of the data field with direct-memory access. + // Use 16-bit (2 byte) chunks as the data size when hashing. + + static_assert( ( sizeof(data.first) == sizeof(data.second)) + && ( sizeof(float_type) >= sizeof(std::uint16_t)) + && ((sizeof(float_type) % sizeof(std::uint16_t)) == std::size_t { UINT8_C(0) }), + "Error: float_type size is inappropriate for hashing routine"); + + auto hash_one + { + [](std::size_t& res, const float_type& val) + { + const std::uint16_t* first { reinterpret_cast(&val) }; + const std::uint16_t* last { first + std::size_t { sizeof(float_type) / sizeof(std::uint16_t) } }; + + while (first != last) + { + boost::multiprecision::detail::hash_combine(res, *first); + + ++first; + } + + return res; + } + }; + + std::size_t result { UINT8_C(0) }; + + static_cast(hash_one(result, data.first)); + static_cast(hash_one(result, data.second)); + + return result; + } + + // The public methods follow. + + constexpr auto isneg_unchecked() const noexcept -> bool { return (data.first < 0); } + + constexpr auto iszero_unchecked() const noexcept -> bool { return (data.first == float_type { 0.0F }); } + + constexpr auto is_one() const noexcept -> bool { return ((data.second == float_type { 0.0F }) && (data.first == float_type { 1.0F })); } + + constexpr auto negate() -> void + { + const bool isinf_u { (cpp_df_qf_detail::ccmath::isinf)(data.first) }; + const bool isnan_u { (cpp_df_qf_detail::ccmath::isnan)(data.first) }; + + if (isnan_u) { } + else if (isinf_u) + { + data.first = -data.first; + } + else + { + data.first = -data.first; + + if (!iszero_unchecked()) + { + data.second = -data.second; + + data = arithmetic::normalize(data.first, data.second); + } + } + } + + // Getters/Setters + constexpr auto my_first () const noexcept -> const float_type& { return data.first; } + constexpr auto my_second() const noexcept -> const float_type& { return data.second; } + + constexpr auto rep() noexcept -> rep_type& { return data; } + + constexpr auto rep() const noexcept -> const rep_type& { return data; } + + constexpr auto crep() const noexcept -> const rep_type& { return data; } + + // Unary add/sub/mul/div follow in the upcoming paragraphs. + + constexpr auto operator+=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend& + { + const int fpc_u { eval_fpclassify(*this) }; + const int fpc_v { eval_fpclassify(v) }; + + if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL)) + { + // Handle special cases like zero, inf and NaN. + + if (fpc_u == FP_NAN) + { + return *this; + } + + const bool isinf_v { (fpc_v == FP_INFINITE) }; + + if (fpc_u == FP_INFINITE) + { + if (isinf_v && (isneg_unchecked() != v.isneg_unchecked())) + { + *this = cpp_double_fp_backend::my_value_nan(); + } + + return *this; + } + + const bool iszero_u { ((fpc_u == FP_ZERO) || (fpc_u == FP_SUBNORMAL)) }; + const bool isnan_v { (fpc_v == FP_NAN) }; + + if (iszero_u || (isnan_v || isinf_v)) + { + if (iszero_u) + { + data.first = float_type { 0.0F }; + data.second = float_type { 0.0F }; + } + + const bool iszero_v { ((fpc_v == FP_ZERO) || (fpc_v == FP_SUBNORMAL)) }; + + return ((!iszero_v) ? operator=(v) : *this); + } + } + + const float_type xlo { data.second }; + + data = arithmetic::two_sum(data.first, v.data.first); + + const rep_type thi_tlo { arithmetic::two_sum(xlo, v.data.second) }; + + data = arithmetic::two_hilo_sum(data.first, data.second + thi_tlo.first); + + data = arithmetic::two_hilo_sum(data.first, thi_tlo.second + data.second); + + return *this; + } + + constexpr auto operator-=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend& + { + const int fpc_u { eval_fpclassify(*this) }; + const int fpc_v { eval_fpclassify(v) }; + + if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL)) + { + // Handle special cases like zero, inf and NaN. + + if (fpc_u == FP_NAN) + { + return *this; + } + + const bool isinf_v { (fpc_v == FP_INFINITE) }; + + if (fpc_u == FP_INFINITE) + { + if (isinf_v && (isneg_unchecked() == v.isneg_unchecked())) + { + *this = cpp_double_fp_backend::my_value_nan(); + } + + return *this; + } + + const bool iszero_u { ((fpc_u == FP_ZERO) || (fpc_u == FP_SUBNORMAL)) }; + const bool isnan_v { (fpc_v == FP_NAN) }; + + if (iszero_u || (isnan_v || isinf_v)) + { + if (iszero_u) + { + data.first = float_type { 0.0F }; + data.second = float_type { 0.0F }; + } + + const bool iszero_v { ((fpc_v == FP_ZERO) || (fpc_v == FP_SUBNORMAL)) }; + + return ((!iszero_v) ? operator=(-v) : *this); + } + } + + if (this == &v) + { + data.first = float_type { 0.0F }; + data.second = float_type { 0.0F }; + + return *this; + } + + const float_type xlo { data.second }; + + data = arithmetic::two_diff(data.first, v.data.first); + + const rep_type thi_tlo { arithmetic::two_diff(xlo, v.data.second) }; + + data = arithmetic::two_hilo_sum(data.first, data.second + thi_tlo.first); + + data = arithmetic::two_hilo_sum(data.first, thi_tlo.second + data.second); + + return *this; + } + + constexpr auto operator*=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend& + { + // Evaluate the sign of the result. + + const int fpc_u { eval_fpclassify(*this) }; + const int fpc_v { eval_fpclassify(v) }; + + if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL)) + { + // Handle special cases like zero, inf and NaN. + const bool isinf_u { (fpc_u == FP_INFINITE) }; + const bool isinf_v { (fpc_v == FP_INFINITE) }; + const bool iszero_u { (fpc_u == FP_ZERO) }; + const bool iszero_v { (fpc_v == FP_ZERO) }; + + if (((fpc_u == FP_NAN) || (fpc_v == FP_NAN)) || (isinf_u && iszero_v) || (isinf_v && iszero_u)) + { + return operator=( cpp_double_fp_backend::my_value_nan()); + } + + if (isinf_u || isinf_v) + { + const bool b_neg { (isneg_unchecked() != v.isneg_unchecked()) }; + + *this = cpp_double_fp_backend::my_value_inf(); + + if (b_neg) + { + negate(); + } + + return *this; + } + + if (iszero_u || iszero_v) + { + return operator=(cpp_double_fp_backend(0)); + } + } + + mul_unchecked(v); + + return *this; + } + + constexpr auto operator/=(const cpp_double_fp_backend& v) -> cpp_double_fp_backend& + { + const int fpc_u { eval_fpclassify(*this) }; + const int fpc_v { eval_fpclassify(v) }; + + if ((fpc_u != FP_NORMAL) || (fpc_v != FP_NORMAL)) + { + // Handle special cases like zero, inf and NaN. + const bool isnan_u { (fpc_u == FP_NAN) }; + const bool isnan_v { (fpc_v == FP_NAN) }; + + if (isnan_u || isnan_v) + { + return operator=(cpp_double_fp_backend::my_value_nan()); + } + + const bool iszero_u { (fpc_u == FP_ZERO) }; + const bool iszero_v { (fpc_v == FP_ZERO) }; + + if (iszero_u) + { + if (iszero_v) + { + return operator=(cpp_double_fp_backend::my_value_nan()); + } + else + { + return operator=(cpp_double_fp_backend(0)); + } + } + + // Handle more special cases like zero, inf and NaN. + if (iszero_v) + { + const bool b_neg = isneg_unchecked(); + + *this = cpp_double_fp_backend::my_value_inf(); + + if (b_neg) + { + negate(); + } + + return *this; + } + + const bool isinf_v { (fpc_v == FP_INFINITE) }; + const bool isinf_u { (fpc_u == FP_INFINITE) }; + + if (isinf_u) + { + if (isinf_v) + { + return operator=(cpp_double_fp_backend::my_value_nan()); + } + else + { + const bool b_neg { isneg_unchecked() }; + + return operator=((!b_neg) ? cpp_double_fp_backend::my_value_inf() : -cpp_double_fp_backend::my_value_inf()); + } + } + + if (isinf_v) + { + return operator=(cpp_double_fp_backend(0)); + } + } + + if (this == &v) + { + data.first = float_type { 1.0F }; + data.second = float_type { 0.0F }; + + return *this; + } + + // The division algorithm has been taken from Victor Shoup, + // package WinNTL-5_3_2. It might originally be related to the + // K. Briggs work. The algorithm has been significantly simplified + // while still atempting to retain proper rounding corrections. + + const float_type C { data.first / v.data.first }; + + float_type c { cpp_df_qf_detail::split_maker::value * C }; + + float_type hc { }; + + if (cpp_df_qf_detail::ccmath::isinf(c)) + { + // Handle overflow by scaling down (and then back up) with the split. + + hc = + cpp_df_qf_detail::ccmath::ldexp + ( + C - float_type { C - cpp_df_qf_detail::ccmath::ldexp(C, -cpp_df_qf_detail::split_maker::n_shl) }, + cpp_df_qf_detail::split_maker::n_shl + ); + } + else + { + hc = c - float_type { c - C }; + } + + float_type u { cpp_df_qf_detail::split_maker::value * v.data.first }; + + float_type hv { }; + + if (cpp_df_qf_detail::ccmath::isinf(u)) + { + // Handle overflow by scaling down (and then back up) with the split. + + hv = + cpp_df_qf_detail::ccmath::ldexp + ( + v.data.first - float_type { v.data.first - cpp_df_qf_detail::ccmath::ldexp(v.data.first, -cpp_df_qf_detail::split_maker::n_shl) }, + cpp_df_qf_detail::split_maker::n_shl + ); + } + else + { + hv = u - float_type { u - v.data.first }; + } + + const float_type U { C * v.data.first }; + + u = cpp_df_qf_detail::ccmath::unsafe::fma(hc, hv, -U); + + const float_type tv { v.data.first - hv }; + + u = cpp_df_qf_detail::ccmath::unsafe::fma(hc, tv, u); + + { + const float_type tc { C - hc }; + + u = cpp_df_qf_detail::ccmath::unsafe::fma(tc, hv, u); + u = cpp_df_qf_detail::ccmath::unsafe::fma(tc, tv, u); + } + + c = float_type { (data.first - U) - u } + data.second; + + c = (c - float_type { C * v.data.second }) / v.data.first; + + // Perform even more simplifications compared to Victor Shoup. + data.first = C + c; + data.second = float_type { C - data.first } + c; + + return *this; + } + + // Unary minus operator. + constexpr auto operator-() const -> cpp_double_fp_backend + { + cpp_double_fp_backend v { *this }; + + v.negate(); + + return v; + } + + // Public helper functions. + static constexpr auto pown(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, int p) -> void + { + using local_float_type = cpp_double_fp_backend; + + if (p == 2) + { + result = x; result.mul_unchecked(x); + } + else if (p == 3) + { + result = x; result.mul_unchecked(x); result.mul_unchecked(x); + } + else if (p == 4) + { + local_float_type x2 { x }; + x2.mul_unchecked(x); + + result = x2; + result.mul_unchecked(x2); + } + else if (p == 1) + { + result = x; + } + else if (p < 0) + { + // The case p == 0 is checked in a higher calling layer. + + pown(result, local_float_type(1U) / x, -p); + } + else + { + result = local_float_type(1U); + + local_float_type y(x); + + auto p_local = static_cast(p); + + for (;;) + { + if (static_cast(p_local & static_cast(UINT8_C(1))) != static_cast(UINT8_C(0))) + { + result.mul_unchecked(y); + } + + p_local >>= 1U; + + if (p_local == static_cast(UINT8_C(0))) + { + break; + } + else + { + y.mul_unchecked(y); + } + } + } + } + + constexpr auto swap(cpp_double_fp_backend& other) -> void + { + if (this != &other) + { + const rep_type tmp { data }; + + data = other.data; + + other.data = tmp; + } + } + + constexpr auto swap(cpp_double_fp_backend&& other) noexcept -> void + { + const rep_type tmp { static_cast(data) }; + + data = other.data; + + other.data = tmp; + } + + constexpr auto compare(const cpp_double_fp_backend& other) const noexcept -> int + { + // Return 1 for *this > other, -1 for *this < other, 0 for *this = other. + + if ((cpp_df_qf_detail::ccmath::isnan)(data.first)) + { + return -1; + } + else + { + return (my_first() > other.my_first()) ? 1 : (my_first() < other.my_first()) + ? -1 : (my_second() > other.my_second()) + ? 1 : (my_second() < other.my_second()) + ? -1 : 0; + } + } + + auto str(std::streamsize number_of_digits, const std::ios::fmtflags format_flags) const -> std::string + { + if (number_of_digits == 0) + { + number_of_digits = cpp_double_fp_backend::my_digits10; + } + + // Use cpp_bin_float when writing to string. This is similar + // to the use of cpp_bin_float when reading from string. + + cpp_bin_float_read_write_type f_bin { data.first }; + + f_bin += data.second; + + return f_bin.str(number_of_digits, format_flags); + } + + constexpr auto order02() const -> int + { + // TBD: Is there another option to get the base-2 log + // that's more unequivocally closer to constexpr? + + // TBD: Either way, integrate this (or something like it) + // into any potential implementation of eval_ilogb(). + + int e2 { }; + + cpp_double_fp_backend dummy { }; + + eval_frexp(dummy, *this, &e2); + + return e2; + } + + static constexpr auto my_value_max() noexcept -> cpp_double_fp_backend + { + // Use the non-normalized sum of two maximum values, where the lower + // value is "shifted" right in the sense of floating-point ldexp. + + constexpr float_type + hi_part + { + (cpp_df_qf_detail::ccmath::numeric_limits::max)() + * ( + static_cast(1.0F) + - static_cast(1.5F) * cpp_df_qf_detail::ccmath::unsafe::sqrt(cpp_df_qf_detail::ccmath::numeric_limits::epsilon()) + ) + }; + + constexpr float_type + lo_part + { + cpp_df_qf_detail::ccmath::unsafe::ldexp + ( + (cpp_df_qf_detail::ccmath::numeric_limits::max)(), + -cpp_df_qf_detail::ccmath::numeric_limits::digits + ) + }; + + constexpr cpp_double_fp_backend my_value_max_constexpr { arithmetic::two_hilo_sum(hi_part, lo_part) }; + + static_assert + ( + eval_gt(my_value_max_constexpr, cpp_double_fp_backend { hi_part }), + "Error: maximum value is too small in relation to the maximum of its constituent type" + ); + + return my_value_max_constexpr; + } + + static constexpr auto my_value_min() noexcept -> cpp_double_fp_backend + { + // Use the non-normalized minimum value, where the lower value + // is "shifted" left in the sense of floating-point ldexp. + + constexpr cpp_double_fp_backend + my_value_min_constexpr + { + cpp_df_qf_detail::ccmath::unsafe::ldexp + ( + (cpp_df_qf_detail::ccmath::numeric_limits::min)(), + cpp_df_qf_detail::ccmath::numeric_limits::digits + ) + }; + + static_assert + ( + eval_gt(my_value_min_constexpr, cpp_double_fp_backend { (cpp_df_qf_detail::ccmath::numeric_limits::min)() }), + "Error: minimum value is too small and must exceed the min of its constituent type" + ); + + return my_value_min_constexpr; + } + + static constexpr auto my_value_eps() noexcept -> cpp_double_fp_backend + { + constexpr cpp_double_fp_backend + my_value_eps_constexpr + { + cpp_df_qf_detail::ccmath::unsafe::ldexp(float_type { 1 }, int { 3 - my_digits }) + }; + + static_assert + ( + eval_lt(cpp_double_fp_backend { 1 } - my_value_eps_constexpr, cpp_double_fp_backend { 1 }), + "Error: epsilon value is too small and must be large enough to differentiate (1 - epsilon) from 1" + ); + + return my_value_eps_constexpr; + } + + static constexpr auto my_value_nan() noexcept -> cpp_double_fp_backend + { + return cpp_double_fp_backend(static_cast(NAN), float_type { 0.0F }); + } + + static constexpr auto my_value_inf() noexcept -> cpp_double_fp_backend + { + return cpp_double_fp_backend(static_cast(HUGE_VAL), float_type { 0.0F }); // conversion from double infinity OK + } + + private: + rep_type data; + + using cpp_bin_float_read_write_backend_type = boost::multiprecision::backends::cpp_bin_float(my_digits), digit_base_2, void, int, cpp_df_qf_detail::ccmath::numeric_limits::min_exponent, cpp_df_qf_detail::ccmath::numeric_limits::max_exponent>; + using cpp_bin_float_read_write_exp_type = typename cpp_bin_float_read_write_backend_type::exponent_type; + + using cpp_bin_float_read_write_type = boost::multiprecision::number; + + auto rd_string(const char* pstr) -> bool; + + constexpr auto mul_unchecked(const cpp_double_fp_backend& v) -> void + { + // The multiplication algorithm has been taken from Victor Shoup, + // package WinNTL-5_3_2. It might originally be related to the + // K. Briggs work. The algorithm has been significantly simplified + // while still atempting to retain proper rounding corrections. + + float_type C { cpp_df_qf_detail::split_maker::value * data.first }; + + float_type hu { }; + + if (cpp_df_qf_detail::ccmath::isinf(C)) + { + // Handle overflow by scaling down (and then back up) with the split. + + C = data.first - cpp_df_qf_detail::ccmath::ldexp(data.first, -cpp_df_qf_detail::split_maker::n_shl); + + hu = cpp_df_qf_detail::ccmath::ldexp(data.first - C, cpp_df_qf_detail::split_maker::n_shl); + } + else + { + hu = C - float_type { C - data.first }; + } + + C = data.first * v.data.first; + + if (cpp_df_qf_detail::ccmath::isinf(C)) + { + // Handle overflow. + const bool b_neg { (isneg_unchecked() != v.isneg_unchecked()) }; + + *this = cpp_double_fp_backend::my_value_inf(); + + if (b_neg) + { + negate(); + } + + return; + } + + float_type c { cpp_df_qf_detail::split_maker::value * v.data.first }; + + float_type hv { }; + + if (cpp_df_qf_detail::ccmath::isinf(c)) + { + // Handle overflow by scaling down (and then back up) with the split. + + c = v.data.first - cpp_df_qf_detail::ccmath::ldexp(v.data.first, -cpp_df_qf_detail::split_maker::n_shl); + + hv = cpp_df_qf_detail::ccmath::ldexp(v.data.first - c, cpp_df_qf_detail::split_maker::n_shl); + } + else + { + hv = c - float_type { c - v.data.first }; + } + + const float_type tv { v.data.first - hv }; + + float_type t1 { cpp_df_qf_detail::ccmath::unsafe::fma(hu, hv, -C) }; + + t1 = cpp_df_qf_detail::ccmath::unsafe::fma(hu, tv, t1); + + const float_type tu { data.first - hu }; + + t1 = cpp_df_qf_detail::ccmath::unsafe::fma(tu, hv, t1); + + c = cpp_df_qf_detail::ccmath::unsafe::fma(tu, tv, t1) + + (data.first * v.data.second) + + (data.second * v.data.first); + + // Perform even more simplifications compared to Victor Shoup. + data.first = C + c; + data.second = float_type { C - data.first } + c; + } + + template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) < 16))>::type const*> + friend constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; + + template ::value && (((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) <= 36)))>::type const*> + friend constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; + + template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) > 36))>::type const*> + friend constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void; +}; + +template +auto cpp_double_fp_backend::rd_string(const char* pstr) -> bool +{ + using local_double_fp_type = cpp_double_fp_backend; + + cpp_bin_float_read_write_type f_bin { pstr }; + + const auto fpc = fpclassify(f_bin); + + const auto is_definitely_nan = (fpc == FP_NAN); + + if (is_definitely_nan) + { + static_cast(operator=(local_double_fp_type::my_value_nan())); + } + else + { + const bool b_neg { (signbit(f_bin) == 1) }; + + if (b_neg) { f_bin = -f_bin; } + + const int + expval_from_f_bin + { + [&f_bin]() + { + int expval; + + frexp(f_bin, &expval); + + return expval; + }() + }; + + const auto is_zero_or_subnormal = + ( + (fpc == FP_ZERO) + || (expval_from_f_bin < static_cast(local_double_fp_type::my_min_exponent)) + ); + + if (is_zero_or_subnormal) + { + data.first = float_type { 0.0F }; + data.second = float_type { 0.0F }; + } + else + { + const auto is_definitely_inf = (f_bin > (::std::numeric_limits::max)()); + + if (is_definitely_inf) + { + static_cast(operator=(local_double_fp_type::my_value_inf())); + + if (b_neg) { negate(); } + } + else + { + data.first = float_type { 0.0F }; + data.second = float_type { 0.0F }; + + using local_builtin_float_type = typename std::conditional<(sizeof(float_type) <= sizeof(double)), double, float_type>::type; + + constexpr int pow2_scaling_for_small_input { cpp_df_qf_detail::ccmath::numeric_limits::digits }; + + const auto has_pow2_scaling_for_small_input = + ( + expval_from_f_bin < static_cast(local_double_fp_type::my_min_exponent + pow2_scaling_for_small_input) + ); + + if (has_pow2_scaling_for_small_input) + { + f_bin = ldexp(f_bin, pow2_scaling_for_small_input); + } + + constexpr auto dig_lim = + static_cast + ( + static_cast + ( + (local_double_fp_type::my_digits / cpp_df_qf_detail::ccmath::numeric_limits::digits) + + (((local_double_fp_type::my_digits % cpp_df_qf_detail::ccmath::numeric_limits::digits) != 0) ? 1 : 0) + ) + * cpp_df_qf_detail::ccmath::numeric_limits::digits + ); + + for(auto i = static_cast(UINT8_C(0)); + i < dig_lim; + i = static_cast(i + static_cast(cpp_df_qf_detail::ccmath::numeric_limits::digits))) + { + const local_builtin_float_type f { static_cast(f_bin) }; + + f_bin -= cpp_bin_float_read_write_type { f }; + + eval_add(*this, local_double_fp_type { f }); + } + + if (has_pow2_scaling_for_small_input) + { + eval_ldexp(*this, *this, -pow2_scaling_for_small_input); + } + + if (b_neg) { negate(); } + } + } + } + + return true; +} + +} } } // namespace boost::multiprecision::backends + +#include + +namespace boost { namespace multiprecision { namespace backends { + +template +constexpr int cpp_double_fp_backend::my_digits; +template +constexpr int cpp_double_fp_backend::my_digits10; +template +constexpr int cpp_double_fp_backend::my_max_digits10; +template +constexpr int cpp_double_fp_backend::my_max_exponent; +template +constexpr int cpp_double_fp_backend::my_min_exponent; +template +constexpr int cpp_double_fp_backend::my_max_exponent10; +template +constexpr int cpp_double_fp_backend::my_min_exponent10; + +template +constexpr cpp_double_fp_backend operator+(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) { return cpp_double_fp_backend(a) += b; } +template +constexpr cpp_double_fp_backend operator-(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) { return cpp_double_fp_backend(a) -= b; } +template +constexpr cpp_double_fp_backend operator*(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) { return cpp_double_fp_backend(a) *= b; } +template +constexpr cpp_double_fp_backend operator/(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) { return cpp_double_fp_backend(a) /= b; } + +template +constexpr auto eval_add(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { result += x; } +template +constexpr auto eval_add(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void { result = cpp_double_fp_backend(a) += b; } +template +constexpr auto eval_subtract(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { result -= x; } +template +constexpr auto eval_subtract(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void { result = cpp_double_fp_backend(a) -= b; } +template +constexpr auto eval_multiply(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { result *= x; } +template +constexpr auto eval_multiply(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void { result = cpp_double_fp_backend(a) *= b; } +template +constexpr auto eval_divide(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void { result /= x; } +template +constexpr auto eval_divide(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> void { result = cpp_double_fp_backend(a) /= b; } +template +constexpr auto eval_eq(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool { return (a.compare(b) == 0); } +template +constexpr auto eval_lt(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool { return (a.compare(b) == -1); } +template +constexpr auto eval_gt(const cpp_double_fp_backend& a, const cpp_double_fp_backend& b) -> bool { return (a.compare(b) == 1); } +template + +constexpr auto eval_is_zero(const cpp_double_fp_backend& x) -> bool +{ + return x.iszero_unchecked(); +} + +template +constexpr auto eval_signbit(const cpp_double_fp_backend& x) -> int +{ + return (x.isneg_unchecked() ? 1 : 0); +} + +template +constexpr auto eval_fabs(cpp_double_fp_backend& result, const cpp_double_fp_backend& a) -> void +{ + result = a; + + if (a.isneg_unchecked()) + { + result.negate(); + } +} + +template +constexpr auto eval_frexp(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, int* v) -> void +{ + using local_backend_type = cpp_double_fp_backend; + using local_float_type = typename local_backend_type::float_type; + + int expptr { }; + + const local_float_type fhi { cpp_df_qf_detail::ccmath::frexp(a.rep().first, &expptr) }; + const local_float_type flo { cpp_df_qf_detail::ccmath::ldexp(a.rep().second, -expptr) }; + + if (v != nullptr) + { + *v = expptr; + } + + result.rep() = local_backend_type::arithmetic::normalize(fhi, flo); +} + +template +constexpr auto eval_ldexp(cpp_double_fp_backend& result, const cpp_double_fp_backend& a, int v) -> void +{ + using local_backend_type = cpp_double_fp_backend; + using local_float_type = typename local_backend_type::float_type; + + const local_float_type fhi { cpp_df_qf_detail::ccmath::ldexp(a.crep().first, v) }; + const local_float_type flo { cpp_df_qf_detail::ccmath::ldexp(a.crep().second, v) }; + + result.rep() = local_backend_type::arithmetic::normalize(fhi, flo); +} + +template +constexpr auto eval_floor(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void +{ + using local_backend_type = cpp_double_fp_backend; + using local_float_type = typename local_backend_type::float_type; + + const local_float_type fhi { cpp_df_qf_detail::ccmath::floor(x.my_first()) }; + + if (fhi != x.my_first()) + { + result.rep().first = fhi; + result.rep().second = local_float_type { 0 }; + } + else + { + const local_float_type flo = { cpp_df_qf_detail::ccmath::floor(x.my_second()) }; + + result.rep() = local_backend_type::arithmetic::normalize(fhi, flo); + } +} + +template +constexpr auto eval_ceil(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void +{ + // Compute -floor(-x); + eval_floor(result, -x); + + result.negate(); +} + +template +constexpr auto eval_fpclassify(const cpp_double_fp_backend& o) -> int +{ + if (cpp_df_qf_detail::ccmath::isnan(o.crep().first)) { return FP_NAN; } + else if (cpp_df_qf_detail::ccmath::isinf(o.crep().first)) { return FP_INFINITE; } + else if (eval_is_zero(o)) { return FP_ZERO; } + else + { + using local_backend_type = cpp_double_fp_backend; + using local_float_type = typename local_backend_type::float_type; + + const local_float_type fabs_x { cpp_df_qf_detail::ccmath::fabs(o.crep().first) }; + + if ((fabs_x > 0) && (fabs_x < (cpp_df_qf_detail::ccmath::numeric_limits::min)())) + { + return FP_SUBNORMAL; + } + else + { + return FP_NORMAL; + } + } +} + +template +constexpr auto eval_sqrt(cpp_double_fp_backend& result, const cpp_double_fp_backend& o) -> void +{ + using double_float_type = cpp_double_fp_backend; + using local_float_type = typename double_float_type::float_type; + + const int fpc { eval_fpclassify(o) }; + + const bool isneg_o { o.isneg_unchecked() }; + + if ((fpc != FP_NORMAL) || isneg_o) + { + if (fpc == FP_ZERO) + { + result = double_float_type(0); + + return; + } + else if ((fpc == FP_NAN) || isneg_o) + { + result = double_float_type::my_value_nan(); + + return; + } + else if (fpc == FP_INFINITE) + { + result = double_float_type::my_value_inf(); + + return; + } + } + + // TBD: Do we need any overflow/underflow guards when multiplying + // by the split or when multiplying (hx * tx) and/or (hx * hx)? + + const local_float_type c { cpp_df_qf_detail::ccmath::sqrt(o.crep().first) }; + + local_float_type p { cpp_df_qf_detail::split_maker::value * c }; + + const local_float_type hx { local_float_type { c - p } + p }; + const local_float_type tx { c - hx }; + + local_float_type q = hx * tx; + + q = q + q; + + p = hx * hx; + + const local_float_type u { p + q }; + + const local_float_type uu { cpp_df_qf_detail::ccmath::unsafe::fma(tx, tx, local_float_type { p - u } + q) }; + + const local_float_type + cc + { + local_float_type { local_float_type { o.crep().first - u } - uu + o.crep().second } / local_float_type { c + c } + }; + + result.rep().first = c + cc; + result.rep().second = local_float_type { c - result.rep().first } + cc; +} + +template +constexpr auto eval_pow(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, const cpp_double_fp_backend& a) -> void +{ + using double_float_type = cpp_double_fp_backend; + + constexpr double_float_type zero { 0 }; + + result = zero; + + signed long long val_nll { }; + + eval_convert_to(&val_nll, a); + + const int na { static_cast(val_nll) }; + + const int fpc_a { eval_fpclassify(a) }; + + if((fpc_a == FP_NORMAL) && (a.compare(double_float_type(na)) == 0)) + { + eval_pow(result, x, na); + } + else + { + constexpr double_float_type one { 1 }; + + const int fpc_x { eval_fpclassify(x) }; + const int fpc_a { eval_fpclassify(a) }; + + if (fpc_a == FP_ZERO) + { + // pow(base, +/-0) returns 1 for any base, even when base is NaN. + + result = one; + } + else if (fpc_a == FP_NAN) + { + result = double_float_type::my_value_nan(); + } + else if (fpc_x == FP_ZERO) + { + if ((fpc_a == FP_NORMAL) || (fpc_a == FP_INFINITE)) + { + // pow(+/-0, exp), where exp is negative and finite, returns +infinity. + // pow(+/-0, exp), where exp is positive non-integer, returns +0. + + // pow(+/-0, -infinity) returns +infinity. + // pow(+/-0, +infinity) returns +0. + + result = (eval_signbit(a) ? double_float_type::my_value_inf() : zero); + } + } + else if (fpc_x == FP_INFINITE) + { + if ((fpc_a == FP_NORMAL) || (fpc_a == FP_INFINITE)) + { + // pow(+infinity, exp) returns +0 for any negative exp. + // pow(-infinity, exp) returns +infinity for any positive exp. + + result = (eval_signbit(a) ? zero : double_float_type::my_value_inf()); + } + } + else if (fpc_x != FP_NORMAL) + { + result = x; + } + else + { + if (fpc_a == FP_INFINITE) + { + result = + ( + (x.compare(one) == -1) ? (eval_signbit(a) ? double_float_type::my_value_inf() : zero) + : (x.compare(one) == +1) ? (eval_signbit(a) ? zero : double_float_type::my_value_inf()) + : one + ); + } + else + { + double_float_type log_x { }; + + eval_log(log_x, x); + + double_float_type a_log_x { }; + + eval_multiply(a_log_x, a, log_x); + + eval_exp(result, a_log_x); + } + } + } +} + +template +constexpr auto eval_pow(cpp_double_fp_backend& result, const cpp_double_fp_backend& x, IntegralType p) -> typename ::std::enable_if::value, void>::type +{ + const int fpc { eval_fpclassify(x) }; + + using double_float_type = cpp_double_fp_backend; + + using local_integral_type = IntegralType; + + const bool p_is_odd { (static_cast(p & 1) != static_cast(0)) }; + + if (p == static_cast(0)) + { + // pow(base, +/-0) returns 1 for any base, even when base is NaN. + + result = double_float_type { unsigned { UINT8_C(1) } }; + } + else if (fpc != FP_NORMAL) + { + if (fpc == FP_ZERO) + { + // pow( +0, exp), where exp is a negative odd integer, returns +infinity. + // pow( -0, exp), where exp is a negative odd integer, returns +infinity. + // pow(+/-0, exp), where exp is a negative even integer, returns +infinity. + + // pow( +0, exp), where exp is a positive odd integer, returns +0. + // pow( -0, exp), where exp is a positive odd integer, returns -0. + // pow(+/-0, exp), where exp is a positive even integer, returns +0. + + result = ((p < static_cast(0)) ? double_float_type::my_value_inf() : double_float_type { unsigned { UINT8_C(0) } }); + } + else if (fpc == FP_INFINITE) + { + if (eval_signbit(x)) + { + if (p < static_cast(0)) + { + // pow(-infinity, exp) returns -0 if exp is a negative odd integer. + // pow(-infinity, exp) returns +0 if exp is a negative even integer. + + result = double_float_type { unsigned { UINT8_C(0) } }; + } + else + { + // pow(-infinity, exp) returns -infinity if exp is a positive odd integer. + // pow(-infinity, exp) returns +infinity if exp is a positive even integer. + + result = (p_is_odd ? -double_float_type::my_value_inf() : double_float_type::my_value_inf()); + } + } + else + { + // pow(+infinity, exp) returns +0 for any negative exp. + // pow(+infinity, exp) returns +infinity for any positive exp. + + result = ((p < static_cast(0)) ? double_float_type { unsigned { UINT8_C(0) } } : double_float_type::my_value_inf()); + } + } + else + { + result = double_float_type::my_value_nan(); + } + } + else + { + double_float_type::pown(result, x, p); + } +} + +template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) < 16))>::type const*> +constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void +{ + const int fpc { eval_fpclassify(x) }; + + using double_float_type = cpp_double_fp_backend; + + if (fpc == FP_ZERO) + { + result = double_float_type { 1.0F }; + } + else if (fpc != FP_NORMAL) + { + if (fpc == FP_INFINITE) + { + result = (x.isneg_unchecked() ? double_float_type { 0.0F } : double_float_type::my_value_inf()); + } + else if (fpc == FP_NAN) + { + result = x; + } + } + else + { + using local_float_type = typename double_float_type::float_type; + + // Get a local copy of the argument and force it to be positive. + const bool b_neg { x.isneg_unchecked() }; + + const double_float_type xx { (!b_neg) ? x : -x }; + + // Check the range of the input. + const double_float_type max_exp_input + { + []() -> local_float_type + { + local_float_type mx { }; + eval_convert_to(&mx, double_float_type::my_value_max()); + + const local_float_type log_of_mx = cpp_df_qf_detail::ccmath::log(mx); + + return log_of_mx; + }() + }; + + const double_float_type min_exp_input + { + []() -> local_float_type + { + local_float_type mn { }; + eval_convert_to(&mn, double_float_type::my_value_min()); + + const local_float_type log_of_mn = cpp_df_qf_detail::ccmath::log(mn); + + return log_of_mn; + }() + }; + + if (eval_lt(x, min_exp_input)) + { + result = double_float_type(0U); + } + else if (eval_gt(xx, max_exp_input)) + { + result = double_float_type::my_value_inf(); + } + else if (xx.is_one()) + { + result = + ((!b_neg) + ? cpp_df_qf_detail::constant_df_exp1() + : double_float_type(1U) / cpp_df_qf_detail::constant_df_exp1()); + } + else + { + // Use an argument reduction algorithm for exp() in classic MPFUN-style. + + double_float_type nf { }; + + eval_floor(nf, xx / cpp_df_qf_detail::constant_df_ln_two()); + + // Prepare the scaled variables. + const auto b_scale = (xx.order02() > -1); + + double_float_type r { }; + + if (b_scale) + { + eval_ldexp(r, xx - (nf * cpp_df_qf_detail::constant_df_ln_two()), -2); + } + else + { + r = xx; + } + + // PadeApproximant[Exp[x] - 1, {x, 0, {6, 6}}] + // FullSimplify[%] + // (84 x (7920 + 240 x^2 + x^4)) + // / (665280 + x (-332640 + x (75600 + x (-10080 + x (840 + (-42 + x) x))))) + + constexpr double_float_type n84(84); + constexpr double_float_type n240(240); + constexpr double_float_type n7920(7920); + + constexpr double_float_type n665280(665280); + constexpr double_float_type n332640(332640); + constexpr double_float_type n75600(75600); + constexpr double_float_type n10080(10080); + constexpr double_float_type n840(840); + constexpr double_float_type n42(42); + + const double_float_type r2 { r * r }; + + // Use the small-argument Pade approximation having coefficients shown above. + const double_float_type top = (n84 * r * (n7920 + (n240 + r2) * r2)); + const double_float_type bot = (n665280 + r * (-n332640 + r * (n75600 + r * (-n10080 + r * (n840 + (-n42 + r) * r))))); + + constexpr double_float_type one { 1 }; + + result = one + (top / bot); + + // Rescale the result. + if (b_scale) + { + result *= result; + result *= result; + + signed long long lln { }; + + eval_convert_to(&lln, nf); + + const int n { static_cast(lln) }; + + if (n > 0) + { + eval_ldexp(result, double_float_type(result), n); + } + } + + if (b_neg) + { + result = double_float_type(1U) / result; + } + } + } +} + +template ::value && (((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) <= 36)))>::type const*> +constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void +{ + const int fpc { eval_fpclassify(x) }; + + using double_float_type = cpp_double_fp_backend; + + if (fpc == FP_ZERO) + { + result = double_float_type(1); + } + else if (fpc != FP_NORMAL) + { + if (fpc == FP_INFINITE) + { + result = (x.isneg_unchecked() ? double_float_type(0) : double_float_type::my_value_inf()); + } + else if (fpc == FP_NAN) + { + result = x; + } + } + else + { + using local_float_type = typename double_float_type::float_type; + + // Get a local copy of the argument and force it to be positive. + const bool b_neg { x.isneg_unchecked() }; + + const double_float_type xx { (!b_neg) ? x : -x }; + + // Check the range of the input. + const double_float_type max_exp_input + { + []() -> local_float_type + { + local_float_type mx { }; + eval_convert_to(&mx, double_float_type::my_value_max()); + + const local_float_type log_of_mx = cpp_df_qf_detail::ccmath::log(mx); + + return log_of_mx; + }() + }; + + const double_float_type min_exp_input + { + []() -> local_float_type + { + local_float_type mn { }; + eval_convert_to(&mn, double_float_type::my_value_min()); + + const local_float_type log_of_mn = cpp_df_qf_detail::ccmath::log(mn); + + return log_of_mn; + }() + }; + + if (eval_lt(x, min_exp_input)) + { + result = double_float_type(0U); + } + else if (eval_gt(xx, max_exp_input)) + { + result = double_float_type::my_value_inf(); + } + else if (xx.is_one()) + { + result = + ((!b_neg) + ? cpp_df_qf_detail::constant_df_exp1() + : double_float_type(1U) / cpp_df_qf_detail::constant_df_exp1()); + } + else + { + // Use an argument reduction algorithm for exp() in classic MPFUN-style. + + double_float_type nf { }; + + eval_floor(nf, xx / cpp_df_qf_detail::constant_df_ln_two()); + + // Prepare the scaled variables. + const auto b_scale = (xx.order02() > -4); + + double_float_type r { }; + + if (b_scale) + { + eval_ldexp(r, xx - (nf * cpp_df_qf_detail::constant_df_ln_two()), -4); + } + else + { + r = xx; + } + + // PadeApproximant[Exp[r], {r, 0, 8, 8}] + // FullSimplify[%] + + constexpr double_float_type n144(144U); + constexpr double_float_type n3603600(3603600UL); + constexpr double_float_type n120120(120120UL); + constexpr double_float_type n770(770U); + + constexpr double_float_type n518918400(518918400UL); + constexpr double_float_type n259459200(259459200UL); + constexpr double_float_type n60540480(60540480UL); + constexpr double_float_type n8648640(8648640UL); + constexpr double_float_type n831600(831600UL); + constexpr double_float_type n55440(55440U); + constexpr double_float_type n2520(2520U); + constexpr double_float_type n72(72U); + + const double_float_type r2 { r * r }; + + const double_float_type top = (n144 * r) * (n3603600 + r2 * (n120120 + r2 * (n770 + r2))); + const double_float_type bot = (n518918400 + r * (-n259459200 + r * (n60540480 + r * (-n8648640 + r * (n831600 + r * (-n55440 + r * (n2520 + r * (-n72 + r)))))))); + + constexpr double_float_type one { 1 }; + + result = one + (top / bot); + + // Rescale the result. + if (b_scale) + { + result *= result; + result *= result; + result *= result; + result *= result; + + signed long long lln { }; + + eval_convert_to(&lln, nf); + + const int n { static_cast(lln) }; + + if (n > 0) + { + eval_ldexp(result, double_float_type(result), n); + } + } + + if (b_neg) + { + result = double_float_type(1U) / result; + } + } + } +} + +template ::value && ((cpp_df_qf_detail::ccmath::numeric_limits::digits10 * 2) > 36))>::type const*> +constexpr auto eval_exp(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void +{ + const int fpc { eval_fpclassify(x) }; + + using double_float_type = cpp_double_fp_backend; + + if (fpc == FP_ZERO) + { + result = double_float_type(1); + } + else if (fpc != FP_NORMAL) + { + if (fpc == FP_INFINITE) + { + result = (x.isneg_unchecked() ? double_float_type(0) : double_float_type::my_value_inf()); + } + else if (fpc == FP_NAN) + { + result = x; + } + } + else + { + using local_float_type = typename double_float_type::float_type; + + // Get a local copy of the argument and force it to be positive. + const bool b_neg { x.isneg_unchecked() }; + + const double_float_type xx { (!b_neg) ? x : -x }; + + // Check the range of the input. + const double_float_type max_exp_input + { + []() -> local_float_type + { + local_float_type mx { }; + eval_convert_to(&mx, double_float_type::my_value_max()); + + const local_float_type log_of_mx = cpp_df_qf_detail::ccmath::log(mx); + + return log_of_mx; + }() + }; + + const double_float_type min_exp_input + { + []() -> local_float_type + { + local_float_type mn { }; + eval_convert_to(&mn, double_float_type::my_value_min()); + + const local_float_type log_of_mn = cpp_df_qf_detail::ccmath::log(mn); + + return log_of_mn; + }() + }; + + if (eval_lt(x, min_exp_input)) + { + result = double_float_type(0U); + } + else if (eval_gt(xx, max_exp_input)) + { + result = double_float_type::my_value_inf(); + } + else if (xx.is_one()) + { + result = + ((!b_neg) + ? cpp_df_qf_detail::constant_df_exp1() + : double_float_type(1U) / cpp_df_qf_detail::constant_df_exp1()); + } + else + { + // Use an argument reduction algorithm for exp() in classic MPFUN-style. + + double_float_type nf { }; + + eval_floor(nf, xx / cpp_df_qf_detail::constant_df_ln_two()); + + // Prepare the scaled variables. + const auto b_scale = (xx.order02() > -4); + + double_float_type xh { }; + + if (b_scale) + { + eval_ldexp(xh, xx - (nf * cpp_df_qf_detail::constant_df_ln_two()), -4); + } + else + { + xh = xx; + } + + double_float_type x_pow_n_div_n_fact(xh); + + result = double_float_type(1U) + x_pow_n_div_n_fact; + + double_float_type dummy { }; + + // Series expansion of hypergeometric_0f0(; ; x). + // For this high(er) digit count, a scaled argument with subsequent + // Taylor series expansion is actually more precise than Pade approximation. + for (unsigned n = 2U; n < 64U; ++n) + { + x_pow_n_div_n_fact *= xh; + x_pow_n_div_n_fact /= typename double_float_type::float_type(n); + + int n_tol { }; + + eval_frexp(dummy, x_pow_n_div_n_fact, &n_tol); + + if ((n > 4U) && (n_tol < -(double_float_type::my_digits - 4))) + { + break; + } + + result += x_pow_n_div_n_fact; + } + + // Rescale the result. + if (b_scale) + { + result *= result; + result *= result; + result *= result; + result *= result; + + signed long long lln { }; + + eval_convert_to(&lln, nf); + + const int n { static_cast(lln) }; + + if (n > 0) + { + eval_ldexp(result, double_float_type(result), n); + } + } + + if (b_neg) + { + result = double_float_type(1U) / result; + } + } + } +} + +template +constexpr auto eval_log(cpp_double_fp_backend& result, const cpp_double_fp_backend& x) -> void +{ + using double_float_type = cpp_double_fp_backend; + + constexpr double_float_type one { 1 }; + + const int result_compare_with_one { x.compare(one) }; + + const int fpc { eval_fpclassify(x) }; + + if (fpc == FP_ZERO) + { + result = -double_float_type::my_value_inf(); + } + else if (x.isneg_unchecked() || (fpc == FP_NAN)) + { + result = double_float_type::my_value_nan(); + } + else if (fpc == FP_INFINITE) + { + result = double_float_type::my_value_inf(); + } + else if (result_compare_with_one == -1) + { + // Use argument inversion and negation of the result. + + double_float_type x_inv { }; + + eval_divide(x_inv, one, x); + eval_log(result, x_inv); + + result.negate(); + } + else if (result_compare_with_one == 1) + { + // Optimize to only use eval_frexp if (and only if) the exponent is + // actually large/small in the sense of above/below a defined cutoff. + + double_float_type x2 { }; + + int n2 { }; + + eval_frexp(x2, x, &n2); + + // Get initial estimate using the (wrapped) standard math function log. + const double_float_type s(cpp_df_qf_detail::ccmath::log(x2.my_first())); + + double_float_type E { }; + + eval_exp(E, s); + + // Do one single step of Newton-Raphson iteration. + result = s + ((x2 - E) / E); + + double_float_type xn2 { n2 }; + + using local_float_type = typename double_float_type::float_type; + + eval_multiply(xn2, cpp_df_qf_detail::constant_df_ln_two()); + + eval_add(result, xn2); + } + else + { + result = double_float_type { 0 }; + } +} + +namespace detail { + +template +constexpr auto extract(cpp_double_fp_backend& source) -> DestType +{ + using double_float_type = cpp_double_fp_backend; + using local_float_type = typename double_float_type::float_type; + + using destination_type = DestType; + + destination_type result { }; + + constexpr double_float_type zero { 0 }; + + result = static_cast(0); + + unsigned fail_safe { UINT32_C(32) }; + + constexpr bool + destination_type_is_longer + { + (std::numeric_limits::digits > cpp_df_qf_detail::ccmath::numeric_limits::digits) + }; + + using float_extract_type = typename std::conditional::type; + + while((source.compare(zero) != 0) && (fail_safe > unsigned { UINT8_C(0) })) + { + const float_extract_type next_flt_val { static_cast(source.my_first()) }; + + result += static_cast(next_flt_val); + + eval_subtract(source, double_float_type(next_flt_val)); + + --fail_safe; + } + + return result; +} + +} // namespace detail + +template +constexpr auto eval_convert_to(signed long long* result, const cpp_double_fp_backend& backend) -> void +{ + const auto fpc = eval_fpclassify(backend); + + if (fpc != FP_NORMAL) + { + *result = static_cast(backend.crep().first); + } + else + { + using double_float_type = cpp_double_fp_backend; + using local_float_type = typename double_float_type::float_type; + + static_assert(std::is_same::value, "Error something went wrong with the limb type"); + + constexpr bool + long_long_is_longer + { + (std::numeric_limits::digits > cpp_df_qf_detail::ccmath::numeric_limits::digits) + }; + + using longer_type = typename ::std::conditional::type; + + constexpr longer_type my_max_val { static_cast((std::numeric_limits::max)()) }; + constexpr longer_type my_min_val { static_cast((std::numeric_limits::min)()) }; + + constexpr double_float_type my_max_val_dd { static_cast(my_max_val) }; + constexpr double_float_type my_min_val_dd { static_cast(my_min_val) }; + + if (backend.compare(my_max_val_dd) >= 0) + { + *result = (std::numeric_limits::max)(); + } + else if (backend.compare(my_min_val_dd) <= 0) + { + *result = (std::numeric_limits::min)(); + } + else + { + double_float_type source { backend }; + + *result = detail::extract(source); + + #if !defined(__x86_64__) && !defined(_M_X64) + + // It has been "empirically found" that non-X64 needs this workaround. + // Even though the same conditions are met for x86_64 on GCC and MSVC, + // this workaround will actually break the long long conversion tests + // on those platforms. + // + // Our assumption is that on x64 there is x87 math (double -> long double) + // being performed in the background. Seemingly these might "aid" the conversion + // of double value to long long. Somehow I get the feeling this issue will arise + // in future evolution of the cpp_double_fp_backend. + // + // This workaround has been tested on: ARM64 (linux and mac), s390x and PPC64LE. + + constexpr bool + needs_workaround + { + (sizeof(signed long long) == 8U) + && (std::is_same::value || (std::is_same::value)) + }; + + BOOST_IF_CONSTEXPR (needs_workaround) + { + // This is the last value stored in a double as 9223372036854775808 + constexpr signed long long upper_bound = 9223372036854775296LL; + + if (!eval_signbit(backend) && *result >= upper_bound) + { + // LONG_MAX is stored with .second = -1, so we compensate for the offset + // We also only need this at the upper end where the values aren't exactly + // representable in double. Below a certain point we are fine + *result += static_cast(1); + } + } + + #endif // non-x64 + } + } +} + +template +constexpr auto eval_convert_to(unsigned long long* result, const cpp_double_fp_backend& backend) -> void +{ + const auto fpc = eval_fpclassify(backend); + + if (fpc != FP_NORMAL) + { + *result = static_cast(backend.crep().first); + } + else + { + using double_float_type = cpp_double_fp_backend; + using local_float_type = typename double_float_type::float_type; + + static_assert(std::is_same::value, "Error something went wrong with the limb type"); + + constexpr bool + ulong_long_is_longer + { + (std::numeric_limits::digits > cpp_df_qf_detail::ccmath::numeric_limits::digits) + }; + + using longer_type = typename ::std::conditional::type; + + constexpr longer_type my_max_val { static_cast((std::numeric_limits::max)()) }; + + constexpr double_float_type my_max_val_dd { static_cast(my_max_val) }; + + if (backend.compare(my_max_val_dd) >= 0) + { + *result = (std::numeric_limits::max)(); + } + else + { + double_float_type source { backend }; + + *result = detail::extract(source); + } + } +} + +#ifdef BOOST_HAS_INT128 +template +constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits::digits > 24), void>::type +{ + const auto fpc = eval_fpclassify(backend); + + if (fpc != FP_NORMAL) + { + *result = static_cast(backend.crep().first); + } + else + { + using double_float_type = cpp_double_fp_backend; + using local_float_type = typename double_float_type::float_type; + + static_assert(std::is_same::value, "Error something went wrong with the limb type"); + + constexpr bool + n128_is_longer + { + ((static_cast(sizeof(boost::int128_type) * static_cast(CHAR_BIT)) - 1) > cpp_df_qf_detail::ccmath::numeric_limits::digits) + }; + + using longer_type = typename ::std::conditional::type; + + constexpr boost::int128_type my_max_val_n128 = (((static_cast(1) << (sizeof(boost::int128_type) * CHAR_BIT - 2)) - 1) << 1) + 1; + constexpr boost::int128_type my_min_val_n128 = static_cast(-my_max_val_n128 - 1); + + constexpr longer_type my_max_val(static_cast(my_max_val_n128)); + constexpr longer_type my_min_val(static_cast(my_min_val_n128)); + + constexpr double_float_type my_max_val_dd(static_cast(my_max_val)); + constexpr double_float_type my_min_val_dd(static_cast(my_min_val)); + + if (backend.compare(my_max_val_dd) >= 0) + { + *result = my_max_val; + } + else if (backend.compare(my_min_val_dd) <= 0) + { + *result = my_min_val; + } + else + { + double_float_type source { backend }; + + *result = detail::extract(source); + } + } +} + +template +constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if::digits > 24), void>::type +{ + const auto fpc = eval_fpclassify(backend); + + if (fpc != FP_NORMAL) + { + *result = static_cast(backend.crep().first); + } + else + { + using double_float_type = cpp_double_fp_backend; + + double_float_type source { backend }; + + *result = detail::extract(source); + } +} + +template +constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits::digits > 24), void>::type +{ + const auto fpc = eval_fpclassify(backend); + + if (fpc != FP_NORMAL) + { + *result = static_cast(backend.crep().first); + } + else + { + using double_float_type = cpp_double_fp_backend; + using local_float_type = typename double_float_type::float_type; + + static_assert(std::is_same::value, "Error something went wrong with the limb type"); + + constexpr bool + u128_is_longer + { + (static_cast(sizeof(boost::uint128_type) * static_cast(CHAR_BIT)) > cpp_df_qf_detail::ccmath::numeric_limits::digits) + }; + + using longer_type = typename ::std::conditional::type; + + constexpr boost::int128_type my_max_val_u128 = static_cast(~static_cast(0)); + + constexpr longer_type my_max_val(static_cast(my_max_val_u128)); + + constexpr double_float_type my_max_val_dd(static_cast(my_max_val)); + + if (backend.compare(my_max_val_dd) >= 0) + { + *result = my_max_val; + } + else + { + double_float_type source { backend }; + + *result = detail::extract(source); + } + } +} + +template +constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend& backend) -> typename std::enable_if::digits > 24), void>::type +{ + const auto fpc = eval_fpclassify(backend); + + if (fpc != FP_NORMAL) + { + *result = static_cast(backend.crep().first); + } + else + { + using double_float_type = cpp_double_fp_backend; + + double_float_type source { backend }; + + *result = detail::extract(source); + } +} +#endif + +template +constexpr auto eval_convert_to(OtherFloatingPointType* result, const cpp_double_fp_backend& backend) -> typename ::std::enable_if::value>::type +{ + const auto fpc = eval_fpclassify(backend); + + // TBD: Implement min/max check for the destination floating-point type result. + + if (fpc != FP_NORMAL) + { + *result = static_cast(backend.my_first()); + } + else + { + BOOST_IF_CONSTEXPR(cpp_df_qf_detail::ccmath::numeric_limits::digits > cpp_df_qf_detail::ccmath::numeric_limits::digits) + { + *result = static_cast(backend.my_first()); + *result += static_cast(backend.my_second()); + } + else + { + cpp_double_fp_backend source = backend; + + *result = 0; + + for(auto digit_count = static_cast(0); + digit_count < cpp_double_fp_backend::my_digits; + digit_count += cpp_df_qf_detail::ccmath::numeric_limits::digits) + { + const auto next = static_cast(source.my_first()); + + *result += next; + + eval_subtract(source, cpp_double_fp_backend(next)); + } + } + } +} + +template +auto hash_value(const cpp_double_fp_backend& a) -> ::std::size_t +{ + return a.hash(); +} + +} // namespace backends + +using backends::cpp_double_fp_backend; + +using cpp_double_float = number, boost::multiprecision::et_off>; +using cpp_double_double = number, boost::multiprecision::et_off>; +using cpp_double_long_double = number, boost::multiprecision::et_off>; +#ifdef BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128 +using cpp_double_float128 = number, boost::multiprecision::et_off>; +#endif + +} } // namespace boost::multiprecision + +namespace std { + +// Specialization of numeric_limits for boost::multiprecision::number> +template +BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE numeric_limits, ExpressionTemplatesOption> > +{ + private: + using local_float_type = FloatingPointType; + using inner_self_type = boost::multiprecision::cpp_double_fp_backend; + + using self_type = boost::multiprecision::number; + + public: + static constexpr bool is_specialized = true; + static constexpr bool is_signed = true; + static constexpr bool is_integer = false; + static constexpr bool is_exact = false; + static constexpr bool is_bounded = true; + static constexpr bool is_modulo = false; + static constexpr bool is_iec559 = false; + static constexpr bool has_infinity = true; + static constexpr bool has_quiet_NaN = true; + static constexpr bool has_signaling_NaN = false; + + // These members were deprecated in C++23, but only MSVC warns (as of June 25) + #ifdef BOOST_MSVC + #pragma warning(push) + #pragma warning(disable:4996) + #endif + + static constexpr float_denorm_style has_denorm = denorm_absent; + static constexpr bool has_denorm_loss = true; + + #ifdef BOOST_MSVC + #pragma warning(pop) + #endif // deprecated members + + static constexpr bool traps = false; + static constexpr bool tinyness_before = false; + static constexpr float_round_style round_style = round_to_nearest; + + static constexpr int radix = 2; + static constexpr int digits = inner_self_type::my_digits; + static constexpr int digits10 = inner_self_type::my_digits10; + static constexpr int max_digits10 = inner_self_type::my_max_digits10; + + static constexpr int max_exponent = inner_self_type::my_max_exponent; + static constexpr int min_exponent = inner_self_type::my_min_exponent; + static constexpr int max_exponent10 = inner_self_type::my_max_exponent10; + static constexpr int min_exponent10 = inner_self_type::my_min_exponent10; + + static constexpr auto(min) () noexcept -> self_type { return static_cast(inner_self_type::my_value_min()); } + static constexpr auto(max) () noexcept -> self_type { return static_cast(inner_self_type::my_value_max()); } + static constexpr auto lowest () noexcept -> self_type { return static_cast(-(max)()); } + static constexpr auto epsilon () noexcept -> self_type { return static_cast(inner_self_type::my_value_eps()); } + static constexpr auto round_error () noexcept -> self_type { return static_cast(static_cast(0.5)); } + static constexpr auto denorm_min () noexcept -> self_type { return static_cast((min)()); } + static constexpr auto infinity () noexcept -> self_type { return static_cast(inner_self_type::my_value_inf()); } + static constexpr auto quiet_NaN () noexcept -> self_type { return static_cast(inner_self_type::my_value_nan()); } + static constexpr auto signaling_NaN() noexcept -> self_type { return static_cast(static_cast(0.0)); } +}; + +} // namespace std + +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::is_specialized; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::is_signed; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::is_integer; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::is_exact; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::is_bounded; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::is_modulo; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::is_iec559; + +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4996) +#endif + +template +constexpr std::float_denorm_style std::numeric_limits, ExpressionTemplatesOption> >::has_denorm; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::has_denorm_loss; + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif // deprecated members + +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::has_infinity; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::has_quiet_NaN; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::has_signaling_NaN; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::traps; +template +constexpr bool std::numeric_limits, ExpressionTemplatesOption> >::tinyness_before; +template +constexpr std::float_round_style std::numeric_limits, ExpressionTemplatesOption> >::round_style; +template +constexpr int std::numeric_limits, ExpressionTemplatesOption> >::radix; +template +constexpr int std::numeric_limits, ExpressionTemplatesOption> >::digits; +template +constexpr int std::numeric_limits, ExpressionTemplatesOption> >::digits10; +template +constexpr int std::numeric_limits, ExpressionTemplatesOption> >::max_digits10; +template +constexpr int std::numeric_limits, ExpressionTemplatesOption> >::max_exponent; +template +constexpr int std::numeric_limits, ExpressionTemplatesOption> >::min_exponent; +template +constexpr int std::numeric_limits, ExpressionTemplatesOption> >::max_exponent10; +template +constexpr int std::numeric_limits, ExpressionTemplatesOption> >::min_exponent10; + +#if defined(BOOST_MP_MATH_AVAILABLE) +namespace boost { namespace math { namespace policies { + +template +struct precision, ExpressionTemplates>, Policy> +{ +private: + using my_multiprecision_backend_type = boost::multiprecision::cpp_double_fp_backend; + + using digits_2 = digits2; + +public: + using precision_type = typename Policy::precision_type; + + using type = + typename std::conditional< + ((digits_2::value <= precision_type::value) || (precision_type::value <= 0)), + digits_2, // Default case: Full precision for RealType. + precision_type>::type; // User customized precision. +}; + +} } } // namespace boost::math::policies +#endif + +#endif // BOOST_MP_CPP_DOUBLE_FP_2021_06_05_HPP diff --git a/include/boost/multiprecision/cpp_int/bitwise.hpp b/include/boost/multiprecision/cpp_int/bitwise.hpp index 99f6358..da15391 100644 --- a/include/boost/multiprecision/cpp_int/bitwise.hpp +++ b/include/boost/multiprecision/cpp_int/bitwise.hpp @@ -563,7 +563,11 @@ inline BOOST_MP_CXX14_CONSTEXPR void right_shift_limb(Int& result, double_limb_t pr[i] = pr[i + offset]; result.resize(rs, rs); } - +// see discussion at https://github.com/boostorg/multiprecision/issues/686 for logic in disabling this. +#if defined(__GNUC__) && !defined(__clang__) && (__GNUC__ >= 12) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Warray-bounds" +#endif template inline BOOST_MP_CXX14_CONSTEXPR void right_shift_generic(Int& result, double_limb_type s) { @@ -598,6 +602,9 @@ inline BOOST_MP_CXX14_CONSTEXPR void right_shift_generic(Int& result, double_lim pr[i] = pr[i + offset] >> shift; result.resize(rs, rs); } +#if defined(__GNUC__) && !defined(__clang__) && (__GNUC__ >= 12) +#pragma GCC diagnostic pop +#endif template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if >::value>::type @@ -700,7 +707,7 @@ BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if& result, T s) noexcept((is_non_throwing_cpp_int >::value)) { is_valid_bitwise_op(result, typename cpp_int_backend::checked_type()); - *result.limbs() = detail::checked_left_shift(*result.limbs(), s, typename cpp_int_backend::checked_type()); + *result.limbs() = detail::checked_left_shift(*result.limbs(), static_cast(s), typename cpp_int_backend::checked_type()); result.normalize(); } @@ -710,7 +717,7 @@ eval_right_shift(cpp_int_backend::checked_type()); - *result.limbs() = (static_cast(s) >= sizeof(*result.limbs()) * CHAR_BIT) ? 0 : (result.sign() ? ((--*result.limbs()) >> s) + 1 : *result.limbs() >> s); + *result.limbs() = (static_cast(s) >= static_cast(sizeof(*result.limbs()) * CHAR_BIT)) ? 0 : (result.sign() ? ((--*result.limbs()) >> s) + 1 : *result.limbs() >> s); if (result.sign() && (*result.limbs() == 0)) result = static_cast(-1); } diff --git a/include/boost/multiprecision/cpp_int/literals.hpp b/include/boost/multiprecision/cpp_int/literals.hpp index aca09c2..0531f0a 100644 --- a/include/boost/multiprecision/cpp_int/literals.hpp +++ b/include/boost/multiprecision/cpp_int/literals.hpp @@ -249,29 +249,30 @@ constexpr typename boost::multiprecision::literals::detail::unsigned_cpp_int_lit return boost::multiprecision::literals::detail::make_backend_from_pack((sizeof...(STR)) - 2u)>::backend_type>::value; } -#define BOOST_MP_DEFINE_SIZED_CPP_INT_LITERAL(Bits) \ - template \ - constexpr boost::multiprecision::number > operator BOOST_JOIN(""_cppi, Bits)() \ - { \ - using pt = typename boost::multiprecision::literals::detail::make_packed_value_from_str::type; \ - return boost::multiprecision::literals::detail::make_backend_from_pack< \ - pt, \ - boost::multiprecision::backends::cpp_int_backend >::value; \ - } \ - template \ - constexpr boost::multiprecision::number > operator BOOST_JOIN(""_cppui, Bits)() \ - { \ - using pt = typename boost::multiprecision::literals::detail::make_packed_value_from_str::type; \ - return boost::multiprecision::literals::detail::make_backend_from_pack< \ - pt, \ - boost::multiprecision::backends::cpp_int_backend >::value; \ +#define BOOST_MP_LIT(P, N) BOOST_JOIN(operator "", BOOST_JOIN(P, N)) +#define BOOST_MP_DEFINE_SIZED_CPP_INT_LITERAL(Bits) \ + template \ + constexpr boost::multiprecision::number > BOOST_MP_LIT(_cppi, Bits)() \ + { \ + using pt = typename boost::multiprecision::literals::detail::make_packed_value_from_str::type; \ + return boost::multiprecision::literals::detail::make_backend_from_pack< \ + pt, \ + boost::multiprecision::backends::cpp_int_backend >::value; \ + } \ + template \ + constexpr boost::multiprecision::number > BOOST_MP_LIT(_cppui, Bits)() \ + { \ + using pt = typename boost::multiprecision::literals::detail::make_packed_value_from_str::type; \ + return boost::multiprecision::literals::detail::make_backend_from_pack< \ + pt, \ + boost::multiprecision::backends::cpp_int_backend >::value; \ } BOOST_MP_DEFINE_SIZED_CPP_INT_LITERAL(128) BOOST_MP_DEFINE_SIZED_CPP_INT_LITERAL(256) BOOST_MP_DEFINE_SIZED_CPP_INT_LITERAL(512) BOOST_MP_DEFINE_SIZED_CPP_INT_LITERAL(1024) - +#undef BOOST_MP_LIT } // namespace literals // diff --git a/include/boost/multiprecision/cpp_int/misc.hpp b/include/boost/multiprecision/cpp_int/misc.hpp index a8a7640..c464cfa 100644 --- a/include/boost/multiprecision/cpp_int/misc.hpp +++ b/include/boost/multiprecision/cpp_int/misc.hpp @@ -11,19 +11,15 @@ #ifndef BOOST_MP_CPP_INT_MISC_HPP #define BOOST_MP_CPP_INT_MISC_HPP -#include -#include #include -#include #include -#include #include // lsb etc +#include +#include #include #include -#include // std::gcd -#include -#include -#include +#include +#include #ifndef BOOST_MP_STANDALONE #include @@ -33,6 +29,15 @@ #include #endif +#if ((defined(__cpp_lib_gcd_lcm) && (__cpp_lib_gcd_lcm >= 201606L)) && (!defined(BOOST_HAS_INT128) || !defined(__STRICT_ANSI__))) +#define BOOST_MP_LIB_GCD_LCM_AVAILABLE +#endif + +#include +#include +#include +#include + #ifdef BOOST_MSVC #pragma warning(push) #pragma warning(disable : 4702) @@ -517,7 +522,7 @@ BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR limb_type eval_gcd(limb_type u, li // boundary cases if (!u || !v) return u | v; -#if (defined(__cpp_lib_gcd_lcm) && (__cpp_lib_gcd_lcm >= 201606L)) +#if defined(BOOST_MP_LIB_GCD_LCM_AVAILABLE) return std::gcd(u, v); #else std::size_t shift = boost::multiprecision::detail::find_lsb(u | v); @@ -535,7 +540,7 @@ BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR limb_type eval_gcd(limb_type u, li inline BOOST_MP_CXX14_CONSTEXPR double_limb_type eval_gcd(double_limb_type u, double_limb_type v) { -#if (defined(__cpp_lib_gcd_lcm) && (__cpp_lib_gcd_lcm >= 201606L)) && (!defined(BOOST_HAS_INT128) || !defined(__STRICT_ANSI__)) +#if defined(BOOST_MP_LIB_GCD_LCM_AVAILABLE) return std::gcd(u, v); #else if (u == 0) diff --git a/include/boost/multiprecision/cpp_int/serialize.hpp b/include/boost/multiprecision/cpp_int/serialize.hpp index e419619..2c20c8b 100644 --- a/include/boost/multiprecision/cpp_int/serialize.hpp +++ b/include/boost/multiprecision/cpp_int/serialize.hpp @@ -19,8 +19,6 @@ class binary_iarchive; namespace serialization { -namespace mp = boost::multiprecision; - namespace cpp_int_detail { using namespace boost::multiprecision; @@ -191,12 +189,12 @@ void do_serialize(Archive& ar, Int& val, std::integral_constant cons } // namespace cpp_int_detail -template -void serialize(Archive& ar, mp::cpp_int_backend& val, const unsigned int /*version*/) +template +void serialize(Archive& ar, multiprecision::cpp_int_backend& val, const unsigned int /*version*/) { using archive_save_tag = typename Archive::is_saving ; using save_tag = std::integral_constant ; - using trivial_tag = std::integral_constant >::value>; + using trivial_tag = std::integral_constant >::value>; using binary_tag = typename cpp_int_detail::is_binary_archive::type ; // Just dispatch to the correct method: diff --git a/include/boost/multiprecision/detail/default_ops.hpp b/include/boost/multiprecision/detail/default_ops.hpp index bd4715b..d66fe4a 100644 --- a/include/boost/multiprecision/detail/default_ops.hpp +++ b/include/boost/multiprecision/detail/default_ops.hpp @@ -2843,7 +2843,7 @@ template inline BOOST_MP_CXX14_CONSTEXPR typename detail::expression::result_type trunc(const detail::expression& v, const Policy& pol) { using number_type = typename detail::expression::result_type; - return std::move(trunc(number_type(v), pol)); + return trunc(number_type(v), pol); } template @@ -2943,7 +2943,7 @@ template inline BOOST_MP_CXX14_CONSTEXPR typename detail::expression::result_type round(const detail::expression& v, const Policy& pol) { using number_type = typename detail::expression::result_type; - return std::move(round(static_cast(v), pol)); + return round(static_cast(v), pol); } template inline BOOST_MP_CXX14_CONSTEXPR number round(const number& v, const Policy&) @@ -3057,7 +3057,7 @@ inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if& v, short* pint) { using number_type = typename detail::expression::result_type; - return std::move(frexp(static_cast(v), pint)); + return frexp(static_cast(v), pint); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_floating_point, number >::type frexp(const number& v, int* pint) @@ -3073,7 +3073,7 @@ inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if& v, int* pint) { using number_type = typename detail::expression::result_type; - return std::move(frexp(static_cast(v), pint)); + return frexp(static_cast(v), pint); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_floating_point, number >::type frexp(const number& v, long* pint) @@ -3089,7 +3089,7 @@ inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if& v, long* pint) { using number_type = typename detail::expression::result_type; - return std::move(frexp(static_cast(v), pint)); + return frexp(static_cast(v), pint); } template inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if::value == number_kind_floating_point, number >::type frexp(const number& v, long long* pint) @@ -3105,7 +3105,7 @@ inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if& v, long long* pint) { using number_type = typename detail::expression::result_type; - return std::move(frexp(static_cast(v), pint)); + return frexp(static_cast(v), pint); } // // modf does not return an expression template since we require the diff --git a/include/boost/multiprecision/detail/float_string_cvt.hpp b/include/boost/multiprecision/detail/float_string_cvt.hpp index 11a0bbd..71c1cbf 100644 --- a/include/boost/multiprecision/detail/float_string_cvt.hpp +++ b/include/boost/multiprecision/detail/float_string_cvt.hpp @@ -139,7 +139,7 @@ std::string convert_to_string(Backend b, std::streamsize digits, std::ios_base:: // // Extract the digits one at a time: // - for (unsigned i = 0; i < digits; ++i) + for (std::streamsize i = 0; i < digits; ++i) { eval_floor(digit, t); eval_convert_to(&cdigit, digit); diff --git a/include/boost/multiprecision/detail/functions/pow.hpp b/include/boost/multiprecision/detail/functions/pow.hpp index 34f4feb..481c7e4 100644 --- a/include/boost/multiprecision/detail/functions/pow.hpp +++ b/include/boost/multiprecision/detail/functions/pow.hpp @@ -1,6 +1,6 @@ -// Copyright Christopher Kormanyos 2002 - 2013. -// Copyright 2011 - 2013 John Maddock. +// Copyright 2011 - 2025 John Maddock. +// Copyright Christopher Kormanyos 2002 - 2025. // Distributed under the Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) @@ -132,7 +132,6 @@ void hyp0F0(T& H0F0, const T& x) BOOST_MP_ASSERT(&H0F0 != &x); long tol = boost::multiprecision::detail::digits2 >::value(); - T t; T x_pow_n_div_n_fact(x); @@ -192,7 +191,7 @@ void hyp1F0(T& H1F0, const T& a, const T& x) lim.negate(); si_type n; - T term, part; + T term; const si_type series_limit = boost::multiprecision::detail::digits2 >::value() < 100 diff --git a/include/boost/multiprecision/detail/precision.hpp b/include/boost/multiprecision/detail/precision.hpp index df58144..a33cc0b 100644 --- a/include/boost/multiprecision/detail/precision.hpp +++ b/include/boost/multiprecision/detail/precision.hpp @@ -136,7 +136,7 @@ inline constexpr unsigned current_precision_of(const expression::value> -struct scoped_default_precision +struct BOOST_ATTRIBUTE_UNUSED scoped_default_precision { template constexpr scoped_default_precision(const T&) {} diff --git a/include/boost/multiprecision/detail/standalone_config.hpp b/include/boost/multiprecision/detail/standalone_config.hpp index 95377fd..e83cbef 100644 --- a/include/boost/multiprecision/detail/standalone_config.hpp +++ b/include/boost/multiprecision/detail/standalone_config.hpp @@ -144,9 +144,13 @@ namespace boost { namespace multiprecision { #define BOOST_MP_NO_CONSTEXPR_DETECTION #endif -#ifdef __has_attribute -# if __has_attribute(fallthrough) -# define BOOST_MP_FALLTHROUGH [[fallthrough]] +// Compilers should ignore unknown attributes, +// but this doesn't stop them from issuing warnings +#if (defined(_MSVC_LANG) && _MSVC_LANG >= 201703L) || (defined(__cplusplus) && __cplusplus >= 201703L) +# ifdef __has_attribute +# if __has_attribute(fallthrough) +# define BOOST_MP_FALLTHROUGH [[fallthrough]] +# endif # endif #endif diff --git a/include/boost/multiprecision/float128.hpp b/include/boost/multiprecision/float128.hpp index 83f5f11..5dfbcdd 100644 --- a/include/boost/multiprecision/float128.hpp +++ b/include/boost/multiprecision/float128.hpp @@ -262,7 +262,7 @@ struct float128_backend int v_max = v; std::unique_ptr buf2; buf2.reset(new char[v + 3]); - v = quadmath_snprintf(&buf2[0], v_max + 3, format.c_str(), digits, m_value); + v = quadmath_snprintf(&buf2[0], static_cast(v_max + 3), format.c_str(), digits, m_value); if (v >= v_max + 3) { BOOST_MP_THROW_EXCEPTION(std::runtime_error("Formatting of float128_type failed.")); diff --git a/include/boost/multiprecision/mpfr.hpp b/include/boost/multiprecision/mpfr.hpp index e518a53..7dbcb3f 100644 --- a/include/boost/multiprecision/mpfr.hpp +++ b/include/boost/multiprecision/mpfr.hpp @@ -2967,6 +2967,8 @@ inline boost::multiprecision::number, ExpressionTemplates> > precision_guard(arg); boost::multiprecision::number, ExpressionTemplates> result; + if (arg < -1) + return policies::raise_domain_error("log1p<%1%>(%1%)", "Argument was %1% but must be greater than -1", arg, pol); mpfr_log1p(result.backend().data(), arg.backend().data(), GMP_RNDN); if (mpfr_inf_p(result.backend().data())) return (arg == -1 ? -1 : 1) * policies::raise_overflow_error, ExpressionTemplates> >("log1p<%1%>(%1%)", nullptr, pol); diff --git a/include/boost/version.hpp b/include/boost/version.hpp index 606b426..7b98764 100644 --- a/include/boost/version.hpp +++ b/include/boost/version.hpp @@ -19,7 +19,7 @@ // BOOST_VERSION / 100 % 1000 is the minor version // BOOST_VERSION / 100000 is the major version -#define BOOST_VERSION 108800 +#define BOOST_VERSION 108900 // // BOOST_LIB_VERSION must be defined to be the same as BOOST_VERSION @@ -27,6 +27,6 @@ // number, y is the minor version number, and z is the patch level if not 0. // This is used by to select which library version to link to. -#define BOOST_LIB_VERSION "1_88" +#define BOOST_LIB_VERSION "1_89" #endif