823 lines · plain
1///////////////////////////////////////////////////////////////////////////////2// Copyright Christopher Kormanyos 2014.3// Copyright John Maddock 2014.4// Copyright Paul Bristow 2014.5// Distributed under the Boost Software License,6// Version 1.0. (See accompanying file LICENSE_1_0.txt7// or copy at http://www.boost.org/LICENSE_1_0.txt)8//9 10// Implement quadruple-precision I/O stream operations.11 12#ifndef BOOST_MATH_CSTDFLOAT_IOSTREAM_2014_02_15_HPP_13 #define BOOST_MATH_CSTDFLOAT_IOSTREAM_2014_02_15_HPP_14 15 #include <boost/math/cstdfloat/cstdfloat_types.hpp>16 #include <boost/math/cstdfloat/cstdfloat_limits.hpp>17 #include <boost/math/cstdfloat/cstdfloat_cmath.hpp>18 19 #if defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_CMATH)20 #error You can not use <boost/math/cstdfloat/cstdfloat_iostream.hpp> with BOOST_CSTDFLOAT_NO_LIBQUADMATH_CMATH defined.21 #endif22 23 #if defined(BOOST_CSTDFLOAT_HAS_INTERNAL_FLOAT128_T) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT)24 25 #include <cstddef>26 #include <istream>27 #include <ostream>28 #include <sstream>29 #include <stdexcept>30 #include <string>31 #include <boost/math/tools/assert.hpp>32 #include <boost/math/tools/nothrow.hpp>33 #include <boost/math/tools/throw_exception.hpp>34 35namespace boost {36 namespace math {37 namespace detail {38 //39 // What follows is the input streaming code: this is not "proper" iostream code at all40 // but that's hard to write.41 // For now just pull in all the characters that could possibly form the number42 // and let libquadmath's string parser make use of it. This fixes most use cases43 // including CSV type formats such as those used by the Random lib.44 //45 inline std::string read_string_while(std::istream& is, std::string const& permitted_chars)46 {47 std::ios_base::iostate state = std::ios_base::goodbit;48 const std::istream::sentry sentry_check(is);49 std::string result;50 51 if (sentry_check)52 {53 int c = is.rdbuf()->sgetc();54 55 for (;; c = is.rdbuf()->snextc())56 if (std::istream::traits_type::eq_int_type(std::istream::traits_type::eof(), c))57 { // end of file:58 state |= std::ios_base::eofbit;59 break;60 }61 else if (permitted_chars.find_first_of(std::istream::traits_type::to_char_type(c)) == std::string::npos)62 {63 // Invalid numeric character, stop reading:64 //is.rdbuf()->sputbackc(static_cast<char>(c));65 break;66 }67 else68 {69 result.append(1, std::istream::traits_type::to_char_type(c));70 }71 }72 73 if (!result.size())74 state |= std::ios_base::failbit;75 is.setstate(state);76 return result;77 }78 79 }80 }81}82 83#if defined(__GNUC__) && !defined(BOOST_MATH_TEST_IO_AS_INTEL_QUAD)84 85 // Forward declarations of quadruple-precision string functions.86 extern "C" int quadmath_snprintf(char *str, size_t size, const char *format, ...) BOOST_MATH_NOTHROW;87 extern "C" boost::math::cstdfloat::detail::float_internal128_t strtoflt128(const char*, char **) BOOST_MATH_NOTHROW;88 89 namespace std90 {91 template<typename char_type, class traits_type>92 inline std::basic_ostream<char_type, traits_type>& operator<<(std::basic_ostream<char_type, traits_type>& os, const boost::math::cstdfloat::detail::float_internal128_t& x)93 {94 std::basic_ostringstream<char_type, traits_type> ostr;95 ostr.flags(os.flags());96 ostr.imbue(os.getloc());97 ostr.precision(os.precision());98 99 char my_buffer[64U];100 101 const int my_prec = static_cast<int>(os.precision());102 const int my_digits = ((my_prec == 0) ? 36 : my_prec);103 104 const std::ios_base::fmtflags my_flags = os.flags();105 106 char my_format_string[8U];107 108 std::size_t my_format_string_index = 0U;109 110 my_format_string[my_format_string_index] = '%';111 ++my_format_string_index;112 113 if(my_flags & std::ios_base::showpos) { my_format_string[my_format_string_index] = '+'; ++my_format_string_index; }114 if(my_flags & std::ios_base::showpoint) { my_format_string[my_format_string_index] = '#'; ++my_format_string_index; }115 116 my_format_string[my_format_string_index + 0U] = '.';117 my_format_string[my_format_string_index + 1U] = '*';118 my_format_string[my_format_string_index + 2U] = 'Q';119 120 my_format_string_index += 3U;121 122 char the_notation_char;123 124 if (my_flags & std::ios_base::scientific) { the_notation_char = 'e'; }125 else if(my_flags & std::ios_base::fixed) { the_notation_char = 'f'; }126 else { the_notation_char = 'g'; }127 128 my_format_string[my_format_string_index + 0U] = the_notation_char;129 my_format_string[my_format_string_index + 1U] = 0;130 131 const int v = ::quadmath_snprintf(my_buffer,132 static_cast<int>(sizeof(my_buffer)),133 my_format_string,134 my_digits,135 x);136 137 if(v < 0) { BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Formatting of boost::float128_t failed internally in quadmath_snprintf().")); }138 139 if(v >= static_cast<int>(sizeof(my_buffer) - 1U))140 {141 // Evidently there is a really long floating-point string here,142 // such as a small decimal representation in non-scientific notation.143 // So we have to use dynamic memory allocation for the output144 // string buffer.145 146 char* my_buffer2 = nullptr;147 148#ifndef BOOST_MATH_NO_EXCEPTIONS149 try150 {151#endif152 my_buffer2 = new char[v + 3];153#ifndef BOOST_MATH_NO_EXCEPTIONS154 }155 catch(const std::bad_alloc&)156 {157 BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Formatting of boost::float128_t failed while allocating memory."));158 }159#endif160 const int v2 = ::quadmath_snprintf(my_buffer2,161 v + 3,162 my_format_string,163 my_digits,164 x);165 166 if(v2 >= v + 3)167 {168 BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Formatting of boost::float128_t failed."));169 }170 171 static_cast<void>(ostr << my_buffer2);172 173 delete [] my_buffer2;174 }175 else176 {177 static_cast<void>(ostr << my_buffer);178 }179 180 return (os << ostr.str());181 }182 183 template<typename char_type, class traits_type>184 inline std::basic_istream<char_type, traits_type>& operator>>(std::basic_istream<char_type, traits_type>& is, boost::math::cstdfloat::detail::float_internal128_t& x)185 {186 std::string str = boost::math::detail::read_string_while(is, "+-eE.0123456789infINFnanNANinfinityINFINITY");187 188 char* p_end;189 190 x = strtoflt128(str.c_str(), &p_end);191 192 if(static_cast<std::ptrdiff_t>(p_end - str.c_str()) != static_cast<std::ptrdiff_t>(str.length()))193 {194 for(std::string::const_reverse_iterator it = str.rbegin(); it != str.rend(); ++it)195 {196 static_cast<void>(is.putback(*it));197 }198 199 is.setstate(ios_base::failbit);200 201 BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Unable to interpret input string as a boost::float128_t"));202 }203 204 return is;205 }206 }207 208#elif defined(__INTEL_COMPILER) || defined(BOOST_MATH_TEST_IO_AS_INTEL_QUAD)209 210 // The section for I/O stream support for the ICC compiler is particularly211 // long, because these functions must be painstakingly synthesized from212 // manually-written routines (ICC does not support I/O stream operations213 // for its _Quad type).214 215 // The following string-extraction routines are based on the methodology216 // used in Boost.Multiprecision by John Maddock and Christopher Kormanyos.217 // This methodology has been slightly modified here for boost::float128_t.218 219 220 #include <cstring>221 #include <cctype>222 223 namespace boost { namespace math { namespace cstdfloat { namespace detail {224 225 template<class string_type>226 void format_float_string(string_type& str,227 int my_exp,228 int digits,229 const std::ios_base::fmtflags f,230 const bool iszero)231 {232 typedef typename string_type::size_type size_type;233 234 const bool scientific = ((f & std::ios_base::scientific) == std::ios_base::scientific);235 const bool fixed = ((f & std::ios_base::fixed) == std::ios_base::fixed);236 const bool showpoint = ((f & std::ios_base::showpoint) == std::ios_base::showpoint);237 const bool showpos = ((f & std::ios_base::showpos) == std::ios_base::showpos);238 239 const bool b_neg = ((str.size() != 0U) && (str[0] == '-'));240 241 if(b_neg)242 {243 str.erase(0, 1);244 }245 246 if(digits == 0)247 {248 digits = static_cast<int>((std::max)(str.size(), size_type(16)));249 }250 251 if(iszero || str.empty() || (str.find_first_not_of('0') == string_type::npos))252 {253 // We will be printing zero, even though the value might not254 // actually be zero (it just may have been rounded to zero).255 str = "0";256 257 if(scientific || fixed)258 {259 str.append(1, '.');260 str.append(size_type(digits), '0');261 262 if(scientific)263 {264 str.append("e+00");265 }266 }267 else268 {269 if(showpoint)270 {271 str.append(1, '.');272 if(digits > 1)273 {274 str.append(size_type(digits - 1), '0');275 }276 }277 }278 279 if(b_neg)280 {281 str.insert(0U, 1U, '-');282 }283 else if(showpos)284 {285 str.insert(0U, 1U, '+');286 }287 288 return;289 }290 291 if(!fixed && !scientific && !showpoint)292 {293 // Suppress trailing zeros.294 typename string_type::iterator pos = str.end();295 296 while(pos != str.begin() && *--pos == '0') { ; }297 298 if(pos != str.end())299 {300 ++pos;301 }302 303 str.erase(pos, str.end());304 305 if(str.empty())306 {307 str = '0';308 }309 }310 else if(!fixed || (my_exp >= 0))311 {312 // Pad out the end with zero's if we need to.313 314 std::ptrdiff_t chars = static_cast<std::ptrdiff_t>(str.size());315 chars = digits - chars;316 317 if(scientific)318 {319 ++chars;320 }321 322 if(chars > 0)323 {324 str.append(static_cast<size_type>(chars), '0');325 }326 }327 328 if(fixed || (!scientific && (my_exp >= -4) && (my_exp < digits)))329 {330 if((1 + my_exp) > static_cast<int>(str.size()))331 {332 // Just pad out the end with zeros.333 str.append(static_cast<size_type>((1 + my_exp) - static_cast<int>(str.size())), '0');334 335 if(showpoint || fixed)336 {337 str.append(".");338 }339 }340 else if(my_exp + 1 < static_cast<int>(str.size()))341 {342 if(my_exp < 0)343 {344 str.insert(0U, static_cast<size_type>(-1 - my_exp), '0');345 str.insert(0U, "0.");346 }347 else348 {349 // Insert the decimal point:350 str.insert(static_cast<size_type>(my_exp + 1), 1, '.');351 }352 }353 else if(showpoint || fixed) // we have exactly the digits we require to left of the point354 {355 str += ".";356 }357 358 if(fixed)359 {360 // We may need to add trailing zeros.361 int l = static_cast<int>(str.find('.') + 1U);362 l = digits - (static_cast<int>(str.size()) - l);363 364 if(l > 0)365 {366 str.append(size_type(l), '0');367 }368 }369 }370 else371 {372 // Scientific format:373 if(showpoint || (str.size() > 1))374 {375 str.insert(1U, 1U, '.');376 }377 378 str.append(1U, 'e');379 380 string_type e = std::to_string(std::abs(my_exp));381 382 if(e.size() < 2U)383 {384 e.insert(0U, 2U - e.size(), '0');385 }386 387 if(my_exp < 0)388 {389 e.insert(0U, 1U, '-');390 }391 else392 {393 e.insert(0U, 1U, '+');394 }395 396 str.append(e);397 }398 399 if(b_neg)400 {401 str.insert(0U, 1U, '-');402 }403 else if(showpos)404 {405 str.insert(0U, 1U, '+');406 }407 }408 409 template<class float_type, class type_a> inline void eval_convert_to(type_a* pa, const float_type& cb) { *pa = static_cast<type_a>(cb); }410 template<class float_type, class type_a> inline void eval_add (float_type& b, const type_a& a) { b += a; }411 template<class float_type, class type_a> inline void eval_subtract (float_type& b, const type_a& a) { b -= a; }412 template<class float_type, class type_a> inline void eval_multiply (float_type& b, const type_a& a) { b *= a; }413 template<class float_type> inline void eval_multiply (float_type& b, const float_type& cb, const float_type& cb2) { b = (cb * cb2); }414 template<class float_type, class type_a> inline void eval_divide (float_type& b, const type_a& a) { b /= a; }415 template<class float_type> inline void eval_log10 (float_type& b, const float_type& cb) { b = std::log10(cb); }416 template<class float_type> inline void eval_floor (float_type& b, const float_type& cb) { b = std::floor(cb); }417 418 inline void round_string_up_at(std::string& s, int pos, int& expon)419 {420 // This subroutine rounds up a string representation of a421 // number at the given position pos.422 423 if(pos < 0)424 {425 s.insert(0U, 1U, '1');426 s.erase(s.size() - 1U);427 ++expon;428 }429 else if(s[pos] == '9')430 {431 s[pos] = '0';432 round_string_up_at(s, pos - 1, expon);433 }434 else435 {436 if((pos == 0) && (s[pos] == '0') && (s.size() == 1))437 {438 ++expon;439 }440 441 ++s[pos];442 }443 }444 445 template<class float_type>446 std::string convert_to_string(float_type& x,447 std::streamsize digits,448 const std::ios_base::fmtflags f)449 {450 const bool isneg = (x < 0);451 const bool iszero = ((!isneg) ? bool(+x < (std::numeric_limits<float_type>::min)())452 : bool(-x < (std::numeric_limits<float_type>::min)()));453 const bool isnan = (x != x);454 const bool isinf = ((!isneg) ? bool(+x > (std::numeric_limits<float_type>::max)())455 : bool(-x > (std::numeric_limits<float_type>::max)()));456 457 int expon = 0;458 459 if(digits <= 0) { digits = std::numeric_limits<float_type>::max_digits10; }460 461 const int org_digits = static_cast<int>(digits);462 463 std::string result;464 465 if(iszero)466 {467 result = "0";468 }469 else if(isinf)470 {471 if(x < 0)472 {473 return "-inf";474 }475 else476 {477 return ((f & std::ios_base::showpos) == std::ios_base::showpos) ? "+inf" : "inf";478 }479 }480 else if(isnan)481 {482 return "nan";483 }484 else485 {486 // Start by figuring out the base-10 exponent.487 if(isneg) { x = -x; }488 489 float_type t;490 constexpr float_type ten = 10;491 492 eval_log10(t, x);493 eval_floor(t, t);494 eval_convert_to(&expon, t);495 496 if(-expon > std::numeric_limits<float_type>::max_exponent10 - 3)497 {498 int e = -expon / 2;499 500 const float_type t2 = boost::math::cstdfloat::detail::pown(ten, e);501 502 eval_multiply(t, t2, x);503 eval_multiply(t, t2);504 505 if((expon & 1) != 0)506 {507 eval_multiply(t, ten);508 }509 }510 else511 {512 t = boost::math::cstdfloat::detail::pown(ten, -expon);513 eval_multiply(t, x);514 }515 516 // Make sure that the value lies between [1, 10), and adjust if not.517 if(t < 1)518 {519 eval_multiply(t, 10);520 521 --expon;522 }523 else if(t >= 10)524 {525 eval_divide(t, 10);526 527 ++expon;528 }529 530 float_type digit;531 int cdigit;532 533 // Adjust the number of digits required based on formatting options.534 if(((f & std::ios_base::fixed) == std::ios_base::fixed) && (expon != -1))535 {536 digits += (expon + 1);537 }538 539 if((f & std::ios_base::scientific) == std::ios_base::scientific)540 {541 ++digits;542 }543 544 // Extract the base-10 digits one at a time.545 for(int i = 0; i < digits; ++i)546 {547 eval_floor(digit, t);548 eval_convert_to(&cdigit, digit);549 550 result += static_cast<char>('0' + cdigit);551 552 eval_subtract(t, digit);553 eval_multiply(t, ten);554 }555 if (result.size() == 0)556 result = "0";557 558 // Possibly round the result.559 if(digits >= 0)560 {561 eval_floor(digit, t);562 eval_convert_to(&cdigit, digit);563 eval_subtract(t, digit);564 565 if((cdigit == 5) && (t == 0))566 {567 // Use simple bankers rounding.568 569 if((static_cast<int>(*result.rbegin() - '0') & 1) != 0)570 {571 round_string_up_at(result, static_cast<int>(result.size() - 1U), expon);572 if (digits == 0) digits = 1;573 }574 }575 else if(cdigit >= 5)576 {577 round_string_up_at(result, static_cast<int>(result.size() - 1u), expon);578 if (digits == 0) digits = 1;579 }580 }581 }582 583 while((result.size() > static_cast<std::string::size_type>(digits)) && result.size())584 {585 // We may get here as a result of rounding.586 587 if(result.size() > 1U)588 {589 result.erase(result.size() - 1U);590 }591 else592 {593 if(expon > 0)594 {595 --expon; // so we put less padding in the result.596 }597 else598 {599 ++expon;600 }601 602 ++digits;603 }604 }605 606 if(isneg)607 {608 result.insert(0U, 1U, '-');609 }610 611 format_float_string(result, expon, org_digits, f, iszero);612 613 return result;614 }615 616 template <class float_type>617 bool convert_from_string(float_type& value, const char* p)618 {619 value = 0;620 621 if((p == nullptr) || (*p == '\0'))622 {623 return false;624 }625 626 bool is_neg = false;627 bool is_neg_expon = false;628 629 constexpr int ten = 10;630 631 int expon = 0;632 int digits_seen = 0;633 634 constexpr int max_digits = std::numeric_limits<float_type>::max_digits10 + 1;635 636 if(*p == '+')637 {638 ++p;639 }640 else if(*p == '-')641 {642 is_neg = true;643 ++p;644 }645 646 const bool isnan = ((std::strcmp(p, "nan") == 0) || (std::strcmp(p, "NaN") == 0) || (std::strcmp(p, "NAN") == 0));647 648 if(isnan)649 {650 eval_divide(value, 0);651 652 if(is_neg)653 {654 value = -value;655 }656 657 return true;658 }659 660 const bool isinf = ((std::strcmp(p, "inf") == 0) || (std::strcmp(p, "Inf") == 0) || (std::strcmp(p, "INF") == 0));661 662 if(isinf)663 {664 value = 1;665 eval_divide(value, 0);666 667 if(is_neg)668 {669 value = -value;670 }671 672 return true;673 }674 675 // Grab all the leading digits before the decimal point.676 while(std::isdigit(*p))677 {678 eval_multiply(value, ten);679 eval_add(value, static_cast<int>(*p - '0'));680 ++p;681 ++digits_seen;682 }683 684 if(*p == '.')685 {686 // Grab everything after the point, stop when we've seen687 // enough digits, even if there are actually more available.688 689 ++p;690 691 while(std::isdigit(*p))692 {693 eval_multiply(value, ten);694 eval_add(value, static_cast<int>(*p - '0'));695 ++p;696 --expon;697 698 if(++digits_seen > max_digits)699 {700 break;701 }702 }703 704 while(std::isdigit(*p))705 {706 ++p;707 }708 }709 710 // Parse the exponent.711 if((*p == 'e') || (*p == 'E'))712 {713 ++p;714 715 if(*p == '+')716 {717 ++p;718 }719 else if(*p == '-')720 {721 is_neg_expon = true;722 ++p;723 }724 725 int e2 = 0;726 727 while(std::isdigit(*p))728 {729 e2 *= 10;730 e2 += (*p - '0');731 ++p;732 }733 734 if(is_neg_expon)735 {736 e2 = -e2;737 }738 739 expon += e2;740 }741 742 if(expon)743 {744 // Scale by 10^expon. Note that 10^expon can be outside the range745 // of our number type, even though the result is within range.746 // If that looks likely, then split the calculation in two parts.747 float_type t;748 t = ten;749 750 if(expon > (std::numeric_limits<float_type>::min_exponent10 + 2))751 {752 t = boost::math::cstdfloat::detail::pown(t, expon);753 eval_multiply(value, t);754 }755 else756 {757 t = boost::math::cstdfloat::detail::pown(t, (expon + digits_seen + 1));758 eval_multiply(value, t);759 t = ten;760 t = boost::math::cstdfloat::detail::pown(t, (-digits_seen - 1));761 eval_multiply(value, t);762 }763 }764 765 if(is_neg)766 {767 value = -value;768 }769 770 return (*p == '\0');771 }772 } } } } // boost::math::cstdfloat::detail773 774 namespace std775 {776 template<typename char_type, class traits_type>777 inline std::basic_ostream<char_type, traits_type>& operator<<(std::basic_ostream<char_type, traits_type>& os, const boost::math::cstdfloat::detail::float_internal128_t& x)778 {779 boost::math::cstdfloat::detail::float_internal128_t non_const_x = x;780 781 const std::string str = boost::math::cstdfloat::detail::convert_to_string(non_const_x,782 os.precision(),783 os.flags());784 785 std::basic_ostringstream<char_type, traits_type> ostr;786 ostr.flags(os.flags());787 ostr.imbue(os.getloc());788 ostr.precision(os.precision());789 790 static_cast<void>(ostr << str);791 792 return (os << ostr.str());793 }794 795 template<typename char_type, class traits_type>796 inline std::basic_istream<char_type, traits_type>& operator>>(std::basic_istream<char_type, traits_type>& is, boost::math::cstdfloat::detail::float_internal128_t& x)797 {798 std::string str = boost::math::detail::read_string_while(is, "+-eE.0123456789infINFnanNANinfinityINFINITY");799 800 const bool conversion_is_ok = boost::math::cstdfloat::detail::convert_from_string(x, str.c_str());801 802 if(false == conversion_is_ok)803 {804 for(std::string::const_reverse_iterator it = str.rbegin(); it != str.rend(); ++it)805 {806 static_cast<void>(is.putback(*it));807 }808 809 is.setstate(ios_base::failbit);810 811 BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Unable to interpret input string as a boost::float128_t"));812 }813 814 return is;815 }816 }817 818 #endif // Use __GNUC__ or __INTEL_COMPILER libquadmath819 820 #endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT (i.e., the user would like to have libquadmath support)821 822#endif // BOOST_MATH_CSTDFLOAT_IOSTREAM_2014_02_15_HPP_823