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_MATH_CPP_BIN_FLOAT_HPP
7#define BOOST_MATH_CPP_BIN_FLOAT_HPP
8
9#include <boost/multiprecision/cpp_int.hpp>
10#include <boost/multiprecision/integer.hpp>
11#include <boost/math/special_functions/trunc.hpp>
12#include <boost/multiprecision/detail/float_string_cvt.hpp>
13
14namespace boost{ namespace multiprecision{ namespace backends{
15
16enum digit_base_type
17{
18 digit_base_2 = 2,
19 digit_base_10 = 10
20};
21
22#ifdef BOOST_MSVC
23#pragma warning(push)
24#pragma warning(disable:4522 6326) // multiple assignment operators specified, comparison of two constants
25#endif
26
27namespace detail{
28
29template <class U>
30inline typename enable_if_c<is_unsigned<U>::value, bool>::type is_negative(U) { return false; }
31template <class S>
32inline typename disable_if_c<is_unsigned<S>::value, bool>::type is_negative(S s) { return s < 0; }
33
34}
35
36template <unsigned Digits, digit_base_type DigitBase = digit_base_10, class Allocator = void, class Exponent = int, Exponent MinExponent = 0, Exponent MaxExponent = 0>
37class cpp_bin_float
38{
39public:
40 static const unsigned bit_count = DigitBase == digit_base_2 ? Digits : (Digits * 1000uL) / 301uL + (((Digits * 1000uL) % 301) ? 2u : 1u);
41 typedef cpp_int_backend<is_void<Allocator>::value ? bit_count : 0, bit_count, is_void<Allocator>::value ? unsigned_magnitude : signed_magnitude, unchecked, Allocator> rep_type;
42 typedef cpp_int_backend<is_void<Allocator>::value ? 2 * bit_count : 0, 2 * bit_count, is_void<Allocator>::value ? unsigned_magnitude : signed_magnitude, unchecked, Allocator> double_rep_type;
43
44 typedef typename rep_type::signed_types signed_types;
45 typedef typename rep_type::unsigned_types unsigned_types;
46 typedef boost::mpl::list<float, double, long double> float_types;
47 typedef Exponent exponent_type;
48
49 static const exponent_type max_exponent_limit = boost::integer_traits<exponent_type>::const_max - 2 * static_cast<exponent_type>(bit_count);
50 static const exponent_type min_exponent_limit = boost::integer_traits<exponent_type>::const_min + 2 * static_cast<exponent_type>(bit_count);
51
52 BOOST_STATIC_ASSERT_MSG(MinExponent >= min_exponent_limit, "Template parameter MinExponent is too negative for our internal logic to function correctly, sorry!");
53 BOOST_STATIC_ASSERT_MSG(MaxExponent <= max_exponent_limit, "Template parameter MaxExponent is too large for our internal logic to function correctly, sorry!");
54 BOOST_STATIC_ASSERT_MSG(MinExponent <= 0, "Template parameter MinExponent can not be positive!");
55 BOOST_STATIC_ASSERT_MSG(MaxExponent >= 0, "Template parameter MaxExponent can not be negative!");
56
57 static const exponent_type max_exponent = MaxExponent == 0 ? max_exponent_limit : MaxExponent;
58 static const exponent_type min_exponent = MinExponent == 0 ? min_exponent_limit : MinExponent;
59
60 static const exponent_type exponent_zero = max_exponent + 1;
61 static const exponent_type exponent_infinity = max_exponent + 2;
62 static const exponent_type exponent_nan = max_exponent + 3;
63
64private:
65
66 rep_type m_data;
67 exponent_type m_exponent;
68 bool m_sign;
69public:
70 cpp_bin_float() BOOST_MP_NOEXCEPT_IF(noexcept(rep_type())) : m_data(), m_exponent(exponent_nan), m_sign(false) {}
71
72 cpp_bin_float(const cpp_bin_float &o) BOOST_MP_NOEXCEPT_IF(noexcept(rep_type(std::declval<const rep_type&>())))
73 : m_data(o.m_data), m_exponent(o.m_exponent), m_sign(o.m_sign) {}
74
75 template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
76 cpp_bin_float(const cpp_bin_float<D, B, A, E, MinE, MaxE> &o, typename boost::enable_if_c<(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = 0)
77 : m_exponent(o.exponent()), m_sign(o.sign())
78 {
79 typename cpp_bin_float<D, B, A, E, MinE, MaxE>::rep_type b(o.bits());
80 this->sign() = o.sign();
81 this->exponent() = o.exponent() + (int)bit_count - (int)cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count;
82 copy_and_round(*this, b);
83 }
84
85 template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
86 explicit cpp_bin_float(const cpp_bin_float<D, B, A, E, MinE, MaxE> &o, typename boost::disable_if_c<(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = 0)
87 : m_exponent(o.exponent()), m_sign(o.sign())
88 {
89 typename cpp_bin_float<D, B, A, E, MinE, MaxE>::rep_type b(o.bits());
90 this->sign() = o.sign();
91 this->exponent() = o.exponent() - (int)(cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count - bit_count);
92 copy_and_round(*this, b);
93 }
94
95 template <class Float>
96 cpp_bin_float(const Float& f,
97 typename boost::enable_if_c<
98 (number_category<Float>::value == number_kind_floating_point)
99 && (std::numeric_limits<Float>::digits <= (int)bit_count)
100 && (std::numeric_limits<Float>::radix == 2)
101 >::type const* = 0)
102 : m_data(), m_exponent(0), m_sign(false)
103 {
104 this->assign_float(f);
105 }
106
107 template <class Float>
108 explicit cpp_bin_float(const Float& f,
109 typename boost::enable_if_c<
110 (number_category<Float>::value == number_kind_floating_point)
111 && (std::numeric_limits<Float>::digits > (int)bit_count)
112 && (std::numeric_limits<Float>::radix == 2)
113 >::type const* = 0)
114 : m_data(), m_exponent(0), m_sign(false)
115 {
116 this->assign_float(f);
117 }
118
119 cpp_bin_float& operator=(const cpp_bin_float &o) BOOST_MP_NOEXCEPT_IF(noexcept(std::declval<rep_type&>() = std::declval<const rep_type&>()))
120 {
121 m_data = o.m_data;
122 m_exponent = o.m_exponent;
123 m_sign = o.m_sign;
124 return *this;
125 }
126
127 template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
128 cpp_bin_float& operator=(const cpp_bin_float<D, B, A, E, MinE, MaxE> &o)
129 {
130 typename cpp_bin_float<D, B, A, E, MinE, MaxE>::rep_type b(o.bits());
131 this->exponent() = o.exponent() + (int)bit_count - (int)cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count;
132 this->sign() = o.sign();
133 copy_and_round(*this, b);
134 return *this;
135 }
136
137 template <class Float>
138 typename boost::enable_if_c<
139 (number_category<Float>::value == number_kind_floating_point)
140 //&& (std::numeric_limits<Float>::digits <= (int)bit_count)
141 && (std::numeric_limits<Float>::radix == 2), cpp_bin_float&>::type operator=(const Float& f)
142 {
143 return assign_float(f);
144 }
145
146 template <class Float>
147 typename boost::enable_if_c<is_floating_point<Float>::value, cpp_bin_float&>::type assign_float(Float f)
148 {
149 BOOST_MATH_STD_USING
150 using default_ops::eval_add;
151 typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type;
152
153 switch((boost::math::fpclassify)(f))
154 {
155 case FP_ZERO:
156 m_data = limb_type(0);
157 m_sign = false;
158 m_exponent = exponent_zero;
159 return *this;
160 case FP_NAN:
161 m_data = limb_type(0);
162 m_sign = false;
163 m_exponent = exponent_nan;
164 return *this;
165 case FP_INFINITE:
166 m_data = limb_type(0);
167 m_sign = false;
168 m_exponent = exponent_infinity;
169 return *this;
170 }
171 if(f < 0)
172 {
173 *this = -f;
174 this->negate();
175 return *this;
176 }
177
178 typedef typename mpl::front<unsigned_types>::type ui_type;
179 m_data = static_cast<ui_type>(0u);
180 m_sign = false;
181 m_exponent = 0;
182
183 static const int bits = sizeof(int) * CHAR_BIT - 1;
184 int e;
185 f = frexp(f, &e);
186 while(f)
187 {
188 f = ldexp(f, bits);
189 e -= bits;
190#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
191 int ipart = itrunc(f);
192#else
193 int ipart = static_cast<int>(f);
194#endif
195 f -= ipart;
196 m_exponent += bits;
197 cpp_bin_float t;
198 t = static_cast<bf_int_type>(ipart);
199 eval_add(*this, t);
200 }
201 m_exponent += static_cast<Exponent>(e);
202 return *this;
203 }
204
205 template <class Float>
206 typename boost::enable_if_c<
207 (number_category<Float>::value == number_kind_floating_point)
208 && !boost::is_floating_point<Float>::value
209 /*&& (std::numeric_limits<number<Float> >::radix == 2)*/,
210 cpp_bin_float&>::type assign_float(Float f)
211 {
212 BOOST_MATH_STD_USING
213 using default_ops::eval_add;
214 using default_ops::eval_get_sign;
215 using default_ops::eval_convert_to;
216 using default_ops::eval_subtract;
217
218 typedef typename boost::multiprecision::detail::canonical<int, Float>::type f_int_type;
219 typedef typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type bf_int_type;
220
221 switch(eval_fpclassify(f))
222 {
223 case FP_ZERO:
224 m_data = limb_type(0);
225 m_sign = false;
226 m_exponent = exponent_zero;
227 return *this;
228 case FP_NAN:
229 m_data = limb_type(0);
230 m_sign = false;
231 m_exponent = exponent_nan;
232 return *this;
233 case FP_INFINITE:
234 m_data = limb_type(0);
235 m_sign = false;
236 m_exponent = exponent_infinity;
237 return *this;
238 }
239 if(eval_get_sign(f) < 0)
240 {
241 f.negate();
242 *this = f;
243 this->negate();
244 return *this;
245 }
246
247 typedef typename mpl::front<unsigned_types>::type ui_type;
248 m_data = static_cast<ui_type>(0u);
249 m_sign = false;
250 m_exponent = 0;
251
252 static const int bits = sizeof(int) * CHAR_BIT - 1;
253 int e;
254 eval_frexp(f, f, &e);
255 while(eval_get_sign(f) != 0)
256 {
257 eval_ldexp(f, f, bits);
258 e -= bits;
259 int ipart;
260 eval_convert_to(&ipart, f);
261 eval_subtract(f, static_cast<f_int_type>(ipart));
262 m_exponent += bits;
263 eval_add(*this, static_cast<bf_int_type>(ipart));
264 }
265 m_exponent += e;
266 if(m_exponent > max_exponent)
267 m_exponent = exponent_infinity;
268 if(m_exponent < min_exponent)
269 {
270 m_data = limb_type(0u);
271 m_exponent = exponent_zero;
272 m_sign = false;
273 }
274 else if(eval_get_sign(m_data) == 0)
275 {
276 m_exponent = exponent_zero;
277 m_sign = false;
278 }
279 return *this;
280 }
281
282 template <class I>
283 typename boost::enable_if<is_integral<I>, cpp_bin_float&>::type operator=(const I& i)
284 {
285 using default_ops::eval_bit_test;
286 if(!i)
287 {
288 m_data = static_cast<limb_type>(0);
289 m_exponent = exponent_zero;
290 m_sign = false;
291 }
292 else
293 {
294 typedef typename make_unsigned<I>::type ui_type;
295 ui_type fi = static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(i));
296 typedef typename boost::multiprecision::detail::canonical<ui_type, rep_type>::type ar_type;
297 m_data = static_cast<ar_type>(fi);
298 unsigned shift = msb(fi);
299 if(shift >= bit_count)
300 {
301 m_exponent = static_cast<Exponent>(shift);
302 m_data = static_cast<ar_type>(fi >> (shift + 1 - bit_count));
303 }
304 else
305 {
306 m_exponent = static_cast<Exponent>(shift);
307 eval_left_shift(m_data, bit_count - shift - 1);
308 }
309 BOOST_ASSERT(eval_bit_test(m_data, bit_count-1));
310 m_sign = detail::is_negative(i);
311 }
312 return *this;
313 }
314
315 cpp_bin_float& operator=(const char *s);
316
317 void swap(cpp_bin_float &o) BOOST_NOEXCEPT
318 {
319 m_data.swap(o.m_data);
320 std::swap(m_exponent, o.m_exponent);
321 std::swap(m_sign, o.m_sign);
322 }
323
324 std::string str(std::streamsize dig, std::ios_base::fmtflags f) const;
325
326 void negate()
327 {
328 if((m_exponent != exponent_zero) && (m_exponent != exponent_nan))
329 m_sign = !m_sign;
330 }
331
332 int compare(const cpp_bin_float &o) const BOOST_NOEXCEPT
333 {
334 if(m_sign != o.m_sign)
335 return m_sign ? -1 : 1;
336 int result;
337 if(m_exponent == exponent_nan)
338 return -1;
339 else if(m_exponent != o.m_exponent)
340 {
341 if(m_exponent == exponent_zero)
342 result = -1;
343 else if(o.m_exponent == exponent_zero)
344 result = 1;
345 else
346 result = m_exponent > o.m_exponent ? 1 : -1;
347 }
348 else
349 result = m_data.compare(o.m_data);
350 if(m_sign)
351 result = -result;
352 return result;
353 }
354 template <class A>
355 int compare(const A& o) const BOOST_NOEXCEPT
356 {
357 cpp_bin_float b;
358 b = o;
359 return compare(b);
360 }
361
362 rep_type& bits() { return m_data; }
363 const rep_type& bits()const { return m_data; }
364 exponent_type& exponent() { return m_exponent; }
365 const exponent_type& exponent()const { return m_exponent; }
366 bool& sign() { return m_sign; }
367 const bool& sign()const { return m_sign; }
368 void check_invariants()
369 {
370 using default_ops::eval_bit_test;
371 using default_ops::eval_is_zero;
372 if((m_exponent <= max_exponent) && (m_exponent >= min_exponent))
373 {
374 BOOST_ASSERT(eval_bit_test(m_data, bit_count - 1));
375 }
376 else
377 {
378 BOOST_ASSERT(m_exponent > max_exponent);
379 BOOST_ASSERT(m_exponent <= exponent_nan);
380 BOOST_ASSERT(eval_is_zero(m_data));
381 }
382 }
383 template<class Archive>
384 void serialize(Archive & ar, const unsigned int /*version*/)
385 {
386 ar & m_data;
387 ar & m_exponent;
388 ar & m_sign;
389 }
390};
391
392#ifdef BOOST_MSVC
393#pragma warning(pop)
394#endif
395
396template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class Int>
397inline void copy_and_round(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, Int &arg)
398{
399 // Precondition: exponent of res must have been set before this function is called
400 // as we may need to adjust it based on how many cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in arg are set.
401 using default_ops::eval_msb;
402 using default_ops::eval_lsb;
403 using default_ops::eval_left_shift;
404 using default_ops::eval_bit_test;
405 using default_ops::eval_right_shift;
406 using default_ops::eval_increment;
407 using default_ops::eval_get_sign;
408
409 // cancellation may have resulted in arg being all zeros:
410 if(eval_get_sign(arg) == 0)
411 {
412 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
413 res.sign() = false;
414 res.bits() = static_cast<limb_type>(0u);
415 return;
416 }
417 int msb = eval_msb(arg);
418 if(static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) > msb + 1)
419 {
420 // Must have had cancellation in subtraction, shift left and copy:
421 res.bits() = arg;
422 eval_left_shift(res.bits(), cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb - 1);
423 res.exponent() -= static_cast<Exponent>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb - 1);
424 }
425 else if(static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) < msb + 1)
426 {
427 // We have more cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count than we need, so round as required,
428 // first get the rounding bit:
429 bool roundup = eval_bit_test(arg, msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count);
430 // Then check for a tie:
431 if(roundup && (msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count == eval_lsb(arg)))
432 {
433 // Ties round towards even:
434 if(!eval_bit_test(arg, msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1))
435 roundup = false;
436 }
437 // Shift off the cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count we don't need:
438 eval_right_shift(arg, msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1);
439 res.exponent() += static_cast<Exponent>(msb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1);
440 if(roundup)
441 {
442 eval_increment(arg);
443 if(eval_bit_test(arg, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
444 {
445 // This happens very very rairly:
446 eval_right_shift(arg, 1u);
447 ++res.exponent();
448 }
449 }
450 res.bits() = arg;
451 }
452 else
453 {
454 res.bits() = arg;
455 }
456 BOOST_ASSERT((eval_msb(res.bits()) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
457
458 if(res.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
459 {
460 // Overflow:
461 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
462 res.bits() = static_cast<limb_type>(0u);
463 }
464 else if(res.exponent() < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
465 {
466 // Underflow:
467 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
468 res.bits() = static_cast<limb_type>(0u);
469 res.sign() = false;
470 }
471}
472
473template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
474inline void do_eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &b)
475{
476 using default_ops::eval_add;
477 using default_ops::eval_bit_test;
478
479 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
480
481 // Special cases first:
482 switch(a.exponent())
483 {
484 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
485 res = b;
486 if(res.sign())
487 res.negate();
488 return;
489 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
490 if(b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)
491 res = b;
492 else
493 res = a;
494 return; // result is still infinite.
495 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
496 res = a;
497 return; // result is still a NaN.
498 }
499 switch(b.exponent())
500 {
501 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
502 res = a;
503 return;
504 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
505 res = b;
506 if(res.sign())
507 res.negate();
508 return; // result is infinite.
509 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
510 res = b;
511 return; // result is a NaN.
512 }
513
514 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e_diff = a.exponent() - b.exponent();
515 bool s = a.sign();
516 if(e_diff >= 0)
517 {
518 dt = a.bits();
519 if(e_diff < (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
520 {
521 eval_left_shift(dt, e_diff);
522 res.exponent() = a.exponent() - e_diff;
523 eval_add(dt, b.bits());
524 }
525 else
526 res.exponent() = a.exponent();
527 }
528 else
529 {
530 dt= b.bits();
531 if(-e_diff < (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
532 {
533 eval_left_shift(dt, -e_diff);
534 res.exponent() = b.exponent() + e_diff;
535 eval_add(dt, a.bits());
536 }
537 else
538 res.exponent() = b.exponent();
539 }
540
541 copy_and_round(res, dt);
542 res.check_invariants();
543 if(res.sign() != s)
544 res.negate();
545}
546
547template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
548inline void do_eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &b)
549{
550 using default_ops::eval_subtract;
551 using default_ops::eval_bit_test;
552
553 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
554
555 // Special cases first:
556 switch(a.exponent())
557 {
558 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
559 if(b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)
560 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
561 else
562 {
563 res = b;
564 if(!res.sign())
565 res.negate();
566 }
567 return;
568 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
569 if((b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan) || (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity))
570 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
571 else
572 res = a;
573 return;
574 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
575 res = a;
576 return; // result is still a NaN.
577 }
578 switch(b.exponent())
579 {
580 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
581 res = a;
582 return;
583 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
584 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan;
585 res.sign() = false;
586 res.bits() = static_cast<limb_type>(0u);
587 return; // result is a NaN.
588 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
589 res = b;
590 return; // result is still a NaN.
591 }
592
593 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e_diff = a.exponent() - b.exponent();
594 bool s = a.sign();
595 if((e_diff > 0) || ((e_diff == 0) && a.bits().compare(b.bits()) >= 0))
596 {
597 dt = a.bits();
598 if(e_diff < (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
599 {
600 eval_left_shift(dt, e_diff);
601 res.exponent() = a.exponent() - e_diff;
602 eval_subtract(dt, b.bits());
603 }
604 else
605 res.exponent() = a.exponent();
606 }
607 else
608 {
609 dt = b.bits();
610 if(-e_diff < (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
611 {
612 eval_left_shift(dt, -e_diff);
613 res.exponent() = b.exponent() + e_diff;
614 eval_subtract(dt, a.bits());
615 }
616 else
617 res.exponent() = b.exponent();
618 s = !s;
619 }
620
621 copy_and_round(res, dt);
622 if(res.sign() != s)
623 res.negate();
624 res.check_invariants();
625}
626
627template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
628inline void eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &b)
629{
630 if(a.sign() == b.sign())
631 do_eval_add(res, a, b);
632 else
633 do_eval_subtract(res, a, b);
634}
635
636template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
637inline void eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a)
638{
639 return eval_add(res, res, a);
640}
641
642template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
643inline void eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &b)
644{
645 if(a.sign() != b.sign())
646 do_eval_add(res, a, b);
647 else
648 do_eval_subtract(res, a, b);
649}
650
651template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
652inline void eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a)
653{
654 return eval_subtract(res, res, a);
655}
656
657template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
658inline void eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &b)
659{
660 using default_ops::eval_bit_test;
661 using default_ops::eval_multiply;
662
663 // Special cases first:
664 switch(a.exponent())
665 {
666 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
667 if(b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)
668 res = b;
669 else
670 res = a;
671 return;
672 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
673 switch(b.exponent())
674 {
675 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
676 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
677 break;
678 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
679 res = b;
680 break;
681 default:
682 res = a;
683 break;
684 }
685 return;
686 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
687 res = a;
688 return;
689 }
690 if(b.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
691 {
692 res = b;
693 return;
694 }
695 if((a.exponent() > 0) && (b.exponent() > 0))
696 {
697 if(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent + 2 - a.exponent() < b.exponent())
698 {
699 // We will certainly overflow:
700 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
701 res.sign() = a.sign() != b.sign();
702 res.bits() = static_cast<limb_type>(0u);
703 return;
704 }
705 }
706 if((a.exponent() < 0) && (b.exponent() < 0))
707 {
708 if(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - 2 - a.exponent() > b.exponent())
709 {
710 // We will certainly underflow:
711 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
712 res.sign() = false;
713 res.bits() = static_cast<limb_type>(0u);
714 return;
715 }
716 }
717
718 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
719 eval_multiply(dt, a.bits(), b.bits());
720 res.exponent() = a.exponent() + b.exponent() - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
721 copy_and_round(res, dt);
722 res.check_invariants();
723 res.sign() = a.sign() != b.sign();
724}
725
726template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
727inline void eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a)
728{
729 eval_multiply(res, res, a);
730}
731
732template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
733inline typename enable_if_c<is_unsigned<U>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a, const U &b)
734{
735 using default_ops::eval_bit_test;
736 using default_ops::eval_multiply;
737
738 // Special cases first:
739 switch(a.exponent())
740 {
741 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
742 res = a;
743 return;
744 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
745 if(b == 0)
746 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
747 else
748 res = a;
749 return;
750 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
751 res = a;
752 return;
753 }
754
755 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
756 typedef typename boost::multiprecision::detail::canonical<U, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::type canon_ui_type;
757 eval_multiply(dt, a.bits(), static_cast<canon_ui_type>(b));
758 res.exponent() = a.exponent();
759 copy_and_round(res, dt);
760 res.check_invariants();
761 res.sign() = a.sign();
762}
763
764template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
765inline typename enable_if_c<is_unsigned<U>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const U &b)
766{
767 eval_multiply(res, res, b);
768}
769
770template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
771inline typename enable_if_c<is_signed<S>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a, const S &b)
772{
773 typedef typename make_unsigned<S>::type ui_type;
774 eval_multiply(res, a, static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(b)));
775 if(b < 0)
776 res.negate();
777}
778
779template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
780inline typename enable_if_c<is_signed<S>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const S &b)
781{
782 eval_multiply(res, res, b);
783}
784
785template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
786inline void eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &u, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &v)
787{
788#ifdef BOOST_MSVC
789#pragma warning(push)
790#pragma warning(disable:6326) // comparison of two constants
791#endif
792 using default_ops::eval_subtract;
793 using default_ops::eval_qr;
794 using default_ops::eval_bit_test;
795 using default_ops::eval_get_sign;
796 using default_ops::eval_increment;
797
798 //
799 // Special cases first:
800 //
801 switch(u.exponent())
802 {
803 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
804 switch(v.exponent())
805 {
806 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
807 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
808 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
809 return;
810 }
811 res = u;
812 return;
813 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
814 switch(v.exponent())
815 {
816 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
817 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
818 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
819 return;
820 }
821 res = u;
822 return;
823 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
824 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
825 return;
826 }
827 switch(v.exponent())
828 {
829 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
830 {
831 bool s = u.sign() != v.sign();
832 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
833 res.sign() = s;
834 return;
835 }
836 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
837 res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
838 res.bits() = limb_type(0);
839 res.sign() = false;
840 return;
841 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
842 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
843 return;
844 }
845
846 // We can scale u and v so that both are integers, then perform integer
847 // division to obtain quotient q and remainder r, such that:
848 //
849 // q * v + r = u
850 //
851 // and hense:
852 //
853 // q + r/v = u/v
854 //
855 // From this, assuming q has cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count
856 // bits we only need to determine whether
857 // r/v is less than, equal to, or greater than 0.5 to determine rounding -
858 // this we can do with a shift and comparison.
859 //
860 // We can set the exponent and sign of the result up front:
861 //
862 res.exponent() = u.exponent() - v.exponent() - 1;
863 res.sign() = u.sign() != v.sign();
864 //
865 // Now get the quotient and remainder:
866 //
867 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(u.bits()), t2(v.bits()), q, r;
868 eval_left_shift(t, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count);
869 eval_qr(t, t2, q, r);
870 //
871 // We now have either "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count"
872 // or "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1" significant
873 // bits in q.
874 //
875 static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
876 if(eval_bit_test(q, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
877 {
878 //
879 // OK we have cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1 bits,
880 // so we already have rounding info,
881 // we just need to changes things if the last bit is 1 and either the
882 // remainder is non-zero (ie we do not have a tie) or the quotient would
883 // be odd if it were shifted to the correct number of bits (ie a tiebreak).
884 //
885 BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count));
886 if((q.limbs()[0] & 1u) && (eval_get_sign(r) || (q.limbs()[0] & 2u)))
887 {
888 eval_increment(q);
889 }
890 }
891 else
892 {
893 //
894 // We have exactly "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" bits in q.
895 // Get rounding info, which we can get by comparing 2r with v.
896 // We want to call copy_and_round to handle rounding and general cleanup,
897 // so we'll left shift q and add some fake digits on the end to represent
898 // how we'll be rounding.
899 //
900 BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
901 static const unsigned lshift = (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count < limb_bits) ? 2 : limb_bits;
902 eval_left_shift(q, lshift);
903 res.exponent() -= lshift;
904 eval_left_shift(r, 1u);
905 int c = r.compare(v.bits());
906 if(c == 0)
907 q.limbs()[0] |= static_cast<limb_type>(1u) << (lshift - 1);
908 else if(c > 0)
909 q.limbs()[0] |= (static_cast<limb_type>(1u) << (lshift - 1)) + static_cast<limb_type>(1u);
910 }
911 copy_and_round(res, q);
912#ifdef BOOST_MSVC
913#pragma warning(pop)
914#endif
915}
916
917template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
918inline void eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
919{
920 eval_divide(res, res, arg);
921}
922
923template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
924inline typename enable_if_c<is_unsigned<U>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &u, const U &v)
925{
926#ifdef BOOST_MSVC
927#pragma warning(push)
928#pragma warning(disable:6326) // comparison of two constants
929#endif
930 using default_ops::eval_subtract;
931 using default_ops::eval_qr;
932 using default_ops::eval_bit_test;
933 using default_ops::eval_get_sign;
934 using default_ops::eval_increment;
935
936 //
937 // Special cases first:
938 //
939 switch(u.exponent())
940 {
941 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
942 if(v == 0)
943 {
944 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
945 return;
946 }
947 res = u;
948 return;
949 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
950 res = u;
951 return;
952 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
953 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
954 return;
955 }
956 if(v == 0)
957 {
958 bool s = u.sign();
959 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
960 res.sign() = s;
961 return;
962 }
963
964 // We can scale u and v so that both are integers, then perform integer
965 // division to obtain quotient q and remainder r, such that:
966 //
967 // q * v + r = u
968 //
969 // and hense:
970 //
971 // q + r/v = u/v
972 //
973 // From this, assuming q has "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count, we only need to determine whether
974 // r/v is less than, equal to, or greater than 0.5 to determine rounding -
975 // this we can do with a shift and comparison.
976 //
977 // We can set the exponent and sign of the result up front:
978 //
979 int gb = msb(v);
980 res.exponent() = u.exponent() - static_cast<Exponent>(gb) - static_cast<Exponent>(1);
981 res.sign() = u.sign();
982 //
983 // Now get the quotient and remainder:
984 //
985 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(u.bits()), q, r;
986 eval_left_shift(t, gb + 1);
987 eval_qr(t, number<typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::canonical_value(v), q, r);
988 //
989 // We now have either "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" or "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1" significant cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in q.
990 //
991 static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
992 if(eval_bit_test(q, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
993 {
994 //
995 // OK we have cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1 cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count, so we already have rounding info,
996 // we just need to changes things if the last bit is 1 and the
997 // remainder is non-zero (ie we do not have a tie).
998 //
999 BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count));
1000 if((q.limbs()[0] & 1u) && eval_get_sign(r))
1001 {
1002 eval_increment(q);
1003 }
1004 }
1005 else
1006 {
1007 //
1008 // We have exactly "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in q.
1009 // Get rounding info, which we can get by comparing 2r with v.
1010 // We want to call copy_and_round to handle rounding and general cleanup,
1011 // so we'll left shift q and add some fake cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count on the end to represent
1012 // how we'll be rounding.
1013 //
1014 BOOST_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
1015 static const unsigned lshift = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count < limb_bits ? 2 : limb_bits;
1016 eval_left_shift(q, lshift);
1017 res.exponent() -= lshift;
1018 eval_left_shift(r, 1u);
1019 int c = r.compare(number<typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::canonical_value(v));
1020 if(c == 0)
1021 q.limbs()[0] |= static_cast<limb_type>(1u) << (lshift - 1);
1022 else if(c > 0)
1023 q.limbs()[0] |= (static_cast<limb_type>(1u) << (lshift - 1)) + static_cast<limb_type>(1u);
1024 }
1025 copy_and_round(res, q);
1026#ifdef BOOST_MSVC
1027#pragma warning(pop)
1028#endif
1029}
1030
1031template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
1032inline typename enable_if_c<is_unsigned<U>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const U &v)
1033{
1034 eval_divide(res, res, v);
1035}
1036
1037template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
1038inline typename enable_if_c<is_signed<S>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &u, const S &v)
1039{
1040 typedef typename make_unsigned<S>::type ui_type;
1041 eval_divide(res, u, static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(v)));
1042 if(v < 0)
1043 res.negate();
1044}
1045
1046template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
1047inline typename enable_if_c<is_signed<S>::value>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const S &v)
1048{
1049 eval_divide(res, res, v);
1050}
1051
1052template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1053inline int eval_get_sign(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1054{
1055 return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero ? 0 : arg.sign() ? -1 : 1;
1056}
1057
1058template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1059inline bool eval_is_zero(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1060{
1061 return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
1062}
1063
1064template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1065inline bool eval_eq(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &b)
1066{
1067 return (a.exponent() == b.exponent())
1068 && (a.sign() == b.sign())
1069 && (a.bits().compare(b.bits()) == 0)
1070 && (a.exponent() != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan);
1071}
1072
1073template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1074inline void eval_convert_to(boost::long_long_type *res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1075{
1076 switch(arg.exponent())
1077 {
1078 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1079 *res = 0;
1080 return;
1081 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1082 BOOST_THROW_EXCEPTION(std::runtime_error("Could not convert NaN to integer."));
1083 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1084 *res = (std::numeric_limits<boost::long_long_type>::max)();
1085 if(arg.sign())
1086 *res = -*res;
1087 return;
1088 }
1089 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type man(arg.bits());
1090 typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift
1091 = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - arg.exponent();
1092 if(shift > (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
1093 {
1094 *res = 0;
1095 return;
1096 }
1097 if(arg.sign() && (arg.compare((std::numeric_limits<boost::long_long_type>::min)()) <= 0))
1098 {
1099 *res = (std::numeric_limits<boost::long_long_type>::min)();
1100 return;
1101 }
1102 else if(!arg.sign() && (arg.compare((std::numeric_limits<boost::long_long_type>::max)()) >= 0))
1103 {
1104 *res = (std::numeric_limits<boost::long_long_type>::max)();
1105 return;
1106 }
1107 eval_right_shift(man, shift);
1108 eval_convert_to(res, man);
1109 if(arg.sign())
1110 {
1111 *res = -*res;
1112 }
1113}
1114
1115template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1116inline void eval_convert_to(boost::ulong_long_type *res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1117{
1118 switch(arg.exponent())
1119 {
1120 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1121 *res = 0;
1122 return;
1123 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1124 BOOST_THROW_EXCEPTION(std::runtime_error("Could not convert NaN to integer."));
1125 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1126 *res = (std::numeric_limits<boost::ulong_long_type>::max)();
1127 return;
1128 }
1129 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type man(arg.bits());
1130 typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift
1131 = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - arg.exponent();
1132 if(shift > (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
1133 {
1134 *res = 0;
1135 return;
1136 }
1137 else if(shift < 0)
1138 {
1139 // TODO: what if we have fewer cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count than a boost::long_long_type?
1140 *res = (std::numeric_limits<boost::long_long_type>::max)();
1141 return;
1142 }
1143 eval_right_shift(man, shift);
1144 eval_convert_to(res, man);
1145}
1146
1147template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1148inline void eval_convert_to(long double *res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1149{
1150 switch(arg.exponent())
1151 {
1152 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1153 *res = 0;
1154 return;
1155 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1156 *res = std::numeric_limits<long double>::quiet_NaN();
1157 return;
1158 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1159 *res = (std::numeric_limits<long double>::infinity)();
1160 if(arg.sign())
1161 *res = -*res;
1162 return;
1163 }
1164 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type e = arg.exponent();
1165 e -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1;
1166 *res = std::ldexp(static_cast<long double>(*arg.bits().limbs()), e);
1167 for(unsigned i = 1; i < arg.bits().size(); ++i)
1168 {
1169 e += sizeof(*arg.bits().limbs()) * CHAR_BIT;
1170 *res += std::ldexp(static_cast<long double>(arg.bits().limbs()[i]), e);
1171 }
1172 if(arg.sign())
1173 *res = -*res;
1174}
1175
1176template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1177inline void eval_frexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg, Exponent *e)
1178{
1179 switch(arg.exponent())
1180 {
1181 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1182 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1183 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1184 *e = 0;
1185 res = arg;
1186 return;
1187 }
1188 res = arg;
1189 *e = arg.exponent() + 1;
1190 res.exponent() = -1;
1191}
1192
1193template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
1194inline void eval_frexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg, I *pe)
1195{
1196 Exponent e;
1197 eval_frexp(res, arg, &e);
1198 if((e > (std::numeric_limits<I>::max)()) || (e < (std::numeric_limits<I>::min)()))
1199 {
1200 BOOST_THROW_EXCEPTION(std::runtime_error("Exponent was outside of the range of the argument type to frexp."));
1201 }
1202 *pe = static_cast<I>(e);
1203}
1204
1205template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1206inline void eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg, Exponent e)
1207{
1208 switch(arg.exponent())
1209 {
1210 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1211 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1212 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1213 res = arg;
1214 return;
1215 }
1216 if((e > 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent - e < arg.exponent()))
1217 {
1218 // Overflow:
1219 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
1220 res.sign() = arg.sign();
1221 }
1222 else if((e < 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - e > arg.exponent()))
1223 {
1224 // Underflow:
1225 res = limb_type(0);
1226 }
1227 else
1228 {
1229 res = arg;
1230 res.exponent() += e;
1231 }
1232}
1233
1234template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
1235inline typename enable_if_c<is_unsigned<I>::value>::type eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg, I e)
1236{
1237 typedef typename make_signed<I>::type si_type;
1238 if(e > static_cast<I>((std::numeric_limits<si_type>::max)()))
1239 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
1240 else
1241 eval_ldexp(res, arg, static_cast<si_type>(e));
1242}
1243
1244template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
1245inline typename enable_if_c<is_signed<I>::value>::type eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg, I e)
1246{
1247 if((e > (std::numeric_limits<Exponent>::max)()) || (e < (std::numeric_limits<Exponent>::min)()))
1248 {
1249 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
1250 if(e < 0)
1251 res.negate();
1252 }
1253 else
1254 eval_ldexp(res, arg, static_cast<Exponent>(e));
1255}
1256
1257/*
1258* Sign manipulation
1259*/
1260
1261template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1262inline void eval_abs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1263{
1264 res = arg;
1265 res.sign() = false;
1266}
1267
1268template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1269inline void eval_fabs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1270{
1271 res = arg;
1272 res.sign() = false;
1273}
1274
1275template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1276inline int eval_fpclassify(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1277{
1278 switch(arg.exponent())
1279 {
1280 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1281 return FP_ZERO;
1282 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1283 return FP_INFINITE;
1284 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1285 return FP_NAN;
1286 }
1287 return FP_NORMAL;
1288}
1289
1290template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1291inline void eval_sqrt(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1292{
1293 using default_ops::eval_integer_sqrt;
1294 using default_ops::eval_bit_test;
1295 using default_ops::eval_increment;
1296 switch(arg.exponent())
1297 {
1298 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1299 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1300 res = arg;
1301 return;
1302 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1303 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1304 return;
1305 }
1306 if(arg.sign())
1307 {
1308 res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
1309 return;
1310 }
1311
1312 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(arg.bits()), r, s;
1313 eval_left_shift(t, arg.exponent() & 1 ? cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count : cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1);
1314 eval_integer_sqrt(s, r, t);
1315
1316 if(!eval_bit_test(s, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
1317 {
1318 // We have exactly the right number of cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in the result, round as required:
1319 if(s.compare(r) < 0)
1320 {
1321 eval_increment(s);
1322 }
1323 }
1324 typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type ae = arg.exponent();
1325 res.exponent() = ae / 2;
1326 if((ae & 1) && (ae < 0))
1327 --res.exponent();
1328 copy_and_round(res, s);
1329}
1330
1331template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1332inline void eval_floor(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1333{
1334 using default_ops::eval_increment;
1335 switch(arg.exponent())
1336 {
1337 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1338 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1339 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1340 res = arg;
1341 return;
1342 }
1343 typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift =
1344 (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - arg.exponent() - 1;
1345 if((arg.exponent() > (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) || (shift <= 0))
1346 {
1347 // Either arg is already an integer, or a special value:
1348 res = arg;
1349 return;
1350 }
1351 if(shift >= (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
1352 {
1353 res = static_cast<signed_limb_type>(arg.sign() ? -1 : 0);
1354 return;
1355 }
1356 bool fractional = (int)eval_lsb(arg.bits()) < shift;
1357 res = arg;
1358 eval_right_shift(res.bits(), shift);
1359 if(fractional && res.sign())
1360 {
1361 eval_increment(res.bits());
1362 if(eval_msb(res.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - shift)
1363 {
1364 // Must have extended result by one bit in the increment:
1365 --shift;
1366 ++res.exponent();
1367 }
1368 }
1369 eval_left_shift(res.bits(), shift);
1370}
1371
1372template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1373inline void eval_ceil(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
1374{
1375 using default_ops::eval_increment;
1376 switch(arg.exponent())
1377 {
1378 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
1379 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
1380 case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
1381 res = arg;
1382 return;
1383 }
1384 typename mpl::if_c<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - arg.exponent() - 1;
1385 if((arg.exponent() > (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) || (shift <= 0))
1386 {
1387 // Either arg is already an integer, or a special value:
1388 res = arg;
1389 return;
1390 }
1391 if(shift >= (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
1392 {
1393 res = static_cast<signed_limb_type>(arg.sign() ? 0 : 1);
1394 return;
1395 }
1396 bool fractional = (int)eval_lsb(arg.bits()) < shift;
1397 res = arg;
1398 eval_right_shift(res.bits(), shift);
1399 if(fractional && !res.sign())
1400 {
1401 eval_increment(res.bits());
1402 if(eval_msb(res.bits()) != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - shift)
1403 {
1404 // Must have extended result by one bit in the increment:
1405 --shift;
1406 ++res.exponent();
1407 }
1408 }
1409 eval_left_shift(res.bits(), shift);
1410}
1411
1412} // namespace backends
1413
1414#ifdef BOOST_NO_SFINAE_EXPR
1415
1416namespace detail{
1417
1418template<unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2, unsigned D2, backends::digit_base_type B2, class A2, class E2, E2 M3, E2 M4>
1419struct is_explicitly_convertible<backends::cpp_bin_float<D1, B1, A1, E1, M1, M2>, backends::cpp_bin_float<D2, B2, A2, E2, M3, M4> > : public mpl::true_ {};
1420template<class FloatT, unsigned D2, backends::digit_base_type B2, class A2, class E2, E2 M3, E2 M4>
1421struct is_explicitly_convertible<FloatT, backends::cpp_bin_float<D2, B2, A2, E2, M3, M4> > : public boost::is_floating_point<FloatT> {};
1422
1423}
1424#endif
1425
1426
1427using backends::cpp_bin_float;
1428using backends::digit_base_2;
1429using backends::digit_base_10;
1430
1431template<unsigned Digits, backends::digit_base_type DigitBase, class Exponent, Exponent MinE, Exponent MaxE, class Allocator>
1432struct number_category<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > : public boost::mpl::int_<boost::multiprecision::number_kind_floating_point>{};
1433
1434template<unsigned Digits, backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
1435struct expression_template_default<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >
1436{
1437 static const expression_template_option value = is_void<Allocator>::value ? et_off : et_on;
1438};
1439
1440typedef number<backends::cpp_bin_float<50> > cpp_bin_float_50;
1441typedef number<backends::cpp_bin_float<100> > cpp_bin_float_100;
1442
1443typedef number<backends::cpp_bin_float<24, backends::digit_base_2, void, boost::int16_t, -126, 127>, et_off> cpp_bin_float_single;
1444typedef number<backends::cpp_bin_float<53, backends::digit_base_2, void, boost::int16_t, -1022, 1023>, et_off> cpp_bin_float_double;
1445typedef number<backends::cpp_bin_float<64, backends::digit_base_2, void, boost::int16_t, -16382, 16383>, et_off> cpp_bin_float_double_extended;
1446typedef number<backends::cpp_bin_float<113, backends::digit_base_2, void, boost::int16_t, -16382, 16383>, et_off> cpp_bin_float_quad;
1447
1448}} // namespaces
1449
1450#include <boost/multiprecision/cpp_bin_float/io.hpp>
1451#include <boost/multiprecision/cpp_bin_float/transcendental.hpp>
1452
1453namespace std{
1454
1455//
1456// numeric_limits [partial] specializations for the types declared in this header:
1457//
1458template<unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1459class numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >
1460{
1461 typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> number_type;
1462public:
1463 BOOST_STATIC_CONSTEXPR bool is_specialized = true;
1464 static number_type (min)()
1465 {
1466 initializer.do_nothing();
1467 static std::pair<bool, number_type> value;
1468 if(!value.first)
1469 {
1470 value.first = true;
1471 value.second = 1u;
1472 value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
1473 }
1474 return value.second;
1475 }
1476 static number_type (max)()
1477 {
1478 initializer.do_nothing();
1479 static std::pair<bool, number_type> value;
1480 if(!value.first)
1481 {
1482 value.first = true;
1483 eval_complement(value.second.backend().bits(), value.second.backend().bits());
1484 value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
1485 }
1486 return value.second;
1487 }
1488 BOOST_STATIC_CONSTEXPR number_type lowest()
1489 {
1490 return -(max)();
1491 }
1492 BOOST_STATIC_CONSTEXPR int digits = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count;
1493 BOOST_STATIC_CONSTEXPR int digits10 = (digits - 1) * 301 / 1000;
1494 // Is this really correct???
1495 BOOST_STATIC_CONSTEXPR int max_digits10 = (digits * 301 / 1000) + 3;
1496 BOOST_STATIC_CONSTEXPR bool is_signed = true;
1497 BOOST_STATIC_CONSTEXPR bool is_integer = false;
1498 BOOST_STATIC_CONSTEXPR bool is_exact = false;
1499 BOOST_STATIC_CONSTEXPR int radix = 2;
1500 static number_type epsilon()
1501 {
1502 initializer.do_nothing();
1503 static std::pair<bool, number_type> value;
1504 if(!value.first)
1505 {
1506 value.first = true;
1507 value.second = 1;
1508 value.second = ldexp(value.second, 1 - (int)digits);
1509 }
1510 return value.second;
1511 }
1512 // What value should this be????
1513 static number_type round_error()
1514 {
1515 // returns 0.5
1516 initializer.do_nothing();
1517 static std::pair<bool, number_type> value;
1518 if(!value.first)
1519 {
1520 value.first = true;
1521 value.second = 1;
1522 value.second = ldexp(value.second, -1);
1523 }
1524 return value.second;
1525 }
1526 BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type min_exponent = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
1527 BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type min_exponent10 = (min_exponent / 1000) * 301L;
1528 BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type max_exponent = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
1529 BOOST_STATIC_CONSTEXPR typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type max_exponent10 = (max_exponent / 1000) * 301L;
1530 BOOST_STATIC_CONSTEXPR bool has_infinity = true;
1531 BOOST_STATIC_CONSTEXPR bool has_quiet_NaN = true;
1532 BOOST_STATIC_CONSTEXPR bool has_signaling_NaN = false;
1533 BOOST_STATIC_CONSTEXPR float_denorm_style has_denorm = denorm_absent;
1534 BOOST_STATIC_CONSTEXPR bool has_denorm_loss = false;
1535 static number_type infinity()
1536 {
1537 // returns epsilon/2
1538 initializer.do_nothing();
1539 static std::pair<bool, number_type> value;
1540 if(!value.first)
1541 {
1542 value.first = true;
1543 value.second.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
1544 }
1545 return value.second;
1546 }
1547 static number_type quiet_NaN()
1548 {
1549 return number_type();
1550 }
1551 BOOST_STATIC_CONSTEXPR number_type signaling_NaN()
1552 {
1553 return number_type(0);
1554 }
1555 BOOST_STATIC_CONSTEXPR number_type denorm_min() { return number_type(0); }
1556 BOOST_STATIC_CONSTEXPR bool is_iec559 = false;
1557 BOOST_STATIC_CONSTEXPR bool is_bounded = true;
1558 BOOST_STATIC_CONSTEXPR bool is_modulo = false;
1559 BOOST_STATIC_CONSTEXPR bool traps = true;
1560 BOOST_STATIC_CONSTEXPR bool tinyness_before = false;
1561 BOOST_STATIC_CONSTEXPR float_round_style round_style = round_to_nearest;
1562private:
1563 struct data_initializer
1564 {
1565 data_initializer()
1566 {
1567 std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::epsilon();
1568 std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::round_error();
1569 (std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::min)();
1570 (std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::max)();
1571 std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity();
1572 std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN();
1573 }
1574 void do_nothing()const{}
1575 };
1576 static const data_initializer initializer;
1577};
1578
1579template<unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1580const typename numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::data_initializer numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::initializer;
1581
1582#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
1583
1584template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1585BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::digits;
1586template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1587BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::digits10;
1588template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1589BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_digits10;
1590template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1591BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_signed;
1592template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1593BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_integer;
1594template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1595BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_exact;
1596template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1597BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::radix;
1598template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1599BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::min_exponent;
1600template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1601BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::min_exponent10;
1602template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1603BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_exponent;
1604template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1605BOOST_CONSTEXPR_OR_CONST typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_exponent10;
1606template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1607BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_infinity;
1608template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1609BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_quiet_NaN;
1610template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1611BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_signaling_NaN;
1612template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1613BOOST_CONSTEXPR_OR_CONST float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_denorm;
1614template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1615BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_denorm_loss;
1616template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1617BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_iec559;
1618template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1619BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_bounded;
1620template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1621BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_modulo;
1622template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1623BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::traps;
1624template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1625BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::tinyness_before;
1626template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
1627BOOST_CONSTEXPR_OR_CONST float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::round_style;
1628
1629#endif
1630
1631} // namespace std
1632
1633#endif
1634

source code of boost/boost/multiprecision/cpp_bin_float.hpp