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