]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 http://www.boost.org/LICENSE_1_ | |
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{ namespace cpp_bf_io_detail{ | |
10 | ||
11 | #ifdef BOOST_MSVC | |
12 | #pragma warning(push) | |
13 | #pragma warning(disable:4127) // conditional expression is constant | |
14 | #endif | |
15 | ||
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> | |
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> | |
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 | ||
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 | ||
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 | |
169 | ||
170 | namespace backends{ | |
171 | ||
172 | template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> | |
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 | } | |
353 | while(true); | |
354 | ||
355 | if(ss != sign()) | |
356 | negate(); | |
357 | } | |
358 | else | |
359 | { | |
360 | // Result is the ratio of two integers: we need to organise the | |
361 | // division so as to produce at least an N-bit result which we can | |
362 | // round according to the remainder. | |
363 | //cpp_int d = pow(cpp_int(5), -decimal_exp); | |
364 | do | |
365 | { | |
366 | cpp_int d; | |
367 | calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -decimal_exp, max_bits, error); | |
368 | int shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb(n) + msb(d); | |
369 | final_exponent = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp - calc_exp; | |
370 | if(shift > 0) | |
371 | { | |
372 | n <<= shift; | |
373 | final_exponent -= static_cast<Exponent>(shift); | |
374 | } | |
375 | cpp_int q, r; | |
376 | divide_qr(n, d, q, r); | |
377 | int gb = msb(q); | |
378 | BOOST_ASSERT((gb >= static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - 1)); | |
379 | // | |
380 | // Check for rounding conditions we have to | |
381 | // handle ourselves: | |
382 | // | |
383 | int roundup = 0; | |
384 | if(gb == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1) | |
385 | { | |
386 | // Exactly the right number of bits, use the remainder to round: | |
387 | roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, q); | |
388 | } | |
389 | 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))) | |
390 | { | |
391 | // Too many bits in q and the bits in q indicate a tie, but we can break that using r, | |
392 | // note that the radius of error in r is error/2 * q: | |
393 | int lshift = gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1; | |
394 | q >>= lshift; | |
395 | final_exponent += static_cast<Exponent>(lshift); | |
396 | BOOST_ASSERT((msb(q) >= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)); | |
397 | if(error && (r < (error / 2) * q)) | |
398 | roundup = -1; | |
399 | else if(error && (r + (error / 2) * q >= d)) | |
400 | roundup = -1; | |
401 | else | |
402 | roundup = r ? 2 : 1; | |
403 | } | |
404 | else if(error && (((error / 2) * q + r >= d) || (r < (error / 2) * q))) | |
405 | { | |
406 | // We might have been rounding up, or got the wrong quotient: can't tell! | |
407 | roundup = -1; | |
408 | } | |
409 | if(roundup < 0) | |
410 | { | |
411 | #ifdef BOOST_MP_STRESS_IO | |
412 | max_bits += 32; | |
413 | #else | |
414 | max_bits *= 2; | |
415 | #endif | |
416 | error = 0; | |
417 | if(shift > 0) | |
418 | { | |
419 | n >>= shift; | |
420 | final_exponent += static_cast<Exponent>(shift); | |
421 | } | |
422 | continue; | |
423 | } | |
424 | else if((roundup == 2) || ((roundup == 1) && q.backend().limbs()[0] & 1)) | |
425 | ++q; | |
426 | if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) | |
427 | { | |
428 | // Overflow: | |
429 | exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent; | |
430 | final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent; | |
431 | } | |
432 | else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent) | |
433 | { | |
434 | // Underflow: | |
435 | exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent; | |
436 | final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent; | |
437 | } | |
438 | else | |
439 | { | |
440 | exponent() = static_cast<Exponent>(final_exponent); | |
441 | final_exponent = 0; | |
442 | } | |
443 | copy_and_round(*this, q.backend()); | |
444 | if(ss != sign()) | |
445 | negate(); | |
446 | break; | |
447 | } | |
448 | while(true); | |
449 | } | |
450 | // | |
451 | // Check for scaling and/or over/under-flow: | |
452 | // | |
453 | final_exponent += exponent(); | |
454 | if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) | |
455 | { | |
456 | // Overflow: | |
457 | exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity; | |
458 | bits() = limb_type(0); | |
459 | } | |
460 | else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent) | |
461 | { | |
462 | // Underflow: | |
463 | exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero; | |
464 | bits() = limb_type(0); | |
465 | sign() = 0; | |
466 | } | |
467 | else | |
468 | { | |
469 | exponent() = static_cast<Exponent>(final_exponent); | |
470 | } | |
471 | return *this; | |
472 | } | |
473 | ||
474 | template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> | |
475 | std::string cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::str(std::streamsize dig, std::ios_base::fmtflags f) const | |
476 | { | |
477 | if(dig == 0) | |
478 | dig = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::max_digits10; | |
479 | ||
480 | bool scientific = (f & std::ios_base::scientific) == std::ios_base::scientific; | |
481 | bool fixed = !scientific && (f & std::ios_base::fixed); | |
482 | ||
483 | std::string s; | |
484 | ||
485 | if(exponent() <= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) | |
486 | { | |
487 | // How far to left-shift in order to demormalise the mantissa: | |
488 | boost::intmax_t shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1; | |
489 | boost::intmax_t digits_wanted = static_cast<int>(dig); | |
490 | 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())); | |
491 | // | |
492 | // For fixed formatting we want /dig/ digits after the decimal point, | |
493 | // so if the exponent is zero, allowing for the one digit before the | |
494 | // decimal point, we want 1 + dig digits etc. | |
495 | // | |
496 | if(fixed) | |
497 | digits_wanted += 1 + base10_exp; | |
498 | if(scientific) | |
499 | digits_wanted += 1; | |
500 | if(digits_wanted < -1) | |
501 | { | |
502 | // Fixed precision, no significant digits, and nothing to round! | |
503 | s = "0"; | |
504 | if(sign()) | |
505 | s.insert(static_cast<std::string::size_type>(0), 1, '-'); | |
506 | boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, true); | |
507 | return s; | |
508 | } | |
509 | // | |
510 | // power10 is the base10 exponent we need to multiply/divide by in order | |
511 | // to convert our denormalised number to an integer with the right number of digits: | |
512 | // | |
513 | boost::intmax_t power10 = digits_wanted - base10_exp - 1; | |
514 | // | |
515 | // If we calculate 5^power10 rather than 10^power10 we need to move | |
516 | // 2^power10 into /shift/ | |
517 | // | |
518 | shift -= power10; | |
519 | cpp_int i; | |
520 | int roundup = 0; // 0=no rounding, 1=tie, 2=up | |
521 | static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT; | |
522 | // | |
523 | // Set our working precision - this is heuristic based, we want | |
524 | // a value as small as possible > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count to avoid large computations | |
525 | // and excessive memory usage, but we also want to avoid having to | |
526 | // up the computation and start again at a higher precision. | |
527 | // So we round cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count up to the nearest whole number of limbs, and add | |
528 | // one limb for good measure. This works very well for small exponents, | |
529 | // but for larger exponents we add a few extra limbs to max_bits: | |
530 | // | |
531 | #ifdef BOOST_MP_STRESS_IO | |
532 | boost::intmax_t max_bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 32; | |
533 | #else | |
534 | 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; | |
535 | if(power10) | |
536 | max_bits += (msb(boost::multiprecision::detail::abs(power10)) / 8) * limb_bits; | |
537 | #endif | |
538 | do | |
539 | { | |
540 | boost::int64_t error = 0; | |
541 | boost::intmax_t calc_exp = 0; | |
542 | // | |
543 | // Our integer result is: bits() * 2^-shift * 5^power10 | |
544 | // | |
545 | i = bits(); | |
546 | if(shift < 0) | |
547 | { | |
548 | if(power10 >= 0) | |
549 | { | |
550 | // We go straight to the answer with all integer arithmetic, | |
551 | // the result is always exact and never needs rounding: | |
552 | BOOST_ASSERT(power10 <= (boost::intmax_t)INT_MAX); | |
553 | i <<= -shift; | |
554 | if(power10) | |
555 | i *= pow(cpp_int(5), static_cast<unsigned>(power10)); | |
556 | } | |
557 | else if(power10 < 0) | |
558 | { | |
559 | cpp_int d; | |
560 | calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -power10, max_bits, error); | |
561 | shift += calc_exp; | |
562 | BOOST_ASSERT(shift < 0); // Must still be true! | |
563 | i <<= -shift; | |
564 | cpp_int r; | |
565 | divide_qr(i, d, i, r); | |
566 | roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(r, d, error, i); | |
567 | if(roundup < 0) | |
568 | { | |
569 | #ifdef BOOST_MP_STRESS_IO | |
570 | max_bits += 32; | |
571 | #else | |
572 | max_bits *= 2; | |
573 | #endif | |
574 | shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10; | |
575 | continue; | |
576 | } | |
577 | } | |
578 | } | |
579 | else | |
580 | { | |
581 | // | |
582 | // Our integer is bits() * 2^-shift * 10^power10 | |
583 | // | |
584 | if(power10 > 0) | |
585 | { | |
586 | if(power10) | |
587 | { | |
588 | cpp_int t; | |
589 | calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(t, cpp_int(5), power10, max_bits, error); | |
590 | calc_exp += boost::multiprecision::cpp_bf_io_detail::restricted_multiply(i, i, t, max_bits, error); | |
591 | shift -= calc_exp; | |
592 | } | |
593 | if((shift < 0) || ((shift == 0) && error)) | |
594 | { | |
595 | // We only get here if we were asked for a crazy number of decimal digits - | |
596 | // more than are present in a 2^max_bits number. | |
597 | #ifdef BOOST_MP_STRESS_IO | |
598 | max_bits += 32; | |
599 | #else | |
600 | max_bits *= 2; | |
601 | #endif | |
602 | shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10; | |
603 | continue; | |
604 | } | |
605 | if(shift) | |
606 | { | |
607 | roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(i, shift - 1, error); | |
608 | if(roundup < 0) | |
609 | { | |
610 | #ifdef BOOST_MP_STRESS_IO | |
611 | max_bits += 32; | |
612 | #else | |
613 | max_bits *= 2; | |
614 | #endif | |
615 | shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10; | |
616 | continue; | |
617 | } | |
618 | i >>= shift; | |
619 | } | |
620 | } | |
621 | else | |
622 | { | |
623 | // We're right shifting, *and* dividing by 5^-power10, | |
624 | // so 5^-power10 can never be that large or we'd simply | |
625 | // get zero as a result, and that case is already handled above: | |
626 | cpp_int r; | |
627 | BOOST_ASSERT(-power10 < INT_MAX); | |
628 | cpp_int d = pow(cpp_int(5), static_cast<unsigned>(-power10)); | |
629 | d <<= shift; | |
630 | divide_qr(i, d, i, r); | |
631 | r <<= 1; | |
632 | int c = r.compare(d); | |
633 | roundup = c < 0 ? 0 : c == 0 ? 1 : 2; | |
634 | } | |
635 | } | |
636 | s = i.str(0, std::ios_base::fmtflags(0)); | |
637 | // | |
638 | // Check if we got the right number of digits, this | |
639 | // is really a test of whether we calculated the | |
640 | // decimal exponent correctly: | |
641 | // | |
642 | boost::intmax_t digits_got = i ? static_cast<boost::intmax_t>(s.size()) : 0; | |
643 | if(digits_got != digits_wanted) | |
644 | { | |
645 | base10_exp += digits_got - digits_wanted; | |
646 | if(fixed) | |
647 | digits_wanted = digits_got; // strange but true. | |
648 | power10 = digits_wanted - base10_exp - 1; | |
649 | shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - exponent() - 1 - power10; | |
650 | if(fixed) | |
651 | break; | |
652 | roundup = 0; | |
653 | } | |
654 | else | |
655 | break; | |
656 | } | |
657 | while(true); | |
658 | // | |
659 | // Check whether we need to round up: note that we could equally round up | |
660 | // the integer /i/ above, but since we need to perform the rounding *after* | |
661 | // the conversion to a string and the digit count check, we might as well | |
662 | // do it here: | |
663 | // | |
664 | if((roundup == 2) || ((roundup == 1) && ((s[s.size() - 1] - '0') & 1))) | |
665 | { | |
666 | boost::multiprecision::detail::round_string_up_at(s, static_cast<int>(s.size() - 1), base10_exp); | |
667 | } | |
668 | ||
669 | if(sign()) | |
670 | s.insert(static_cast<std::string::size_type>(0), 1, '-'); | |
671 | ||
672 | boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, false); | |
673 | } | |
674 | else | |
675 | { | |
676 | switch(exponent()) | |
677 | { | |
678 | case exponent_zero: | |
679 | s = sign() ? "-0" : f & std::ios_base::showpos ? "+0" : "0"; | |
680 | boost::multiprecision::detail::format_float_string(s, 0, dig, f, true); | |
681 | break; | |
682 | case exponent_nan: | |
683 | s = "nan"; | |
684 | break; | |
685 | case exponent_infinity: | |
686 | s = sign() ? "-inf" : f & std::ios_base::showpos ? "+inf" : "inf"; | |
687 | break; | |
688 | } | |
689 | } | |
690 | return s; | |
691 | } | |
692 | ||
693 | #ifdef BOOST_MSVC | |
694 | #pragma warning(pop) | |
695 | #endif | |
696 | ||
697 | }}} // namespaces | |
698 | ||
699 | #endif | |
700 |