/////////////////////////////////////////////////////////////////////////////// // Copyright 2011 John Maddock. Distributed under the Boost // Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_MATH_BN_MPFI_HPP #define BOOST_MATH_BN_MPFI_HPP #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifndef BOOST_MULTIPRECISION_MPFI_DEFAULT_PRECISION #define BOOST_MULTIPRECISION_MPFI_DEFAULT_PRECISION 20 #endif namespace boost { namespace multiprecision { namespace backends { template struct mpfi_float_backend; template struct debug_adaptor; } // namespace backends template struct number_category > : public std::integral_constant {}; struct interval_error : public std::runtime_error { interval_error(const std::string& s) : std::runtime_error(s) {} }; namespace detail { template <> struct is_variable_precision > : public std::integral_constant {}; } // namespace detail namespace backends { namespace detail { inline int mpfi_sgn(mpfi_srcptr p) { if (mpfi_is_zero(p)) return 0; if (mpfi_is_strictly_pos(p)) return 1; if (mpfi_is_strictly_neg(p)) return -1; BOOST_THROW_EXCEPTION(interval_error("Sign of interval is ambiguous.")); } template struct mpfi_float_imp; template struct mpfi_float_imp { #ifdef BOOST_HAS_LONG_LONG using signed_types = std::tuple ; using unsigned_types = std::tuple; #else using signed_types = std::tuple ; using unsigned_types = std::tuple; #endif using float_types = std::tuple; using exponent_type = long ; mpfi_float_imp() { mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); mpfi_set_ui(m_data, 0u); } mpfi_float_imp(unsigned prec) { mpfi_init2(m_data, prec); mpfi_set_ui(m_data, 0u); } mpfi_float_imp(const mpfi_float_imp& o) { mpfi_init2(m_data, preserve_source_precision() ? mpfi_get_prec(o.data()) : boost::multiprecision::detail::digits10_2_2(get_default_precision())); if (o.m_data[0].left._mpfr_d) mpfi_set(m_data, o.m_data); } template mpfi_float_imp(const mpfr_float_imp& o) { mpfi_init2(m_data, (D == 0 ? preserve_component_precision() : preserve_related_precision()) ? mpfr_get_prec(o.data()) : boost::multiprecision::detail::digits10_2_2(get_default_precision())); if (o.data()[0]._mpfr_d) mpfi_set_fr(m_data, o.data()); } // rvalue copy mpfi_float_imp(mpfi_float_imp&& o) noexcept { mpfr_prec_t binary_default_precision = boost::multiprecision::detail::digits10_2_2(get_default_precision()); if ((this->get_default_options() != variable_precision_options::preserve_target_precision) || (mpfi_get_prec(o.data()) == binary_default_precision)) { m_data[0] = o.m_data[0]; o.m_data[0].left._mpfr_d = 0; } else { // NOTE: C allocation interface must not throw: mpfi_init2(m_data, binary_default_precision); mpfi_set(m_data, o.m_data); } } mpfi_float_imp& operator=(const mpfi_float_imp& o) { if (this != &o) { if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, preserve_source_precision() ? mpfi_get_prec(o.m_data) : boost::multiprecision::detail::digits10_2_2(get_default_precision())); else if (preserve_source_precision() && (mpfi_get_prec(o.data()) != mpfi_get_prec(data()))) { mpfi_set_prec(m_data, mpfi_get_prec(o.m_data)); } mpfi_set(m_data, o.m_data); } return *this; } // rvalue assign mpfi_float_imp& operator=(mpfi_float_imp&& o) noexcept { if ((this->get_default_options() != variable_precision_options::preserve_target_precision) || (mpfi_get_prec(o.data()) == mpfi_get_prec(data()))) mpfi_swap(m_data, o.m_data); else *this = static_cast(o); return *this; } #ifdef BOOST_HAS_LONG_LONG #ifdef _MPFR_H_HAVE_INTMAX_T mpfi_float_imp& operator=(boost::ulong_long_type i) { if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); mpfr_set_uj(left_data(), i, GMP_RNDD); mpfr_set_uj(right_data(), i, GMP_RNDU); return *this; } mpfi_float_imp& operator=(boost::long_long_type i) { if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); mpfr_set_sj(left_data(), i, GMP_RNDD); mpfr_set_sj(right_data(), i, GMP_RNDU); return *this; } #else mpfi_float_imp& operator=(boost::ulong_long_type i) { if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); boost::ulong_long_type mask = ((((1uLL << (std::numeric_limits::digits - 1)) - 1) << 1) | 1u); unsigned shift = 0; mpfi_t t; mpfi_init2(t, (std::max)(static_cast(std::numeric_limits::digits), static_cast(multiprecision::detail::digits10_2_2(digits10)))); mpfi_set_ui(m_data, 0); while (i) { mpfi_set_ui(t, static_cast(i & mask)); if (shift) mpfi_mul_2exp(t, t, shift); mpfi_add(m_data, m_data, t); shift += std::numeric_limits::digits; i >>= std::numeric_limits::digits; } mpfi_clear(t); return *this; } mpfi_float_imp& operator=(boost::long_long_type i) { if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); bool neg = i < 0; *this = boost::multiprecision::detail::unsigned_abs(i); if (neg) mpfi_neg(m_data, m_data); return *this; } #endif #endif mpfi_float_imp& operator=(unsigned long i) { if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); mpfi_set_ui(m_data, i); return *this; } mpfi_float_imp& operator=(long i) { if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); mpfi_set_si(m_data, i); return *this; } mpfi_float_imp& operator=(double d) { if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); mpfi_set_d(m_data, d); return *this; } mpfi_float_imp& operator=(long double a) { if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); mpfr_set_ld(left_data(), a, GMP_RNDD); mpfr_set_ld(right_data(), a, GMP_RNDU); return *this; } mpfi_float_imp& operator=(const char* s) { using default_ops::eval_fpclassify; if (m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); if (s && (*s == '{')) { mpfr_float_backend a, b; std::string part; const char* p = ++s; while (*p && (*p != ',') && (*p != '}')) ++p; part.assign(s + 1, p); a = part.c_str(); s = p; if (*p && (*p != '}')) { ++p; while (*p && (*p != ',') && (*p != '}')) ++p; part.assign(s + 1, p); } else part.erase(); b = part.c_str(); if (eval_fpclassify(a) == (int)FP_NAN) { mpfi_set_fr(this->data(), a.data()); } else if (eval_fpclassify(b) == (int)FP_NAN) { mpfi_set_fr(this->data(), b.data()); } else { if (a.compare(b) > 0) { BOOST_THROW_EXCEPTION(std::runtime_error("Attempt to create interval with invalid range (start is greater than end).")); } mpfi_interv_fr(m_data, a.data(), b.data()); } } else if (mpfi_set_str(m_data, s, 10) != 0) { BOOST_THROW_EXCEPTION(std::runtime_error(std::string("Unable to parse string \"") + s + std::string("\"as a valid floating point number."))); } return *this; } void swap(mpfi_float_imp& o) noexcept { mpfi_swap(m_data, o.m_data); } std::string str(std::streamsize digits, std::ios_base::fmtflags f) const { BOOST_ASSERT(m_data[0].left._mpfr_d); mpfr_float_backend a, b; mpfi_get_left(a.data(), m_data); mpfi_get_right(b.data(), m_data); if (a.compare(b) == 0) return a.str(digits, f); return "{" + a.str(digits, f) + "," + b.str(digits, f) + "}"; } ~mpfi_float_imp() noexcept { if (m_data[0].left._mpfr_d) mpfi_clear(m_data); } void negate() noexcept { BOOST_ASSERT(m_data[0].left._mpfr_d); mpfi_neg(m_data, m_data); } int compare(const mpfi_float_imp& o) const noexcept { BOOST_ASSERT(m_data[0].left._mpfr_d && o.m_data[0].left._mpfr_d); if (mpfr_cmp(right_data(), o.left_data()) < 0) return -1; if (mpfr_cmp(left_data(), o.right_data()) > 0) return 1; if ((mpfr_cmp(left_data(), o.left_data()) == 0) && (mpfr_cmp(right_data(), o.right_data()) == 0)) return 0; BOOST_THROW_EXCEPTION(interval_error("Ambiguous comparison between two values.")); return 0; } template int compare(V v) const noexcept { mpfi_float_imp d; d = v; return compare(d); } mpfi_t& data() noexcept { BOOST_ASSERT(m_data[0].left._mpfr_d); return m_data; } const mpfi_t& data() const noexcept { BOOST_ASSERT(m_data[0].left._mpfr_d); return m_data; } mpfr_ptr left_data() noexcept { BOOST_ASSERT(m_data[0].left._mpfr_d); return &(m_data[0].left); } mpfr_srcptr left_data() const noexcept { BOOST_ASSERT(m_data[0].left._mpfr_d); return &(m_data[0].left); } mpfr_ptr right_data() noexcept { BOOST_ASSERT(m_data[0].left._mpfr_d); return &(m_data[0].right); } mpfr_srcptr right_data() const noexcept { BOOST_ASSERT(m_data[0].left._mpfr_d); return &(m_data[0].right); } protected: mpfi_t m_data; static boost::multiprecision::detail::precision_type& get_global_default_precision() noexcept { static boost::multiprecision::detail::precision_type val(BOOST_MULTIPRECISION_MPFI_DEFAULT_PRECISION); return val; } static unsigned& get_default_precision() noexcept { static BOOST_MP_THREAD_LOCAL unsigned val(get_global_default_precision()); return val; } #ifndef BOOST_MT_NO_ATOMIC_INT static std::atomic& get_global_default_options() noexcept #else static variable_precision_options& get_global_default_options() noexcept #endif { #ifndef BOOST_MT_NO_ATOMIC_INT static std::atomic val{variable_precision_options::preserve_related_precision}; #else static variable_precision_optionss val{variable_precision_options::preserve_related_precision}; #endif return val; } static variable_precision_options& get_default_options() noexcept { static BOOST_MP_THREAD_LOCAL variable_precision_options val(get_global_default_options()); return val; } static bool preserve_source_precision() noexcept { return get_default_options() >= variable_precision_options::preserve_source_precision; } static bool preserve_component_precision() noexcept { return get_default_options() >= variable_precision_options::preserve_component_precision; } static bool preserve_related_precision() noexcept { return get_default_options() >= variable_precision_options::preserve_related_precision; } static bool preserve_all_precision() noexcept { return get_default_options() >= variable_precision_options::preserve_all_precision; } }; } // namespace detail template struct mpfi_float_backend : public detail::mpfi_float_imp { mpfi_float_backend() : detail::mpfi_float_imp() {} mpfi_float_backend(const mpfi_float_backend& o) : detail::mpfi_float_imp(o) {} // rvalue copy mpfi_float_backend(mpfi_float_backend&& o) : detail::mpfi_float_imp(static_cast&&>(o)) {} template mpfi_float_backend(const mpfi_float_backend& val, typename std::enable_if::type* = 0) : detail::mpfi_float_imp() { mpfi_set(this->m_data, val.data()); } template mpfi_float_backend(const mpfr_float_backend& val, typename std::enable_if::type* = 0) : detail::mpfi_float_imp(val) {} template explicit mpfi_float_backend(const mpfi_float_backend& val, typename std::enable_if::type* = 0) : detail::mpfi_float_imp() { mpfi_set(this->m_data, val.data()); } mpfi_float_backend(const mpfi_t val) : detail::mpfi_float_imp() { mpfi_set(this->m_data, val); } mpfi_float_backend& operator=(const mpfi_float_backend& o) { *static_cast*>(this) = static_cast const&>(o); return *this; } template mpfi_float_backend(const mpfr_float_backend& val, typename std::enable_if::type* = 0) : detail::mpfi_float_imp() { mpfi_set_fr(this->m_data, val.data()); } template mpfi_float_backend& operator=(const mpfr_float_backend& val) { mpfi_set_fr(this->m_data, val.data()); return *this; } template explicit mpfi_float_backend(const mpfr_float_backend& val, typename std::enable_if::type* = 0) : detail::mpfi_float_imp() { mpfi_set_fr(this->m_data, val.data()); } // rvalue copy mpfi_float_backend& operator=(mpfi_float_backend&& o) noexcept { *static_cast*>(this) = static_cast&&>(o); return *this; } template mpfi_float_backend& operator=(const V& v) { *static_cast*>(this) = v; return *this; } mpfi_float_backend& operator=(const mpfi_t val) { mpfi_set(this->m_data, val); return *this; } // We don't change our precision here, this is a fixed precision type: template mpfi_float_backend& operator=(const mpfi_float_backend& val) { mpfi_set(this->m_data, val.data()); return *this; } }; template void assign_components(mpfi_float_backend& result, const mpfr_float_backend& a, const mpfr_float_backend& b); template typename std::enable_if, et_on>, V>::value || std::is_convertible::value>::type assign_components(mpfi_float_backend& result, const V& a, const V& b); template <> struct mpfi_float_backend<0> : public detail::mpfi_float_imp<0> { mpfi_float_backend() : detail::mpfi_float_imp<0>() {} mpfi_float_backend(const mpfi_t val) : detail::mpfi_float_imp<0>(mpfi_get_prec(val)) { mpfi_set(this->m_data, val); } mpfi_float_backend(const mpfi_float_backend& o) : detail::mpfi_float_imp<0>(o) {} template mpfi_float_backend(const mpfr_float_backend& val) : detail::mpfi_float_imp<0>(val) {} // rvalue copy mpfi_float_backend(mpfi_float_backend&& o) noexcept : detail::mpfi_float_imp<0>(static_cast&&>(o)) {} mpfi_float_backend(const mpfi_float_backend& o, unsigned digits10) : detail::mpfi_float_imp<0>(multiprecision::detail::digits10_2_2(digits10)) { mpfi_set(this->m_data, o.data()); } template mpfi_float_backend(const V& a, const V& b, unsigned digits10) : detail::mpfi_float_imp<0>(multiprecision::detail::digits10_2_2(digits10)) { boost::multiprecision::detail::scoped_target_precision > opts; assign_components(*this, a, b); } template mpfi_float_backend(const mpfi_float_backend& val) : detail::mpfi_float_imp<0>(mpfi_get_prec(val.data())) { mpfi_set(this->m_data, val.data()); } mpfi_float_backend& operator=(const mpfi_float_backend& o) = default; // rvalue assign mpfi_float_backend& operator=(mpfi_float_backend&& o) noexcept = default; template mpfi_float_backend& operator=(const V& v) { constexpr unsigned d10 = std::is_floating_point::value ? std::numeric_limits::digits10 : std::numeric_limits::digits10 ? 1 + std::numeric_limits::digits10 : 1 + boost::multiprecision::detail::digits2_2_10(std::numeric_limits::digits); if (thread_default_variable_precision_options() >= variable_precision_options::preserve_all_precision) { BOOST_IF_CONSTEXPR(std::is_floating_point::value) { if (std::numeric_limits::digits > mpfi_get_prec(this->data())) mpfi_set_prec(this->data(), std::numeric_limits::digits); } else { if (precision() < d10) this->precision(d10); } } *static_cast*>(this) = v; return *this; } mpfi_float_backend& operator=(const mpfi_t val) { mpfi_set_prec(this->m_data, mpfi_get_prec(val)); mpfi_set(this->m_data, val); return *this; } template mpfi_float_backend& operator=(const mpfi_float_backend& val) { mpfi_set_prec(this->m_data, mpfi_get_prec(val.data())); mpfi_set(this->m_data, val.data()); return *this; } static unsigned thread_default_precision() noexcept { return get_default_precision(); } static void thread_default_precision(unsigned v) noexcept { get_default_precision() = v; } static unsigned default_precision() noexcept { return get_global_default_precision(); } static void default_precision(unsigned v) noexcept { get_global_default_precision() = v; } unsigned precision() const noexcept { return multiprecision::detail::digits2_2_10(mpfi_get_prec(this->m_data)); } void precision(unsigned digits10) noexcept { mpfi_float_backend t(*this, digits10); this->swap(t); } // // Variable precision options: // static variable_precision_options default_variable_precision_options() noexcept { return get_global_default_options(); } static variable_precision_options thread_default_variable_precision_options() noexcept { return get_default_options(); } static void default_variable_precision_options(variable_precision_options opts) { get_global_default_options() = opts; } static void thread_default_variable_precision_options(variable_precision_options opts) { get_default_options() = opts; } }; template inline typename std::enable_if::value, bool>::type eval_eq(const mpfi_float_backend& a, const T& b) noexcept { return a.compare(b) == 0; } template inline typename std::enable_if::value, bool>::type eval_lt(const mpfi_float_backend& a, const T& b) noexcept { return a.compare(b) < 0; } template inline typename std::enable_if::value, bool>::type eval_gt(const mpfi_float_backend& a, const T& b) noexcept { return a.compare(b) > 0; } template inline void eval_add(mpfi_float_backend& result, const mpfi_float_backend& o) { mpfi_add(result.data(), result.data(), o.data()); } template inline void eval_subtract(mpfi_float_backend& result, const mpfi_float_backend& o) { mpfi_sub(result.data(), result.data(), o.data()); } template inline void eval_multiply(mpfi_float_backend& result, const mpfi_float_backend& o) { if ((void*)&result == (void*)&o) mpfi_sqr(result.data(), o.data()); else mpfi_mul(result.data(), result.data(), o.data()); } template inline void eval_divide(mpfi_float_backend& result, const mpfi_float_backend& o) { mpfi_div(result.data(), result.data(), o.data()); } template inline void eval_add(mpfi_float_backend& result, unsigned long i) { mpfi_add_ui(result.data(), result.data(), i); } template inline void eval_subtract(mpfi_float_backend& result, unsigned long i) { mpfi_sub_ui(result.data(), result.data(), i); } template inline void eval_multiply(mpfi_float_backend& result, unsigned long i) { mpfi_mul_ui(result.data(), result.data(), i); } template inline void eval_divide(mpfi_float_backend& result, unsigned long i) { mpfi_div_ui(result.data(), result.data(), i); } template inline void eval_add(mpfi_float_backend& result, long i) { if (i > 0) mpfi_add_ui(result.data(), result.data(), i); else mpfi_sub_ui(result.data(), result.data(), boost::multiprecision::detail::unsigned_abs(i)); } template inline void eval_subtract(mpfi_float_backend& result, long i) { if (i > 0) mpfi_sub_ui(result.data(), result.data(), i); else mpfi_add_ui(result.data(), result.data(), boost::multiprecision::detail::unsigned_abs(i)); } template inline void eval_multiply(mpfi_float_backend& result, long i) { mpfi_mul_ui(result.data(), result.data(), boost::multiprecision::detail::unsigned_abs(i)); if (i < 0) mpfi_neg(result.data(), result.data()); } template inline void eval_divide(mpfi_float_backend& result, long i) { mpfi_div_ui(result.data(), result.data(), boost::multiprecision::detail::unsigned_abs(i)); if (i < 0) mpfi_neg(result.data(), result.data()); } // // Specialised 3 arg versions of the basic operators: // template inline void eval_add(mpfi_float_backend& a, const mpfi_float_backend& x, const mpfi_float_backend& y) { mpfi_add(a.data(), x.data(), y.data()); } template inline void eval_add(mpfi_float_backend& a, const mpfi_float_backend& x, unsigned long y) { mpfi_add_ui(a.data(), x.data(), y); } template inline void eval_add(mpfi_float_backend& a, const mpfi_float_backend& x, long y) { if (y < 0) mpfi_sub_ui(a.data(), x.data(), boost::multiprecision::detail::unsigned_abs(y)); else mpfi_add_ui(a.data(), x.data(), y); } template inline void eval_add(mpfi_float_backend& a, unsigned long x, const mpfi_float_backend& y) { mpfi_add_ui(a.data(), y.data(), x); } template inline void eval_add(mpfi_float_backend& a, long x, const mpfi_float_backend& y) { if (x < 0) { mpfi_ui_sub(a.data(), boost::multiprecision::detail::unsigned_abs(x), y.data()); mpfi_neg(a.data(), a.data()); } else mpfi_add_ui(a.data(), y.data(), x); } template inline void eval_subtract(mpfi_float_backend& a, const mpfi_float_backend& x, const mpfi_float_backend& y) { mpfi_sub(a.data(), x.data(), y.data()); } template inline void eval_subtract(mpfi_float_backend& a, const mpfi_float_backend& x, unsigned long y) { mpfi_sub_ui(a.data(), x.data(), y); } template inline void eval_subtract(mpfi_float_backend& a, const mpfi_float_backend& x, long y) { if (y < 0) mpfi_add_ui(a.data(), x.data(), boost::multiprecision::detail::unsigned_abs(y)); else mpfi_sub_ui(a.data(), x.data(), y); } template inline void eval_subtract(mpfi_float_backend& a, unsigned long x, const mpfi_float_backend& y) { mpfi_ui_sub(a.data(), x, y.data()); } template inline void eval_subtract(mpfi_float_backend& a, long x, const mpfi_float_backend& y) { if (x < 0) { mpfi_add_ui(a.data(), y.data(), boost::multiprecision::detail::unsigned_abs(x)); mpfi_neg(a.data(), a.data()); } else mpfi_ui_sub(a.data(), x, y.data()); } template inline void eval_multiply(mpfi_float_backend& a, const mpfi_float_backend& x, const mpfi_float_backend& y) { if ((void*)&x == (void*)&y) mpfi_sqr(a.data(), x.data()); else mpfi_mul(a.data(), x.data(), y.data()); } template inline void eval_multiply(mpfi_float_backend& a, const mpfi_float_backend& x, unsigned long y) { mpfi_mul_ui(a.data(), x.data(), y); } template inline void eval_multiply(mpfi_float_backend& a, const mpfi_float_backend& x, long y) { if (y < 0) { mpfi_mul_ui(a.data(), x.data(), boost::multiprecision::detail::unsigned_abs(y)); a.negate(); } else mpfi_mul_ui(a.data(), x.data(), y); } template inline void eval_multiply(mpfi_float_backend& a, unsigned long x, const mpfi_float_backend& y) { mpfi_mul_ui(a.data(), y.data(), x); } template inline void eval_multiply(mpfi_float_backend& a, long x, const mpfi_float_backend& y) { if (x < 0) { mpfi_mul_ui(a.data(), y.data(), boost::multiprecision::detail::unsigned_abs(x)); mpfi_neg(a.data(), a.data()); } else mpfi_mul_ui(a.data(), y.data(), x); } template inline void eval_divide(mpfi_float_backend& a, const mpfi_float_backend& x, const mpfi_float_backend& y) { mpfi_div(a.data(), x.data(), y.data()); } template inline void eval_divide(mpfi_float_backend& a, const mpfi_float_backend& x, unsigned long y) { mpfi_div_ui(a.data(), x.data(), y); } template inline void eval_divide(mpfi_float_backend& a, const mpfi_float_backend& x, long y) { if (y < 0) { mpfi_div_ui(a.data(), x.data(), boost::multiprecision::detail::unsigned_abs(y)); a.negate(); } else mpfi_div_ui(a.data(), x.data(), y); } template inline void eval_divide(mpfi_float_backend& a, unsigned long x, const mpfi_float_backend& y) { mpfi_ui_div(a.data(), x, y.data()); } template inline void eval_divide(mpfi_float_backend& a, long x, const mpfi_float_backend& y) { if (x < 0) { mpfi_ui_div(a.data(), boost::multiprecision::detail::unsigned_abs(x), y.data()); mpfi_neg(a.data(), a.data()); } else mpfi_ui_div(a.data(), x, y.data()); } template inline bool eval_is_zero(const mpfi_float_backend& val) noexcept { return 0 != mpfi_is_zero(val.data()); } template inline int eval_get_sign(const mpfi_float_backend& val) { return detail::mpfi_sgn(val.data()); } template inline void eval_convert_to(unsigned long* result, const mpfi_float_backend& val) { mpfr_float_backend t; mpfi_mid(t.data(), val.data()); eval_convert_to(result, t); } template inline void eval_convert_to(long* result, const mpfi_float_backend& val) { mpfr_float_backend t; mpfi_mid(t.data(), val.data()); eval_convert_to(result, t); } #ifdef _MPFR_H_HAVE_INTMAX_T template inline void eval_convert_to(boost::ulong_long_type* result, const mpfi_float_backend& val) { mpfr_float_backend t; mpfi_mid(t.data(), val.data()); eval_convert_to(result, t); } template inline void eval_convert_to(boost::long_long_type* result, const mpfi_float_backend& val) { mpfr_float_backend t; mpfi_mid(t.data(), val.data()); eval_convert_to(result, t); } #endif template inline void eval_convert_to(double* result, const mpfi_float_backend& val) noexcept { *result = mpfi_get_d(val.data()); } template inline void eval_convert_to(long double* result, const mpfi_float_backend& val) noexcept { mpfr_float_backend t; mpfi_mid(t.data(), val.data()); eval_convert_to(result, t); } template inline void assign_components_set_precision(mpfi_float_backend<0>& result, const mpfr_float_backend<0, AllocationType>& a, const mpfr_float_backend<0, AllocationType>& b) { if (result.thread_default_variable_precision_options() >= variable_precision_options::preserve_component_precision) { unsigned long prec = (std::max)(mpfr_get_prec(a.data()), mpfr_get_prec(b.data())); mpfi_set_prec(result.data(), prec); } } template inline void assign_components_set_precision(mpfi_float_backend<0>& result, const mpfr_float_backend& a, const mpfr_float_backend& b) { if (result.thread_default_variable_precision_options() >= variable_precision_options::preserve_related_precision) { unsigned long prec = (std::max)(mpfr_get_prec(a.data()), mpfr_get_prec(b.data())); mpfi_set_prec(result.data(), prec); } } template inline void assign_components_set_precision(mpfi_float_backend&, const mpfr_float_backend&, const mpfr_float_backend&) { } template inline void assign_components(mpfi_float_backend& result, const mpfr_float_backend& a, const mpfr_float_backend& b) { // // This is called from class number's constructors, so if we have variable // precision, then copy the precision of the source variables. // assign_components_set_precision(result, a, b); using default_ops::eval_fpclassify; if (eval_fpclassify(a) == (int)FP_NAN) { mpfi_set_fr(result.data(), a.data()); } else if (eval_fpclassify(b) == (int)FP_NAN) { mpfi_set_fr(result.data(), b.data()); } else { if (a.compare(b) > 0) { BOOST_THROW_EXCEPTION(std::runtime_error("Attempt to create interval with invalid range (start is greater than end).")); } mpfi_interv_fr(result.data(), a.data(), b.data()); } } template inline typename std::enable_if, et_on>, V>::value || std::is_convertible::value>::type assign_components(mpfi_float_backend& result, const V& a, const V& b) { number, et_on> x(a), y(b); assign_components(result, x.backend(), y.backend()); } // // Native non-member operations: // template inline void eval_sqrt(mpfi_float_backend& result, const mpfi_float_backend& val) { mpfi_sqrt(result.data(), val.data()); } template inline void eval_abs(mpfi_float_backend& result, const mpfi_float_backend& val) { mpfi_abs(result.data(), val.data()); } template inline void eval_fabs(mpfi_float_backend& result, const mpfi_float_backend& val) { mpfi_abs(result.data(), val.data()); } template inline void eval_ceil(mpfi_float_backend& result, const mpfi_float_backend& val) { mpfr_float_backend a, b; mpfr_set(a.data(), val.left_data(), GMP_RNDN); mpfr_set(b.data(), val.right_data(), GMP_RNDN); eval_ceil(a, a); eval_ceil(b, b); if (a.compare(b) != 0) { BOOST_THROW_EXCEPTION(interval_error("Attempt to take the ceil of a value that straddles an integer boundary.")); } mpfi_set_fr(result.data(), a.data()); } template inline void eval_floor(mpfi_float_backend& result, const mpfi_float_backend& val) { mpfr_float_backend a, b; mpfr_set(a.data(), val.left_data(), GMP_RNDN); mpfr_set(b.data(), val.right_data(), GMP_RNDN); eval_floor(a, a); eval_floor(b, b); if (a.compare(b) != 0) { BOOST_THROW_EXCEPTION(interval_error("Attempt to take the floor of a value that straddles an integer boundary.")); } mpfi_set_fr(result.data(), a.data()); } template inline void eval_ldexp(mpfi_float_backend& result, const mpfi_float_backend& val, long e) { if (e > 0) mpfi_mul_2exp(result.data(), val.data(), e); else if (e < 0) mpfi_div_2exp(result.data(), val.data(), -e); else result = val; } template inline void eval_frexp(mpfi_float_backend& result, const mpfi_float_backend& val, int* e) { mpfr_float_backend t, rt; mpfi_mid(t.data(), val.data()); eval_frexp(rt, t, e); eval_ldexp(result, val, -*e); } template inline void eval_frexp(mpfi_float_backend& result, const mpfi_float_backend& val, long* e) { mpfr_float_backend t, rt; mpfi_mid(t.data(), val.data()); eval_frexp(rt, t, e); eval_ldexp(result, val, -*e); } template inline int eval_fpclassify(const mpfi_float_backend& val) noexcept { return mpfi_inf_p(val.data()) ? FP_INFINITE : mpfi_nan_p(val.data()) ? FP_NAN : mpfi_is_zero(val.data()) ? FP_ZERO : FP_NORMAL; } template inline void eval_pow(mpfi_float_backend& result, const mpfi_float_backend& b, const mpfi_float_backend& e) { using ui_type = typename boost::multiprecision::detail::canonical >::type; using default_ops::eval_get_sign; int s = eval_get_sign(b); if (s == 0) { if (eval_get_sign(e) == 0) { result = ui_type(1); } else { result = ui_type(0); } return; } if (s < 0) { if (eval_get_sign(e) < 0) { mpfi_float_backend t1, t2; t1 = e; t1.negate(); eval_pow(t2, b, t1); t1 = ui_type(1); eval_divide(result, t1, t2); return; } typename boost::multiprecision::detail::canonical >::type an; #ifndef BOOST_NO_EXCEPTIONS try { #endif using default_ops::eval_convert_to; eval_convert_to(&an, e); if (e.compare(an) == 0) { mpfi_float_backend pb(b); pb.negate(); eval_pow(result, pb, e); if (an & 1u) result.negate(); return; } #ifndef BOOST_NO_EXCEPTIONS } catch (const std::exception&) { // conversion failed, just fall through, value is not an integer. } #endif result = std::numeric_limits, et_on> >::quiet_NaN().backend(); return; } mpfi_log(result.data(), b.data()); mpfi_mul(result.data(), result.data(), e.data()); mpfi_exp(result.data(), result.data()); } template inline void eval_exp(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_exp(result.data(), arg.data()); } template inline void eval_exp2(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_exp2(result.data(), arg.data()); } template inline void eval_log(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_log(result.data(), arg.data()); } template inline void eval_log10(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_log10(result.data(), arg.data()); } template inline void eval_sin(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_sin(result.data(), arg.data()); } template inline void eval_cos(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_cos(result.data(), arg.data()); } template inline void eval_tan(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_tan(result.data(), arg.data()); } template inline void eval_asin(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_asin(result.data(), arg.data()); } template inline void eval_acos(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_acos(result.data(), arg.data()); } template inline void eval_atan(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_atan(result.data(), arg.data()); } template inline void eval_atan2(mpfi_float_backend& result, const mpfi_float_backend& arg1, const mpfi_float_backend& arg2) { mpfi_atan2(result.data(), arg1.data(), arg2.data()); } template inline void eval_sinh(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_sinh(result.data(), arg.data()); } template inline void eval_cosh(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_cosh(result.data(), arg.data()); } template inline void eval_tanh(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_tanh(result.data(), arg.data()); } template inline void eval_log2(mpfi_float_backend& result, const mpfi_float_backend& arg) { mpfi_log2(result.data(), arg.data()); } template inline std::size_t hash_value(const mpfi_float_backend& val) { std::size_t result = 0; std::size_t len = val.left_data()[0]._mpfr_prec / mp_bits_per_limb; if (val.left_data()[0]._mpfr_prec % mp_bits_per_limb) ++len; for (std::size_t i = 0; i < len; ++i) boost::multiprecision::detail::hash_combine(result, val.left_data()[0]._mpfr_d[i]); boost::multiprecision::detail::hash_combine(result, val.left_data()[0]._mpfr_exp, val.left_data()[0]._mpfr_sign); len = val.right_data()[0]._mpfr_prec / mp_bits_per_limb; if (val.right_data()[0]._mpfr_prec % mp_bits_per_limb) ++len; for (std::size_t i = 0; i < len; ++i) boost::multiprecision::detail::hash_combine(result, val.right_data()[0]._mpfr_d[i]); boost::multiprecision::detail::hash_combine(result, val.right_data()[0]._mpfr_exp, val.right_data()[0]._mpfr_sign); return result; } template void generic_interconvert(To& to, const mpfi_float_backend& from, const std::integral_constant& to_type, const std::integral_constant& from_type) { using boost::multiprecision::detail::generic_interconvert; boost::multiprecision::detail::scoped_precision_options>> scoped(from); mpfr_float_backend t; mpfi_mid(t.data(), from.data()); generic_interconvert(to, t, to_type, from_type); } template void generic_interconvert(To& to, const mpfi_float_backend& from, const std::integral_constant& to_type, const std::integral_constant& from_type) { using boost::multiprecision::detail::generic_interconvert; boost::multiprecision::detail::scoped_precision_options>> scoped(from); mpfr_float_backend t; mpfi_mid(t.data(), from.data()); generic_interconvert(to, t, to_type, from_type); } template void generic_interconvert(To& to, const mpfi_float_backend& from, const std::integral_constant& to_type, const std::integral_constant& from_type) { using boost::multiprecision::detail::generic_interconvert; boost::multiprecision::detail::scoped_precision_options>> scoped(from); mpfr_float_backend t; mpfi_mid(t.data(), from.data()); generic_interconvert(to, t, to_type, from_type); } } // namespace backends template <> struct number_category >::type> : public std::integral_constant {}; template struct is_interval_number > : public std::integral_constant {}; using boost::multiprecision::backends::mpfi_float_backend; using mpfi_float_50 = number > ; using mpfi_float_100 = number > ; using mpfi_float_500 = number > ; using mpfi_float_1000 = number >; using mpfi_float = number > ; // // Special interval specific functions: // template inline number, ExpressionTemplates> lower(const number, ExpressionTemplates>& val) { boost::multiprecision::detail::scoped_precision_options, ExpressionTemplates> > precision_guard(val); number > result; mpfr_set(result.backend().data(), val.backend().left_data(), GMP_RNDN); return result; } template inline number, ExpressionTemplates> upper(const number, ExpressionTemplates>& val) { boost::multiprecision::detail::scoped_precision_options, ExpressionTemplates> > precision_guard(val); number > result; mpfr_set(result.backend().data(), val.backend().right_data(), GMP_RNDN); return result; } template inline number, ExpressionTemplates> median(const number, ExpressionTemplates>& val) { boost::multiprecision::detail::scoped_precision_options, ExpressionTemplates> > precision_guard(val); number > result; mpfi_mid(result.backend().data(), val.backend().data()); return result; } template inline number, ExpressionTemplates> width(const number, ExpressionTemplates>& val) { boost::multiprecision::detail::scoped_precision_options, ExpressionTemplates> > precision_guard(val); number > result; mpfi_diam_abs(result.backend().data(), val.backend().data()); return result; } template inline number, ExpressionTemplates> intersect(const number, ExpressionTemplates>& a, const number, ExpressionTemplates>& b) { boost::multiprecision::detail::scoped_default_precision, ExpressionTemplates> > precision_guard(a, b); number, ExpressionTemplates> result; mpfi_intersect(result.backend().data(), a.backend().data(), b.backend().data()); return result; } template inline number, ExpressionTemplates> hull(const number, ExpressionTemplates>& a, const number, ExpressionTemplates>& b) { boost::multiprecision::detail::scoped_default_precision, ExpressionTemplates> > precision_guard(a, b); number, ExpressionTemplates> result; mpfi_union(result.backend().data(), a.backend().data(), b.backend().data()); return result; } template inline bool overlap(const number, ExpressionTemplates>& a, const number, ExpressionTemplates>& b) { return (lower(a) <= lower(b) && lower(b) <= upper(a)) || (lower(b) <= lower(a) && lower(a) <= upper(b)); } template inline bool in(const number, ExpressionTemplates1>& a, const number, ExpressionTemplates2>& b) { return mpfi_is_inside_fr(a.backend().data(), b.backend().data()) != 0; } template inline bool zero_in(const number, ExpressionTemplates>& a) { return mpfi_has_zero(a.backend().data()) != 0; } template