/////////////////////////////////////////////////////////////////////////////// // 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 namespace boost{ namespace multiprecision{ namespace backends{ template struct mpfi_float_backend; } // namespace backends template struct number_category > : public mpl::int_{}; struct interval_error : public std::runtime_error { interval_error(const std::string& s) : std::runtime_error(s) {} }; 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 { typedef mpl::list signed_types; typedef mpl::list unsigned_types; typedef mpl::list float_types; typedef long exponent_type; mpfi_float_imp() { mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : get_default_precision())); } mpfi_float_imp(unsigned prec) { mpfi_init2(m_data, prec); } mpfi_float_imp(const mpfi_float_imp& o) { mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : get_default_precision())); if(o.m_data[0].left._mpfr_d) mpfi_set(m_data, o.m_data); } #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES mpfi_float_imp(mpfi_float_imp&& o) BOOST_NOEXCEPT { m_data[0] = o.m_data[0]; o.m_data[0].left._mpfr_d = 0; } #endif mpfi_float_imp& operator = (const mpfi_float_imp& o) { if(m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : get_default_precision())); if(o.m_data[0].left._mpfr_d) mpfi_set(m_data, o.m_data); return *this; } #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES mpfi_float_imp& operator = (mpfi_float_imp&& o) BOOST_NOEXCEPT { mpfi_swap(m_data, o.m_data); return *this; } #endif #ifdef _MPFR_H_HAVE_INTMAX_T mpfi_float_imp& operator = (unsigned long long i) { if(m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : 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 = (long long i) { if(m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : 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 = (unsigned long long i) { if(m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : get_default_precision())); unsigned long long mask = ((1uLL << std::numeric_limits::digits) - 1); 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 = (long long i) { if(m_data[0].left._mpfr_d == 0) mpfi_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : 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 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 : 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 : 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 : 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 : 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 : 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) BOOST_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() BOOST_NOEXCEPT { if(m_data[0].left._mpfr_d) mpfi_clear(m_data); } void negate() BOOST_NOEXCEPT { BOOST_ASSERT(m_data[0].left._mpfr_d); mpfi_neg(m_data, m_data); } int compare(const mpfi_float_imp& o)const BOOST_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 BOOST_NOEXCEPT { mpfi_float_imp d; d = v; return compare(d); } mpfi_t& data() BOOST_NOEXCEPT { BOOST_ASSERT(m_data[0].left._mpfr_d); return m_data; } const mpfi_t& data()const BOOST_NOEXCEPT { BOOST_ASSERT(m_data[0].left._mpfr_d); return m_data; } mpfr_ptr left_data() BOOST_NOEXCEPT { BOOST_ASSERT(m_data[0].left._mpfr_d); return &(m_data[0].left); } mpfr_srcptr left_data()const BOOST_NOEXCEPT { BOOST_ASSERT(m_data[0].left._mpfr_d); return &(m_data[0].left); } mpfr_ptr right_data() BOOST_NOEXCEPT { BOOST_ASSERT(m_data[0].left._mpfr_d); return &(m_data[0].right); } mpfr_srcptr right_data()const BOOST_NOEXCEPT { BOOST_ASSERT(m_data[0].left._mpfr_d); return &(m_data[0].right); } protected: mpfi_t m_data; static unsigned& get_default_precision() BOOST_NOEXCEPT { static unsigned val = 50; return val; } }; } // 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) {} #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES mpfi_float_backend(mpfi_float_backend&& o) : detail::mpfi_float_imp(static_cast&&>(o)) {} #endif template mpfi_float_backend(const mpfi_float_backend& val, typename enable_if_c::type* = 0) : detail::mpfi_float_imp() { mpfi_set(this->m_data, val.data()); } template explicit mpfi_float_backend(const mpfi_float_backend& val, typename disable_if_c::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 enable_if_c::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 disable_if_c::type* = 0) : detail::mpfi_float_imp() { mpfi_set_fr(this->m_data, val.data()); } #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES mpfi_float_backend& operator=(mpfi_float_backend&& o) BOOST_NOEXCEPT { *static_cast*>(this) = static_cast&&>(o); return *this; } #endif 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 <> 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) {} #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES mpfi_float_backend(mpfi_float_backend&& o) BOOST_NOEXCEPT : detail::mpfi_float_imp<0>(static_cast&&>(o)) {} #endif mpfi_float_backend(const mpfi_float_backend& o, unsigned digits10) : detail::mpfi_float_imp<0>(digits10) { *this = o; } 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) { mpfi_set_prec(this->m_data, mpfi_get_prec(o.data())); mpfi_set(this->m_data, o.data()); return *this; } #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES mpfi_float_backend& operator=(mpfi_float_backend&& o) BOOST_NOEXCEPT { *static_cast*>(this) = static_cast &&>(o); return *this; } #endif template mpfi_float_backend& operator=(const V& v) { *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 default_precision() BOOST_NOEXCEPT { return get_default_precision(); } static void default_precision(unsigned v) BOOST_NOEXCEPT { get_default_precision() = v; } unsigned precision()const BOOST_NOEXCEPT { return multiprecision::detail::digits2_2_10(mpfi_get_prec(this->m_data)); } void precision(unsigned digits10) BOOST_NOEXCEPT { mpfi_set_prec(this->m_data, multiprecision::detail::digits2_2_10((digits10))); } }; template inline typename enable_if, bool>::type eval_eq(const mpfi_float_backend& a, const T& b) BOOST_NOEXCEPT { return a.compare(b) == 0; } template inline typename enable_if, bool>::type eval_lt(const mpfi_float_backend& a, const T& b) BOOST_NOEXCEPT { return a.compare(b) < 0; } template inline typename enable_if, bool>::type eval_gt(const mpfi_float_backend& a, const T& b) BOOST_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) BOOST_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(unsigned long 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 long* 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) BOOST_NOEXCEPT { *result = mpfi_get_d(val.data()); } template inline void eval_convert_to(long double* result, const mpfi_float_backend& val) BOOST_NOEXCEPT { mpfr_float_backend t; mpfi_mid(t.data(), val.data()); eval_convert_to(result, t); } template inline void assign_components(mpfi_float_backend& result, const mpfr_float_backend& a, const mpfr_float_backend& 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 enable_if_c, et_on> >::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) BOOST_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) { typedef typename boost::multiprecision::detail::canonical >::type ui_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; try { 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; } } catch(const std::exception&) { // conversion failed, just fall through, value is not an integer. } 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_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()); } } // namespace backends #ifdef BOOST_NO_SFINAE_EXPR namespace detail{ template struct is_explicitly_convertible, backends::mpfi_float_backend > : public mpl::true_ {}; } #endif template<> struct number_category >::type> : public mpl::int_{}; template struct is_interval_number > : public mpl::true_ {}; using boost::multiprecision::backends::mpfi_float_backend; typedef number > mpfi_float_50; typedef number > mpfi_float_100; typedef number > mpfi_float_500; typedef number > mpfi_float_1000; typedef number > mpfi_float; // // Special interval specific functions: // template inline number, ExpressionTemplates> lower(const number, ExpressionTemplates>& 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) { number > result; mpfr_set(result.backend().data(), val.backend().right_data(), GMP_RNDN); return result; } template inline number, ExpressionTemplates> median(const number, ExpressionTemplates>& val) { number > result; mpfi_mid(result.backend().data(), val.backend().data()); return result; } template inline number, ExpressionTemplates> width(const number, ExpressionTemplates>& 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) { 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) { 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 inline bool subset(const number, ExpressionTemplates>& a, const number, ExpressionTemplates>& b) { return mpfi_is_inside(a.backend().data(), b.backend().data()) != 0; } template inline bool proper_subset(const number, ExpressionTemplates>& a, const number, ExpressionTemplates>& b) { return mpfi_is_strictly_inside(a.backend().data(), b.backend().data()) != 0; } template inline bool empty(const number, ExpressionTemplates>& a) { return mpfi_is_empty(a.backend().data()) != 0; } template inline bool singleton(const number, ExpressionTemplates>& a) { return mpfr_cmp(a.backend().left_data(), a.backend().right_data()) == 0; } template struct component_type, ExpressionTemplates> > { typedef number, ExpressionTemplates> type; }; } // namespace multiprecision namespace math{ namespace tools{ template <> inline int digits() { return boost::multiprecision::backends::detail::get_default_precision(); } template <> inline int digits, boost::multiprecision::et_off> >() { return boost::multiprecision::backends::detail::get_default_precision(); } } // namespace tools namespace constants{ namespace detail{ template struct constant_pi; template struct constant_ln_two; template struct constant_euler; template struct constant_catalan; // // Initializer: ensure all our constants are initialized prior to the first call of main: // template struct mpfi_initializer { struct init { init() { boost::math::constants::pi(); boost::math::constants::ln_two(); boost::math::constants::euler(); boost::math::constants::catalan(); } void force_instantiate()const{} }; static const init initializer; static void force_instantiate() { initializer.force_instantiate(); } }; template const typename mpfi_initializer::init mpfi_initializer::initializer; template struct constant_pi, ExpressionTemplates> > { typedef boost::multiprecision::number, ExpressionTemplates> result_type; template static inline result_type const& get(const mpl::int_&) { mpfi_initializer::force_instantiate(); static result_type result; static bool has_init = false; if(!has_init) { has_init = true; mpfi_const_pi(result.backend().data()); } return result; } }; template struct constant_ln_two, ExpressionTemplates> > { typedef boost::multiprecision::number, ExpressionTemplates> result_type; template static inline result_type const& get(const mpl::int_&) { mpfi_initializer::force_instantiate(); static result_type result; static bool has_init = false; if(!has_init) { has_init = true; mpfi_const_log2(result.backend().data()); } return result; } }; template struct constant_euler, ExpressionTemplates> > { typedef boost::multiprecision::number, ExpressionTemplates> result_type; template static inline result_type const& get(const mpl::int_&) { mpfi_initializer::force_instantiate(); static result_type result; static bool has_init = false; if(!has_init) { has_init = true; mpfi_const_euler(result.backend().data()); } return result; } }; template struct constant_catalan, ExpressionTemplates> > { typedef boost::multiprecision::number, ExpressionTemplates> result_type; template static inline result_type const& get(const mpl::int_&) { mpfi_initializer::force_instantiate(); static result_type result; static bool has_init = false; if(!has_init) { has_init = true; mpfi_const_catalan(result.backend().data()); } return result; } }; }} // namespaces }} // namespaces namespace std{ // // numeric_limits [partial] specializations for the types declared in this header: // template class numeric_limits, ExpressionTemplates> > { typedef boost::multiprecision::number, ExpressionTemplates> number_type; public: BOOST_STATIC_CONSTEXPR bool is_specialized = true; static number_type (min)() { initializer.do_nothing(); static std::pair value; if(!value.first) { value.first = true; value.second = 0.5; mpfi_div_2exp(value.second.backend().data(), value.second.backend().data(), -mpfr_get_emin()); } return value.second; } static number_type (max)() { initializer.do_nothing(); static std::pair value; if(!value.first) { value.first = true; value.second = 0.5; mpfi_mul_2exp(value.second.backend().data(), value.second.backend().data(), mpfr_get_emax()); } return value.second; } BOOST_STATIC_CONSTEXPR number_type lowest() { return -(max)(); } BOOST_STATIC_CONSTEXPR int digits = static_cast((Digits10 * 1000L) / 301L + ((Digits10 * 1000L) % 301 ? 2 : 1)); BOOST_STATIC_CONSTEXPR int digits10 = Digits10; // Is this really correct??? BOOST_STATIC_CONSTEXPR int max_digits10 = Digits10 + 2; BOOST_STATIC_CONSTEXPR bool is_signed = true; BOOST_STATIC_CONSTEXPR bool is_integer = false; BOOST_STATIC_CONSTEXPR bool is_exact = false; BOOST_STATIC_CONSTEXPR int radix = 2; static number_type epsilon() { initializer.do_nothing(); static std::pair value; if(!value.first) { value.first = true; value.second = 1; mpfi_div_2exp(value.second.backend().data(), value.second.backend().data(), std::numeric_limits::digits - 1); } return value.second; } // What value should this be???? static number_type round_error() { // returns epsilon/2 initializer.do_nothing(); static std::pair value; if(!value.first) { value.first = true; value.second = 1; mpfi_div_2exp(value.second.backend().data(), value.second.backend().data(), 1); } return value.second; } BOOST_STATIC_CONSTEXPR long min_exponent = MPFR_EMIN_DEFAULT; BOOST_STATIC_CONSTEXPR long min_exponent10 = (MPFR_EMIN_DEFAULT / 1000) * 301L; BOOST_STATIC_CONSTEXPR long max_exponent = MPFR_EMAX_DEFAULT; BOOST_STATIC_CONSTEXPR long max_exponent10 = (MPFR_EMAX_DEFAULT / 1000) * 301L; BOOST_STATIC_CONSTEXPR bool has_infinity = true; BOOST_STATIC_CONSTEXPR bool has_quiet_NaN = true; BOOST_STATIC_CONSTEXPR bool has_signaling_NaN = false; BOOST_STATIC_CONSTEXPR float_denorm_style has_denorm = denorm_absent; BOOST_STATIC_CONSTEXPR bool has_denorm_loss = false; static number_type infinity() { initializer.do_nothing(); static std::pair value; if(!value.first) { boost::multiprecision::mpfr_float_backend t; mpfr_set_inf(t.data(), 1); value.first = true; mpfi_set_fr(value.second.backend().data(), t.data()); } return value.second; } static number_type quiet_NaN() { initializer.do_nothing(); static std::pair value; if(!value.first) { boost::multiprecision::mpfr_float_backend t; mpfr_set_nan(t.data()); value.first = true; mpfi_set_fr(value.second.backend().data(), t.data()); } return value.second; } BOOST_STATIC_CONSTEXPR number_type signaling_NaN() { return number_type(0); } BOOST_STATIC_CONSTEXPR number_type denorm_min() { return number_type(0); } BOOST_STATIC_CONSTEXPR bool is_iec559 = false; BOOST_STATIC_CONSTEXPR bool is_bounded = true; BOOST_STATIC_CONSTEXPR bool is_modulo = false; BOOST_STATIC_CONSTEXPR bool traps = true; BOOST_STATIC_CONSTEXPR bool tinyness_before = false; BOOST_STATIC_CONSTEXPR float_round_style round_style = round_to_nearest; private: struct data_initializer { data_initializer() { std::numeric_limits > >::epsilon(); std::numeric_limits > >::round_error(); (std::numeric_limits > >::min)(); (std::numeric_limits > >::max)(); std::numeric_limits > >::infinity(); std::numeric_limits > >::quiet_NaN(); } void do_nothing()const{} }; static const data_initializer initializer; }; template const typename numeric_limits, ExpressionTemplates> >::data_initializer numeric_limits, ExpressionTemplates> >::initializer; #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION template BOOST_CONSTEXPR_OR_CONST int numeric_limits