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_MP_CPP_BIN_FLOAT_IO_HPP
7 #define BOOST_MP_CPP_BIN_FLOAT_IO_HPP
8
9 namespace boost { namespace multiprecision {
10 namespace cpp_bf_io_detail {
11
12 #ifdef BOOST_MSVC
13 #pragma warning(push)
14 #pragma warning(disable : 4127) // conditional expression is constant
15 #endif
16
17 //
18 // Multiplies a by b and shifts the result so it fits inside max_bits bits,
19 // returns by how much the result was shifted.
20 //
21 template <class I>
restricted_multiply(cpp_int & result,const cpp_int & a,const cpp_int & b,I max_bits,boost::int64_t & error)22 inline I restricted_multiply(cpp_int& result, const cpp_int& a, const cpp_int& b, I max_bits, boost::int64_t& error)
23 {
24 result = a * b;
25 I gb = msb(result);
26 I rshift = 0;
27 if (gb > max_bits)
28 {
29 rshift = gb - max_bits;
30 I lb = lsb(result);
31 int roundup = 0;
32 // The error rate increases by the error of both a and b,
33 // this may be overly pessimistic in many case as we're assuming
34 // that a and b have the same level of uncertainty...
35 if (lb < rshift)
36 error = error ? error * 2 : 1;
37 if (rshift)
38 {
39 BOOST_ASSERT(rshift < INT_MAX);
40 if (bit_test(result, static_cast<unsigned>(rshift - 1)))
41 {
42 if (lb == rshift - 1)
43 roundup = 1;
44 else
45 roundup = 2;
46 }
47 result >>= rshift;
48 }
49 if ((roundup == 2) || ((roundup == 1) && (result.backend().limbs()[0] & 1)))
50 ++result;
51 }
52 return rshift;
53 }
54 //
55 // Computes a^e shifted to the right so it fits in max_bits, returns how far
56 // to the right we are shifted.
57 //
58 template <class I>
restricted_pow(cpp_int & result,const cpp_int & a,I e,I max_bits,boost::int64_t & error)59 inline I restricted_pow(cpp_int& result, const cpp_int& a, I e, I max_bits, boost::int64_t& error)
60 {
61 BOOST_ASSERT(&result != &a);
62 I exp = 0;
63 if (e == 1)
64 {
65 result = a;
66 return exp;
67 }
68 else if (e == 2)
69 {
70 return restricted_multiply(result, a, a, max_bits, error);
71 }
72 else if (e == 3)
73 {
74 exp = restricted_multiply(result, a, a, max_bits, error);
75 exp += restricted_multiply(result, result, a, max_bits, error);
76 return exp;
77 }
78 I p = e / 2;
79 exp = restricted_pow(result, a, p, max_bits, error);
80 exp *= 2;
81 exp += restricted_multiply(result, result, result, max_bits, error);
82 if (e & 1)
83 exp += restricted_multiply(result, result, a, max_bits, error);
84 return exp;
85 }
86
get_round_mode(const cpp_int & what,boost::int64_t location,boost::int64_t error)87 inline int get_round_mode(const cpp_int& what, boost::int64_t location, boost::int64_t error)
88 {
89 //
90 // Can we round what at /location/, if the error in what is /error/ in
91 // units of 0.5ulp. Return:
92 //
93 // -1: Can't round.
94 // 0: leave as is.
95 // 1: tie.
96 // 2: round up.
97 //
98 BOOST_ASSERT(location >= 0);
99 BOOST_ASSERT(location < INT_MAX);
100 boost::int64_t error_radius = error & 1 ? (1 + error) / 2 : error / 2;
101 if (error_radius && ((int)msb(error_radius) >= location))
102 return -1;
103 if (bit_test(what, static_cast<unsigned>(location)))
104 {
105 if ((int)lsb(what) == location)
106 return error ? -1 : 1; // Either a tie or can't round depending on whether we have any error
107 if (!error)
108 return 2; // no error, round up.
109 cpp_int t = what - error_radius;
110 if ((int)lsb(t) >= location)
111 return -1;
112 return 2;
113 }
114 else if (error)
115 {
116 cpp_int t = what + error_radius;
117 return bit_test(t, static_cast<unsigned>(location)) ? -1 : 0;
118 }
119 return 0;
120 }
121
get_round_mode(cpp_int & r,cpp_int & d,boost::int64_t error,const cpp_int & q)122 inline int get_round_mode(cpp_int& r, cpp_int& d, boost::int64_t error, const cpp_int& q)
123 {
124 //
125 // Lets suppose we have an inexact division by d+delta, where the true
126 // value for the divisor is d, and with |delta| <= error/2, then
127 // we have calculated q and r such that:
128 //
129 // n r
130 // --- = q + -----------
131 // d + error d + error
132 //
133 // Rearranging for n / d we get:
134 //
135 // n delta*q + r
136 // --- = q + -------------
137 // d d
138 //
139 // So rounding depends on whether 2r + error * q > d.
140 //
141 // We return:
142 // 0 = down down.
143 // 1 = tie.
144 // 2 = round up.
145 // -1 = couldn't decide.
146 //
147 r <<= 1;
148 int c = r.compare(d);
149 if (c == 0)
150 return error ? -1 : 1;
151 if (c > 0)
152 {
153 if (error)
154 {
155 r -= error * q;
156 return r.compare(d) > 0 ? 2 : -1;
157 }
158 return 2;
159 }
160 if (error)
161 {
162 r += error * q;
163 return r.compare(d) < 0 ? 0 : -1;
164 }
165 return 0;
166 }
167
168 } // namespace cpp_bf_io_detail
169
170 namespace backends {
171
172 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
operator =(const char * s)173 cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::operator=(const char* s)
174 {
175 cpp_int n;
176 boost::intmax_t decimal_exp = 0;
177 boost::intmax_t digits_seen = 0;
178 static const boost::intmax_t max_digits_seen = 4 + (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count * 301L) / 1000;
179 bool ss = false;
180 //
181 // Extract the sign:
182 //
183 if (*s == '-')
184 {
185 ss = true;
186 ++s;
187 }
188 else if (*s == '+')
189 ++s;
190 //
191 // Special cases first:
192 //
193 if ((std::strcmp(s, "nan") == 0) || (std::strcmp(s, "NaN") == 0) || (std::strcmp(s, "NAN") == 0))
194 {
195 return *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
196 }
197 if ((std::strcmp(s, "inf") == 0) || (std::strcmp(s, "Inf") == 0) || (std::strcmp(s, "INF") == 0) || (std::strcmp(s, "infinity") == 0) || (std::strcmp(s, "Infinity") == 0) || (std::strcmp(s, "INFINITY") == 0))
198 {
199 *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
200 if (ss)
201 negate();
202 return *this;
203 }
204 //
205 // Digits before the point:
206 //
207 while (*s && (*s >= '0') && (*s <= '9'))
208 {
209 n *= 10u;
210 n += *s - '0';
211 if (digits_seen || (*s != '0'))
212 ++digits_seen;
213 ++s;
214 }
215 // The decimal point (we really should localise this!!)
216 if (*s && (*s == '.'))
217 ++s;
218 //
219 // Digits after the point:
220 //
221 while (*s && (*s >= '0') && (*s <= '9'))
222 {
223 n *= 10u;
224 n += *s - '0';
225 --decimal_exp;
226 if (digits_seen || (*s != '0'))
227 ++digits_seen;
228 ++s;
229 if (digits_seen > max_digits_seen)
230 break;
231 }
232 //
233 // Digits we're skipping:
234 //
235 while (*s && (*s >= '0') && (*s <= '9'))
236 ++s;
237 //
238 // See if there's an exponent:
239 //
240 if (*s && ((*s == 'e') || (*s == 'E')))
241 {
242 ++s;
243 boost::intmax_t e = 0;
244 bool es = false;
245 if (*s && (*s == '-'))
246 {
247 es = true;
248 ++s;
249 }
250 else if (*s && (*s == '+'))
251 ++s;
252 while (*s && (*s >= '0') && (*s <= '9'))
253 {
254 e *= 10u;
255 e += *s - '0';
256 ++s;
257 }
258 if (es)
259 e = -e;
260 decimal_exp += e;
261 }
262 if (*s)
263 {
264 //
265 // Oops unexpected input at the end of the number:
266 //
267 BOOST_THROW_EXCEPTION(std::runtime_error("Unable to parse string as a valid floating point number."));
268 }
269 if (n == 0)
270 {
271 // Result is necessarily zero:
272 *this = static_cast<limb_type>(0u);
273 return *this;
274 }
275
276 static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
277 //
278 // Set our working precision - this is heuristic based, we want
279 // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations
280 // and excessive memory usage, but we also want to avoid having to
281 // up the computation and start again at a higher precision.
282 // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add
283 // one limb for good measure. This works very well for small exponents,
284 // but for larger exponents we may may need to restart, we could add some
285 // extra precision right from the start for larger exponents, but this
286 // seems to be slightly slower in the *average* case:
287 //
288 #ifdef BOOST_MP_STRESS_IO
289 boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32;
290 #else
291 boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + ((cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count % limb_bits) ? (limb_bits - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count % limb_bits) : 0) + limb_bits;
292 #endif
293 boost::int64_t error = 0;
294 boost::intmax_t calc_exp = 0;
295 boost::intmax_t final_exponent = 0;
296
297 if (decimal_exp >= 0)
298 {
299 // Nice and simple, the result is an integer...
300 do
301 {
302 cpp_int t;
303 if (decimal_exp)
304 {
305 calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), decimal_exp, max_bits, error);
306 calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(t, t, n, max_bits, error);
307 }
308 else
309 t = n;
310 final_exponent = (boost::int64_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp + calc_exp;
311 int rshift = msb(t) - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
312 if (rshift > 0)
313 {
314 final_exponent += rshift;
315 int roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(t, rshift - 1, error);
316 t >>= rshift;
317 if ((roundup == 2) || ((roundup == 1) && t.backend().limbs()[0] & 1))
318 ++t;
319 else if (roundup < 0)
320 {
321 #ifdef BOOST_MP_STRESS_IO
322 max_bits += 32;
323 #else
324 max_bits *= 2;
325 #endif
326 error = 0;
327 continue;
328 }
329 }
330 else
331 {
332 BOOST_ASSERT(!error);
333 }
334 if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
335 {
336 exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
337 final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
338 }
339 else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
340 {
341 // Underflow:
342 exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
343 final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
344 }
345 else
346 {
347 exponent() = static_cast<Exponent>(final_exponent);
348 final_exponent = 0;
349 }
350 copy_and_round(*this, t.backend());
351 break;
352 } while (true);
353
354 if (ss != sign())
355 negate();
356 }
357 else
358 {
359 // Result is the ratio of two integers: we need to organise the
360 // division so as to produce at least an N-bit result which we can
361 // round according to the remainder.
362 //cpp_int d = pow(cpp_int(5), -decimal_exp);
363 do
364 {
365 cpp_int d;
366 calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -decimal_exp, max_bits, error);
367 int shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb(n) + msb(d);
368 final_exponent = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp - calc_exp;
369 if (shift > 0)
370 {
371 n <<= shift;
372 final_exponent -= static_cast<Exponent>(shift);
373 }
374 cpp_int q, r;
375 divide_qr(n, d, q, r);
376 int gb = msb(q);
377 BOOST_ASSERT((gb >= static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - 1));
378 //
379 // Check for rounding conditions we have to
380 // handle ourselves:
381 //
382 int roundup = 0;
383 if (gb == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
384 {
385 // Exactly the right number of bits, use the remainder to round:
386 roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, q);
387 }
388 else if (bit_test(q, gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) && ((int)lsb(q) == (gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)))
389 {
390 // Too many bits in q and the bits in q indicate a tie, but we can break that using r,
391 // note that the radius of error in r is error/2 * q:
392 int lshift = gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
393 q >>= lshift;
394 final_exponent += static_cast<Exponent>(lshift);
395 BOOST_ASSERT((msb(q) >= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
396 if (error && (r < (error / 2) * q))
397 roundup = -1;
398 else if (error && (r + (error / 2) * q >= d))
399 roundup = -1;
400 else
401 roundup = r ? 2 : 1;
402 }
403 else if (error && (((error / 2) * q + r >= d) || (r < (error / 2) * q)))
404 {
405 // We might have been rounding up, or got the wrong quotient: can't tell!
406 roundup = -1;
407 }
408 if (roundup < 0)
409 {
410 #ifdef BOOST_MP_STRESS_IO
411 max_bits += 32;
412 #else
413 max_bits *= 2;
414 #endif
415 error = 0;
416 if (shift > 0)
417 {
418 n >>= shift;
419 final_exponent += static_cast<Exponent>(shift);
420 }
421 continue;
422 }
423 else if ((roundup == 2) || ((roundup == 1) && q.backend().limbs()[0] & 1))
424 ++q;
425 if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
426 {
427 // Overflow:
428 exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
429 final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
430 }
431 else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
432 {
433 // Underflow:
434 exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
435 final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
436 }
437 else
438 {
439 exponent() = static_cast<Exponent>(final_exponent);
440 final_exponent = 0;
441 }
442 copy_and_round(*this, q.backend());
443 if (ss != sign())
444 negate();
445 break;
446 } while (true);
447 }
448 //
449 // Check for scaling and/or over/under-flow:
450 //
451 final_exponent += exponent();
452 if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
453 {
454 // Overflow:
455 exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
456 bits() = limb_type(0);
457 }
458 else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
459 {
460 // Underflow:
461 exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
462 bits() = limb_type(0);
463 sign() = 0;
464 }
465 else
466 {
467 exponent() = static_cast<Exponent>(final_exponent);
468 }
469 return *this;
470 }
471
472 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
str(std::streamsize dig,std::ios_base::fmtflags f) const473 std::string cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::str(std::streamsize dig, std::ios_base::fmtflags f) const
474 {
475 if (dig == 0)
476 dig = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::max_digits10;
477
478 bool scientific = (f & std::ios_base::scientific) == std::ios_base::scientific;
479 bool fixed = !scientific && (f & std::ios_base::fixed);
480
481 std::string s;
482
483 if (exponent() <= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
484 {
485 // How far to left-shift in order to demormalise the mantissa:
486 boost::intmax_t shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - (boost::intmax_t)exponent() - 1;
487 boost::intmax_t digits_wanted = static_cast<int>(dig);
488 boost::intmax_t base10_exp = exponent() >= 0 ? static_cast<boost::intmax_t>(std::floor(0.30103 * exponent())) : static_cast<boost::intmax_t>(std::ceil(0.30103 * exponent()));
489 //
490 // For fixed formatting we want /dig/ digits after the decimal point,
491 // so if the exponent is zero, allowing for the one digit before the
492 // decimal point, we want 1 + dig digits etc.
493 //
494 if (fixed)
495 digits_wanted += 1 + base10_exp;
496 if (scientific)
497 digits_wanted += 1;
498 if (digits_wanted < -1)
499 {
500 // Fixed precision, no significant digits, and nothing to round!
501 s = "0";
502 if (sign())
503 s.insert(static_cast<std::string::size_type>(0), 1, '-');
504 boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, true);
505 return s;
506 }
507 //
508 // power10 is the base10 exponent we need to multiply/divide by in order
509 // to convert our denormalised number to an integer with the right number of digits:
510 //
511 boost::intmax_t power10 = digits_wanted - base10_exp - 1;
512 //
513 // If we calculate 5^power10 rather than 10^power10 we need to move
514 // 2^power10 into /shift/
515 //
516 shift -= power10;
517 cpp_int i;
518 int roundup = 0; // 0=no rounding, 1=tie, 2=up
519 static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
520 //
521 // Set our working precision - this is heuristic based, we want
522 // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations
523 // and excessive memory usage, but we also want to avoid having to
524 // up the computation and start again at a higher precision.
525 // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add
526 // one limb for good measure. This works very well for small exponents,
527 // but for larger exponents we add a few extra limbs to max_bits:
528 //
529 #ifdef BOOST_MP_STRESS_IO
530 boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32;
531 #else
532 boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + ((cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count % limb_bits) ? (limb_bits - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count % limb_bits) : 0) + limb_bits;
533 if (power10)
534 max_bits += (msb(boost::multiprecision::detail::abs(power10)) / 8) * limb_bits;
535 #endif
536 do
537 {
538 boost::int64_t error = 0;
539 boost::intmax_t calc_exp = 0;
540 //
541 // Our integer result is: bits() * 2^-shift * 5^power10
542 //
543 i = bits();
544 if (shift < 0)
545 {
546 if (power10 >= 0)
547 {
548 // We go straight to the answer with all integer arithmetic,
549 // the result is always exact and never needs rounding:
550 BOOST_ASSERT(power10 <= (boost::intmax_t)INT_MAX);
551 i <<= -shift;
552 if (power10)
553 i *= pow(cpp_int(5), static_cast<unsigned>(power10));
554 }
555 else if (power10 < 0)
556 {
557 cpp_int d;
558 calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -power10, max_bits, error);
559 shift += calc_exp;
560 BOOST_ASSERT(shift < 0); // Must still be true!
561 i <<= -shift;
562 cpp_int r;
563 divide_qr(i, d, i, r);
564 roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, i);
565 if (roundup < 0)
566 {
567 #ifdef BOOST_MP_STRESS_IO
568 max_bits += 32;
569 #else
570 max_bits *= 2;
571 #endif
572 shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
573 continue;
574 }
575 }
576 }
577 else
578 {
579 //
580 // Our integer is bits() * 2^-shift * 10^power10
581 //
582 if (power10 > 0)
583 {
584 if (power10)
585 {
586 cpp_int t;
587 calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), power10, max_bits, error);
588 calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(i, i, t, max_bits, error);
589 shift -= calc_exp;
590 }
591 if ((shift < 0) || ((shift == 0) && error))
592 {
593 // We only get here if we were asked for a crazy number of decimal digits -
594 // more than are present in a 2^max_bits number.
595 #ifdef BOOST_MP_STRESS_IO
596 max_bits += 32;
597 #else
598 max_bits *= 2;
599 #endif
600 shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
601 continue;
602 }
603 if (shift)
604 {
605 roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(i, shift - 1, error);
606 if (roundup < 0)
607 {
608 #ifdef BOOST_MP_STRESS_IO
609 max_bits += 32;
610 #else
611 max_bits *= 2;
612 #endif
613 shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
614 continue;
615 }
616 i >>= shift;
617 }
618 }
619 else
620 {
621 // We're right shifting, *and* dividing by 5^-power10,
622 // so 5^-power10 can never be that large or we'd simply
623 // get zero as a result, and that case is already handled above:
624 cpp_int r;
625 BOOST_ASSERT(-power10 < INT_MAX);
626 cpp_int d = pow(cpp_int(5), static_cast<unsigned>(-power10));
627 d <<= shift;
628 divide_qr(i, d, i, r);
629 r <<= 1;
630 int c = r.compare(d);
631 roundup = c < 0 ? 0 : c == 0 ? 1 : 2;
632 }
633 }
634 s = i.str(0, std::ios_base::fmtflags(0));
635 //
636 // Check if we got the right number of digits, this
637 // is really a test of whether we calculated the
638 // decimal exponent correctly:
639 //
640 boost::intmax_t digits_got = i ? static_cast<boost::intmax_t>(s.size()) : 0;
641 if (digits_got != digits_wanted)
642 {
643 base10_exp += digits_got - digits_wanted;
644 if (fixed)
645 digits_wanted = digits_got; // strange but true.
646 power10 = digits_wanted - base10_exp - 1;
647 shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
648 if (fixed)
649 break;
650 roundup = 0;
651 }
652 else
653 break;
654 } while (true);
655 //
656 // Check whether we need to round up: note that we could equally round up
657 // the integer /i/ above, but since we need to perform the rounding *after*
658 // the conversion to a string and the digit count check, we might as well
659 // do it here:
660 //
661 if ((roundup == 2) || ((roundup == 1) && ((s[s.size() - 1] - '0') & 1)))
662 {
663 boost::multiprecision::detail::round_string_up_at(s, static_cast<int>(s.size() - 1), base10_exp);
664 }
665
666 if (sign())
667 s.insert(static_cast<std::string::size_type>(0), 1, '-');
668
669 boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, false);
670 }
671 else
672 {
673 switch (exponent())
674 {
675 case exponent_zero:
676 s = sign() ? "-0" : f & std::ios_base::showpos ? "+0" : "0";
677 boost::multiprecision::detail::format_float_string(s, 0, dig, f, true);
678 break;
679 case exponent_nan:
680 s = "nan";
681 break;
682 case exponent_infinity:
683 s = sign() ? "-inf" : f & std::ios_base::showpos ? "+inf" : "inf";
684 break;
685 }
686 }
687 return s;
688 }
689
690 #ifdef BOOST_MSVC
691 #pragma warning(pop)
692 #endif
693
694 } // namespace backends
695 }} // namespace boost::multiprecision
696
697 #endif
698