2634 lines
100 KiB
C++
2634 lines
100 KiB
C++
///////////////////////////////////////////////////////////////////////////////
|
|
// 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 <boost/multiprecision/cpp_bin_float.hpp>
|
|
#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail.hpp>
|
|
#include <boost/multiprecision/detail/hash.hpp>
|
|
#include <boost/multiprecision/traits/max_digits10.hpp>
|
|
#include <boost/multiprecision/traits/std_integer_traits.hpp>
|
|
|
|
#ifdef BOOST_MP_MATH_AVAILABLE
|
|
//
|
|
// Headers required for Boost.Math integration:
|
|
//
|
|
#include <boost/math/policies/policy.hpp>
|
|
//
|
|
// Some includes we need from Boost.Math, since we rely on that library to provide these functions:
|
|
//
|
|
#include <boost/math/special_functions/acosh.hpp>
|
|
#include <boost/math/special_functions/asinh.hpp>
|
|
#include <boost/math/special_functions/atanh.hpp>
|
|
#include <boost/math/special_functions/cbrt.hpp>
|
|
#include <boost/math/special_functions/expm1.hpp>
|
|
#include <boost/math/special_functions/gamma.hpp>
|
|
#endif
|
|
|
|
#include <limits>
|
|
#include <string>
|
|
#include <type_traits>
|
|
|
|
#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 <typename FloatingPointType>
|
|
class cpp_double_fp_backend;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_add(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_add(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_subtract(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_subtract(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_multiply(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_multiply(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_divide(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_divide(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_eq(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_lt(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_gt(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_is_zero(const cpp_double_fp_backend<FloatingPointType>& x) -> bool;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_signbit(const cpp_double_fp_backend<FloatingPointType>& x) -> int;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_fabs(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_frexp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, int* v) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_ldexp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, int v) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_floor(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_ceil(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_fpclassify(const cpp_double_fp_backend<FloatingPointType>& o) -> int;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_sqrt(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& o) -> void;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_pow(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x, const cpp_double_fp_backend<FloatingPointType>& a) -> void;
|
|
|
|
template <typename FloatingPointType,
|
|
typename IntegralType>
|
|
constexpr auto eval_pow(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x, IntegralType n) -> typename ::std::enable_if<boost::multiprecision::detail::is_integral<IntegralType>::value, void>::type;
|
|
|
|
template <typename FloatingPointType,
|
|
typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) < 16))>::type const* = nullptr>
|
|
constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
|
|
template <typename FloatingPointType,
|
|
typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && (((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) <= 36)))>::type const* = nullptr>
|
|
constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
|
|
template <typename FloatingPointType,
|
|
typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) > 36))>::type const* = nullptr>
|
|
constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_log(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(signed long long* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> void;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(unsigned long long* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> void;
|
|
|
|
#ifdef BOOST_HAS_INT128
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<!(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<!(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type;
|
|
#endif
|
|
|
|
template <typename FloatingPointType,
|
|
typename OtherFloatingPointType>
|
|
constexpr auto eval_convert_to(OtherFloatingPointType* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename ::std::enable_if<cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value>::type;
|
|
|
|
// TBD: constexpr on hash_value?
|
|
template <typename FloatingPointType>
|
|
auto hash_value(const cpp_double_fp_backend<FloatingPointType>& a) -> ::std::size_t;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto fabs(const cpp_double_fp_backend<FloatingPointType>& a) -> cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
} } } // namespace boost::multiprecision::backends
|
|
|
|
namespace std {
|
|
|
|
// Foward declarations of various specializations of std::numeric_limits
|
|
|
|
template <typename FloatingPointType,
|
|
const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >;
|
|
|
|
} // namespace std
|
|
|
|
namespace boost { namespace multiprecision {
|
|
|
|
template <typename FloatingPointType>
|
|
struct number_category<backends::cpp_double_fp_backend<FloatingPointType>> : public std::integral_constant<int, number_kind_floating_point> { };
|
|
|
|
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 <typename FloatingPointType>
|
|
class cpp_double_fp_backend
|
|
{
|
|
public:
|
|
using float_type = FloatingPointType;
|
|
|
|
static_assert
|
|
(
|
|
cpp_df_qf_detail::is_floating_point<float_type>::value
|
|
&& bool
|
|
{
|
|
(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits == 24)
|
|
|| (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits == 53)
|
|
|| (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits == 64)
|
|
|| (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits == 113)
|
|
}, "Error: float_type does not fulfil the backend requirements of cpp_double_fp_backend"
|
|
);
|
|
|
|
using rep_type = cpp_df_qf_detail::pair<float_type, float_type>;
|
|
using arithmetic = cpp_df_qf_detail::exact_arithmetic<float_type>;
|
|
|
|
using signed_types = std::tuple<signed char, signed short, signed int, signed long, signed long long, std::intmax_t>;
|
|
using unsigned_types = std::tuple<unsigned char, unsigned short, unsigned int, unsigned long, unsigned long long, std::uintmax_t>;
|
|
using float_types = std::tuple<float, double, long double>;
|
|
using exponent_type = int;
|
|
|
|
static constexpr int my_digits = static_cast<int>(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits * static_cast<int>(INT8_C(2)));
|
|
static constexpr int my_digits10 = static_cast<int>(boost::multiprecision::detail::calc_digits10 <static_cast<unsigned>(my_digits)>::value);
|
|
static constexpr int my_max_digits10 = static_cast<int>(boost::multiprecision::detail::calc_max_digits10<static_cast<unsigned>(my_digits)>::value);
|
|
static constexpr int my_max_exponent = cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max_exponent;
|
|
static constexpr int my_max_exponent10 = cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max_exponent10;
|
|
|
|
static constexpr int my_min_exponent =
|
|
static_cast<int>
|
|
(
|
|
cpp_df_qf_detail::ccmath::numeric_limits<float_type>::min_exponent
|
|
+ cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits
|
|
);
|
|
|
|
static constexpr int my_min_exponent10 =
|
|
static_cast<int>
|
|
(
|
|
-static_cast<int>
|
|
(
|
|
boost::multiprecision::detail::calc_digits10<static_cast<unsigned>(-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<rep_type&&>(other.data)) { }
|
|
|
|
// Constructors from other floating-point types.
|
|
template <typename OtherFloatType,
|
|
typename ::std::enable_if<( cpp_df_qf_detail::is_floating_point<OtherFloatType>::value
|
|
&& (cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatType>::digits <= cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
|
|
constexpr cpp_double_fp_backend(const OtherFloatType& f)
|
|
: data(f, static_cast<float_type>(0.0F)) { }
|
|
|
|
template <typename OtherFloatType,
|
|
typename ::std::enable_if<( cpp_df_qf_detail::is_floating_point<OtherFloatType>::value
|
|
&& (cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatType>::digits > cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
|
|
constexpr cpp_double_fp_backend(const OtherFloatType& f)
|
|
: data(static_cast<float_type>(f),
|
|
static_cast<float_type>(f - static_cast<OtherFloatType>(static_cast<float_type>(f)))) { }
|
|
|
|
// Construtor from another kind of cpp_double_fp_backend<> object.
|
|
|
|
template <typename OtherFloatType,
|
|
typename ::std::enable_if<( cpp_df_qf_detail::is_floating_point<OtherFloatType>::value
|
|
&& (!std::is_same<FloatingPointType, OtherFloatType>::value))>::type const* = nullptr>
|
|
constexpr cpp_double_fp_backend(const cpp_double_fp_backend<OtherFloatType>& 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 <typename SignedIntegralType,
|
|
typename ::std::enable_if<( boost::multiprecision::detail::is_integral<SignedIntegralType>::value
|
|
&& (!boost::multiprecision::detail::is_unsigned<SignedIntegralType>::value)
|
|
&& ((static_cast<int>(sizeof(SignedIntegralType) * 8) - 1) <= cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
|
|
constexpr cpp_double_fp_backend(const SignedIntegralType& n)
|
|
: data(static_cast<float_type>(n), static_cast<float_type>(0.0F)) { }
|
|
|
|
template <typename UnsignedIntegralType,
|
|
typename ::std::enable_if<( boost::multiprecision::detail::is_integral<UnsignedIntegralType>::value
|
|
&& boost::multiprecision::detail::is_unsigned<UnsignedIntegralType>::value
|
|
&& (static_cast<int>(sizeof(UnsignedIntegralType) * 8u) <= cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits))>::type const* = nullptr>
|
|
constexpr cpp_double_fp_backend(const UnsignedIntegralType& u)
|
|
: data(static_cast<float_type>(u), static_cast<float_type>(0.0F)) { }
|
|
|
|
// Constructors from integers which hold more information than *this can contain.
|
|
template <typename UnsignedIntegralType,
|
|
typename ::std::enable_if<( boost::multiprecision::detail::is_integral<UnsignedIntegralType>::value
|
|
&& boost::multiprecision::detail::is_unsigned<UnsignedIntegralType>::value
|
|
&& (static_cast<int>(sizeof(UnsignedIntegralType) * 8u) > cpp_df_qf_detail::ccmath::numeric_limits<float_type>::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<local_unsigned_integral_type>
|
|
(
|
|
static_cast<local_unsigned_integral_type>
|
|
(
|
|
static_cast<local_unsigned_integral_type>(UINT8_C(1)) << static_cast<unsigned>(std::numeric_limits<float_type>::digits)
|
|
)
|
|
- static_cast<local_unsigned_integral_type>(UINT8_C(1))
|
|
)
|
|
};
|
|
|
|
data.second = static_cast<float_type>(0.0F);
|
|
data.first = static_cast<float_type>(u & flt_mask);
|
|
|
|
int u_index { 1 };
|
|
|
|
while (u > static_cast<local_unsigned_integral_type>(UINT8_C(0)))
|
|
{
|
|
u >>= static_cast<unsigned>(std::numeric_limits<float_type>::digits);
|
|
|
|
const float_type
|
|
xhi
|
|
{
|
|
static_cast<float_type>(u & flt_mask)
|
|
* cpp_df_qf_detail::pow2_maker<float_type>::value(std::numeric_limits<float_type>::digits * u_index)
|
|
};
|
|
|
|
eval_add(*this, cpp_double_fp_backend(xhi));
|
|
|
|
++u_index;
|
|
}
|
|
}
|
|
|
|
template <typename SignedIntegralType,
|
|
typename ::std::enable_if<( boost::multiprecision::detail::is_integral<SignedIntegralType>::value
|
|
&& (!boost::multiprecision::detail::is_unsigned<SignedIntegralType>::value)
|
|
&& ((static_cast<int>(sizeof(SignedIntegralType) * 8) - 1) > cpp_df_qf_detail::ccmath::numeric_limits<float_type>::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<SignedIntegralType>::type;
|
|
|
|
const local_unsigned_integral_type
|
|
u_val
|
|
{
|
|
(!is_neg)
|
|
? static_cast<local_unsigned_integral_type>(n)
|
|
: static_cast<local_unsigned_integral_type>
|
|
(
|
|
static_cast<local_unsigned_integral_type>(~n)
|
|
+ static_cast<local_unsigned_integral_type>(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<float_type, float_type>& 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<rep_type&&>(other.data);
|
|
|
|
return *this;
|
|
}
|
|
|
|
// Assignment operator from another kind of cpp_double_fp_backend<> object.
|
|
template <typename OtherFloatType,
|
|
typename ::std::enable_if<( cpp_df_qf_detail::is_floating_point<OtherFloatType>::value
|
|
&& (!std::is_same<FloatingPointType, OtherFloatType>::value))>::type const* = nullptr>
|
|
constexpr auto operator=(const cpp_double_fp_backend<OtherFloatType>& other) -> cpp_double_fp_backend&
|
|
{
|
|
return operator=(cpp_double_fp_backend(other));
|
|
}
|
|
|
|
template <typename OtherFloatType>
|
|
constexpr auto operator=(const OtherFloatType f) -> typename ::std::enable_if<cpp_df_qf_detail::is_floating_point<OtherFloatType>::value, cpp_double_fp_backend&>::type
|
|
{
|
|
return operator=(cpp_double_fp_backend(f));
|
|
}
|
|
|
|
template <typename IntegralType>
|
|
constexpr auto operator=(const IntegralType n) -> typename ::std::enable_if<boost::multiprecision::detail::is_integral<IntegralType>::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<const std::uint16_t*>(&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<void>(hash_one(result, data.first));
|
|
static_cast<void>(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<float_type>::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<float_type>::n_shl) },
|
|
cpp_df_qf_detail::split_maker<float_type>::n_shl
|
|
);
|
|
}
|
|
else
|
|
{
|
|
hc = c - float_type { c - C };
|
|
}
|
|
|
|
float_type u { cpp_df_qf_detail::split_maker<float_type>::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<float_type>::n_shl) },
|
|
cpp_df_qf_detail::split_maker<float_type>::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<std::uint32_t>(p);
|
|
|
|
for (;;)
|
|
{
|
|
if (static_cast<std::uint_fast8_t>(p_local & static_cast<std::uint32_t>(UINT8_C(1))) != static_cast<std::uint_fast8_t>(UINT8_C(0)))
|
|
{
|
|
result.mul_unchecked(y);
|
|
}
|
|
|
|
p_local >>= 1U;
|
|
|
|
if (p_local == static_cast<std::uint32_t>(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<typename cpp_double_fp_backend::rep_type&&>(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<float_type>::max)()
|
|
* (
|
|
static_cast<float_type>(1.0F)
|
|
- static_cast<float_type>(1.5F) * cpp_df_qf_detail::ccmath::unsafe::sqrt(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::epsilon())
|
|
)
|
|
};
|
|
|
|
constexpr float_type
|
|
lo_part
|
|
{
|
|
cpp_df_qf_detail::ccmath::unsafe::ldexp
|
|
(
|
|
(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)(),
|
|
-cpp_df_qf_detail::ccmath::numeric_limits<float_type>::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<float_type>::min)(),
|
|
cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits
|
|
)
|
|
};
|
|
|
|
static_assert
|
|
(
|
|
eval_gt(my_value_min_constexpr, cpp_double_fp_backend { (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::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<float_type>(NAN), float_type { 0.0F });
|
|
}
|
|
|
|
static constexpr auto my_value_inf() noexcept -> cpp_double_fp_backend
|
|
{
|
|
return cpp_double_fp_backend(static_cast<float_type>(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<static_cast<unsigned>(my_digits), digit_base_2, void, int, cpp_df_qf_detail::ccmath::numeric_limits<float_type>::min_exponent, cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max_exponent>;
|
|
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<cpp_bin_float_read_write_backend_type, boost::multiprecision::et_off>;
|
|
|
|
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<float_type>::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<float_type>::n_shl);
|
|
|
|
hu = cpp_df_qf_detail::ccmath::ldexp(data.first - C, cpp_df_qf_detail::split_maker<float_type>::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<float_type>::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<float_type>::n_shl);
|
|
|
|
hv = cpp_df_qf_detail::ccmath::ldexp(v.data.first - c, cpp_df_qf_detail::split_maker<float_type>::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 <typename OtherFloatingPointType,
|
|
typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits10 * 2) < 16))>::type const*>
|
|
friend constexpr auto eval_exp(cpp_double_fp_backend<OtherFloatingPointType>& result, const cpp_double_fp_backend<OtherFloatingPointType>& x) -> void;
|
|
|
|
template <typename OtherFloatingPointType,
|
|
typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value && (((cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits10 * 2) <= 36)))>::type const*>
|
|
friend constexpr auto eval_exp(cpp_double_fp_backend<OtherFloatingPointType>& result, const cpp_double_fp_backend<OtherFloatingPointType>& x) -> void;
|
|
|
|
template <typename OtherFloatingPointType,
|
|
typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits10 * 2) > 36))>::type const*>
|
|
friend constexpr auto eval_exp(cpp_double_fp_backend<OtherFloatingPointType>& result, const cpp_double_fp_backend<OtherFloatingPointType>& x) -> void;
|
|
};
|
|
|
|
template <typename FloatingPointType>
|
|
auto cpp_double_fp_backend<FloatingPointType>::rd_string(const char* pstr) -> bool
|
|
{
|
|
using local_double_fp_type = cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
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<void>(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<cpp_bin_float_read_write_exp_type>(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<cpp_bin_float_read_write_type>::max)());
|
|
|
|
if (is_definitely_inf)
|
|
{
|
|
static_cast<void>(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<float_type>::digits };
|
|
|
|
const auto has_pow2_scaling_for_small_input =
|
|
(
|
|
expval_from_f_bin < static_cast<int>(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<unsigned>
|
|
(
|
|
static_cast<int>
|
|
(
|
|
(local_double_fp_type::my_digits / cpp_df_qf_detail::ccmath::numeric_limits<local_builtin_float_type>::digits)
|
|
+ (((local_double_fp_type::my_digits % cpp_df_qf_detail::ccmath::numeric_limits<local_builtin_float_type>::digits) != 0) ? 1 : 0)
|
|
)
|
|
* cpp_df_qf_detail::ccmath::numeric_limits<local_builtin_float_type>::digits
|
|
);
|
|
|
|
for(auto i = static_cast<unsigned>(UINT8_C(0));
|
|
i < dig_lim;
|
|
i = static_cast<unsigned>(i + static_cast<unsigned>(cpp_df_qf_detail::ccmath::numeric_limits<local_builtin_float_type>::digits)))
|
|
{
|
|
const local_builtin_float_type f { static_cast<local_builtin_float_type>(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 <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_constants.hpp>
|
|
|
|
namespace boost { namespace multiprecision { namespace backends {
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr int cpp_double_fp_backend<FloatingPointType>::my_digits;
|
|
template <typename FloatingPointType>
|
|
constexpr int cpp_double_fp_backend<FloatingPointType>::my_digits10;
|
|
template <typename FloatingPointType>
|
|
constexpr int cpp_double_fp_backend<FloatingPointType>::my_max_digits10;
|
|
template <typename FloatingPointType>
|
|
constexpr int cpp_double_fp_backend<FloatingPointType>::my_max_exponent;
|
|
template <typename FloatingPointType>
|
|
constexpr int cpp_double_fp_backend<FloatingPointType>::my_min_exponent;
|
|
template <typename FloatingPointType>
|
|
constexpr int cpp_double_fp_backend<FloatingPointType>::my_max_exponent10;
|
|
template <typename FloatingPointType>
|
|
constexpr int cpp_double_fp_backend<FloatingPointType>::my_min_exponent10;
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr cpp_double_fp_backend<FloatingPointType> operator+(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) { return cpp_double_fp_backend<FloatingPointType>(a) += b; }
|
|
template <typename FloatingPointType>
|
|
constexpr cpp_double_fp_backend<FloatingPointType> operator-(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) { return cpp_double_fp_backend<FloatingPointType>(a) -= b; }
|
|
template <typename FloatingPointType>
|
|
constexpr cpp_double_fp_backend<FloatingPointType> operator*(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) { return cpp_double_fp_backend<FloatingPointType>(a) *= b; }
|
|
template <typename FloatingPointType>
|
|
constexpr cpp_double_fp_backend<FloatingPointType> operator/(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) { return cpp_double_fp_backend<FloatingPointType>(a) /= b; }
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_add(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void { result += x; }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_add(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void { result = cpp_double_fp_backend<FloatingPointType>(a) += b; }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_subtract(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void { result -= x; }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_subtract(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void { result = cpp_double_fp_backend<FloatingPointType>(a) -= b; }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_multiply(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void { result *= x; }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_multiply(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void { result = cpp_double_fp_backend<FloatingPointType>(a) *= b; }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_divide(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void { result /= x; }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_divide(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> void { result = cpp_double_fp_backend<FloatingPointType>(a) /= b; }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_eq(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool { return (a.compare(b) == 0); }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_lt(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool { return (a.compare(b) == -1); }
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_gt(const cpp_double_fp_backend<FloatingPointType>& a, const cpp_double_fp_backend<FloatingPointType>& b) -> bool { return (a.compare(b) == 1); }
|
|
template <typename FloatingPointType>
|
|
|
|
constexpr auto eval_is_zero(const cpp_double_fp_backend<FloatingPointType>& x) -> bool
|
|
{
|
|
return x.iszero_unchecked();
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_signbit(const cpp_double_fp_backend<FloatingPointType>& x) -> int
|
|
{
|
|
return (x.isneg_unchecked() ? 1 : 0);
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_fabs(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a) -> void
|
|
{
|
|
result = a;
|
|
|
|
if (a.isneg_unchecked())
|
|
{
|
|
result.negate();
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_frexp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, int* v) -> void
|
|
{
|
|
using local_backend_type = cpp_double_fp_backend<FloatingPointType>;
|
|
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 <typename FloatingPointType>
|
|
constexpr auto eval_ldexp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& a, int v) -> void
|
|
{
|
|
using local_backend_type = cpp_double_fp_backend<FloatingPointType>;
|
|
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 <typename FloatingPointType>
|
|
constexpr auto eval_floor(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
|
|
{
|
|
using local_backend_type = cpp_double_fp_backend<FloatingPointType>;
|
|
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 <typename FloatingPointType>
|
|
constexpr auto eval_ceil(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
|
|
{
|
|
// Compute -floor(-x);
|
|
eval_floor(result, -x);
|
|
|
|
result.negate();
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_fpclassify(const cpp_double_fp_backend<FloatingPointType>& 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<FloatingPointType>;
|
|
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<local_float_type>::min)()))
|
|
{
|
|
return FP_SUBNORMAL;
|
|
}
|
|
else
|
|
{
|
|
return FP_NORMAL;
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_sqrt(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& o) -> void
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
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<local_float_type>::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 <typename FloatingPointType>
|
|
constexpr auto eval_pow(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x, const cpp_double_fp_backend<FloatingPointType>& a) -> void
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
constexpr double_float_type zero { 0 };
|
|
|
|
result = zero;
|
|
|
|
signed long long val_nll { };
|
|
|
|
eval_convert_to(&val_nll, a);
|
|
|
|
const int na { static_cast<int>(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 <typename FloatingPointType,
|
|
typename IntegralType>
|
|
constexpr auto eval_pow(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x, IntegralType p) -> typename ::std::enable_if<boost::multiprecision::detail::is_integral<IntegralType>::value, void>::type
|
|
{
|
|
const int fpc { eval_fpclassify(x) };
|
|
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
using local_integral_type = IntegralType;
|
|
|
|
const bool p_is_odd { (static_cast<local_integral_type>(p & 1) != static_cast<local_integral_type>(0)) };
|
|
|
|
if (p == static_cast<local_integral_type>(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<local_integral_type>(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<local_integral_type>(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<local_integral_type>(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 <typename FloatingPointType,
|
|
typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) < 16))>::type const*>
|
|
constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
|
|
{
|
|
const int fpc { eval_fpclassify(x) };
|
|
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
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<local_float_type>()
|
|
: double_float_type(1U) / cpp_df_qf_detail::constant_df_exp1<local_float_type>());
|
|
}
|
|
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<local_float_type>());
|
|
|
|
// 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<local_float_type>()), -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<int>(lln) };
|
|
|
|
if (n > 0)
|
|
{
|
|
eval_ldexp(result, double_float_type(result), n);
|
|
}
|
|
}
|
|
|
|
if (b_neg)
|
|
{
|
|
result = double_float_type(1U) / result;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType,
|
|
typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && (((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) <= 36)))>::type const*>
|
|
constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
|
|
{
|
|
const int fpc { eval_fpclassify(x) };
|
|
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
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<local_float_type>()
|
|
: double_float_type(1U) / cpp_df_qf_detail::constant_df_exp1<local_float_type>());
|
|
}
|
|
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<local_float_type>());
|
|
|
|
// 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<local_float_type>()), -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<int>(lln) };
|
|
|
|
if (n > 0)
|
|
{
|
|
eval_ldexp(result, double_float_type(result), n);
|
|
}
|
|
}
|
|
|
|
if (b_neg)
|
|
{
|
|
result = double_float_type(1U) / result;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType,
|
|
typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) > 36))>::type const*>
|
|
constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
|
|
{
|
|
const int fpc { eval_fpclassify(x) };
|
|
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
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<local_float_type>()
|
|
: double_float_type(1U) / cpp_df_qf_detail::constant_df_exp1<local_float_type>());
|
|
}
|
|
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<local_float_type>());
|
|
|
|
// 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<local_float_type>()), -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<int>(lln) };
|
|
|
|
if (n > 0)
|
|
{
|
|
eval_ldexp(result, double_float_type(result), n);
|
|
}
|
|
}
|
|
|
|
if (b_neg)
|
|
{
|
|
result = double_float_type(1U) / result;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_log(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
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<local_float_type>());
|
|
|
|
eval_add(result, xn2);
|
|
}
|
|
else
|
|
{
|
|
result = double_float_type { 0 };
|
|
}
|
|
}
|
|
|
|
namespace detail {
|
|
|
|
template<typename DestType, typename FloatingPointType>
|
|
constexpr auto extract(cpp_double_fp_backend<FloatingPointType>& source) -> DestType
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
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<destination_type>(0);
|
|
|
|
unsigned fail_safe { UINT32_C(32) };
|
|
|
|
constexpr bool
|
|
destination_type_is_longer
|
|
{
|
|
(std::numeric_limits<signed long long>::digits > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
|
|
};
|
|
|
|
using float_extract_type = typename std::conditional<destination_type_is_longer, local_float_type, float>::type;
|
|
|
|
while((source.compare(zero) != 0) && (fail_safe > unsigned { UINT8_C(0) }))
|
|
{
|
|
const float_extract_type next_flt_val { static_cast<float_extract_type>(source.my_first()) };
|
|
|
|
result += static_cast<destination_type>(next_flt_val);
|
|
|
|
eval_subtract(source, double_float_type(next_flt_val));
|
|
|
|
--fail_safe;
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(signed long long* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> void
|
|
{
|
|
const auto fpc = eval_fpclassify(backend);
|
|
|
|
if (fpc != FP_NORMAL)
|
|
{
|
|
*result = static_cast<signed long long>(backend.crep().first);
|
|
}
|
|
else
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
using local_float_type = typename double_float_type::float_type;
|
|
|
|
static_assert(std::is_same<local_float_type, FloatingPointType>::value, "Error something went wrong with the limb type");
|
|
|
|
constexpr bool
|
|
long_long_is_longer
|
|
{
|
|
(std::numeric_limits<signed long long>::digits > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
|
|
};
|
|
|
|
using longer_type = typename ::std::conditional<long_long_is_longer, signed long long, double_float_type>::type;
|
|
|
|
constexpr longer_type my_max_val { static_cast<longer_type>((std::numeric_limits<signed long long>::max)()) };
|
|
constexpr longer_type my_min_val { static_cast<longer_type>((std::numeric_limits<signed long long>::min)()) };
|
|
|
|
constexpr double_float_type my_max_val_dd { static_cast<double_float_type>(my_max_val) };
|
|
constexpr double_float_type my_min_val_dd { static_cast<double_float_type>(my_min_val) };
|
|
|
|
if (backend.compare(my_max_val_dd) >= 0)
|
|
{
|
|
*result = (std::numeric_limits<signed long long>::max)();
|
|
}
|
|
else if (backend.compare(my_min_val_dd) <= 0)
|
|
{
|
|
*result = (std::numeric_limits<signed long long>::min)();
|
|
}
|
|
else
|
|
{
|
|
double_float_type source { backend };
|
|
|
|
*result = detail::extract<signed long long>(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<local_float_type, double>::value || (std::is_same<local_float_type, long double>::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<signed long long>(1);
|
|
}
|
|
}
|
|
|
|
#endif // non-x64
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(unsigned long long* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> void
|
|
{
|
|
const auto fpc = eval_fpclassify(backend);
|
|
|
|
if (fpc != FP_NORMAL)
|
|
{
|
|
*result = static_cast<signed long long>(backend.crep().first);
|
|
}
|
|
else
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
using local_float_type = typename double_float_type::float_type;
|
|
|
|
static_assert(std::is_same<local_float_type, FloatingPointType>::value, "Error something went wrong with the limb type");
|
|
|
|
constexpr bool
|
|
ulong_long_is_longer
|
|
{
|
|
(std::numeric_limits<unsigned long long>::digits > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
|
|
};
|
|
|
|
using longer_type = typename ::std::conditional<ulong_long_is_longer, unsigned long long, double_float_type>::type;
|
|
|
|
constexpr longer_type my_max_val { static_cast<longer_type>((std::numeric_limits<unsigned long long>::max)()) };
|
|
|
|
constexpr double_float_type my_max_val_dd { static_cast<double_float_type>(my_max_val) };
|
|
|
|
if (backend.compare(my_max_val_dd) >= 0)
|
|
{
|
|
*result = (std::numeric_limits<unsigned long long>::max)();
|
|
}
|
|
else
|
|
{
|
|
double_float_type source { backend };
|
|
|
|
*result = detail::extract<unsigned long long>(source);
|
|
}
|
|
}
|
|
}
|
|
|
|
#ifdef BOOST_HAS_INT128
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type
|
|
{
|
|
const auto fpc = eval_fpclassify(backend);
|
|
|
|
if (fpc != FP_NORMAL)
|
|
{
|
|
*result = static_cast<boost::int128_type>(backend.crep().first);
|
|
}
|
|
else
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
using local_float_type = typename double_float_type::float_type;
|
|
|
|
static_assert(std::is_same<local_float_type, FloatingPointType>::value, "Error something went wrong with the limb type");
|
|
|
|
constexpr bool
|
|
n128_is_longer
|
|
{
|
|
((static_cast<int>(sizeof(boost::int128_type) * static_cast<std::size_t>(CHAR_BIT)) - 1) > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
|
|
};
|
|
|
|
using longer_type = typename ::std::conditional<n128_is_longer, boost::int128_type, double_float_type>::type;
|
|
|
|
constexpr boost::int128_type my_max_val_n128 = (((static_cast<boost::int128_type>(1) << (sizeof(boost::int128_type) * CHAR_BIT - 2)) - 1) << 1) + 1;
|
|
constexpr boost::int128_type my_min_val_n128 = static_cast<boost::int128_type>(-my_max_val_n128 - 1);
|
|
|
|
constexpr longer_type my_max_val(static_cast<longer_type>(my_max_val_n128));
|
|
constexpr longer_type my_min_val(static_cast<longer_type>(my_min_val_n128));
|
|
|
|
constexpr double_float_type my_max_val_dd(static_cast<double_float_type>(my_max_val));
|
|
constexpr double_float_type my_min_val_dd(static_cast<double_float_type>(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<boost::int128_type>(source);
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(boost::int128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<!(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type
|
|
{
|
|
const auto fpc = eval_fpclassify(backend);
|
|
|
|
if (fpc != FP_NORMAL)
|
|
{
|
|
*result = static_cast<boost::int128_type>(backend.crep().first);
|
|
}
|
|
else
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
double_float_type source { backend };
|
|
|
|
*result = detail::extract<boost::int128_type>(source);
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type
|
|
{
|
|
const auto fpc = eval_fpclassify(backend);
|
|
|
|
if (fpc != FP_NORMAL)
|
|
{
|
|
*result = static_cast<boost::int128_type>(backend.crep().first);
|
|
}
|
|
else
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
using local_float_type = typename double_float_type::float_type;
|
|
|
|
static_assert(std::is_same<local_float_type, FloatingPointType>::value, "Error something went wrong with the limb type");
|
|
|
|
constexpr bool
|
|
u128_is_longer
|
|
{
|
|
(static_cast<int>(sizeof(boost::uint128_type) * static_cast<std::size_t>(CHAR_BIT)) > cpp_df_qf_detail::ccmath::numeric_limits<local_float_type>::digits)
|
|
};
|
|
|
|
using longer_type = typename ::std::conditional<u128_is_longer, boost::uint128_type, double_float_type>::type;
|
|
|
|
constexpr boost::int128_type my_max_val_u128 = static_cast<boost::uint128_type>(~static_cast<boost::uint128_type>(0));
|
|
|
|
constexpr longer_type my_max_val(static_cast<longer_type>(my_max_val_u128));
|
|
|
|
constexpr double_float_type my_max_val_dd(static_cast<double_float_type>(my_max_val));
|
|
|
|
if (backend.compare(my_max_val_dd) >= 0)
|
|
{
|
|
*result = my_max_val;
|
|
}
|
|
else
|
|
{
|
|
double_float_type source { backend };
|
|
|
|
*result = detail::extract<boost::uint128_type>(source);
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
constexpr auto eval_convert_to(boost::uint128_type* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename std::enable_if<!(cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits > 24), void>::type
|
|
{
|
|
const auto fpc = eval_fpclassify(backend);
|
|
|
|
if (fpc != FP_NORMAL)
|
|
{
|
|
*result = static_cast<boost::uint128_type>(backend.crep().first);
|
|
}
|
|
else
|
|
{
|
|
using double_float_type = cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
double_float_type source { backend };
|
|
|
|
*result = detail::extract<boost::uint128_type>(source);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
template <typename FloatingPointType,
|
|
typename OtherFloatingPointType>
|
|
constexpr auto eval_convert_to(OtherFloatingPointType* result, const cpp_double_fp_backend<FloatingPointType>& backend) -> typename ::std::enable_if<cpp_df_qf_detail::is_floating_point<OtherFloatingPointType>::value>::type
|
|
{
|
|
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<OtherFloatingPointType>(backend.my_first());
|
|
}
|
|
else
|
|
{
|
|
BOOST_IF_CONSTEXPR(cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits > cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits)
|
|
{
|
|
*result = static_cast<OtherFloatingPointType>(backend.my_first());
|
|
*result += static_cast<OtherFloatingPointType>(backend.my_second());
|
|
}
|
|
else
|
|
{
|
|
cpp_double_fp_backend<FloatingPointType> source = backend;
|
|
|
|
*result = 0;
|
|
|
|
for(auto digit_count = static_cast<int>(0);
|
|
digit_count < cpp_double_fp_backend<FloatingPointType>::my_digits;
|
|
digit_count += cpp_df_qf_detail::ccmath::numeric_limits<OtherFloatingPointType>::digits)
|
|
{
|
|
const auto next = static_cast<OtherFloatingPointType>(source.my_first());
|
|
|
|
*result += next;
|
|
|
|
eval_subtract(source, cpp_double_fp_backend<FloatingPointType>(next));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename FloatingPointType>
|
|
auto hash_value(const cpp_double_fp_backend<FloatingPointType>& a) -> ::std::size_t
|
|
{
|
|
return a.hash();
|
|
}
|
|
|
|
} // namespace backends
|
|
|
|
using backends::cpp_double_fp_backend;
|
|
|
|
using cpp_double_float = number<cpp_double_fp_backend<float>, boost::multiprecision::et_off>;
|
|
using cpp_double_double = number<cpp_double_fp_backend<double>, boost::multiprecision::et_off>;
|
|
using cpp_double_long_double = number<cpp_double_fp_backend<long double>, boost::multiprecision::et_off>;
|
|
#ifdef BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128
|
|
using cpp_double_float128 = number<cpp_double_fp_backend<::boost::float128_type>, boost::multiprecision::et_off>;
|
|
#endif
|
|
|
|
} } // namespace boost::multiprecision
|
|
|
|
namespace std {
|
|
|
|
// Specialization of numeric_limits for boost::multiprecision::number<cpp_double_fp_backend<>>
|
|
template <typename FloatingPointType,
|
|
const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
BOOST_MP_DF_QF_NUM_LIMITS_CLASS_TYPE numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >
|
|
{
|
|
private:
|
|
using local_float_type = FloatingPointType;
|
|
using inner_self_type = boost::multiprecision::cpp_double_fp_backend<local_float_type>;
|
|
|
|
using self_type = boost::multiprecision::number<inner_self_type, ExpressionTemplatesOption>;
|
|
|
|
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<self_type>(inner_self_type::my_value_min()); }
|
|
static constexpr auto(max) () noexcept -> self_type { return static_cast<self_type>(inner_self_type::my_value_max()); }
|
|
static constexpr auto lowest () noexcept -> self_type { return static_cast<self_type>(-(max)()); }
|
|
static constexpr auto epsilon () noexcept -> self_type { return static_cast<self_type>(inner_self_type::my_value_eps()); }
|
|
static constexpr auto round_error () noexcept -> self_type { return static_cast<self_type>(static_cast<local_float_type>(0.5)); }
|
|
static constexpr auto denorm_min () noexcept -> self_type { return static_cast<self_type>((min)()); }
|
|
static constexpr auto infinity () noexcept -> self_type { return static_cast<self_type>(inner_self_type::my_value_inf()); }
|
|
static constexpr auto quiet_NaN () noexcept -> self_type { return static_cast<self_type>(inner_self_type::my_value_nan()); }
|
|
static constexpr auto signaling_NaN() noexcept -> self_type { return static_cast<self_type>(static_cast<local_float_type>(0.0)); }
|
|
};
|
|
|
|
} // namespace std
|
|
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_specialized;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_signed;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_integer;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_exact;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_bounded;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_modulo;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::is_iec559;
|
|
|
|
#ifdef BOOST_MSVC
|
|
#pragma warning(push)
|
|
#pragma warning(disable:4996)
|
|
#endif
|
|
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr std::float_denorm_style std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_denorm;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_denorm_loss;
|
|
|
|
#ifdef BOOST_MSVC
|
|
#pragma warning(pop)
|
|
#endif // deprecated members
|
|
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_infinity;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_quiet_NaN;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::has_signaling_NaN;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::traps;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr bool std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::tinyness_before;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr std::float_round_style std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::round_style;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr int std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::radix;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr int std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::digits;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr int std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::digits10;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr int std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::max_digits10;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr int std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::max_exponent;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr int std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::min_exponent;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr int std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::max_exponent10;
|
|
template <typename FloatingPointType, const boost::multiprecision::expression_template_option ExpressionTemplatesOption>
|
|
constexpr int std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplatesOption> >::min_exponent10;
|
|
|
|
#if defined(BOOST_MP_MATH_AVAILABLE)
|
|
namespace boost { namespace math { namespace policies {
|
|
|
|
template <class FloatingPointType, class Policy, boost::multiprecision::expression_template_option ExpressionTemplates>
|
|
struct precision<boost::multiprecision::number<boost::multiprecision::cpp_double_fp_backend<FloatingPointType>, ExpressionTemplates>, Policy>
|
|
{
|
|
private:
|
|
using my_multiprecision_backend_type = boost::multiprecision::cpp_double_fp_backend<FloatingPointType>;
|
|
|
|
using digits_2 = digits2<my_multiprecision_backend_type::my_digits>;
|
|
|
|
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
|