]>
Commit | Line | Data |
---|---|---|
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 |
9 | namespace boost { namespace multiprecision { |
10 | namespace 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 | // | |
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 | { | |
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 | // | |
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; | |
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 | ||
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; | |
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 | ||
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); | |
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 | 170 | namespace backends { |
7c673cae FG |
171 | |
172 | template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> | |
92f5a8d4 | 173 | cpp_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 | ||
472 | template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE> | |
473 | std::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 |