]> git.proxmox.com Git - ceph.git/blame - ceph/src/boost/boost/multiprecision/cpp_bin_float/io.hpp
import new upstream nautilus stable release 14.2.8
[ceph.git] / ceph / src / boost / boost / multiprecision / cpp_bin_float / io.hpp
CommitLineData
7c673cae
FG
1///////////////////////////////////////////////////////////////
2// Copyright 2013 John Maddock. Distributed under the Boost
3// Software License, Version 1.0. (See accompanying file
92f5a8d4 4// LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
7c673cae
FG
5
6#ifndef BOOST_MP_CPP_BIN_FLOAT_IO_HPP
7#define BOOST_MP_CPP_BIN_FLOAT_IO_HPP
8
92f5a8d4
TL
9namespace boost { namespace multiprecision {
10namespace cpp_bf_io_detail {
7c673cae
FG
11
12#ifdef BOOST_MSVC
13#pragma warning(push)
92f5a8d4 14#pragma warning(disable : 4127) // conditional expression is constant
7c673cae
FG
15#endif
16
7c673cae 17//
92f5a8d4 18// Multiplies a by b and shifts the result so it fits inside max_bits bits,
7c673cae
FG
19// returns by how much the result was shifted.
20//
21template <class I>
22inline I restricted_multiply(cpp_int& result, const cpp_int& a, const cpp_int& b, I max_bits, boost::int64_t& error)
23{
92f5a8d4
TL
24 result = a * b;
25 I gb = msb(result);
7c673cae 26 I rshift = 0;
92f5a8d4 27 if (gb > max_bits)
7c673cae 28 {
92f5a8d4
TL
29 rshift = gb - max_bits;
30 I lb = lsb(result);
7c673cae 31 int roundup = 0;
92f5a8d4 32 // The error rate increases by the error of both a and b,
7c673cae
FG
33 // this may be overly pessimistic in many case as we're assuming
34 // that a and b have the same level of uncertainty...
92f5a8d4 35 if (lb < rshift)
7c673cae 36 error = error ? error * 2 : 1;
92f5a8d4 37 if (rshift)
7c673cae
FG
38 {
39 BOOST_ASSERT(rshift < INT_MAX);
92f5a8d4 40 if (bit_test(result, static_cast<unsigned>(rshift - 1)))
7c673cae 41 {
92f5a8d4 42 if (lb == rshift - 1)
7c673cae
FG
43 roundup = 1;
44 else
45 roundup = 2;
46 }
47 result >>= rshift;
48 }
92f5a8d4 49 if ((roundup == 2) || ((roundup == 1) && (result.backend().limbs()[0] & 1)))
7c673cae
FG
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//
58template <class I>
59inline 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;
92f5a8d4 63 if (e == 1)
7c673cae
FG
64 {
65 result = a;
66 return exp;
67 }
92f5a8d4 68 else if (e == 2)
7c673cae
FG
69 {
70 return restricted_multiply(result, a, a, max_bits, error);
71 }
92f5a8d4 72 else if (e == 3)
7c673cae
FG
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);
92f5a8d4 82 if (e & 1)
7c673cae
FG
83 exp += restricted_multiply(result, result, a, max_bits, error);
84 return exp;
85}
86
87inline 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;
92f5a8d4 101 if (error_radius && ((int)msb(error_radius) >= location))
7c673cae 102 return -1;
92f5a8d4 103 if (bit_test(what, static_cast<unsigned>(location)))
7c673cae 104 {
92f5a8d4
TL
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.
7c673cae 109 cpp_int t = what - error_radius;
92f5a8d4 110 if ((int)lsb(t) >= location)
7c673cae
FG
111 return -1;
112 return 2;
113 }
92f5a8d4 114 else if (error)
7c673cae
FG
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
122inline 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);
92f5a8d4 149 if (c == 0)
7c673cae 150 return error ? -1 : 1;
92f5a8d4 151 if (c > 0)
7c673cae 152 {
92f5a8d4 153 if (error)
7c673cae
FG
154 {
155 r -= error * q;
156 return r.compare(d) > 0 ? 2 : -1;
157 }
158 return 2;
159 }
92f5a8d4 160 if (error)
7c673cae
FG
161 {
162 r += error * q;
163 return r.compare(d) < 0 ? 0 : -1;
164 }
165 return 0;
166}
167
92f5a8d4 168} // namespace cpp_bf_io_detail
7c673cae 169
92f5a8d4 170namespace backends {
7c673cae
FG
171
172template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
92f5a8d4 173cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::operator=(const char* s)
7c673cae 174{
92f5a8d4
TL
175 cpp_int n;
176 boost::intmax_t decimal_exp = 0;
177 boost::intmax_t digits_seen = 0;
7c673cae 178 static const boost::intmax_t max_digits_seen = 4 + (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count * 301L) / 1000;
92f5a8d4 179 bool ss = false;
7c673cae
FG
180 //
181 // Extract the sign:
182 //
92f5a8d4 183 if (*s == '-')
7c673cae
FG
184 {
185 ss = true;
186 ++s;
187 }
92f5a8d4 188 else if (*s == '+')
7c673cae
FG
189 ++s;
190 //
191 // Special cases first:
192 //
92f5a8d4 193 if ((std::strcmp(s, "nan") == 0) || (std::strcmp(s, "NaN") == 0) || (std::strcmp(s, "NAN") == 0))
7c673cae
FG
194 {
195 return *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
196 }
92f5a8d4 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))
7c673cae
FG
198 {
199 *this = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
92f5a8d4 200 if (ss)
7c673cae
FG
201 negate();
202 return *this;
203 }
204 //
205 // Digits before the point:
206 //
92f5a8d4 207 while (*s && (*s >= '0') && (*s <= '9'))
7c673cae
FG
208 {
209 n *= 10u;
210 n += *s - '0';
92f5a8d4 211 if (digits_seen || (*s != '0'))
7c673cae
FG
212 ++digits_seen;
213 ++s;
214 }
215 // The decimal point (we really should localise this!!)
92f5a8d4 216 if (*s && (*s == '.'))
7c673cae
FG
217 ++s;
218 //
219 // Digits after the point:
220 //
92f5a8d4 221 while (*s && (*s >= '0') && (*s <= '9'))
7c673cae
FG
222 {
223 n *= 10u;
224 n += *s - '0';
225 --decimal_exp;
92f5a8d4 226 if (digits_seen || (*s != '0'))
7c673cae
FG
227 ++digits_seen;
228 ++s;
92f5a8d4 229 if (digits_seen > max_digits_seen)
7c673cae
FG
230 break;
231 }
232 //
233 // Digits we're skipping:
234 //
92f5a8d4 235 while (*s && (*s >= '0') && (*s <= '9'))
7c673cae
FG
236 ++s;
237 //
238 // See if there's an exponent:
239 //
92f5a8d4 240 if (*s && ((*s == 'e') || (*s == 'E')))
7c673cae
FG
241 {
242 ++s;
92f5a8d4
TL
243 boost::intmax_t e = 0;
244 bool es = false;
245 if (*s && (*s == '-'))
7c673cae
FG
246 {
247 es = true;
248 ++s;
249 }
92f5a8d4 250 else if (*s && (*s == '+'))
7c673cae 251 ++s;
92f5a8d4 252 while (*s && (*s >= '0') && (*s <= '9'))
7c673cae
FG
253 {
254 e *= 10u;
255 e += *s - '0';
256 ++s;
257 }
92f5a8d4 258 if (es)
7c673cae
FG
259 e = -e;
260 decimal_exp += e;
261 }
92f5a8d4 262 if (*s)
7c673cae
FG
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 }
92f5a8d4 269 if (n == 0)
7c673cae
FG
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
92f5a8d4
TL
293 boost::int64_t error = 0;
294 boost::intmax_t calc_exp = 0;
7c673cae
FG
295 boost::intmax_t final_exponent = 0;
296
92f5a8d4 297 if (decimal_exp >= 0)
7c673cae
FG
298 {
299 // Nice and simple, the result is an integer...
300 do
301 {
302 cpp_int t;
92f5a8d4 303 if (decimal_exp)
7c673cae
FG
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;
92f5a8d4
TL
311 int rshift = msb(t) - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
312 if (rshift > 0)
7c673cae
FG
313 {
314 final_exponent += rshift;
315 int roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(t, rshift - 1, error);
316 t >>= rshift;
92f5a8d4 317 if ((roundup == 2) || ((roundup == 1) && t.backend().limbs()[0] & 1))
7c673cae 318 ++t;
92f5a8d4 319 else if (roundup < 0)
7c673cae
FG
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 }
92f5a8d4 334 if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
7c673cae
FG
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 }
92f5a8d4 339 else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
7c673cae
FG
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 {
92f5a8d4 347 exponent() = static_cast<Exponent>(final_exponent);
7c673cae
FG
348 final_exponent = 0;
349 }
350 copy_and_round(*this, t.backend());
351 break;
92f5a8d4 352 } while (true);
7c673cae 353
92f5a8d4 354 if (ss != sign())
7c673cae
FG
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;
92f5a8d4
TL
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);
7c673cae 368 final_exponent = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp - calc_exp;
92f5a8d4 369 if (shift > 0)
7c673cae
FG
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;
92f5a8d4 383 if (gb == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
7c673cae
FG
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 }
92f5a8d4 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)))
7c673cae
FG
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));
92f5a8d4 396 if (error && (r < (error / 2) * q))
7c673cae 397 roundup = -1;
92f5a8d4 398 else if (error && (r + (error / 2) * q >= d))
7c673cae
FG
399 roundup = -1;
400 else
401 roundup = r ? 2 : 1;
402 }
92f5a8d4 403 else if (error && (((error / 2) * q + r >= d) || (r < (error / 2) * q)))
7c673cae
FG
404 {
405 // We might have been rounding up, or got the wrong quotient: can't tell!
406 roundup = -1;
407 }
92f5a8d4 408 if (roundup < 0)
7c673cae
FG
409 {
410#ifdef BOOST_MP_STRESS_IO
411 max_bits += 32;
412#else
413 max_bits *= 2;
414#endif
415 error = 0;
92f5a8d4 416 if (shift > 0)
7c673cae
FG
417 {
418 n >>= shift;
419 final_exponent += static_cast<Exponent>(shift);
420 }
421 continue;
422 }
92f5a8d4 423 else if ((roundup == 2) || ((roundup == 1) && q.backend().limbs()[0] & 1))
7c673cae 424 ++q;
92f5a8d4 425 if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
7c673cae
FG
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 }
92f5a8d4 431 else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
7c673cae
FG
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 {
92f5a8d4 439 exponent() = static_cast<Exponent>(final_exponent);
7c673cae
FG
440 final_exponent = 0;
441 }
442 copy_and_round(*this, q.backend());
92f5a8d4 443 if (ss != sign())
7c673cae
FG
444 negate();
445 break;
92f5a8d4 446 } while (true);
7c673cae
FG
447 }
448 //
449 // Check for scaling and/or over/under-flow:
450 //
451 final_exponent += exponent();
92f5a8d4 452 if (final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
7c673cae
FG
453 {
454 // Overflow:
455 exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
92f5a8d4 456 bits() = limb_type(0);
7c673cae 457 }
92f5a8d4 458 else if (final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
7c673cae
FG
459 {
460 // Underflow:
461 exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
92f5a8d4
TL
462 bits() = limb_type(0);
463 sign() = 0;
7c673cae
FG
464 }
465 else
466 {
467 exponent() = static_cast<Exponent>(final_exponent);
468 }
469 return *this;
470}
471
472template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
473std::string cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::str(std::streamsize dig, std::ios_base::fmtflags f) const
474{
92f5a8d4 475 if (dig == 0)
7c673cae
FG
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;
92f5a8d4 479 bool fixed = !scientific && (f & std::ios_base::fixed);
7c673cae
FG
480
481 std::string s;
482
92f5a8d4 483 if (exponent() <= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
7c673cae
FG
484 {
485 // How far to left-shift in order to demormalise the mantissa:
92f5a8d4 486 boost::intmax_t shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - (boost::intmax_t)exponent() - 1;
7c673cae 487 boost::intmax_t digits_wanted = static_cast<int>(dig);
92f5a8d4 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()));
7c673cae
FG
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 //
92f5a8d4 494 if (fixed)
7c673cae 495 digits_wanted += 1 + base10_exp;
92f5a8d4 496 if (scientific)
7c673cae 497 digits_wanted += 1;
92f5a8d4 498 if (digits_wanted < -1)
7c673cae
FG
499 {
500 // Fixed precision, no significant digits, and nothing to round!
501 s = "0";
92f5a8d4 502 if (sign())
7c673cae
FG
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;
92f5a8d4
TL
517 cpp_int i;
518 int roundup = 0; // 0=no rounding, 1=tie, 2=up
7c673cae
FG
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;
92f5a8d4 533 if (power10)
7c673cae
FG
534 max_bits += (msb(boost::multiprecision::detail::abs(power10)) / 8) * limb_bits;
535#endif
536 do
537 {
92f5a8d4 538 boost::int64_t error = 0;
7c673cae
FG
539 boost::intmax_t calc_exp = 0;
540 //
541 // Our integer result is: bits() * 2^-shift * 5^power10
542 //
543 i = bits();
92f5a8d4 544 if (shift < 0)
7c673cae 545 {
92f5a8d4 546 if (power10 >= 0)
7c673cae
FG
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;
92f5a8d4 552 if (power10)
7c673cae
FG
553 i *= pow(cpp_int(5), static_cast<unsigned>(power10));
554 }
92f5a8d4 555 else if (power10 < 0)
7c673cae
FG
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);
92f5a8d4 565 if (roundup < 0)
7c673cae
FG
566 {
567#ifdef BOOST_MP_STRESS_IO
568 max_bits += 32;
569#else
570 max_bits *= 2;
571#endif
11fdf7f2 572 shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
7c673cae
FG
573 continue;
574 }
575 }
576 }
577 else
578 {
579 //
580 // Our integer is bits() * 2^-shift * 10^power10
581 //
92f5a8d4 582 if (power10 > 0)
7c673cae 583 {
92f5a8d4 584 if (power10)
7c673cae
FG
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 }
92f5a8d4 591 if ((shift < 0) || ((shift == 0) && error))
7c673cae
FG
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
11fdf7f2 600 shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
7c673cae
FG
601 continue;
602 }
92f5a8d4 603 if (shift)
7c673cae
FG
604 {
605 roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(i, shift - 1, error);
92f5a8d4 606 if (roundup < 0)
7c673cae
FG
607 {
608#ifdef BOOST_MP_STRESS_IO
609 max_bits += 32;
610#else
611 max_bits *= 2;
612#endif
11fdf7f2 613 shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
7c673cae
FG
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;
92f5a8d4 630 int c = r.compare(d);
7c673cae
FG
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;
92f5a8d4 641 if (digits_got != digits_wanted)
7c673cae
FG
642 {
643 base10_exp += digits_got - digits_wanted;
92f5a8d4
TL
644 if (fixed)
645 digits_wanted = digits_got; // strange but true.
7c673cae 646 power10 = digits_wanted - base10_exp - 1;
92f5a8d4
TL
647 shift = (boost::intmax_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10;
648 if (fixed)
7c673cae
FG
649 break;
650 roundup = 0;
651 }
652 else
653 break;
92f5a8d4 654 } while (true);
7c673cae
FG
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 //
92f5a8d4 661 if ((roundup == 2) || ((roundup == 1) && ((s[s.size() - 1] - '0') & 1)))
7c673cae
FG
662 {
663 boost::multiprecision::detail::round_string_up_at(s, static_cast<int>(s.size() - 1), base10_exp);
664 }
665
92f5a8d4 666 if (sign())
7c673cae
FG
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 {
92f5a8d4 673 switch (exponent())
7c673cae
FG
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
92f5a8d4
TL
694} // namespace backends
695}} // namespace boost::multiprecision
7c673cae
FG
696
697#endif