1 /////////////////////////////////////////////////////////////// 2 // Copyright 2013 John Maddock. Distributed under the Boost 3 // Software License, Version 1.0. (See accompanying file 4 // LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt 5 6 #ifndef BOOST_MATH_CPP_BIN_FLOAT_HPP 7 #define BOOST_MATH_CPP_BIN_FLOAT_HPP 8 9 #include <boost/multiprecision/cpp_int.hpp> 10 #include <boost/multiprecision/integer.hpp> 11 #include <boost/math/special_functions/trunc.hpp> 12 #include <boost/multiprecision/detail/float_string_cvt.hpp> 13 #include <boost/multiprecision/traits/max_digits10.hpp> 14 15 // 16 // Some includes we need from Boost.Math, since we rely on that library to provide these functions: 17 // 18 #include <boost/math/special_functions/asinh.hpp> 19 #include <boost/math/special_functions/acosh.hpp> 20 #include <boost/math/special_functions/atanh.hpp> 21 #include <boost/math/special_functions/cbrt.hpp> 22 #include <boost/math/special_functions/expm1.hpp> 23 #include <boost/math/special_functions/gamma.hpp> 24 25 #ifdef BOOST_HAS_FLOAT128 26 #include <quadmath.h> 27 #endif 28 29 namespace boost { 30 namespace multiprecision { 31 namespace backends { 32 33 enum digit_base_type 34 { 35 digit_base_2 = 2, 36 digit_base_10 = 10 37 }; 38 39 #ifdef BOOST_MSVC 40 #pragma warning(push) 41 #pragma warning(disable : 4522 6326) // multiple assignment operators specified, comparison of two constants 42 #endif 43 44 namespace detail { 45 46 template <class U> 47 inline typename enable_if_c<is_unsigned<U>::value, bool>::type is_negative(U) { return false; } 48 template <class S> 49 inline typename disable_if_c<is_unsigned<S>::value, bool>::type is_negative(S s) { return s < 0; } 50 51 template <class Float, int, bool = number_category<Float>::value == number_kind_floating_point> 52 struct is_cpp_bin_float_implicitly_constructible_from_type 53 { 54 static const bool value = false; 55 }; 56 57 template <class Float, int bit_count> 58 struct is_cpp_bin_float_implicitly_constructible_from_type<Float, bit_count, true> 59 { 60 static const bool value = (std::numeric_limits<Float>::digits <= (int)bit_count) && (std::numeric_limits<Float>::radix == 2) && std::numeric_limits<Float>::is_specialized 61 #ifdef BOOST_HAS_FLOAT128 62 && !boost::is_same<Float, __float128>::value 63 #endif 64 && (is_floating_point<Float>::value || is_number<Float>::value); 65 }; 66 67 template <class Float, int, bool = number_category<Float>::value == number_kind_floating_point> 68 struct is_cpp_bin_float_explicitly_constructible_from_type 69 { 70 static const bool value = false; 71 }; 72 73 template <class Float, int bit_count> 74 struct is_cpp_bin_float_explicitly_constructible_from_type<Float, bit_count, true> 75 { 76 static const bool value = (std::numeric_limits<Float>::digits > (int)bit_count) && (std::numeric_limits<Float>::radix == 2) && std::numeric_limits<Float>::is_specialized 77 #ifdef BOOST_HAS_FLOAT128 78 && !boost::is_same<Float, __float128>::value 79 #endif 80 ; 81 }; 82 83 } // namespace detail 84 85 template <unsigned Digits, digit_base_type DigitBase = digit_base_10, class Allocator = void, class Exponent = int, Exponent MinExponent = 0, Exponent MaxExponent = 0> 86 class cpp_bin_float 87 { 88 public: 89 static const unsigned bit_count = DigitBase == digit_base_2 ? Digits : (Digits * 1000uL) / 301uL + (((Digits * 1000uL) % 301) ? 2u : 1u); 90 typedef cpp_int_backend<is_void<Allocator>::value ? bit_count : 0, bit_count, is_void<Allocator>::value ? unsigned_magnitude : signed_magnitude, unchecked, Allocator> rep_type; 91 typedef cpp_int_backend<is_void<Allocator>::value ? 2 * bit_count : 0, 2 * bit_count, is_void<Allocator>::value ? unsigned_magnitude : signed_magnitude, unchecked, Allocator> double_rep_type; 92 93 typedef typename rep_type::signed_types signed_types; 94 typedef typename rep_type::unsigned_types unsigned_types; 95 typedef boost::mpl::list<float, double, long double> float_types; 96 typedef Exponent exponent_type; 97 98 static const exponent_type max_exponent_limit = boost::integer_traits<exponent_type>::const_max - 2 * static_cast<exponent_type>(bit_count); 99 static const exponent_type min_exponent_limit = boost::integer_traits<exponent_type>::const_min + 2 * static_cast<exponent_type>(bit_count); 100 101 BOOST_STATIC_ASSERT_MSG(MinExponent >= min_exponent_limit, "Template parameter MinExponent is too negative for our internal logic to function correctly, sorry!"); 102 BOOST_STATIC_ASSERT_MSG(MaxExponent <= max_exponent_limit, "Template parameter MaxExponent is too large for our internal logic to function correctly, sorry!"); 103 BOOST_STATIC_ASSERT_MSG(MinExponent <= 0, "Template parameter MinExponent can not be positive!"); 104 BOOST_STATIC_ASSERT_MSG(MaxExponent >= 0, "Template parameter MaxExponent can not be negative!"); 105 106 static const exponent_type max_exponent = MaxExponent == 0 ? max_exponent_limit : MaxExponent; 107 static const exponent_type min_exponent = MinExponent == 0 ? min_exponent_limit : MinExponent; 108 109 static const exponent_type exponent_zero = max_exponent + 1; 110 static const exponent_type exponent_infinity = max_exponent + 2; 111 static const exponent_type exponent_nan = max_exponent + 3; 112 113 private: 114 rep_type m_data; 115 exponent_type m_exponent; 116 bool m_sign; 117 118 public: 119 cpp_bin_float() BOOST_MP_NOEXCEPT_IF(noexcept(rep_type())) : m_data(), m_exponent(exponent_zero), m_sign(false) {} 120 121 cpp_bin_float(const cpp_bin_float& o) BOOST_MP_NOEXCEPT_IF(noexcept(rep_type(std::declval<const rep_type&>()))) 122 : m_data(o.m_data), m_exponent(o.m_exponent), m_sign(o.m_sign) {} 123 124 template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE> 125 cpp_bin_float(const cpp_bin_float<D, B, A, E, MinE, MaxE>& o, typename boost::enable_if_c<(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = 0) 126 { 127 *this = o; 128 } 129 template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE> 130 explicit cpp_bin_float(const cpp_bin_float<D, B, A, E, MinE, MaxE>& o, typename boost::disable_if_c<(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = 0) 131 : m_exponent(o.exponent()), m_sign(o.sign()) 132 { 133 *this = o; 134 } 135 template <class Float> 136 cpp_bin_float(const Float& f, 137 typename boost::enable_if_c<detail::is_cpp_bin_float_implicitly_constructible_from_type<Float, bit_count>::value>::type const* = 0) 138 : m_data(), m_exponent(0), m_sign(false) 139 { 140 this->assign_float(f); 141 } 142 143 template <class Float> 144 explicit cpp_bin_float(const Float& f, 145 typename boost::enable_if_c<detail::is_cpp_bin_float_explicitly_constructible_from_type<Float, bit_count>::value>::type const* = 0) 146 : m_data(), m_exponent(0), m_sign(false) 147 { 148 this->assign_float(f); 149 } 150 #ifdef BOOST_HAS_FLOAT128 151 template <class Float> 152 cpp_bin_float(const Float& f, 153 typename boost::enable_if_c< 154 boost::is_same<Float, __float128>::value && ((int)bit_count >= 113)>::type const* = 0) 155 : m_data(), m_exponent(0), m_sign(false) 156 { 157 this->assign_float(f); 158 } 159 template <class Float> 160 explicit cpp_bin_float(const Float& f, 161 typename boost::enable_if_c< 162 boost::is_same<Float, __float128>::value && ((int)bit_count < 113)>::type const* = 0) 163 : m_data(), m_exponent(0), m_sign(false) 164 { 165 this->assign_float(f); 166 } 167 #endif 168 cpp_bin_float& operator=(const cpp_bin_float& o) BOOST_MP_NOEXCEPT_IF(noexcept(std::declval<rep_type&>() = std::declval<const rep_type&>())) 169 { 170 m_data = o.m_data; 171 m_exponent = o.m_exponent; 172 m_sign = o.m_sign; 173 return *this; 174 } 175 176 template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE> 177 cpp_bin_float& operator=(const cpp_bin_float<D, B, A, E, MinE, MaxE>& f) 178 { 179 switch (eval_fpclassify(f)) 180 { 181 case FP_ZERO: 182 m_data = limb_type(0); 183 m_sign = f.sign(); 184 m_exponent = exponent_zero; 185 break; 186 case FP_NAN: 187 m_data = limb_type(0); 188 m_sign = false; 189 m_exponent = exponent_nan; 190 break; 191 ; 192 case FP_INFINITE: 193 m_data = limb_type(0); 194 m_sign = f.sign(); 195 m_exponent = exponent_infinity; 196 break; 197 default: 198 typename cpp_bin_float<D, B, A, E, MinE, MaxE>::rep_type b(f.bits()); 199 this->exponent() = f.exponent() + (E)bit_count - (E)cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count; 200 this->sign() = f.sign(); 201 copy_and_round(*this, b); 202 } 203 return *this; 204 } 205 #ifdef BOOST_HAS_FLOAT128 206 template <class Float> 207 typename boost::enable_if_c< 208 (number_category<Float>::value == number_kind_floating_point) 209 //&& (std::numeric_limits<Float>::digits <= (int)bit_count) 210 && ((std::numeric_limits<Float>::radix == 2) || (boost::is_same<Float, __float128>::value)), 211 cpp_bin_float&>::type 212 operator=(const Float& f) 213 #else 214 template <class Float> 215 typename boost::enable_if_c< 216 (number_category<Float>::value == number_kind_floating_point) 217 //&& (std::numeric_limits<Float>::digits <= (int)bit_count) 218 && (std::numeric_limits<Float>::radix == 2), 219 cpp_bin_float&>::type 220 operator=(const Float& f) 221 #endif 222 { 223 return assign_float(f); 224 } 225 226 #ifdef BOOST_HAS_FLOAT128 227 template <class Float> 228 typename boost::enable_if_c<boost::is_same<Float, __float128>::value, cpp_bin_float&>::type assign_float(Float f) 229 { 230 using default_ops::eval_add; 231 typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type; 232 if (f == 0) 233 { 234 m_data = limb_type(0); 235 m_sign = (signbitq(f) > 0); 236 m_exponent = exponent_zero; 237 return *this; 238 } 239 else if (isnanq(f)) 240 { 241 m_data = limb_type(0); 242 m_sign = false; 243 m_exponent = exponent_nan; 244 return *this; 245 } 246 else if (isinfq(f)) 247 { 248 m_data = limb_type(0); 249 m_sign = (f < 0); 250 m_exponent = exponent_infinity; 251 return *this; 252 } 253 if (f < 0) 254 { 255 *this = -f; 256 this->negate(); 257 return *this; 258 } 259 260 typedef typename mpl::front<unsigned_types>::type ui_type; 261 m_data = static_cast<ui_type>(0u); 262 m_sign = false; 263 m_exponent = 0; 264 265 static const int bits = sizeof(int) * CHAR_BIT - 1; 266 int e; 267 f = frexpq(f, &e); 268 while (f) 269 { 270 f = ldexpq(f, bits); 271 e -= bits; 272 int ipart = (int)truncq(f); 273 f -= ipart; 274 m_exponent += bits; 275 cpp_bin_float t; 276 t = static_cast<bf_int_type>(ipart); 277 eval_add(*this, t); 278 } 279 m_exponent += static_cast<Exponent>(e); 280 return *this; 281 } 282 #endif 283 #ifdef BOOST_HAS_FLOAT128 284 template <class Float> 285 typename boost::enable_if_c<is_floating_point<Float>::value && !is_same<Float, __float128>::value, cpp_bin_float&>::type assign_float(Float f) 286 #else 287 template <class Float> 288 typename boost::enable_if_c<is_floating_point<Float>::value, cpp_bin_float&>::type assign_float(Float f) 289 #endif 290 { 291 BOOST_MATH_STD_USING 292 using default_ops::eval_add; 293 typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type; 294 295 switch ((boost::math::fpclassify)(f)) 296 { 297 case FP_ZERO: 298 m_data = limb_type(0); 299 m_sign = ((boost::math::signbit)(f) > 0); 300 m_exponent = exponent_zero; 301 return *this; 302 case FP_NAN: 303 m_data = limb_type(0); 304 m_sign = false; 305 m_exponent = exponent_nan; 306 return *this; 307 case FP_INFINITE: 308 m_data = limb_type(0); 309 m_sign = (f < 0); 310 m_exponent = exponent_infinity; 311 return *this; 312 } 313 if (f < 0) 314 { 315 *this = -f; 316 this->negate(); 317 return *this; 318 } 319 320 typedef typename mpl::front<unsigned_types>::type ui_type; 321 m_data = static_cast<ui_type>(0u); 322 m_sign = false; 323 m_exponent = 0; 324 325 static const int bits = sizeof(int) * CHAR_BIT - 1; 326 int e; 327 f = frexp(f, &e); 328 while (f) 329 { 330 f = ldexp(f, bits); 331 e -= bits; 332 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS 333 int ipart = itrunc(f); 334 #else 335 int ipart = static_cast<int>(f); 336 #endif 337 f -= ipart; 338 m_exponent += bits; 339 cpp_bin_float t; 340 t = static_cast<bf_int_type>(ipart); 341 eval_add(*this, t); 342 } 343 m_exponent += static_cast<Exponent>(e); 344 return *this; 345 } 346 347 template <class Float> 348 typename boost::enable_if_c< 349 (number_category<Float>::value == number_kind_floating_point) && !boost::is_floating_point<Float>::value && (number_category<Float>::value == number_kind_floating_point), 350 cpp_bin_float&>::type 351 assign_float(Float f) 352 { 353 BOOST_MATH_STD_USING 354 using default_ops::eval_add; 355 using default_ops::eval_convert_to; 356 using default_ops::eval_get_sign; 357 using default_ops::eval_subtract; 358 359 typedef typename boost::multiprecision::detail::canonical<int, Float>::type f_int_type; 360 typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type; 361 362 switch (eval_fpclassify(f)) 363 { 364 case FP_ZERO: 365 m_data = limb_type(0); 366 m_sign = (eval_get_sign(f) > 0); 367 m_exponent = exponent_zero; 368 return *this; 369 case FP_NAN: 370 m_data = limb_type(0); 371 m_sign = false; 372 m_exponent = exponent_nan; 373 return *this; 374 case FP_INFINITE: 375 m_data = limb_type(0); 376 m_sign = eval_get_sign(f) < 0; 377 m_exponent = exponent_infinity; 378 return *this; 379 } 380 if (eval_get_sign(f) < 0) 381 { 382 f.negate(); 383 assign_float(f); 384 this->negate(); 385 return *this; 386 } 387 388 typedef typename mpl::front<unsigned_types>::type ui_type; 389 m_data = static_cast<ui_type>(0u); 390 m_sign = false; 391 m_exponent = 0; 392 393 static const int bits = sizeof(int) * CHAR_BIT - 1; 394 int e; 395 eval_frexp(f, f, &e); 396 while (eval_get_sign(f) != 0) 397 { 398 eval_ldexp(f, f, bits); 399 e -= bits; 400 int ipart; 401 eval_convert_to(&ipart, f); 402 eval_subtract(f, static_cast<f_int_type>(ipart)); 403 m_exponent += bits; 404 eval_add(*this, static_cast<bf_int_type>(ipart)); 405 } 406 m_exponent += e; 407 if (m_exponent > max_exponent) 408 m_exponent = exponent_infinity; 409 if (m_exponent < min_exponent) 410 { 411 m_data = limb_type(0u); 412 m_exponent = exponent_zero; 413 m_sign = (eval_get_sign(f) > 0); 414 } 415 else if (eval_get_sign(m_data) == 0) 416 { 417 m_exponent = exponent_zero; 418 m_sign = (eval_get_sign(f) > 0); 419 } 420 return *this; 421 } 422 template <class B, expression_template_option et> 423 cpp_bin_float& assign_float(const number<B, et>& f) 424 { 425 return assign_float(f.backend()); 426 } 427 428 template <class I> 429 typename boost::enable_if<is_integral<I>, cpp_bin_float&>::type operator=(const I& i) 430 { 431 using default_ops::eval_bit_test; 432 if (!i) 433 { 434 m_data = static_cast<limb_type>(0); 435 m_exponent = exponent_zero; 436 m_sign = false; 437 } 438 else 439 { 440 typedef typename make_unsigned<I>::type ui_type; 441 ui_type fi = static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(i)); 442 typedef typename boost::multiprecision::detail::canonical<ui_type, rep_type>::type ar_type; 443 m_data = static_cast<ar_type>(fi); 444 unsigned shift = msb(fi); 445 if (shift >= bit_count) 446 { 447 m_exponent = static_cast<Exponent>(shift); 448 m_data = static_cast<ar_type>(fi >> (shift + 1 - bit_count)); 449 } 450 else 451 { 452 m_exponent = static_cast<Exponent>(shift); 453 eval_left_shift(m_data, bit_count - shift - 1); 454 } 455 BOOST_ASSERT(eval_bit_test(m_data, bit_count - 1)); 456 m_sign = detail::is_negative(i); 457 } 458 return *this; 459 } 460 461 cpp_bin_float& operator=(const char* s); 462 463 void swap(cpp_bin_float& o) BOOST_NOEXCEPT 464 { 465 m_data.swap(o.m_data); 466 std::swap(m_exponent, o.m_exponent); 467 std::swap(m_sign, o.m_sign); 468 } 469 470 std::string str(std::streamsize dig, std::ios_base::fmtflags f) const; 471 472 void negate() 473 { 474 if (m_exponent != exponent_nan) 475 m_sign = !m_sign; 476 } 477 478 int compare(const cpp_bin_float& o) const BOOST_NOEXCEPT 479 { 480 if (m_sign != o.m_sign) 481 return (m_exponent == exponent_zero) && (m_exponent == o.m_exponent) ? 0 : m_sign ? -1 : 1; 482 int result; 483 if (m_exponent == exponent_nan) 484 return -1; 485 else if (m_exponent != o.m_exponent) 486 { 487 if (m_exponent == exponent_zero) 488 result = -1; 489 else if (o.m_exponent == exponent_zero) 490 result = 1; 491 else 492 result = m_exponent > o.m_exponent ? 1 : -1; 493 } 494 else 495 result = m_data.compare(o.m_data); 496 if (m_sign) 497 result = -result; 498 return result; 499 } 500 template <class A> 501 int compare(const A& o) const BOOST_NOEXCEPT 502 { 503 cpp_bin_float b; 504 b = o; 505 return compare(b); 506 } 507 508 rep_type& bits() { return m_data; } 509 const rep_type& bits() const { return m_data; } 510 exponent_type& exponent() { return m_exponent; } 511 const exponent_type& exponent() const { return m_exponent; } 512 bool& sign() { return m_sign; } 513 const bool& sign() const { return m_sign; } 514 void check_invariants() 515 { 516 using default_ops::eval_bit_test; 517 using default_ops::eval_is_zero; 518 if ((m_exponent <= max_exponent) && (m_exponent >= min_exponent)) 519 { 520 BOOST_ASSERT(eval_bit_test(m_data, bit_count - 1)); 521 } 522 else 523 { 524 BOOST_ASSERT(m_exponent > max_exponent); 525 BOOST_ASSERT(m_exponent <= exponent_nan); 526 BOOST_ASSERT(eval_is_zero(m_data)); 527 } 528 } 529 template <class Archive> 530 void serialize(Archive& ar, const unsigned int /*version*/) 531 { 532 ar& boost::make_nvp("data", m_data); 533 ar& boost::make_nvp("exponent", m_exponent); 534 ar& boost::make_nvp("sign", m_sign); 535 } 536 }; 537 538 #ifdef BOOST_MSVC 539 #pragma warning(pop) 540 #endif 541 542 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class Int> 543 inline void copy_and_round(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, Int& arg, int bits_to_keep = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) 544 { 545 // Precondition: exponent of res must have been set before this function is called 546 // as we may need to adjust it based on how many bits_to_keep in arg are set. 547 using default_ops::eval_bit_test; 548 using default_ops::eval_get_sign; 549 using default_ops::eval_increment; 550 using default_ops::eval_left_shift; 551 using default_ops::eval_lsb; 552 using default_ops::eval_msb; 553 using default_ops::eval_right_shift; 554 555 // cancellation may have resulted in arg being all zeros: 556 if (eval_get_sign(arg) == 0) 557 { 558 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero; 559 res.sign() = false; 560 res.bits() = static_cast<limb_type>(0u); 561 return; 562 } 563 int msb = eval_msb(arg); 564 if (static_cast<int>(bits_to_keep) > msb + 1) 565 { 566 // Must have had cancellation in subtraction, 567 // or be converting from a narrower type, so shift left: 568 res.bits() = arg; 569 eval_left_shift(res.bits(), bits_to_keep - msb - 1); 570 res.exponent() -= static_cast<Exponent>(bits_to_keep - msb - 1); 571 } 572 else if (static_cast<int>(bits_to_keep) < msb + 1) 573 { 574 // We have more bits_to_keep than we need, so round as required, 575 // first get the rounding bit: 576 bool roundup = eval_bit_test(arg, msb - bits_to_keep); 577 // Then check for a tie: 578 if (roundup && (msb - bits_to_keep == (int)eval_lsb(arg))) 579 { 580 // Ties round towards even: 581 if (!eval_bit_test(arg, msb - bits_to_keep + 1)) 582 roundup = false; 583 } 584 // Shift off the bits_to_keep we don't need: 585 eval_right_shift(arg, msb - bits_to_keep + 1); 586 res.exponent() += static_cast<Exponent>(msb - bits_to_keep + 1); 587 if (roundup) 588 { 589 eval_increment(arg); 590 if (bits_to_keep) 591 { 592 if (eval_bit_test(arg, bits_to_keep)) 593 { 594 // This happens very very rairly, all the bits left after 595 // truncation must be 1's and we're rounding up an order of magnitude: 596 eval_right_shift(arg, 1u); 597 ++res.exponent(); 598 } 599 } 600 else 601 { 602 // We get here when bits_to_keep is zero but we're rounding up, 603 // as a result we end up with a single digit that is a 1: 604 ++bits_to_keep; 605 } 606 } 607 if (bits_to_keep != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) 608 { 609 // Normalize result when we're rounding to fewer bits than we can hold, only happens in conversions 610 // to narrower types: 611 eval_left_shift(arg, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - bits_to_keep); 612 res.exponent() -= static_cast<Exponent>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - bits_to_keep); 613 } 614 res.bits() = arg; 615 } 616 else 617 { 618 res.bits() = arg; 619 } 620 if (!bits_to_keep && !res.bits().limbs()[0]) 621 { 622 // We're keeping zero bits and did not round up, so result is zero: 623 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero; 624 return; 625 } 626 // Result must be normalized: 627 BOOST_ASSERT(((int)eval_msb(res.bits()) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)); 628 629 if (res.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) 630 { 631 // Overflow: 632 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity; 633 res.bits() = static_cast<limb_type>(0u); 634 } 635 else if (res.exponent() < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent) 636 { 637 // Underflow: 638 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero; 639 res.bits() = static_cast<limb_type>(0u); 640 } 641 } 642 643 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 644 inline void do_eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b) 645 { 646 if (a.exponent() < b.exponent()) 647 { 648 bool s = a.sign(); 649 do_eval_add(res, b, a); 650 if (res.sign() != s) 651 res.negate(); 652 return; 653 } 654 655 using default_ops::eval_add; 656 using default_ops::eval_bit_test; 657 658 typedef typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type exponent_type; 659 660 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt; 661 662 // Special cases first: 663 switch (a.exponent()) 664 { 665 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 666 { 667 bool s = a.sign(); 668 res = b; 669 res.sign() = s; 670 return; 671 } 672 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 673 if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan) 674 res = b; 675 else 676 res = a; 677 return; // result is still infinite. 678 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 679 res = a; 680 return; // result is still a NaN. 681 } 682 switch (b.exponent()) 683 { 684 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 685 res = a; 686 return; 687 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 688 res = b; 689 if (res.sign()) 690 res.negate(); 691 return; // result is infinite. 692 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 693 res = b; 694 return; // result is a NaN. 695 } 696 697 BOOST_STATIC_ASSERT(boost::integer_traits<exponent_type>::const_max - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent); 698 699 bool s = a.sign(); 700 dt = a.bits(); 701 if (a.exponent() > (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent()) 702 { 703 res.exponent() = a.exponent(); 704 } 705 else 706 { 707 exponent_type e_diff = a.exponent() - b.exponent(); 708 BOOST_ASSERT(e_diff >= 0); 709 eval_left_shift(dt, e_diff); 710 res.exponent() = a.exponent() - e_diff; 711 eval_add(dt, b.bits()); 712 } 713 714 copy_and_round(res, dt); 715 res.check_invariants(); 716 if (res.sign() != s) 717 res.negate(); 718 } 719 720 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 721 inline void do_eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b) 722 { 723 using default_ops::eval_bit_test; 724 using default_ops::eval_decrement; 725 using default_ops::eval_subtract; 726 727 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt; 728 729 // Special cases first: 730 switch (a.exponent()) 731 { 732 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 733 if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan) 734 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 735 else 736 { 737 bool s = a.sign(); 738 res = b; 739 if (res.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero) 740 res.sign() = false; 741 else if (res.sign() == s) 742 res.negate(); 743 } 744 return; 745 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 746 if ((b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan) || (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity)) 747 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 748 else 749 res = a; 750 return; 751 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 752 res = a; 753 return; // result is still a NaN. 754 } 755 switch (b.exponent()) 756 { 757 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 758 res = a; 759 return; 760 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 761 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity; 762 res.sign() = !a.sign(); 763 res.bits() = static_cast<limb_type>(0u); 764 return; // result is a NaN. 765 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 766 res = b; 767 return; // result is still a NaN. 768 } 769 770 bool s = a.sign(); 771 if ((a.exponent() > b.exponent()) || ((a.exponent() == b.exponent()) && a.bits().compare(b.bits()) >= 0)) 772 { 773 dt = a.bits(); 774 if (a.exponent() <= (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent()) 775 { 776 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e_diff = a.exponent() - b.exponent(); 777 eval_left_shift(dt, e_diff); 778 res.exponent() = a.exponent() - e_diff; 779 eval_subtract(dt, b.bits()); 780 } 781 else if (a.exponent() == (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent() + 1) 782 { 783 if ((eval_lsb(a.bits()) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1) 784 && (eval_lsb(b.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)) 785 { 786 eval_left_shift(dt, 1); 787 eval_decrement(dt); 788 res.exponent() = a.exponent() - 1; 789 } 790 else 791 res.exponent() = a.exponent(); 792 } 793 else 794 res.exponent() = a.exponent(); 795 } 796 else 797 { 798 dt = b.bits(); 799 if (b.exponent() <= (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + a.exponent()) 800 { 801 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e_diff = a.exponent() - b.exponent(); 802 eval_left_shift(dt, -e_diff); 803 res.exponent() = b.exponent() + e_diff; 804 eval_subtract(dt, a.bits()); 805 } 806 else if (b.exponent() == (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + a.exponent() + 1) 807 { 808 if ((eval_lsb(a.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1) 809 && eval_lsb(b.bits())) 810 { 811 eval_left_shift(dt, 1); 812 eval_decrement(dt); 813 res.exponent() = b.exponent() - 1; 814 } 815 else 816 res.exponent() = b.exponent(); 817 } 818 else 819 res.exponent() = b.exponent(); 820 s = !s; 821 } 822 823 copy_and_round(res, dt); 824 if (res.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero) 825 res.sign() = false; 826 else if (res.sign() != s) 827 res.negate(); 828 res.check_invariants(); 829 } 830 831 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 832 inline void eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b) 833 { 834 if (a.sign() == b.sign()) 835 do_eval_add(res, a, b); 836 else 837 do_eval_subtract(res, a, b); 838 } 839 840 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 841 inline void eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a) 842 { 843 return eval_add(res, res, a); 844 } 845 846 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 847 inline void eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b) 848 { 849 if (a.sign() != b.sign()) 850 do_eval_add(res, a, b); 851 else 852 do_eval_subtract(res, a, b); 853 } 854 855 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 856 inline void eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a) 857 { 858 return eval_subtract(res, res, a); 859 } 860 861 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 862 inline void eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b) 863 { 864 using default_ops::eval_bit_test; 865 using default_ops::eval_multiply; 866 867 // Special cases first: 868 switch (a.exponent()) 869 { 870 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 871 { 872 if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan) 873 res = b; 874 else if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity) 875 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 876 else 877 { 878 bool s = a.sign() != b.sign(); 879 res = a; 880 res.sign() = s; 881 } 882 return; 883 } 884 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 885 switch (b.exponent()) 886 { 887 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 888 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 889 break; 890 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 891 res = b; 892 break; 893 default: 894 bool s = a.sign() != b.sign(); 895 res = a; 896 res.sign() = s; 897 break; 898 } 899 return; 900 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 901 res = a; 902 return; 903 } 904 if (b.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) 905 { 906 bool s = a.sign() != b.sign(); 907 res = b; 908 res.sign() = s; 909 return; 910 } 911 if ((a.exponent() > 0) && (b.exponent() > 0)) 912 { 913 if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent + 2 - a.exponent() < b.exponent()) 914 { 915 // We will certainly overflow: 916 bool s = a.sign() != b.sign(); 917 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity; 918 res.sign() = s; 919 res.bits() = static_cast<limb_type>(0u); 920 return; 921 } 922 } 923 if ((a.exponent() < 0) && (b.exponent() < 0)) 924 { 925 if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - 2 - a.exponent() > b.exponent()) 926 { 927 // We will certainly underflow: 928 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero; 929 res.sign() = a.sign() != b.sign(); 930 res.bits() = static_cast<limb_type>(0u); 931 return; 932 } 933 } 934 935 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt; 936 eval_multiply(dt, a.bits(), b.bits()); 937 res.exponent() = a.exponent() + b.exponent() - (Exponent)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1; 938 copy_and_round(res, dt); 939 res.check_invariants(); 940 res.sign() = a.sign() != b.sign(); 941 } 942 943 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 944 inline void eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a) 945 { 946 eval_multiply(res, res, a); 947 } 948 949 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U> 950 inline typename enable_if_c<is_unsigned<U>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const U& b) 951 { 952 using default_ops::eval_bit_test; 953 using default_ops::eval_multiply; 954 955 // Special cases first: 956 switch (a.exponent()) 957 { 958 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 959 { 960 bool s = a.sign(); 961 res = a; 962 res.sign() = s; 963 return; 964 } 965 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 966 if (b == 0) 967 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 968 else 969 res = a; 970 return; 971 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 972 res = a; 973 return; 974 } 975 976 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt; 977 typedef typename boost::multiprecision::detail::canonical<U, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::type canon_ui_type; 978 eval_multiply(dt, a.bits(), static_cast<canon_ui_type>(b)); 979 res.exponent() = a.exponent(); 980 copy_and_round(res, dt); 981 res.check_invariants(); 982 res.sign() = a.sign(); 983 } 984 985 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U> 986 inline typename enable_if_c<is_unsigned<U>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const U& b) 987 { 988 eval_multiply(res, res, b); 989 } 990 991 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S> 992 inline typename enable_if_c<is_signed<S>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, const S& b) 993 { 994 typedef typename make_unsigned<S>::type ui_type; 995 eval_multiply(res, a, static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(b))); 996 if (b < 0) 997 res.negate(); 998 } 999 1000 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S> 1001 inline typename enable_if_c<is_signed<S>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const S& b) 1002 { 1003 eval_multiply(res, res, b); 1004 } 1005 1006 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1007 inline void eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& u, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& v) 1008 { 1009 #ifdef BOOST_MSVC 1010 #pragma warning(push) 1011 #pragma warning(disable : 6326) // comparison of two constants 1012 #endif 1013 using default_ops::eval_bit_test; 1014 using default_ops::eval_get_sign; 1015 using default_ops::eval_increment; 1016 using default_ops::eval_qr; 1017 using default_ops::eval_subtract; 1018 1019 // 1020 // Special cases first: 1021 // 1022 switch (u.exponent()) 1023 { 1024 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1025 { 1026 switch (v.exponent()) 1027 { 1028 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1029 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1030 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 1031 return; 1032 } 1033 bool s = u.sign() != v.sign(); 1034 res = u; 1035 res.sign() = s; 1036 return; 1037 } 1038 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1039 { 1040 switch (v.exponent()) 1041 { 1042 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1043 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1044 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 1045 return; 1046 } 1047 bool s = u.sign() != v.sign(); 1048 res = u; 1049 res.sign() = s; 1050 return; 1051 } 1052 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1053 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 1054 return; 1055 } 1056 switch (v.exponent()) 1057 { 1058 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1059 { 1060 bool s = u.sign() != v.sign(); 1061 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend(); 1062 res.sign() = s; 1063 return; 1064 } 1065 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1066 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero; 1067 res.bits() = limb_type(0); 1068 res.sign() = u.sign() != v.sign(); 1069 return; 1070 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1071 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 1072 return; 1073 } 1074 1075 // We can scale u and v so that both are integers, then perform integer 1076 // division to obtain quotient q and remainder r, such that: 1077 // 1078 // q * v + r = u 1079 // 1080 // and hense: 1081 // 1082 // q + r/v = u/v 1083 // 1084 // From this, assuming q has cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count 1085 // bits we only need to determine whether 1086 // r/v is less than, equal to, or greater than 0.5 to determine rounding - 1087 // this we can do with a shift and comparison. 1088 // 1089 // We can set the exponent and sign of the result up front: 1090 // 1091 if ((v.exponent() < 0) && (u.exponent() > 0)) 1092 { 1093 // Check for overflow: 1094 if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent + v.exponent() < u.exponent() - 1) 1095 { 1096 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity; 1097 res.sign() = u.sign() != v.sign(); 1098 res.bits() = static_cast<limb_type>(0u); 1099 return; 1100 } 1101 } 1102 else if ((v.exponent() > 0) && (u.exponent() < 0)) 1103 { 1104 // Check for underflow: 1105 if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent + v.exponent() > u.exponent()) 1106 { 1107 // We will certainly underflow: 1108 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero; 1109 res.sign() = u.sign() != v.sign(); 1110 res.bits() = static_cast<limb_type>(0u); 1111 return; 1112 } 1113 } 1114 res.exponent() = u.exponent() - v.exponent() - 1; 1115 res.sign() = u.sign() != v.sign(); 1116 // 1117 // Now get the quotient and remainder: 1118 // 1119 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(u.bits()), t2(v.bits()), q, r; 1120 eval_left_shift(t, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count); 1121 eval_qr(t, t2, q, r); 1122 // 1123 // We now have either "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" 1124 // or "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1" significant 1125 // bits in q. 1126 // 1127 static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT; 1128 if (eval_bit_test(q, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)) 1129 { 1130 // 1131 // OK we have cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1 bits, 1132 // so we already have rounding info, 1133 // we just need to changes things if the last bit is 1 and either the 1134 // remainder is non-zero (ie we do not have a tie) or the quotient would 1135 // be odd if it were shifted to the correct number of bits (ie a tiebreak). 1136 // 1137 BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)); 1138 if ((q.limbs()[0] & 1u) && (eval_get_sign(r) || (q.limbs()[0] & 2u))) 1139 { 1140 eval_increment(q); 1141 } 1142 } 1143 else 1144 { 1145 // 1146 // We have exactly "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" bits in q. 1147 // Get rounding info, which we can get by comparing 2r with v. 1148 // We want to call copy_and_round to handle rounding and general cleanup, 1149 // so we'll left shift q and add some fake digits on the end to represent 1150 // how we'll be rounding. 1151 // 1152 BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)); 1153 static const unsigned lshift = (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count < limb_bits) ? 2 : limb_bits; 1154 eval_left_shift(q, lshift); 1155 res.exponent() -= lshift; 1156 eval_left_shift(r, 1u); 1157 int c = r.compare(v.bits()); 1158 if (c == 0) 1159 q.limbs()[0] |= static_cast<limb_type>(1u) << (lshift - 1); 1160 else if (c > 0) 1161 q.limbs()[0] |= (static_cast<limb_type>(1u) << (lshift - 1)) + static_cast<limb_type>(1u); 1162 } 1163 copy_and_round(res, q); 1164 #ifdef BOOST_MSVC 1165 #pragma warning(pop) 1166 #endif 1167 } 1168 1169 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1170 inline void eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1171 { 1172 eval_divide(res, res, arg); 1173 } 1174 1175 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U> 1176 inline typename enable_if_c<is_unsigned<U>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& u, const U& v) 1177 { 1178 #ifdef BOOST_MSVC 1179 #pragma warning(push) 1180 #pragma warning(disable : 6326) // comparison of two constants 1181 #endif 1182 using default_ops::eval_bit_test; 1183 using default_ops::eval_get_sign; 1184 using default_ops::eval_increment; 1185 using default_ops::eval_qr; 1186 using default_ops::eval_subtract; 1187 1188 // 1189 // Special cases first: 1190 // 1191 switch (u.exponent()) 1192 { 1193 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1194 { 1195 if (v == 0) 1196 { 1197 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 1198 return; 1199 } 1200 bool s = u.sign() != (v < 0); 1201 res = u; 1202 res.sign() = s; 1203 return; 1204 } 1205 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1206 res = u; 1207 return; 1208 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1209 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 1210 return; 1211 } 1212 if (v == 0) 1213 { 1214 bool s = u.sign(); 1215 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend(); 1216 res.sign() = s; 1217 return; 1218 } 1219 1220 // We can scale u and v so that both are integers, then perform integer 1221 // division to obtain quotient q and remainder r, such that: 1222 // 1223 // q * v + r = u 1224 // 1225 // and hense: 1226 // 1227 // q + r/v = u/v 1228 // 1229 // From this, assuming q has "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count, we only need to determine whether 1230 // r/v is less than, equal to, or greater than 0.5 to determine rounding - 1231 // this we can do with a shift and comparison. 1232 // 1233 // We can set the exponent and sign of the result up front: 1234 // 1235 int gb = msb(v); 1236 res.exponent() = u.exponent() - static_cast<Exponent>(gb) - static_cast<Exponent>(1); 1237 res.sign() = u.sign(); 1238 // 1239 // Now get the quotient and remainder: 1240 // 1241 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(u.bits()), q, r; 1242 eval_left_shift(t, gb + 1); 1243 eval_qr(t, number<typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::canonical_value(v), q, r); 1244 // 1245 // We now have either "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" or "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1" significant cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in q. 1246 // 1247 static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT; 1248 if (eval_bit_test(q, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)) 1249 { 1250 // 1251 // OK we have cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1 cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count, so we already have rounding info, 1252 // we just need to changes things if the last bit is 1 and the 1253 // remainder is non-zero (ie we do not have a tie). 1254 // 1255 BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)); 1256 if ((q.limbs()[0] & 1u) && eval_get_sign(r)) 1257 { 1258 eval_increment(q); 1259 } 1260 } 1261 else 1262 { 1263 // 1264 // We have exactly "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in q. 1265 // Get rounding info, which we can get by comparing 2r with v. 1266 // We want to call copy_and_round to handle rounding and general cleanup, 1267 // so we'll left shift q and add some fake cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count on the end to represent 1268 // how we'll be rounding. 1269 // 1270 BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)); 1271 static const unsigned lshift = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count < limb_bits ? 2 : limb_bits; 1272 eval_left_shift(q, lshift); 1273 res.exponent() -= lshift; 1274 eval_left_shift(r, 1u); 1275 int c = r.compare(number<typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::canonical_value(v)); 1276 if (c == 0) 1277 q.limbs()[0] |= static_cast<limb_type>(1u) << (lshift - 1); 1278 else if (c > 0) 1279 q.limbs()[0] |= (static_cast<limb_type>(1u) << (lshift - 1)) + static_cast<limb_type>(1u); 1280 } 1281 copy_and_round(res, q); 1282 #ifdef BOOST_MSVC 1283 #pragma warning(pop) 1284 #endif 1285 } 1286 1287 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U> 1288 inline typename enable_if_c<is_unsigned<U>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const U& v) 1289 { 1290 eval_divide(res, res, v); 1291 } 1292 1293 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S> 1294 inline typename enable_if_c<is_signed<S>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& u, const S& v) 1295 { 1296 typedef typename make_unsigned<S>::type ui_type; 1297 eval_divide(res, u, static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(v))); 1298 if (v < 0) 1299 res.negate(); 1300 } 1301 1302 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S> 1303 inline typename enable_if_c<is_signed<S>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const S& v) 1304 { 1305 eval_divide(res, res, v); 1306 } 1307 1308 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1309 inline int eval_get_sign(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1310 { 1311 return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero ? 0 : arg.sign() ? -1 : 1; 1312 } 1313 1314 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1315 inline bool eval_is_zero(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1316 { 1317 return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero; 1318 } 1319 1320 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1321 inline bool eval_eq(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b) 1322 { 1323 if (a.exponent() == b.exponent()) 1324 { 1325 if (a.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero) 1326 return true; 1327 return (a.sign() == b.sign()) && (a.bits().compare(b.bits()) == 0) && (a.exponent() != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan); 1328 } 1329 return false; 1330 } 1331 1332 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1333 inline void eval_convert_to(boost::long_long_type* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1334 { 1335 switch (arg.exponent()) 1336 { 1337 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1338 *res = 0; 1339 return; 1340 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1341 BOOST_THROW_EXCEPTION(std::runtime_error("Could not convert NaN to integer.")); 1342 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1343 *res = (std::numeric_limits<boost::long_long_type>::max)(); 1344 if (arg.sign()) 1345 *res = -*res; 1346 return; 1347 } 1348 typedef typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift_type; 1349 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type man(arg.bits()); 1350 shift_type shift = (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - arg.exponent(); 1351 if (shift > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1) 1352 { 1353 *res = 0; 1354 return; 1355 } 1356 if (arg.sign() && (arg.compare((std::numeric_limits<boost::long_long_type>::min)()) <= 0)) 1357 { 1358 *res = (std::numeric_limits<boost::long_long_type>::min)(); 1359 return; 1360 } 1361 else if (!arg.sign() && (arg.compare((std::numeric_limits<boost::long_long_type>::max)()) >= 0)) 1362 { 1363 *res = (std::numeric_limits<boost::long_long_type>::max)(); 1364 return; 1365 } 1366 1367 if (shift < 0) 1368 { 1369 if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - shift <= std::numeric_limits<boost::long_long_type>::digits) 1370 { 1371 // We have more bits in long_long_type than the float, so it's OK to left shift: 1372 eval_convert_to(res, man); 1373 *res <<= -shift; 1374 } 1375 else 1376 { 1377 *res = (std::numeric_limits<boost::long_long_type>::max)(); 1378 return; 1379 } 1380 } 1381 else 1382 { 1383 eval_right_shift(man, shift); 1384 eval_convert_to(res, man); 1385 } 1386 if (arg.sign()) 1387 { 1388 *res = -*res; 1389 } 1390 } 1391 1392 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1393 inline void eval_convert_to(boost::ulong_long_type* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1394 { 1395 switch (arg.exponent()) 1396 { 1397 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1398 *res = 0; 1399 return; 1400 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1401 BOOST_THROW_EXCEPTION(std::runtime_error("Could not convert NaN to integer.")); 1402 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1403 *res = (std::numeric_limits<boost::ulong_long_type>::max)(); 1404 return; 1405 } 1406 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type man(arg.bits()); 1407 typedef typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift_type; 1408 shift_type shift = (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - arg.exponent(); 1409 if (shift > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1) 1410 { 1411 *res = 0; 1412 return; 1413 } 1414 else if (shift < 0) 1415 { 1416 if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - shift <= std::numeric_limits<boost::ulong_long_type>::digits) 1417 { 1418 // We have more bits in ulong_long_type than the float, so it's OK to left shift: 1419 eval_convert_to(res, man); 1420 *res <<= -shift; 1421 return; 1422 } 1423 *res = (std::numeric_limits<boost::ulong_long_type>::max)(); 1424 return; 1425 } 1426 eval_right_shift(man, shift); 1427 eval_convert_to(res, man); 1428 } 1429 1430 template <class Float, unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1431 inline typename boost::enable_if_c<boost::is_float<Float>::value>::type eval_convert_to(Float* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& original_arg) 1432 { 1433 typedef cpp_bin_float<std::numeric_limits<Float>::digits, digit_base_2, void, Exponent, MinE, MaxE> conv_type; 1434 typedef typename common_type<typename conv_type::exponent_type, int>::type common_exp_type; 1435 // 1436 // Special cases first: 1437 // 1438 switch (original_arg.exponent()) 1439 { 1440 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1441 *res = 0; 1442 if (original_arg.sign()) 1443 *res = -*res; 1444 return; 1445 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1446 *res = std::numeric_limits<Float>::quiet_NaN(); 1447 return; 1448 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1449 *res = (std::numeric_limits<Float>::infinity)(); 1450 if (original_arg.sign()) 1451 *res = -*res; 1452 return; 1453 } 1454 // 1455 // Check for super large exponent that must be converted to infinity: 1456 // 1457 if (original_arg.exponent() > std::numeric_limits<Float>::max_exponent) 1458 { 1459 *res = std::numeric_limits<Float>::has_infinity ? std::numeric_limits<Float>::infinity() : (std::numeric_limits<Float>::max)(); 1460 if (original_arg.sign()) 1461 *res = -*res; 1462 return; 1463 } 1464 // 1465 // Figure out how many digits we will have in our result, 1466 // allowing for a possibly denormalized result: 1467 // 1468 common_exp_type digits_to_round_to = std::numeric_limits<Float>::digits; 1469 if (original_arg.exponent() < std::numeric_limits<Float>::min_exponent - 1) 1470 { 1471 common_exp_type diff = original_arg.exponent(); 1472 diff -= std::numeric_limits<Float>::min_exponent - 1; 1473 digits_to_round_to += diff; 1474 } 1475 if (digits_to_round_to < 0) 1476 { 1477 // Result must be zero: 1478 *res = 0; 1479 if (original_arg.sign()) 1480 *res = -*res; 1481 return; 1482 } 1483 // 1484 // Perform rounding first, then afterwards extract the digits: 1485 // 1486 cpp_bin_float<std::numeric_limits<Float>::digits, digit_base_2, Allocator, Exponent, MinE, MaxE> arg; 1487 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type bits(original_arg.bits()); 1488 arg.exponent() = original_arg.exponent(); 1489 copy_and_round(arg, bits, (int)digits_to_round_to); 1490 common_exp_type e = arg.exponent(); 1491 e -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1; 1492 static const unsigned limbs_needed = std::numeric_limits<Float>::digits / (sizeof(*arg.bits().limbs()) * CHAR_BIT) + (std::numeric_limits<Float>::digits % (sizeof(*arg.bits().limbs()) * CHAR_BIT) ? 1 : 0); 1493 unsigned first_limb_needed = arg.bits().size() - limbs_needed; 1494 *res = 0; 1495 e += first_limb_needed * sizeof(*arg.bits().limbs()) * CHAR_BIT; 1496 while (first_limb_needed < arg.bits().size()) 1497 { 1498 *res += std::ldexp(static_cast<Float>(arg.bits().limbs()[first_limb_needed]), static_cast<int>(e)); 1499 ++first_limb_needed; 1500 e += sizeof(*arg.bits().limbs()) * CHAR_BIT; 1501 } 1502 if (original_arg.sign()) 1503 *res = -*res; 1504 } 1505 1506 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1507 inline void eval_frexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, Exponent* e) 1508 { 1509 switch (arg.exponent()) 1510 { 1511 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1512 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1513 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1514 *e = 0; 1515 res = arg; 1516 return; 1517 } 1518 res = arg; 1519 *e = arg.exponent() + 1; 1520 res.exponent() = -1; 1521 } 1522 1523 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I> 1524 inline void eval_frexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, I* pe) 1525 { 1526 Exponent e; 1527 eval_frexp(res, arg, &e); 1528 if ((e > (std::numeric_limits<I>::max)()) || (e < (std::numeric_limits<I>::min)())) 1529 { 1530 BOOST_THROW_EXCEPTION(std::runtime_error("Exponent was outside of the range of the argument type to frexp.")); 1531 } 1532 *pe = static_cast<I>(e); 1533 } 1534 1535 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1536 inline void eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, Exponent e) 1537 { 1538 switch (arg.exponent()) 1539 { 1540 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1541 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1542 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1543 res = arg; 1544 return; 1545 } 1546 if ((e > 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent - e < arg.exponent())) 1547 { 1548 // Overflow: 1549 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend(); 1550 res.sign() = arg.sign(); 1551 } 1552 else if ((e < 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - e > arg.exponent())) 1553 { 1554 // Underflow: 1555 res = limb_type(0); 1556 } 1557 else 1558 { 1559 res = arg; 1560 res.exponent() += e; 1561 } 1562 } 1563 1564 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I> 1565 inline typename enable_if_c<is_unsigned<I>::value>::type eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, I e) 1566 { 1567 typedef typename make_signed<I>::type si_type; 1568 if (e > static_cast<I>((std::numeric_limits<si_type>::max)())) 1569 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend(); 1570 else 1571 eval_ldexp(res, arg, static_cast<si_type>(e)); 1572 } 1573 1574 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I> 1575 inline typename enable_if_c<is_signed<I>::value>::type eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, I e) 1576 { 1577 if ((e > (std::numeric_limits<Exponent>::max)()) || (e < (std::numeric_limits<Exponent>::min)())) 1578 { 1579 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend(); 1580 if (e < 0) 1581 res.negate(); 1582 } 1583 else 1584 eval_ldexp(res, arg, static_cast<Exponent>(e)); 1585 } 1586 1587 /* 1588 * Sign manipulation 1589 */ 1590 1591 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1592 inline void eval_abs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1593 { 1594 res = arg; 1595 res.sign() = false; 1596 } 1597 1598 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1599 inline void eval_fabs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1600 { 1601 res = arg; 1602 res.sign() = false; 1603 } 1604 1605 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1606 inline int eval_fpclassify(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1607 { 1608 switch (arg.exponent()) 1609 { 1610 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1611 return FP_ZERO; 1612 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1613 return FP_INFINITE; 1614 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1615 return FP_NAN; 1616 } 1617 return FP_NORMAL; 1618 } 1619 1620 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1621 inline void eval_sqrt(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1622 { 1623 using default_ops::eval_bit_test; 1624 using default_ops::eval_increment; 1625 using default_ops::eval_integer_sqrt; 1626 switch (arg.exponent()) 1627 { 1628 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1629 errno = EDOM; 1630 // fallthrough... 1631 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1632 res = arg; 1633 return; 1634 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1635 if (arg.sign()) 1636 { 1637 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 1638 errno = EDOM; 1639 } 1640 else 1641 res = arg; 1642 return; 1643 } 1644 if (arg.sign()) 1645 { 1646 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend(); 1647 errno = EDOM; 1648 return; 1649 } 1650 1651 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(arg.bits()), r, s; 1652 eval_left_shift(t, arg.exponent() & 1 ? cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count : cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1); 1653 eval_integer_sqrt(s, r, t); 1654 1655 if (!eval_bit_test(s, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)) 1656 { 1657 // We have exactly the right number of cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in the result, round as required: 1658 if (s.compare(r) < 0) 1659 { 1660 eval_increment(s); 1661 } 1662 } 1663 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type ae = arg.exponent(); 1664 res.exponent() = ae / 2; 1665 res.sign() = false; 1666 if ((ae & 1) && (ae < 0)) 1667 --res.exponent(); 1668 copy_and_round(res, s); 1669 } 1670 1671 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1672 inline void eval_floor(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1673 { 1674 using default_ops::eval_increment; 1675 switch (arg.exponent()) 1676 { 1677 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1678 errno = EDOM; 1679 // fallthrough... 1680 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1681 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1682 res = arg; 1683 return; 1684 } 1685 typedef typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift_type; 1686 shift_type shift = 1687 (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - arg.exponent() - 1; 1688 if ((arg.exponent() > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) || (shift <= 0)) 1689 { 1690 // Either arg is already an integer, or a special value: 1691 res = arg; 1692 return; 1693 } 1694 if (shift >= (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) 1695 { 1696 res = static_cast<signed_limb_type>(arg.sign() ? -1 : 0); 1697 return; 1698 } 1699 bool fractional = (shift_type)eval_lsb(arg.bits()) < shift; 1700 res = arg; 1701 eval_right_shift(res.bits(), shift); 1702 if (fractional && res.sign()) 1703 { 1704 eval_increment(res.bits()); 1705 if (eval_msb(res.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - shift) 1706 { 1707 // Must have extended result by one bit in the increment: 1708 --shift; 1709 ++res.exponent(); 1710 } 1711 } 1712 eval_left_shift(res.bits(), shift); 1713 } 1714 1715 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1716 inline void eval_ceil(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg) 1717 { 1718 using default_ops::eval_increment; 1719 switch (arg.exponent()) 1720 { 1721 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity: 1722 errno = EDOM; 1723 // fallthrough... 1724 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero: 1725 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan: 1726 res = arg; 1727 return; 1728 } 1729 typedef typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift_type; 1730 shift_type shift = (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - arg.exponent() - 1; 1731 if ((arg.exponent() > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) || (shift <= 0)) 1732 { 1733 // Either arg is already an integer, or a special value: 1734 res = arg; 1735 return; 1736 } 1737 if (shift >= (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) 1738 { 1739 bool s = arg.sign(); // takes care of signed zeros 1740 res = static_cast<signed_limb_type>(arg.sign() ? 0 : 1); 1741 res.sign() = s; 1742 return; 1743 } 1744 bool fractional = (shift_type)eval_lsb(arg.bits()) < shift; 1745 res = arg; 1746 eval_right_shift(res.bits(), shift); 1747 if (fractional && !res.sign()) 1748 { 1749 eval_increment(res.bits()); 1750 if (eval_msb(res.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - shift) 1751 { 1752 // Must have extended result by one bit in the increment: 1753 --shift; 1754 ++res.exponent(); 1755 } 1756 } 1757 eval_left_shift(res.bits(), shift); 1758 } 1759 1760 template <unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2> 1761 int eval_signbit(const cpp_bin_float<D1, B1, A1, E1, M1, M2>& val) 1762 { 1763 return val.sign(); 1764 } 1765 1766 template <unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2> 1767 inline std::size_t hash_value(const cpp_bin_float<D1, B1, A1, E1, M1, M2>& val) 1768 { 1769 std::size_t result = hash_value(val.bits()); 1770 boost::hash_combine(result, val.exponent()); 1771 boost::hash_combine(result, val.sign()); 1772 return result; 1773 } 1774 1775 } // namespace backends 1776 1777 namespace detail { 1778 1779 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinExponent, Exponent MaxExponent> 1780 struct transcendental_reduction_type<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent> > 1781 { 1782 // 1783 // The type used for trigonometric reduction needs 3 times the precision of the base type. 1784 // This is double the precision of the original type, plus the largest exponent supported. 1785 // As a practical measure the largest argument supported is 1/eps, as supporting larger 1786 // arguments requires the division of argument by PI/2 to also be done at higher precision, 1787 // otherwise the result (an integer) can not be represented exactly. 1788 // 1789 // See ARGUMENT REDUCTION FOR HUGE ARGUMENTS. K C Ng. 1790 // 1791 typedef boost::multiprecision::backends::cpp_bin_float< 1792 boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent>::bit_count * 3, 1793 boost::multiprecision::backends::digit_base_2, 1794 Allocator, Exponent, MinExponent, MaxExponent> type; 1795 }; 1796 1797 #ifdef BOOST_NO_SFINAE_EXPR 1798 1799 template <unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2, unsigned D2, backends::digit_base_type B2, class A2, class E2, E2 M3, E2 M4> 1800 struct is_explicitly_convertible<backends::cpp_bin_float<D1, B1, A1, E1, M1, M2>, backends::cpp_bin_float<D2, B2, A2, E2, M3, M4> > : public mpl::true_ 1801 {}; 1802 template <class FloatT, unsigned D2, backends::digit_base_type B2, class A2, class E2, E2 M3, E2 M4> 1803 struct is_explicitly_convertible<FloatT, backends::cpp_bin_float<D2, B2, A2, E2, M3, M4> > : public boost::is_floating_point<FloatT> 1804 {}; 1805 1806 #endif 1807 1808 } // namespace detail 1809 1810 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Exponent, Exponent MinE, Exponent MaxE, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates> 1811 inline boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> 1812 copysign BOOST_PREVENT_MACRO_SUBSTITUTION( 1813 const boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>& a, 1814 const boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>& b) 1815 { 1816 boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> res(a); 1817 res.backend().sign() = b.backend().sign(); 1818 return res; 1819 } 1820 1821 using backends::cpp_bin_float; 1822 using backends::digit_base_10; 1823 using backends::digit_base_2; 1824 1825 template <unsigned Digits, backends::digit_base_type DigitBase, class Exponent, Exponent MinE, Exponent MaxE, class Allocator> 1826 struct number_category<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > : public boost::mpl::int_<boost::multiprecision::number_kind_floating_point> 1827 {}; 1828 1829 template <unsigned Digits, backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> 1830 struct expression_template_default<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > 1831 { 1832 static const expression_template_option value = is_void<Allocator>::value ? et_off : et_on; 1833 }; 1834 1835 typedef number<backends::cpp_bin_float<50> > cpp_bin_float_50; 1836 typedef number<backends::cpp_bin_float<100> > cpp_bin_float_100; 1837 1838 typedef number<backends::cpp_bin_float<24, backends::digit_base_2, void, boost::int16_t, -126, 127>, et_off> cpp_bin_float_single; 1839 typedef number<backends::cpp_bin_float<53, backends::digit_base_2, void, boost::int16_t, -1022, 1023>, et_off> cpp_bin_float_double; 1840 typedef number<backends::cpp_bin_float<64, backends::digit_base_2, void, boost::int16_t, -16382, 16383>, et_off> cpp_bin_float_double_extended; 1841 typedef number<backends::cpp_bin_float<113, backends::digit_base_2, void, boost::int16_t, -16382, 16383>, et_off> cpp_bin_float_quad; 1842 typedef number<backends::cpp_bin_float<237, backends::digit_base_2, void, boost::int32_t, -262142, 262143>, et_off> cpp_bin_float_oct; 1843 1844 } // namespace multiprecision 1845 1846 namespace math { 1847 1848 using boost::multiprecision::copysign; 1849 using boost::multiprecision::signbit; 1850 1851 } // namespace math 1852 1853 } // namespace boost 1854 1855 #include <boost/multiprecision/cpp_bin_float/io.hpp> 1856 #include <boost/multiprecision/cpp_bin_float/transcendental.hpp> 1857 1858 namespace std { 1859 1860 // 1861 // numeric_limits [partial] specializations for the types declared in this header: 1862 // 1863 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 1864 class numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> > 1865 { 1866 typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> number_type; 1867 1868 public: 1869 BOOST_STATIC_CONSTEXPR bool is_specialized = true; 1870 static number_type(min)() 1871 { 1872 initializer.do_nothing(); 1873 static std::pair<bool, number_type> value; 1874 if (!value.first) 1875 { 1876 value.first = true; 1877 typedef typename boost::mpl::front<typename number_type::backend_type::unsigned_types>::type ui_type; 1878 value.second.backend() = ui_type(1u); 1879 value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent; 1880 } 1881 return value.second; 1882 } 1883 static number_type(max)() 1884 { 1885 initializer.do_nothing(); 1886 static std::pair<bool, number_type> value; 1887 if (!value.first) 1888 { 1889 value.first = true; 1890 if (boost::is_void<Allocator>::value) 1891 eval_complement(value.second.backend().bits(), value.second.backend().bits()); 1892 else 1893 { 1894 // We jump through hoops here using the backend type directly just to keep VC12 happy 1895 // (ie compiler workaround, for very strange compiler bug): 1896 using boost::multiprecision::default_ops::eval_add; 1897 using boost::multiprecision::default_ops::eval_decrement; 1898 using boost::multiprecision::default_ops::eval_left_shift; 1899 typedef typename number_type::backend_type::rep_type int_backend_type; 1900 typedef typename boost::mpl::front<typename int_backend_type::unsigned_types>::type ui_type; 1901 int_backend_type i; 1902 i = ui_type(1u); 1903 eval_left_shift(i, boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1); 1904 int_backend_type j(i); 1905 eval_decrement(i); 1906 eval_add(j, i); 1907 value.second.backend().bits() = j; 1908 } 1909 value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent; 1910 } 1911 return value.second; 1912 } 1913 BOOST_STATIC_CONSTEXPR number_type lowest() 1914 { 1915 return -(max)(); 1916 } 1917 BOOST_STATIC_CONSTEXPR int digits = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count; 1918 BOOST_STATIC_CONSTEXPR int digits10 = boost::multiprecision::detail::calc_digits10<digits>::value; 1919 // Is this really correct??? 1920 BOOST_STATIC_CONSTEXPR int max_digits10 = boost::multiprecision::detail::calc_max_digits10<digits>::value; 1921 BOOST_STATIC_CONSTEXPR bool is_signed = true; 1922 BOOST_STATIC_CONSTEXPR bool is_integer = false; 1923 BOOST_STATIC_CONSTEXPR bool is_exact = false; 1924 BOOST_STATIC_CONSTEXPR int radix = 2; 1925 static number_type epsilon() 1926 { 1927 initializer.do_nothing(); 1928 static std::pair<bool, number_type> value; 1929 if (!value.first) 1930 { 1931 // We jump through hoops here just to keep VC12 happy (ie compiler workaround, for very strange compiler bug): 1932 typedef typename boost::mpl::front<typename number_type::backend_type::unsigned_types>::type ui_type; 1933 value.first = true; 1934 value.second.backend() = ui_type(1u); 1935 value.second = ldexp(value.second, 1 - (int)digits); 1936 } 1937 return value.second; 1938 } 1939 // What value should this be???? 1940 static number_type round_error() 1941 { 1942 // returns 0.5 1943 initializer.do_nothing(); 1944 static std::pair<bool, number_type> value; 1945 if (!value.first) 1946 { 1947 value.first = true; 1948 // We jump through hoops here just to keep VC12 happy (ie compiler workaround, for very strange compiler bug): 1949 typedef typename boost::mpl::front<typename number_type::backend_type::unsigned_types>::type ui_type; 1950 value.second.backend() = ui_type(1u); 1951 value.second = ldexp(value.second, -1); 1952 } 1953 return value.second; 1954 } 1955 BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type min_exponent = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent; 1956 BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type min_exponent10 = (min_exponent / 1000) * 301L; 1957 BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type max_exponent = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent; 1958 BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type max_exponent10 = (max_exponent / 1000) * 301L; 1959 BOOST_STATIC_CONSTEXPR bool has_infinity = true; 1960 BOOST_STATIC_CONSTEXPR bool has_quiet_NaN = true; 1961 BOOST_STATIC_CONSTEXPR bool has_signaling_NaN = false; 1962 BOOST_STATIC_CONSTEXPR float_denorm_style has_denorm = denorm_absent; 1963 BOOST_STATIC_CONSTEXPR bool has_denorm_loss = false; 1964 static number_type infinity() 1965 { 1966 initializer.do_nothing(); 1967 static std::pair<bool, number_type> value; 1968 if (!value.first) 1969 { 1970 value.first = true; 1971 value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity; 1972 } 1973 return value.second; 1974 } 1975 static number_type quiet_NaN() 1976 { 1977 initializer.do_nothing(); 1978 static std::pair<bool, number_type> value; 1979 if (!value.first) 1980 { 1981 value.first = true; 1982 value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan; 1983 } 1984 return value.second; 1985 } 1986 BOOST_STATIC_CONSTEXPR number_type signaling_NaN() 1987 { 1988 return number_type(0); 1989 } 1990 BOOST_STATIC_CONSTEXPR number_type denorm_min() { return number_type(0); } 1991 BOOST_STATIC_CONSTEXPR bool is_iec559 = false; 1992 BOOST_STATIC_CONSTEXPR bool is_bounded = true; 1993 BOOST_STATIC_CONSTEXPR bool is_modulo = false; 1994 BOOST_STATIC_CONSTEXPR bool traps = true; 1995 BOOST_STATIC_CONSTEXPR bool tinyness_before = false; 1996 BOOST_STATIC_CONSTEXPR float_round_style round_style = round_to_nearest; 1997 1998 private: 1999 struct data_initializer 2000 { 2001 data_initializer() 2002 { 2003 std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::epsilon(); 2004 std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::round_error(); 2005 (std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::min)(); 2006 (std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::max)(); 2007 std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity(); 2008 std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN(); 2009 } 2010 void do_nothing() const {} 2011 }; 2012 static const data_initializer initializer; 2013 }; 2014 2015 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2016 const typename numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::data_initializer numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::initializer; 2017 2018 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION 2019 2020 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2021 BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::digits; 2022 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2023 BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::digits10; 2024 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2025 BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_digits10; 2026 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2027 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_signed; 2028 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2029 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_integer; 2030 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2031 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_exact; 2032 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2033 BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::radix; 2034 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2035 BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::min_exponent; 2036 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2037 BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::min_exponent10; 2038 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2039 BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_exponent; 2040 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2041 BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_exponent10; 2042 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2043 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_infinity; 2044 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2045 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_quiet_NaN; 2046 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2047 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_signaling_NaN; 2048 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2049 BOOST_CONSTEXPR_OR_CONST float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_denorm; 2050 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2051 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_denorm_loss; 2052 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2053 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_iec559; 2054 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2055 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_bounded; 2056 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2057 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_modulo; 2058 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2059 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::traps; 2060 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2061 BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::tinyness_before; 2062 template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates> 2063 BOOST_CONSTEXPR_OR_CONST float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::round_style; 2064 2065 #endif 2066 2067 } // namespace std 2068 2069 #endif 2070