1 ///////////////////////////////////////////////////////////////////////////////
2 // Copyright 2011 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_0.txt)
5
6 #ifndef BOOST_MATH_BIG_NUM_DEF_OPS
7 #define BOOST_MATH_BIG_NUM_DEF_OPS
8
9 #include <boost/math/policies/error_handling.hpp>
10 #include <boost/multiprecision/detail/number_base.hpp>
11 #include <boost/math/special_functions/fpclassify.hpp>
12 #include <boost/math/special_functions/next.hpp>
13 #include <boost/utility/enable_if.hpp>
14 #include <boost/mpl/front.hpp>
15 #include <boost/mpl/fold.hpp>
16 #include <boost/cstdint.hpp>
17 #include <boost/type_traits/make_unsigned.hpp>
18
19 #ifndef INSTRUMENT_BACKEND
20 #ifndef BOOST_MP_INSTRUMENT
21 #define INSTRUMENT_BACKEND(x)
22 #else
23 #define INSTRUMENT_BACKEND(x)\
24 std::cout << BOOST_STRINGIZE(x) << " = " << x.str(0, std::ios_base::scientific) << std::endl;
25 #endif
26 #endif
27
28
29 namespace boost{ namespace multiprecision{
30
31 namespace detail {
32
33 template <class T>
34 struct is_backend;
35
36 template <class To, class From>
37 void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_floating_point>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/);
38 template <class To, class From>
39 void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_integer>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/);
40 template <class To, class From>
41 void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_floating_point>& /*to_type*/, const mpl::int_<number_kind_floating_point>& /*from_type*/);
42 template <class To, class From>
43 void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_rational>& /*from_type*/);
44 template <class To, class From>
45 void generic_interconvert(To& to, const From& from, const mpl::int_<number_kind_rational>& /*to_type*/, const mpl::int_<number_kind_integer>& /*from_type*/);
46
47 }
48
49 namespace default_ops{
50
51 #ifdef BOOST_MSVC
52 // warning C4127: conditional expression is constant
53 #pragma warning(push)
54 #pragma warning(disable:4127)
55 #endif
56 //
57 // Default versions of mixed arithmetic, these just construct a temporary
58 // from the arithmetic value and then do the arithmetic on that, two versions
59 // of each depending on whether the backend can be directly constructed from type V.
60 //
61 // Note that we have to provide *all* the template parameters to class number when used in
62 // enable_if as MSVC-10 won't compile the code if we rely on a computed-default parameter.
63 // Since the result of the test doesn't depend on whether expression templates are on or off
64 // we just use et_on everywhere. We could use a BOOST_WORKAROUND but that just obfuscates the
65 // code even more....
66 //
67 template <class T, class V>
68 inline typename disable_if_c<is_convertible<V, T>::value >::type
eval_add(T & result,V const & v)69 eval_add(T& result, V const& v)
70 {
71 T t;
72 t = v;
73 eval_add(result, t);
74 }
75 template <class T, class V>
76 inline typename enable_if_c<is_convertible<V, T>::value >::type
eval_add(T & result,V const & v)77 eval_add(T& result, V const& v)
78 {
79 T t(v);
80 eval_add(result, t);
81 }
82 template <class T, class V>
83 inline typename disable_if_c<is_convertible<V, T>::value>::type
eval_subtract(T & result,V const & v)84 eval_subtract(T& result, V const& v)
85 {
86 T t;
87 t = v;
88 eval_subtract(result, t);
89 }
90 template <class T, class V>
91 inline typename enable_if_c<is_convertible<V, T>::value>::type
eval_subtract(T & result,V const & v)92 eval_subtract(T& result, V const& v)
93 {
94 T t(v);
95 eval_subtract(result, t);
96 }
97 template <class T, class V>
98 inline typename disable_if_c<is_convertible<V, T>::value>::type
eval_multiply(T & result,V const & v)99 eval_multiply(T& result, V const& v)
100 {
101 T t;
102 t = v;
103 eval_multiply(result, t);
104 }
105 template <class T, class V>
106 inline typename enable_if_c<is_convertible<V, T>::value>::type
eval_multiply(T & result,V const & v)107 eval_multiply(T& result, V const& v)
108 {
109 T t(v);
110 eval_multiply(result, t);
111 }
112
113 template <class T, class U, class V>
114 void eval_multiply(T& t, const U& u, const V& v);
115
116 template <class T, class U, class V>
eval_multiply_add(T & t,const U & u,const V & v)117 inline typename disable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_add(T& t, const U& u, const V& v)
118 {
119 T z;
120 eval_multiply(z, u, v);
121 eval_add(t, z);
122 }
123 template <class T, class U, class V>
eval_multiply_add(T & t,const U & u,const V & v)124 inline typename enable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_add(T& t, const U& u, const V& v)
125 {
126 eval_multiply_add(t, v, u);
127 }
128 template <class T, class U, class V>
eval_multiply_subtract(T & t,const U & u,const V & v)129 inline typename disable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_subtract(T& t, const U& u, const V& v)
130 {
131 T z;
132 eval_multiply(z, u, v);
133 eval_subtract(t, z);
134 }
135 template <class T, class U, class V>
eval_multiply_subtract(T & t,const U & u,const V & v)136 inline typename enable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_subtract(T& t, const U& u, const V& v)
137 {
138 eval_multiply_subtract(t, v, u);
139 }
140 template <class T, class V>
141 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && !is_convertible<V, T>::value>::type
eval_divide(T & result,V const & v)142 eval_divide(T& result, V const& v)
143 {
144 T t;
145 t = v;
146 eval_divide(result, t);
147 }
148 template <class T, class V>
149 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && is_convertible<V, T>::value>::type
eval_divide(T & result,V const & v)150 eval_divide(T& result, V const& v)
151 {
152 T t(v);
153 eval_divide(result, t);
154 }
155 template <class T, class V>
156 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && !is_convertible<V, T>::value>::type
eval_modulus(T & result,V const & v)157 eval_modulus(T& result, V const& v)
158 {
159 T t;
160 t = v;
161 eval_modulus(result, t);
162 }
163 template <class T, class V>
164 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value&& is_convertible<V, T>::value>::type
eval_modulus(T & result,V const & v)165 eval_modulus(T& result, V const& v)
166 {
167 T t(v);
168 eval_modulus(result, t);
169 }
170 template <class T, class V>
171 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && !is_convertible<V, T>::value>::type
eval_bitwise_and(T & result,V const & v)172 eval_bitwise_and(T& result, V const& v)
173 {
174 T t;
175 t = v;
176 eval_bitwise_and(result, t);
177 }
178 template <class T, class V>
179 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && is_convertible<V, T>::value>::type
eval_bitwise_and(T & result,V const & v)180 eval_bitwise_and(T& result, V const& v)
181 {
182 T t(v);
183 eval_bitwise_and(result, t);
184 }
185 template <class T, class V>
186 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && !is_convertible<V, T>::value>::type
eval_bitwise_or(T & result,V const & v)187 eval_bitwise_or(T& result, V const& v)
188 {
189 T t;
190 t = v;
191 eval_bitwise_or(result, t);
192 }
193 template <class T, class V>
194 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && is_convertible<V, T>::value>::type
eval_bitwise_or(T & result,V const & v)195 eval_bitwise_or(T& result, V const& v)
196 {
197 T t(v);
198 eval_bitwise_or(result, t);
199 }
200 template <class T, class V>
201 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && !is_convertible<V, T>::value>::type
eval_bitwise_xor(T & result,V const & v)202 eval_bitwise_xor(T& result, V const& v)
203 {
204 T t;
205 t = v;
206 eval_bitwise_xor(result, t);
207 }
208 template <class T, class V>
209 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && is_convertible<V, T>::value>::type
eval_bitwise_xor(T & result,V const & v)210 eval_bitwise_xor(T& result, V const& v)
211 {
212 T t(v);
213 eval_bitwise_xor(result, t);
214 }
215
216 template <class T, class V>
217 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && !is_convertible<V, T>::value>::type
eval_complement(T & result,V const & v)218 eval_complement(T& result, V const& v)
219 {
220 T t;
221 t = v;
222 eval_complement(result, t);
223 }
224 template <class T, class V>
225 inline typename enable_if_c<is_convertible<V, number<T, et_on> >::value && is_convertible<V, T>::value>::type
eval_complement(T & result,V const & v)226 eval_complement(T& result, V const& v)
227 {
228 T t(v);
229 eval_complement(result, t);
230 }
231
232 //
233 // Default versions of 3-arg arithmetic functions, these mostly just forward to the 2 arg versions:
234 //
235 template <class T, class U, class V>
236 void eval_add(T& t, const U& u, const V& v);
237
238 template <class T>
eval_add_default(T & t,const T & u,const T & v)239 inline void eval_add_default(T& t, const T& u, const T& v)
240 {
241 if(&t == &v)
242 {
243 eval_add(t, u);
244 }
245 else if(&t == &u)
246 {
247 eval_add(t, v);
248 }
249 else
250 {
251 t = u;
252 eval_add(t, v);
253 }
254 }
255 template <class T, class U>
eval_add_default(T & t,const T & u,const U & v)256 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value>::type eval_add_default(T& t, const T& u, const U& v)
257 {
258 T vv;
259 vv = v;
260 eval_add(t, u, vv);
261 }
262 template <class T, class U>
eval_add_default(T & t,const T & u,const U & v)263 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value>::type eval_add_default(T& t, const T& u, const U& v)
264 {
265 T vv(v);
266 eval_add(t, u, vv);
267 }
268 template <class T, class U>
eval_add_default(T & t,const U & u,const T & v)269 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value>::type eval_add_default(T& t, const U& u, const T& v)
270 {
271 eval_add(t, v, u);
272 }
273 template <class T, class U, class V>
eval_add_default(T & t,const U & u,const V & v)274 inline void eval_add_default(T& t, const U& u, const V& v)
275 {
276 if(is_same<T, V>::value && ((void*)&t == (void*)&v))
277 {
278 eval_add(t, u);
279 }
280 else
281 {
282 t = u;
283 eval_add(t, v);
284 }
285 }
286 template <class T, class U, class V>
eval_add(T & t,const U & u,const V & v)287 inline void eval_add(T& t, const U& u, const V& v)
288 {
289 eval_add_default(t, u, v);
290 }
291
292 template <class T, class U, class V>
293 void eval_subtract(T& t, const U& u, const V& v);
294
295 template <class T>
eval_subtract_default(T & t,const T & u,const T & v)296 inline void eval_subtract_default(T& t, const T& u, const T& v)
297 {
298 if((&t == &v) && is_signed_number<T>::value)
299 {
300 eval_subtract(t, u);
301 t.negate();
302 }
303 else if(&t == &u)
304 {
305 eval_subtract(t, v);
306 }
307 else
308 {
309 t = u;
310 eval_subtract(t, v);
311 }
312 }
313 template <class T, class U>
eval_subtract_default(T & t,const T & u,const U & v)314 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value>::type eval_subtract_default(T& t, const T& u, const U& v)
315 {
316 T vv;
317 vv = v;
318 eval_subtract(t, u, vv);
319 }
320 template <class T, class U>
eval_subtract_default(T & t,const T & u,const U & v)321 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value>::type eval_subtract_default(T& t, const T& u, const U& v)
322 {
323 T vv(v);
324 eval_subtract(t, u, vv);
325 }
326 template <class T, class U>
eval_subtract_default(T & t,const U & u,const T & v)327 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_signed_number<T>::value>::type eval_subtract_default(T& t, const U& u, const T& v)
328 {
329 eval_subtract(t, v, u);
330 t.negate();
331 }
332 template <class T, class U>
eval_subtract_default(T & t,const U & u,const T & v)333 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value && is_unsigned_number<T>::value>::type eval_subtract_default(T& t, const U& u, const T& v)
334 {
335 T temp;
336 temp = u;
337 eval_subtract(t, temp, v);
338 }
339 template <class T, class U>
eval_subtract_default(T & t,const U & u,const T & v)340 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value && is_unsigned_number<T>::value>::type eval_subtract_default(T& t, const U& u, const T& v)
341 {
342 T temp(u);
343 eval_subtract(t, temp, v);
344 }
345 template <class T, class U, class V>
eval_subtract_default(T & t,const U & u,const V & v)346 inline void eval_subtract_default(T& t, const U& u, const V& v)
347 {
348 if(is_same<T, V>::value && ((void*)&t == (void*)&v))
349 {
350 eval_subtract(t, u);
351 t.negate();
352 }
353 else
354 {
355 t = u;
356 eval_subtract(t, v);
357 }
358 }
359 template <class T, class U, class V>
eval_subtract(T & t,const U & u,const V & v)360 inline void eval_subtract(T& t, const U& u, const V& v)
361 {
362 eval_subtract_default(t, u, v);
363 }
364
365 template <class T>
eval_multiply_default(T & t,const T & u,const T & v)366 inline void eval_multiply_default(T& t, const T& u, const T& v)
367 {
368 if(&t == &v)
369 {
370 eval_multiply(t, u);
371 }
372 else if(&t == &u)
373 {
374 eval_multiply(t, v);
375 }
376 else
377 {
378 t = u;
379 eval_multiply(t, v);
380 }
381 }
382 #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
383 template <class T, class U>
eval_multiply_default(T & t,const T & u,const U & v)384 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value>::type eval_multiply_default(T& t, const T& u, const U& v)
385 {
386 T vv;
387 vv = v;
388 eval_multiply(t, u, vv);
389 }
390 template <class T, class U>
eval_multiply_default(T & t,const T & u,const U & v)391 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value>::type eval_multiply_default(T& t, const T& u, const U& v)
392 {
393 T vv(v);
394 eval_multiply(t, u, vv);
395 }
396 template <class T, class U>
eval_multiply_default(T & t,const U & u,const T & v)397 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value>::type eval_multiply_default(T& t, const U& u, const T& v)
398 {
399 eval_multiply(t, v, u);
400 }
401 #endif
402 template <class T, class U, class V>
eval_multiply_default(T & t,const U & u,const V & v)403 inline void eval_multiply_default(T& t, const U& u, const V& v)
404 {
405 if(is_same<T, V>::value && ((void*)&t == (void*)&v))
406 {
407 eval_multiply(t, u);
408 }
409 else
410 {
411 t = number<T>::canonical_value(u);
412 eval_multiply(t, v);
413 }
414 }
415 template <class T, class U, class V>
eval_multiply(T & t,const U & u,const V & v)416 inline void eval_multiply(T& t, const U& u, const V& v)
417 {
418 eval_multiply_default(t, u, v);
419 }
420
421 template <class T>
eval_multiply_add(T & t,const T & u,const T & v,const T & x)422 inline void eval_multiply_add(T& t, const T& u, const T& v, const T& x)
423 {
424 if((void*)&x == (void*)&t)
425 {
426 T z;
427 z = number<T>::canonical_value(x);
428 eval_multiply_add(t, u, v, z);
429 }
430 else
431 {
432 eval_multiply(t, u, v);
433 eval_add(t, x);
434 }
435 }
436
437 template <class T, class U>
make_T(const U & u)438 inline typename boost::disable_if_c<boost::is_same<T, U>::value, T>::type make_T(const U& u)
439 {
440 T t;
441 t = number<T>::canonical_value(u);
442 return BOOST_MP_MOVE(t);
443 }
444 template <class T>
make_T(const T & t)445 inline const T& make_T(const T& t)
446 {
447 return t;
448 }
449
450 template <class T, class U, class V, class X>
eval_multiply_add(T & t,const U & u,const V & v,const X & x)451 inline typename disable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_add(T& t, const U& u, const V& v, const X& x)
452 {
453 eval_multiply_add(t, make_T<T>(u), make_T<T>(v), make_T<T>(x));
454 }
455 template <class T, class U, class V, class X>
eval_multiply_add(T & t,const U & u,const V & v,const X & x)456 inline typename enable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_add(T& t, const U& u, const V& v, const X& x)
457 {
458 eval_multiply_add(t, v, u, x);
459 }
460 template <class T, class U, class V, class X>
eval_multiply_subtract(T & t,const U & u,const V & v,const X & x)461 inline typename disable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_subtract(T& t, const U& u, const V& v, const X& x)
462 {
463 if((void*)&x == (void*)&t)
464 {
465 T z;
466 z = x;
467 eval_multiply_subtract(t, u, v, z);
468 }
469 else
470 {
471 eval_multiply(t, u, v);
472 eval_subtract(t, x);
473 }
474 }
475 template <class T, class U, class V, class X>
eval_multiply_subtract(T & t,const U & u,const V & v,const X & x)476 inline typename enable_if_c<!is_same<T, U>::value && is_same<T, V>::value>::type eval_multiply_subtract(T& t, const U& u, const V& v, const X& x)
477 {
478 eval_multiply_subtract(t, v, u, x);
479 }
480
481 template <class T, class U, class V>
482 void eval_divide(T& t, const U& u, const V& v);
483
484 template <class T>
eval_divide_default(T & t,const T & u,const T & v)485 inline void eval_divide_default(T& t, const T& u, const T& v)
486 {
487 if(&t == &u)
488 eval_divide(t, v);
489 else if(&t == &v)
490 {
491 T temp;
492 eval_divide(temp, u, v);
493 temp.swap(t);
494 }
495 else
496 {
497 t = u;
498 eval_divide(t, v);
499 }
500 }
501 #if !BOOST_WORKAROUND(BOOST_MSVC, < 1900)
502 template <class T, class U>
eval_divide_default(T & t,const T & u,const U & v)503 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value>::type eval_divide_default(T& t, const T& u, const U& v)
504 {
505 T vv;
506 vv = v;
507 eval_divide(t, u, vv);
508 }
509 template <class T, class U>
eval_divide_default(T & t,const T & u,const U & v)510 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value>::type eval_divide_default(T& t, const T& u, const U& v)
511 {
512 T vv(v);
513 eval_divide(t, u, vv);
514 }
515 template <class T, class U>
eval_divide_default(T & t,const U & u,const T & v)516 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value>::type eval_divide_default(T& t, const U& u, const T& v)
517 {
518 T uu;
519 uu = u;
520 eval_divide(t, uu, v);
521 }
522 template <class T, class U>
eval_divide_default(T & t,const U & u,const T & v)523 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value>::type eval_divide_default(T& t, const U& u, const T& v)
524 {
525 T uu(u);
526 eval_divide(t, uu, v);
527 }
528 #endif
529 template <class T, class U, class V>
eval_divide_default(T & t,const U & u,const V & v)530 inline void eval_divide_default(T& t, const U& u, const V& v)
531 {
532 if(is_same<T, V>::value && ((void*)&t == (void*)&v))
533 {
534 T temp;
535 temp = u;
536 eval_divide(temp, v);
537 t = temp;
538 }
539 else
540 {
541 t = u;
542 eval_divide(t, v);
543 }
544 }
545 template <class T, class U, class V>
eval_divide(T & t,const U & u,const V & v)546 inline void eval_divide(T& t, const U& u, const V& v)
547 {
548 eval_divide_default(t, u, v);
549 }
550
551 template <class T, class U, class V>
552 void eval_modulus(T& t, const U& u, const V& v);
553
554 template <class T>
eval_modulus_default(T & t,const T & u,const T & v)555 inline void eval_modulus_default(T& t, const T& u, const T& v)
556 {
557 if(&t == &u)
558 eval_modulus(t, v);
559 else if(&t == &v)
560 {
561 T temp;
562 eval_modulus(temp, u, v);
563 temp.swap(t);
564 }
565 else
566 {
567 t = u;
568 eval_modulus(t, v);
569 }
570 }
571 template <class T, class U>
eval_modulus_default(T & t,const T & u,const U & v)572 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value>::type eval_modulus_default(T& t, const T& u, const U& v)
573 {
574 T vv;
575 vv = v;
576 eval_modulus(t, u, vv);
577 }
578 template <class T, class U>
eval_modulus_default(T & t,const T & u,const U & v)579 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value>::type eval_modulus_default(T& t, const T& u, const U& v)
580 {
581 T vv(v);
582 eval_modulus(t, u, vv);
583 }
584 template <class T, class U>
eval_modulus_default(T & t,const U & u,const T & v)585 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value>::type eval_modulus_default(T& t, const U& u, const T& v)
586 {
587 T uu;
588 uu = u;
589 eval_modulus(t, uu, v);
590 }
591 template <class T, class U>
eval_modulus_default(T & t,const U & u,const T & v)592 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value>::type eval_modulus_default(T& t, const U& u, const T& v)
593 {
594 T uu(u);
595 eval_modulus(t, uu, v);
596 }
597 template <class T, class U, class V>
eval_modulus_default(T & t,const U & u,const V & v)598 inline void eval_modulus_default(T& t, const U& u, const V& v)
599 {
600 if(is_same<T, V>::value && ((void*)&t == (void*)&v))
601 {
602 T temp(u);
603 eval_modulus(temp, v);
604 t = temp;
605 }
606 else
607 {
608 t = u;
609 eval_modulus(t, v);
610 }
611 }
612 template <class T, class U, class V>
eval_modulus(T & t,const U & u,const V & v)613 inline void eval_modulus(T& t, const U& u, const V& v)
614 {
615 eval_modulus_default(t, u, v);
616 }
617
618 template <class T, class U, class V>
619 void eval_bitwise_and(T& t, const U& u, const V& v);
620
621 template <class T>
eval_bitwise_and_default(T & t,const T & u,const T & v)622 inline void eval_bitwise_and_default(T& t, const T& u, const T& v)
623 {
624 if(&t == &v)
625 {
626 eval_bitwise_and(t, u);
627 }
628 else if(&t == &u)
629 {
630 eval_bitwise_and(t, v);
631 }
632 else
633 {
634 t = u;
635 eval_bitwise_and(t, v);
636 }
637 }
638 template <class T, class U>
eval_bitwise_and_default(T & t,const T & u,const U & v)639 inline typename disable_if_c<is_convertible<U, T>::value>::type eval_bitwise_and_default(T& t, const T& u, const U& v)
640 {
641 T vv;
642 vv = v;
643 eval_bitwise_and(t, u, vv);
644 }
645 template <class T, class U>
eval_bitwise_and_default(T & t,const T & u,const U & v)646 inline typename enable_if_c<is_convertible<U, T>::value>::type eval_bitwise_and_default(T& t, const T& u, const U& v)
647 {
648 T vv(v);
649 eval_bitwise_and(t, u, vv);
650 }
651 template <class T, class U>
eval_bitwise_and_default(T & t,const U & u,const T & v)652 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value>::type eval_bitwise_and_default(T& t, const U& u, const T& v)
653 {
654 eval_bitwise_and(t, v, u);
655 }
656 template <class T, class U, class V>
eval_bitwise_and_default(T & t,const U & u,const V & v)657 inline typename disable_if_c<is_same<T, U>::value || is_same<T, V>::value>::type eval_bitwise_and_default(T& t, const U& u, const V& v)
658 {
659 t = u;
660 eval_bitwise_and(t, v);
661 }
662 template <class T, class U, class V>
eval_bitwise_and(T & t,const U & u,const V & v)663 inline void eval_bitwise_and(T& t, const U& u, const V& v)
664 {
665 eval_bitwise_and_default(t, u, v);
666 }
667
668 template <class T, class U, class V>
669 void eval_bitwise_or(T& t, const U& u, const V& v);
670
671 template <class T>
eval_bitwise_or_default(T & t,const T & u,const T & v)672 inline void eval_bitwise_or_default(T& t, const T& u, const T& v)
673 {
674 if(&t == &v)
675 {
676 eval_bitwise_or(t, u);
677 }
678 else if(&t == &u)
679 {
680 eval_bitwise_or(t, v);
681 }
682 else
683 {
684 t = u;
685 eval_bitwise_or(t, v);
686 }
687 }
688 template <class T, class U>
eval_bitwise_or_default(T & t,const T & u,const U & v)689 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value>::type eval_bitwise_or_default(T& t, const T& u, const U& v)
690 {
691 T vv;
692 vv = v;
693 eval_bitwise_or(t, u, vv);
694 }
695 template <class T, class U>
eval_bitwise_or_default(T & t,const T & u,const U & v)696 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value>::type eval_bitwise_or_default(T& t, const T& u, const U& v)
697 {
698 T vv(v);
699 eval_bitwise_or(t, u, vv);
700 }
701 template <class T, class U>
eval_bitwise_or_default(T & t,const U & u,const T & v)702 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value>::type eval_bitwise_or_default(T& t, const U& u, const T& v)
703 {
704 eval_bitwise_or(t, v, u);
705 }
706 template <class T, class U, class V>
eval_bitwise_or_default(T & t,const U & u,const V & v)707 inline void eval_bitwise_or_default(T& t, const U& u, const V& v)
708 {
709 if(is_same<T, V>::value && ((void*)&t == (void*)&v))
710 {
711 eval_bitwise_or(t, u);
712 }
713 else
714 {
715 t = u;
716 eval_bitwise_or(t, v);
717 }
718 }
719 template <class T, class U, class V>
eval_bitwise_or(T & t,const U & u,const V & v)720 inline void eval_bitwise_or(T& t, const U& u, const V& v)
721 {
722 eval_bitwise_or_default(t, u, v);
723 }
724
725 template <class T, class U, class V>
726 void eval_bitwise_xor(T& t, const U& u, const V& v);
727
728 template <class T>
eval_bitwise_xor_default(T & t,const T & u,const T & v)729 inline void eval_bitwise_xor_default(T& t, const T& u, const T& v)
730 {
731 if(&t == &v)
732 {
733 eval_bitwise_xor(t, u);
734 }
735 else if(&t == &u)
736 {
737 eval_bitwise_xor(t, v);
738 }
739 else
740 {
741 t = u;
742 eval_bitwise_xor(t, v);
743 }
744 }
745 template <class T, class U>
eval_bitwise_xor_default(T & t,const T & u,const U & v)746 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && !is_convertible<U, T>::value>::type eval_bitwise_xor_default(T& t, const T& u, const U& v)
747 {
748 T vv;
749 vv = v;
750 eval_bitwise_xor(t, u, vv);
751 }
752 template <class T, class U>
eval_bitwise_xor_default(T & t,const T & u,const U & v)753 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value && is_convertible<U, T>::value>::type eval_bitwise_xor_default(T& t, const T& u, const U& v)
754 {
755 T vv(v);
756 eval_bitwise_xor(t, u, vv);
757 }
758 template <class T, class U>
eval_bitwise_xor_default(T & t,const U & u,const T & v)759 inline typename enable_if_c<is_convertible<U, number<T, et_on> >::value>::type eval_bitwise_xor_default(T& t, const U& u, const T& v)
760 {
761 eval_bitwise_xor(t, v, u);
762 }
763 template <class T, class U, class V>
eval_bitwise_xor_default(T & t,const U & u,const V & v)764 inline void eval_bitwise_xor_default(T& t, const U& u, const V& v)
765 {
766 if(is_same<T, V>::value && ((void*)&t == (void*)&v))
767 {
768 eval_bitwise_xor(t, u);
769 }
770 else
771 {
772 t = u;
773 eval_bitwise_xor(t, v);
774 }
775 }
776 template <class T, class U, class V>
eval_bitwise_xor(T & t,const U & u,const V & v)777 inline void eval_bitwise_xor(T& t, const U& u, const V& v)
778 {
779 eval_bitwise_xor_default(t, u, v);
780 }
781
782 template <class T>
eval_increment(T & val)783 inline void eval_increment(T& val)
784 {
785 typedef typename mpl::front<typename T::unsigned_types>::type ui_type;
786 eval_add(val, static_cast<ui_type>(1u));
787 }
788 template <class T>
eval_decrement(T & val)789 inline void eval_decrement(T& val)
790 {
791 typedef typename mpl::front<typename T::unsigned_types>::type ui_type;
792 eval_subtract(val, static_cast<ui_type>(1u));
793 }
794
795 template <class T, class V>
eval_left_shift(T & result,const T & arg,const V val)796 inline void eval_left_shift(T& result, const T& arg, const V val)
797 {
798 result = arg;
799 eval_left_shift(result, val);
800 }
801
802 template <class T, class V>
eval_right_shift(T & result,const T & arg,const V val)803 inline void eval_right_shift(T& result, const T& arg, const V val)
804 {
805 result = arg;
806 eval_right_shift(result, val);
807 }
808
809 template <class T>
eval_is_zero(const T & val)810 inline bool eval_is_zero(const T& val)
811 {
812 typedef typename mpl::front<typename T::unsigned_types>::type ui_type;
813 return val.compare(static_cast<ui_type>(0)) == 0;
814 }
815 template <class T>
eval_get_sign(const T & val)816 inline int eval_get_sign(const T& val)
817 {
818 typedef typename mpl::front<typename T::unsigned_types>::type ui_type;
819 return val.compare(static_cast<ui_type>(0));
820 }
821
822 template <class T, class V>
assign_components_imp(T & result,const V & v1,const V & v2,const mpl::int_<number_kind_rational> &)823 inline void assign_components_imp(T& result, const V& v1, const V& v2, const mpl::int_<number_kind_rational>&)
824 {
825 result = v1;
826 T t;
827 t = v2;
828 eval_divide(result, t);
829 }
830
831 template <class T, class V>
assign_components(T & result,const V & v1,const V & v2)832 inline void assign_components(T& result, const V& v1, const V& v2)
833 {
834 return assign_components_imp(result, v1, v2, typename number_category<T>::type());
835 }
836
837 template <class R, int b>
838 struct has_enough_bits
839 {
840 template <class T>
841 struct type : public mpl::and_<mpl::not_<is_same<R, T> >, mpl::bool_<std::numeric_limits<T>::digits >= b> >{};
842 };
843
844 template <class R>
845 struct terminal
846 {
terminalboost::multiprecision::default_ops::terminal847 terminal(const R& v) : value(v){}
terminalboost::multiprecision::default_ops::terminal848 terminal(){}
operator =boost::multiprecision::default_ops::terminal849 terminal& operator = (R val) { value = val; return *this; }
850 R value;
operator Rboost::multiprecision::default_ops::terminal851 operator R()const { return value; }
852 };
853
854 template<class R, class B>
855 struct calculate_next_larger_type
856 {
857 // Find which list we're looking through:
858 typedef typename mpl::if_<
859 is_signed<R>,
860 typename B::signed_types,
861 typename mpl::if_<
862 is_unsigned<R>,
863 typename B::unsigned_types,
864 typename B::float_types
865 >::type
866 >::type list_type;
867 // A predicate to find a type with enough bits:
868 typedef typename has_enough_bits<R, std::numeric_limits<R>::digits>::template type<mpl::_> pred_type;
869 // See if the last type is in the list, if so we have to start after this:
870 typedef typename mpl::find_if<
871 list_type,
872 is_same<R, mpl::_>
873 >::type start_last;
874 // Where we're starting from, either the start of the sequence or the last type found:
875 typedef typename mpl::if_<is_same<start_last, typename mpl::end<list_type>::type>, typename mpl::begin<list_type>::type, start_last>::type start_seq;
876 // The range we're searching:
877 typedef mpl::iterator_range<start_seq, typename mpl::end<list_type>::type> range;
878 // Find the next type:
879 typedef typename mpl::find_if<
880 range,
881 pred_type
882 >::type iter_type;
883 // Either the next type, or a "terminal" to indicate we've run out of types to search:
884 typedef typename mpl::eval_if<
885 is_same<typename mpl::end<list_type>::type, iter_type>,
886 mpl::identity<terminal<R> >,
887 mpl::deref<iter_type>
888 >::type type;
889 };
890
891 template <class R, class T>
check_in_range(const T & t)892 inline bool check_in_range(const T& t)
893 {
894 // Can t fit in an R?
895 if(std::numeric_limits<R>::is_specialized && std::numeric_limits<R>::is_bounded && (t > (std::numeric_limits<R>::max)()))
896 return true;
897 return false;
898 }
899
900 template <class R, class T>
check_in_range(const terminal<T> &)901 inline bool check_in_range(const terminal<T>&)
902 {
903 return false;
904 }
905
906 template <class R, class B>
eval_convert_to(R * result,const B & backend)907 inline void eval_convert_to(R* result, const B& backend)
908 {
909 typedef typename calculate_next_larger_type<R, B>::type next_type;
910 next_type n;
911 eval_convert_to(&n, backend);
912 if(check_in_range<R>(n))
913 {
914 *result = (std::numeric_limits<R>::max)();
915 }
916 else
917 *result = static_cast<R>(n);
918 }
919
920 template <class R, class B>
eval_convert_to(terminal<R> * result,const B & backend)921 inline void eval_convert_to(terminal<R>* result, const B& backend)
922 {
923 //
924 // We ran out of types to try for the conversion, try
925 // a lexical_cast and hope for the best:
926 //
927 result->value = boost::lexical_cast<R>(backend.str(0, std::ios_base::fmtflags(0)));
928 }
929
930 template <class B1, class B2, expression_template_option et>
eval_convert_to(terminal<number<B1,et>> * result,const B2 & backend)931 inline void eval_convert_to(terminal<number<B1, et> >* result, const B2& backend)
932 {
933 //
934 // We ran out of types to try for the conversion, try
935 // a generic conversion and hope for the best:
936 //
937 boost::multiprecision::detail::generic_interconvert(result->value.backend(), backend, number_category<B1>(), number_category<B2>());
938 }
939
940 template <class B>
eval_convert_to(std::string * result,const B & backend)941 inline void eval_convert_to(std::string* result, const B& backend)
942 {
943 *result = backend.str(0, std::ios_base::fmtflags(0));
944 }
945 //
946 // Functions:
947 //
948 template <class T>
eval_abs(T & result,const T & arg)949 void eval_abs(T& result, const T& arg)
950 {
951 typedef typename T::signed_types type_list;
952 typedef typename mpl::front<type_list>::type front;
953 result = arg;
954 if(arg.compare(front(0)) < 0)
955 result.negate();
956 }
957 template <class T>
eval_fabs(T & result,const T & arg)958 void eval_fabs(T& result, const T& arg)
959 {
960 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The fabs function is only valid for floating point types.");
961 typedef typename T::signed_types type_list;
962 typedef typename mpl::front<type_list>::type front;
963 result = arg;
964 if(arg.compare(front(0)) < 0)
965 result.negate();
966 }
967
968 template <class Backend>
eval_fpclassify(const Backend & arg)969 inline int eval_fpclassify(const Backend& arg)
970 {
971 BOOST_STATIC_ASSERT_MSG(number_category<Backend>::value == number_kind_floating_point, "The fpclassify function is only valid for floating point types.");
972 return eval_is_zero(arg) ? FP_ZERO : FP_NORMAL;
973 }
974
975 template <class T>
eval_fmod(T & result,const T & a,const T & b)976 inline void eval_fmod(T& result, const T& a, const T& b)
977 {
978 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The fmod function is only valid for floating point types.");
979 if((&result == &a) || (&result == &b))
980 {
981 T temp;
982 eval_fmod(temp, a, b);
983 result = temp;
984 return;
985 }
986 switch(eval_fpclassify(a))
987 {
988 case FP_ZERO:
989 result = a;
990 return;
991 case FP_INFINITE:
992 case FP_NAN:
993 result = std::numeric_limits<number<T> >::quiet_NaN().backend();
994 errno = EDOM;
995 return;
996 }
997 switch(eval_fpclassify(b))
998 {
999 case FP_ZERO:
1000 case FP_NAN:
1001 result = std::numeric_limits<number<T> >::quiet_NaN().backend();
1002 errno = EDOM;
1003 return;
1004 }
1005 T n;
1006 eval_divide(result, a, b);
1007 if(eval_get_sign(result) < 0)
1008 eval_ceil(n, result);
1009 else
1010 eval_floor(n, result);
1011 eval_multiply(n, b);
1012 eval_subtract(result, a, n);
1013 }
1014 template<class T, class A>
eval_fmod(T & result,const T & x,const A & a)1015 inline typename enable_if<is_arithmetic<A>, void>::type eval_fmod(T& result, const T& x, const A& a)
1016 {
1017 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
1018 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1019 cast_type c;
1020 c = a;
1021 eval_fmod(result, x, c);
1022 }
1023
1024 template<class T, class A>
eval_fmod(T & result,const A & x,const T & a)1025 inline typename enable_if<is_arithmetic<A>, void>::type eval_fmod(T& result, const A& x, const T& a)
1026 {
1027 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
1028 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1029 cast_type c;
1030 c = x;
1031 eval_fmod(result, c, a);
1032 }
1033
1034 template <class T>
1035 void eval_round(T& result, const T& a);
1036
1037 template <class T>
eval_remquo(T & result,const T & a,const T & b,int * pi)1038 inline void eval_remquo(T& result, const T& a, const T& b, int* pi)
1039 {
1040 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The remquo function is only valid for floating point types.");
1041 if((&result == &a) || (&result == &b))
1042 {
1043 T temp;
1044 eval_remquo(temp, a, b, pi);
1045 result = temp;
1046 return;
1047 }
1048 T n;
1049 eval_divide(result, a, b);
1050 eval_round(n, result);
1051 eval_convert_to(pi, n);
1052 eval_multiply(n, b);
1053 eval_subtract(result, a, n);
1054 }
1055 template<class T, class A>
eval_remquo(T & result,const T & x,const A & a,int * pi)1056 inline typename enable_if<is_arithmetic<A>, void>::type eval_remquo(T& result, const T& x, const A& a, int* pi)
1057 {
1058 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
1059 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1060 cast_type c;
1061 c = a;
1062 eval_remquo(result, x, c, pi);
1063 }
1064 template<class T, class A>
eval_remquo(T & result,const A & x,const T & a,int * pi)1065 inline typename enable_if<is_arithmetic<A>, void>::type eval_remquo(T& result, const A& x, const T& a, int* pi)
1066 {
1067 typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
1068 typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1069 cast_type c;
1070 c = x;
1071 eval_remquo(result, c, a, pi);
1072 }
1073 template <class T, class U, class V>
eval_remainder(T & result,const U & a,const V & b)1074 inline void eval_remainder(T& result, const U& a, const V& b)
1075 {
1076 int i;
1077 eval_remquo(result, a, b, &i);
1078 }
1079
1080 template <class B>
1081 bool eval_gt(const B& a, const B& b);
1082 template <class T, class U>
1083 bool eval_gt(const T& a, const U& b);
1084 template <class B>
1085 bool eval_lt(const B& a, const B& b);
1086 template <class T, class U>
1087 bool eval_lt(const T& a, const U& b);
1088
1089 template<class T>
eval_fdim(T & result,const T & a,const T & b)1090 inline void eval_fdim(T& result, const T& a, const T& b)
1091 {
1092 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1093 static const ui_type zero = 0u;
1094 switch(eval_fpclassify(b))
1095 {
1096 case FP_NAN:
1097 case FP_INFINITE:
1098 result = zero;
1099 return;
1100 }
1101 switch(eval_fpclassify(a))
1102 {
1103 case FP_NAN:
1104 result = zero;
1105 return;
1106 case FP_INFINITE:
1107 result = a;
1108 return;
1109 }
1110 if(eval_gt(a, b))
1111 {
1112 eval_subtract(result, a, b);
1113 }
1114 else
1115 result = zero;
1116 }
1117
1118 template<class T, class A>
eval_fdim(T & result,const T & a,const A & b)1119 inline typename boost::enable_if_c<boost::is_arithmetic<A>::value>::type eval_fdim(T& result, const T& a, const A& b)
1120 {
1121 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1122 typedef typename boost::multiprecision::detail::canonical<A, T>::type arithmetic_type;
1123 static const ui_type zero = 0u;
1124 arithmetic_type canonical_b = b;
1125 switch((::boost::math::fpclassify)(b))
1126 {
1127 case FP_NAN:
1128 case FP_INFINITE:
1129 result = zero;
1130 return;
1131 }
1132 switch(eval_fpclassify(a))
1133 {
1134 case FP_NAN:
1135 result = zero;
1136 return;
1137 case FP_INFINITE:
1138 result = a;
1139 return;
1140 }
1141 if(eval_gt(a, canonical_b))
1142 {
1143 eval_subtract(result, a, canonical_b);
1144 }
1145 else
1146 result = zero;
1147 }
1148
1149 template<class T, class A>
eval_fdim(T & result,const A & a,const T & b)1150 inline typename boost::enable_if_c<boost::is_arithmetic<A>::value>::type eval_fdim(T& result, const A& a, const T& b)
1151 {
1152 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1153 typedef typename boost::multiprecision::detail::canonical<A, T>::type arithmetic_type;
1154 static const ui_type zero = 0u;
1155 arithmetic_type canonical_a = a;
1156 switch(eval_fpclassify(b))
1157 {
1158 case FP_NAN:
1159 case FP_INFINITE:
1160 result = zero;
1161 return;
1162 }
1163 switch((::boost::math::fpclassify)(a))
1164 {
1165 case FP_NAN:
1166 result = zero;
1167 return;
1168 case FP_INFINITE:
1169 result = std::numeric_limits<number<T> >::infinity().backend();
1170 return;
1171 }
1172 if(eval_gt(canonical_a, b))
1173 {
1174 eval_subtract(result, canonical_a, b);
1175 }
1176 else
1177 result = zero;
1178 }
1179
1180 template <class T>
eval_trunc(T & result,const T & a)1181 inline void eval_trunc(T& result, const T& a)
1182 {
1183 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The trunc function is only valid for floating point types.");
1184 switch(eval_fpclassify(a))
1185 {
1186 case FP_NAN:
1187 errno = EDOM;
1188 // fallthrough...
1189 case FP_ZERO:
1190 case FP_INFINITE:
1191 result = a;
1192 return;
1193 }
1194 if(eval_get_sign(a) < 0)
1195 eval_ceil(result, a);
1196 else
1197 eval_floor(result, a);
1198 }
1199
1200 template <class T>
eval_modf(T & result,T const & arg,T * pipart)1201 inline void eval_modf(T& result, T const& arg, T* pipart)
1202 {
1203 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1204 int c = eval_fpclassify(arg);
1205 if(c == (int)FP_NAN)
1206 {
1207 if(pipart)
1208 *pipart = arg;
1209 result = arg;
1210 return;
1211 }
1212 else if(c == (int)FP_INFINITE)
1213 {
1214 if(pipart)
1215 *pipart = arg;
1216 result = ui_type(0u);
1217 return;
1218 }
1219 if(pipart)
1220 {
1221 eval_trunc(*pipart, arg);
1222 eval_subtract(result, arg, *pipart);
1223 }
1224 else
1225 {
1226 T ipart;
1227 eval_trunc(ipart, arg);
1228 eval_subtract(result, arg, ipart);
1229 }
1230 }
1231
1232 template <class T>
eval_round(T & result,const T & a)1233 inline void eval_round(T& result, const T& a)
1234 {
1235 BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The round function is only valid for floating point types.");
1236 typedef typename boost::multiprecision::detail::canonical<float, T>::type fp_type;
1237 int c = eval_fpclassify(a);
1238 if(c == (int)FP_NAN)
1239 {
1240 result = a;
1241 errno = EDOM;
1242 return;
1243 }
1244 if((c == FP_ZERO) || (c == (int)FP_INFINITE))
1245 {
1246 result = a;
1247 }
1248 else if(eval_get_sign(a) < 0)
1249 {
1250 eval_subtract(result, a, fp_type(0.5f));
1251 eval_ceil(result, result);
1252 }
1253 else
1254 {
1255 eval_add(result, a, fp_type(0.5f));
1256 eval_floor(result, result);
1257 }
1258 }
1259
1260 template <class B>
1261 void eval_lcm(B& result, const B& a, const B& b);
1262 template <class B>
1263 void eval_gcd(B& result, const B& a, const B& b);
1264
1265 template <class T, class Arithmetic>
eval_gcd(T & result,const T & a,const Arithmetic & b)1266 inline typename enable_if<is_integral<Arithmetic> >::type eval_gcd(T& result, const T& a, const Arithmetic& b)
1267 {
1268 typedef typename boost::multiprecision::detail::canonical<Arithmetic, T>::type si_type;
1269 using default_ops::eval_gcd;
1270 T t;
1271 t = static_cast<si_type>(b);
1272 eval_gcd(result, a, t);
1273 }
1274 template <class T, class Arithmetic>
eval_gcd(T & result,const Arithmetic & a,const T & b)1275 inline typename enable_if<is_integral<Arithmetic> >::type eval_gcd(T& result, const Arithmetic& a, const T& b)
1276 {
1277 eval_gcd(result, b, a);
1278 }
1279 template <class T, class Arithmetic>
eval_lcm(T & result,const T & a,const Arithmetic & b)1280 inline typename enable_if<is_integral<Arithmetic> >::type eval_lcm(T& result, const T& a, const Arithmetic& b)
1281 {
1282 typedef typename boost::multiprecision::detail::canonical<Arithmetic, T>::type si_type;
1283 using default_ops::eval_lcm;
1284 T t;
1285 t = static_cast<si_type>(b);
1286 eval_lcm(result, a, t);
1287 }
1288 template <class T, class Arithmetic>
eval_lcm(T & result,const Arithmetic & a,const T & b)1289 inline typename enable_if<is_integral<Arithmetic> >::type eval_lcm(T& result, const Arithmetic& a, const T& b)
1290 {
1291 eval_lcm(result, b, a);
1292 }
1293
1294 template <class T>
eval_lsb(const T & val)1295 inline unsigned eval_lsb(const T& val)
1296 {
1297 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1298 int c = eval_get_sign(val);
1299 if(c == 0)
1300 {
1301 BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
1302 }
1303 if(c < 0)
1304 {
1305 BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
1306 }
1307 unsigned result = 0;
1308 T mask, t;
1309 mask = ui_type(1);
1310 do
1311 {
1312 eval_bitwise_and(t, mask, val);
1313 ++result;
1314 eval_left_shift(mask, 1);
1315 }
1316 while(eval_is_zero(t));
1317
1318 return --result;
1319 }
1320
1321 template <class T>
eval_msb(const T & val)1322 inline int eval_msb(const T& val)
1323 {
1324 int c = eval_get_sign(val);
1325 if(c == 0)
1326 {
1327 BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
1328 }
1329 if(c < 0)
1330 {
1331 BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
1332 }
1333 //
1334 // This implementation is really really rubbish - it does
1335 // a linear scan for the most-significant-bit. We should really
1336 // do a binary search, but as none of our backends actually needs
1337 // this implementation, we'll leave it for now. In fact for most
1338 // backends it's likely that there will always be a more efficient
1339 // native implementation possible.
1340 //
1341 unsigned result = 0;
1342 T t(val);
1343 while(!eval_is_zero(t))
1344 {
1345 eval_right_shift(t, 1);
1346 ++result;
1347 }
1348 return --result;
1349 }
1350
1351 template <class T>
eval_bit_test(const T & val,unsigned index)1352 inline bool eval_bit_test(const T& val, unsigned index)
1353 {
1354 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1355 T mask, t;
1356 mask = ui_type(1);
1357 eval_left_shift(mask, index);
1358 eval_bitwise_and(t, mask, val);
1359 return !eval_is_zero(t);
1360 }
1361
1362 template <class T>
eval_bit_set(T & val,unsigned index)1363 inline void eval_bit_set(T& val, unsigned index)
1364 {
1365 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1366 T mask;
1367 mask = ui_type(1);
1368 eval_left_shift(mask, index);
1369 eval_bitwise_or(val, mask);
1370 }
1371
1372 template <class T>
eval_bit_flip(T & val,unsigned index)1373 inline void eval_bit_flip(T& val, unsigned index)
1374 {
1375 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1376 T mask;
1377 mask = ui_type(1);
1378 eval_left_shift(mask, index);
1379 eval_bitwise_xor(val, mask);
1380 }
1381
1382 template <class T>
eval_bit_unset(T & val,unsigned index)1383 inline void eval_bit_unset(T& val, unsigned index)
1384 {
1385 typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1386 T mask, t;
1387 mask = ui_type(1);
1388 eval_left_shift(mask, index);
1389 eval_bitwise_and(t, mask, val);
1390 if(!eval_is_zero(t))
1391 eval_bitwise_xor(val, mask);
1392 }
1393
1394 template <class B>
eval_integer_sqrt(B & s,B & r,const B & x)1395 void eval_integer_sqrt(B& s, B& r, const B& x)
1396 {
1397 //
1398 // This is slow bit-by-bit integer square root, see for example
1399 // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
1400 // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
1401 // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
1402 // at some point.
1403 //
1404 typedef typename boost::multiprecision::detail::canonical<unsigned char, B>::type ui_type;
1405
1406 s = ui_type(0u);
1407 if(eval_get_sign(x) == 0)
1408 {
1409 r = ui_type(0u);
1410 return;
1411 }
1412 int g = eval_msb(x);
1413 if(g <= 1)
1414 {
1415 s = ui_type(1);
1416 eval_subtract(r, x, s);
1417 return;
1418 }
1419
1420 B t;
1421 r = x;
1422 g /= 2;
1423 int org_g = g;
1424 eval_bit_set(s, g);
1425 eval_bit_set(t, 2 * g);
1426 eval_subtract(r, x, t);
1427 --g;
1428 if(eval_get_sign(r) == 0)
1429 return;
1430 int msbr = eval_msb(r);
1431 do
1432 {
1433 if(msbr >= org_g + g + 1)
1434 {
1435 t = s;
1436 eval_left_shift(t, g + 1);
1437 eval_bit_set(t, 2 * g);
1438 if(t.compare(r) <= 0)
1439 {
1440 BOOST_ASSERT(g >= 0);
1441 eval_bit_set(s, g);
1442 eval_subtract(r, t);
1443 if(eval_get_sign(r) == 0)
1444 return;
1445 msbr = eval_msb(r);
1446 }
1447 }
1448 --g;
1449 }
1450 while(g >= 0);
1451 }
1452
1453 //
1454 // These have to implemented by the backend, declared here so that our macro generated code compiles OK.
1455 //
1456 template <class T>
1457 typename enable_if_c<sizeof(T) == 0>::type eval_floor();
1458 template <class T>
1459 typename enable_if_c<sizeof(T) == 0>::type eval_ceil();
1460 template <class T>
1461 typename enable_if_c<sizeof(T) == 0>::type eval_trunc();
1462 template <class T>
1463 typename enable_if_c<sizeof(T) == 0>::type eval_sqrt();
1464 template <class T>
1465 typename enable_if_c<sizeof(T) == 0>::type eval_ldexp();
1466 template <class T>
1467 typename enable_if_c<sizeof(T) == 0>::type eval_frexp();
1468
1469 //
1470 // eval_logb and eval_scalbn simply assume base 2 and forward to
1471 // eval_ldexp and eval_frexp:
1472 //
1473 template <class B>
eval_ilogb(const B & val)1474 inline typename B::exponent_type eval_ilogb(const B& val)
1475 {
1476 BOOST_STATIC_ASSERT_MSG(!std::numeric_limits<number<B> >::is_specialized || (std::numeric_limits<number<B> >::radix == 2), "The default implementation of ilogb requires a base 2 number type");
1477 typename B::exponent_type e;
1478 switch(eval_fpclassify(val))
1479 {
1480 case FP_NAN:
1481 #ifdef FP_ILOGBNAN
1482 return FP_ILOGBNAN > 0 ? (std::numeric_limits<typename B::exponent_type>::max)() : (std::numeric_limits<typename B::exponent_type>::min)();
1483 #else
1484 return (std::numeric_limits<typename B::exponent_type>::max)();
1485 #endif
1486 case FP_INFINITE:
1487 return (std::numeric_limits<typename B::exponent_type>::max)();
1488 case FP_ZERO:
1489 return (std::numeric_limits<typename B::exponent_type>::min)();
1490 }
1491 B result;
1492 eval_frexp(result, val, &e);
1493 return e - 1;
1494 }
1495
1496 template <class T>
1497 int eval_signbit(const T& val);
1498
1499 template <class B>
eval_logb(B & result,const B & val)1500 inline void eval_logb(B& result, const B& val)
1501 {
1502 switch(eval_fpclassify(val))
1503 {
1504 case FP_NAN:
1505 result = val;
1506 errno = EDOM;
1507 return;
1508 case FP_ZERO:
1509 result = std::numeric_limits<number<B> >::infinity().backend();
1510 result.negate();
1511 errno = ERANGE;
1512 return;
1513 case FP_INFINITE:
1514 result = val;
1515 if(eval_signbit(val))
1516 result.negate();
1517 return;
1518 }
1519 typedef typename boost::mpl::if_c<boost::is_same<boost::intmax_t, long>::value, boost::long_long_type, boost::intmax_t>::type max_t;
1520 result = static_cast<max_t>(eval_ilogb(val));
1521 }
1522 template <class B, class A>
eval_scalbn(B & result,const B & val,A e)1523 inline void eval_scalbn(B& result, const B& val, A e)
1524 {
1525 BOOST_STATIC_ASSERT_MSG(!std::numeric_limits<number<B> >::is_specialized || (std::numeric_limits<number<B> >::radix == 2), "The default implementation of scalbn requires a base 2 number type");
1526 eval_ldexp(result, val, static_cast<typename B::exponent_type>(e));
1527 }
1528 template <class B, class A>
eval_scalbln(B & result,const B & val,A e)1529 inline void eval_scalbln(B& result, const B& val, A e)
1530 {
1531 eval_scalbn(result, val, e);
1532 }
1533
1534 template <class T>
is_arg_nan(const T & val,mpl::true_ const &,const mpl::false_ &)1535 inline bool is_arg_nan(const T& val, mpl::true_ const&, const mpl::false_&)
1536 {
1537 return eval_fpclassify(val) == FP_NAN;
1538 }
1539 template <class T>
is_arg_nan(const T & val,mpl::false_ const &,const mpl::true_ &)1540 inline bool is_arg_nan(const T& val, mpl::false_ const&, const mpl::true_&)
1541 {
1542 return (boost::math::isnan)(val);
1543 }
1544 template <class T>
is_arg_nan(const T &,mpl::false_ const &,const mpl::false_ &)1545 inline bool is_arg_nan(const T&, mpl::false_ const&, const mpl::false_&)
1546 {
1547 return false;
1548 }
1549
1550 template <class T>
is_arg_nan(const T & val)1551 inline bool is_arg_nan(const T& val)
1552 {
1553 return is_arg_nan(val, mpl::bool_<boost::multiprecision::detail::is_backend<T>::value>(), is_floating_point<T>());
1554 }
1555
1556 template <class T, class U, class V>
eval_fmax(T & result,const U & a,const V & b)1557 inline void eval_fmax(T& result, const U& a, const V& b)
1558 {
1559 if(is_arg_nan(a))
1560 result = number<T>::canonical_value(b);
1561 else if(is_arg_nan(b))
1562 result = number<T>::canonical_value(a);
1563 else if(eval_lt(number<T>::canonical_value(a), number<T>::canonical_value(b)))
1564 result = number<T>::canonical_value(b);
1565 else
1566 result = number<T>::canonical_value(a);
1567 }
1568 template <class T, class U, class V>
eval_fmin(T & result,const U & a,const V & b)1569 inline void eval_fmin(T& result, const U& a, const V& b)
1570 {
1571 if(is_arg_nan(a))
1572 result = number<T>::canonical_value(b);
1573 else if(is_arg_nan(b))
1574 result = number<T>::canonical_value(a);
1575 else if(eval_lt(number<T>::canonical_value(a), number<T>::canonical_value(b)))
1576 result = number<T>::canonical_value(a);
1577 else
1578 result = number<T>::canonical_value(b);
1579 }
1580
1581 template <class R, class T, class U>
eval_hypot(R & result,const T & a,const U & b)1582 inline void eval_hypot(R& result, const T& a, const U& b)
1583 {
1584 //
1585 // Normalize x and y, so that both are positive and x >= y:
1586 //
1587 R x, y;
1588 x = number<R>::canonical_value(a);
1589 y = number<R>::canonical_value(b);
1590 if(eval_get_sign(x) < 0)
1591 x.negate();
1592 if(eval_get_sign(y) < 0)
1593 y.negate();
1594
1595 // Special case, see C99 Annex F.
1596 // The order of the if's is important: do not change!
1597 int c1 = eval_fpclassify(x);
1598 int c2 = eval_fpclassify(y);
1599
1600 if(c1 == FP_ZERO)
1601 {
1602 result = y;
1603 return;
1604 }
1605 if(c2 == FP_ZERO)
1606 {
1607 result = x;
1608 return;
1609 }
1610 if(c1 == FP_INFINITE)
1611 {
1612 result = x;
1613 return;
1614 }
1615 if((c2 == FP_INFINITE) || (c2 == FP_NAN))
1616 {
1617 result = y;
1618 return;
1619 }
1620 if(c1 == FP_NAN)
1621 {
1622 result = x;
1623 return;
1624 }
1625
1626 if(eval_gt(y, x))
1627 x.swap(y);
1628
1629 eval_multiply(result, x, std::numeric_limits<number<R> >::epsilon().backend());
1630
1631 if(eval_gt(result, y))
1632 {
1633 result = x;
1634 return;
1635 }
1636
1637 R rat;
1638 eval_divide(rat, y, x);
1639 eval_multiply(result, rat, rat);
1640 eval_increment(result);
1641 eval_sqrt(rat, result);
1642 eval_multiply(result, rat, x);
1643 }
1644
1645 template <class R, class T>
eval_nearbyint(R & result,const T & a)1646 inline void eval_nearbyint(R& result, const T& a)
1647 {
1648 eval_round(result, a);
1649 }
1650 template <class R, class T>
eval_rint(R & result,const T & a)1651 inline void eval_rint(R& result, const T& a)
1652 {
1653 eval_nearbyint(result, a);
1654 }
1655
1656 template <class T>
eval_signbit(const T & val)1657 inline int eval_signbit(const T& val)
1658 {
1659 return eval_get_sign(val) < 0 ? 1 : 0;
1660 }
1661
1662 //
1663 // These functions are implemented in separate files, but expanded inline here,
1664 // DO NOT CHANGE THE ORDER OF THESE INCLUDES:
1665 //
1666 #include <boost/multiprecision/detail/functions/constants.hpp>
1667 #include <boost/multiprecision/detail/functions/pow.hpp>
1668 #include <boost/multiprecision/detail/functions/trig.hpp>
1669
1670 }
1671
1672 //
1673 // Default versions of floating point classification routines:
1674 //
1675 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1676 inline int fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1677 {
1678 using multiprecision::default_ops::eval_fpclassify;
1679 return eval_fpclassify(arg.backend());
1680 }
1681 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1682 inline int fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1683 {
1684 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1685 return fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1686 }
1687 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1688 inline bool isfinite BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1689 {
1690 int v = fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(arg);
1691 return (v != (int)FP_INFINITE) && (v != (int)FP_NAN);
1692 }
1693 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1694 inline bool isfinite BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1695 {
1696 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1697 return isfinite BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1698 }
1699 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1700 inline bool isnan BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1701 {
1702 return fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(arg) == (int)FP_NAN;
1703 }
1704 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1705 inline bool isnan BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1706 {
1707 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1708 return isnan BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1709 }
1710 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1711 inline bool isinf BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1712 {
1713 return fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(arg) == (int)FP_INFINITE;
1714 }
1715 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1716 inline bool isinf BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1717 {
1718 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1719 return isinf BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1720 }
1721 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1722 inline bool isnormal BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1723 {
1724 return fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(arg) == (int)FP_NORMAL;
1725 }
1726 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1727 inline bool isnormal BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1728 {
1729 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1730 return isnormal BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1731 }
1732
1733 // Default versions of sign manipulation functions, if individual backends can do better than this
1734 // (for example with signed zero), then they should overload these functions further:
1735
1736 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1737 inline int sign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1738 {
1739 return arg.sign();
1740 }
1741 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1742 inline int sign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1743 {
1744 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1745 return sign BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1746 }
1747
1748 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1749 inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1750 {
1751 using default_ops::eval_signbit;
1752 return eval_signbit(arg.backend());
1753 }
1754 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1755 inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1756 {
1757 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1758 return signbit BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1759 }
1760 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1761 inline multiprecision::number<Backend, ExpressionTemplates> changesign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1762 {
1763 return -arg;
1764 }
1765 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1766 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type changesign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1767 {
1768 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1769 return changesign BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1770 }
1771 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & a,const multiprecision::number<Backend,ExpressionTemplates> & b)1772 inline multiprecision::number<Backend, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& a, const multiprecision::number<Backend, ExpressionTemplates>& b)
1773 {
1774 return (boost::multiprecision::signbit)(a) != (boost::multiprecision::signbit)(b) ? (boost::multiprecision::changesign)(a) : a;
1775 }
1776 template <class Backend, multiprecision::expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & a,const multiprecision::detail::expression<tag,A1,A2,A3,A4> & b)1777 inline multiprecision::number<Backend, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& a, const multiprecision::detail::expression<tag, A1, A2, A3, A4>& b)
1778 {
1779 return copysign BOOST_PREVENT_MACRO_SUBSTITUTION(a, multiprecision::number<Backend, ExpressionTemplates>(b));
1780 }
1781 template <class tag, class A1, class A2, class A3, class A4, class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & a,const multiprecision::number<Backend,ExpressionTemplates> & b)1782 inline multiprecision::number<Backend, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& a, const multiprecision::number<Backend, ExpressionTemplates>& b)
1783 {
1784 return copysign BOOST_PREVENT_MACRO_SUBSTITUTION(multiprecision::number<Backend, ExpressionTemplates>(a), b);
1785 }
1786 template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & a,const multiprecision::detail::expression<tagb,A1b,A2b,A3b,A4b> & b)1787 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& a, const multiprecision::detail::expression<tagb, A1b, A2b, A3b, A4b>& b)
1788 {
1789 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1790 return copysign BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(a), value_type(b));
1791 }
1792
1793 } // namespace multiprecision
1794
1795 namespace math {
1796
1797 //
1798 // Import Math functions here, so they can be found by Boost.Math:
1799 //
1800 using boost::multiprecision::signbit;
1801 using boost::multiprecision::sign;
1802 using boost::multiprecision::copysign;
1803 using boost::multiprecision::changesign;
1804 using boost::multiprecision::fpclassify;
1805 using boost::multiprecision::isinf;
1806 using boost::multiprecision::isnan;
1807 using boost::multiprecision::isnormal;
1808 using boost::multiprecision::isfinite;
1809
1810 }
1811
1812 namespace multiprecision{
1813
1814 typedef ::boost::math::policies::policy<
1815 ::boost::math::policies::domain_error< ::boost::math::policies::errno_on_error>,
1816 ::boost::math::policies::pole_error< ::boost::math::policies::errno_on_error>,
1817 ::boost::math::policies::overflow_error< ::boost::math::policies::errno_on_error>,
1818 ::boost::math::policies::evaluation_error< ::boost::math::policies::errno_on_error>,
1819 ::boost::math::policies::rounding_error< ::boost::math::policies::errno_on_error>
1820 > c99_error_policy;
1821
1822 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1823 inline multiprecision::number<Backend, ExpressionTemplates> asinh BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1824 {
1825 return boost::math::asinh(arg, c99_error_policy());
1826 }
1827 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1828 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type asinh BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1829 {
1830 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1831 return asinh(value_type(arg));
1832 }
1833 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1834 inline multiprecision::number<Backend, ExpressionTemplates> acosh BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1835 {
1836 return boost::math::acosh(arg, c99_error_policy());
1837 }
1838 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1839 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type acosh BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1840 {
1841 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1842 return acosh(value_type(arg));
1843 }
1844 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1845 inline multiprecision::number<Backend, ExpressionTemplates> atanh BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1846 {
1847 return boost::math::atanh(arg, c99_error_policy());
1848 }
1849 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1850 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type atanh BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1851 {
1852 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1853 return atanh(value_type(arg));
1854 }
1855 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1856 inline multiprecision::number<Backend, ExpressionTemplates> cbrt BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1857 {
1858 return boost::math::cbrt(arg, c99_error_policy());
1859 }
1860 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1861 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type cbrt BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1862 {
1863 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1864 return cbrt(value_type(arg));
1865 }
1866 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1867 inline multiprecision::number<Backend, ExpressionTemplates> erf BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1868 {
1869 return boost::math::erf(arg, c99_error_policy());
1870 }
1871 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1872 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type erf BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1873 {
1874 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1875 return erf(value_type(arg));
1876 }
1877 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1878 inline multiprecision::number<Backend, ExpressionTemplates> erfc BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1879 {
1880 return boost::math::erfc(arg, c99_error_policy());
1881 }
1882 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1883 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type erfc BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1884 {
1885 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1886 return erfc(value_type(arg));
1887 }
1888 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1889 inline multiprecision::number<Backend, ExpressionTemplates> expm1 BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1890 {
1891 return boost::math::expm1(arg, c99_error_policy());
1892 }
1893 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1894 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type expm1 BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1895 {
1896 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1897 return expm1(value_type(arg));
1898 }
1899 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1900 inline multiprecision::number<Backend, ExpressionTemplates> lgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1901 {
1902 multiprecision::number<Backend, ExpressionTemplates> result;
1903 result = boost::math::lgamma(arg, c99_error_policy());
1904 if((boost::multiprecision::isnan)(result) && !(boost::multiprecision::isnan)(arg))
1905 {
1906 result = std::numeric_limits<multiprecision::number<Backend, ExpressionTemplates> >::infinity();
1907 errno = ERANGE;
1908 }
1909 return result;
1910 }
1911 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1912 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type lgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1913 {
1914 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1915 return lgamma(value_type(arg));
1916 }
1917 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1918 inline multiprecision::number<Backend, ExpressionTemplates> tgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1919 {
1920 if((arg == 0) && std::numeric_limits<multiprecision::number<Backend, ExpressionTemplates> >::has_infinity)
1921 {
1922 errno = ERANGE;
1923 return 1 / arg;
1924 }
1925 return boost::math::tgamma(arg, c99_error_policy());
1926 }
1927 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1928 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type tgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1929 {
1930 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1931 return tgamma(value_type(arg));
1932 }
1933
1934 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1935 inline long lrint BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1936 {
1937 return lround(arg);
1938 }
1939 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1940 inline long lrint BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1941 {
1942 return lround(arg);
1943 }
1944 #ifndef BOOST_NO_LONG_LONG
1945 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1946 inline boost::long_long_type llrint BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1947 {
1948 return llround(arg);
1949 }
1950 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1951 inline boost::long_long_type llrint BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1952 {
1953 return llround(arg);
1954 }
1955 #endif
1956 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1957 inline multiprecision::number<Backend, ExpressionTemplates> log1p BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1958 {
1959 return boost::math::log1p(arg, c99_error_policy());
1960 }
1961 template <class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & arg)1962 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type log1p BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1963 {
1964 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1965 return log1p(value_type(arg));
1966 }
1967
1968 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & a,const multiprecision::number<Backend,ExpressionTemplates> & b)1969 inline multiprecision::number<Backend, ExpressionTemplates> nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& a, const multiprecision::number<Backend, ExpressionTemplates>& b)
1970 {
1971 return boost::math::nextafter(a, b, c99_error_policy());
1972 }
1973 template <class Backend, multiprecision::expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & a,const multiprecision::detail::expression<tag,A1,A2,A3,A4> & b)1974 inline multiprecision::number<Backend, ExpressionTemplates> nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& a, const multiprecision::detail::expression<tag, A1, A2, A3, A4>& b)
1975 {
1976 return nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(a, multiprecision::number<Backend, ExpressionTemplates>(b));
1977 }
1978 template <class tag, class A1, class A2, class A3, class A4, class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & a,const multiprecision::number<Backend,ExpressionTemplates> & b)1979 inline multiprecision::number<Backend, ExpressionTemplates> nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& a, const multiprecision::number<Backend, ExpressionTemplates>& b)
1980 {
1981 return nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(multiprecision::number<Backend, ExpressionTemplates>(a), b);
1982 }
1983 template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & a,const multiprecision::detail::expression<tagb,A1b,A2b,A3b,A4b> & b)1984 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& a, const multiprecision::detail::expression<tagb, A1b, A2b, A3b, A4b>& b)
1985 {
1986 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1987 return nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(a), value_type(b));
1988 }
1989 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & a,const multiprecision::number<Backend,ExpressionTemplates> & b)1990 inline multiprecision::number<Backend, ExpressionTemplates> nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& a, const multiprecision::number<Backend, ExpressionTemplates>& b)
1991 {
1992 return boost::math::nextafter(a, b, c99_error_policy());
1993 }
1994 template <class Backend, multiprecision::expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & a,const multiprecision::detail::expression<tag,A1,A2,A3,A4> & b)1995 inline multiprecision::number<Backend, ExpressionTemplates> nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& a, const multiprecision::detail::expression<tag, A1, A2, A3, A4>& b)
1996 {
1997 return nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(a, multiprecision::number<Backend, ExpressionTemplates>(b));
1998 }
1999 template <class tag, class A1, class A2, class A3, class A4, class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & a,const multiprecision::number<Backend,ExpressionTemplates> & b)2000 inline multiprecision::number<Backend, ExpressionTemplates> nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& a, const multiprecision::number<Backend, ExpressionTemplates>& b)
2001 {
2002 return nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(multiprecision::number<Backend, ExpressionTemplates>(a), b);
2003 }
2004 template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag,A1,A2,A3,A4> & a,const multiprecision::detail::expression<tagb,A1b,A2b,A3b,A4b> & b)2005 inline typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& a, const multiprecision::detail::expression<tagb, A1b, A2b, A3b, A4b>& b)
2006 {
2007 typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
2008 return nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(a), value_type(b));
2009 }
2010
2011 template <class B1, class B2, class B3, expression_template_option ET1, expression_template_option ET2, expression_template_option ET3>
add(number<B1,ET1> & result,const number<B2,ET2> & a,const number<B3,ET3> & b)2012 inline number<B1, ET1>& add(number<B1, ET1>& result, const number<B2, ET2>& a, const number<B3, ET3>& b)
2013 {
2014 BOOST_STATIC_ASSERT_MSG((is_convertible<B2, B1>::value), "No conversion to the target of a mixed precision addition exists");
2015 BOOST_STATIC_ASSERT_MSG((is_convertible<B3, B1>::value), "No conversion to the target of a mixed precision addition exists");
2016 using default_ops::eval_add;
2017 eval_add(result.backend(), a.backend(), b.backend());
2018 return result;
2019 }
2020
2021 template <class B1, class B2, class B3, expression_template_option ET1, expression_template_option ET2, expression_template_option ET3>
subtract(number<B1,ET1> & result,const number<B2,ET2> & a,const number<B3,ET3> & b)2022 inline number<B1, ET1>& subtract(number<B1, ET1>& result, const number<B2, ET2>& a, const number<B3, ET3>& b)
2023 {
2024 BOOST_STATIC_ASSERT_MSG((is_convertible<B2, B1>::value), "No conversion to the target of a mixed precision addition exists");
2025 BOOST_STATIC_ASSERT_MSG((is_convertible<B3, B1>::value), "No conversion to the target of a mixed precision addition exists");
2026 using default_ops::eval_subtract;
2027 eval_subtract(result.backend(), a.backend(), b.backend());
2028 return result;
2029 }
2030
2031 template <class B1, class B2, class B3, expression_template_option ET1, expression_template_option ET2, expression_template_option ET3>
multiply(number<B1,ET1> & result,const number<B2,ET2> & a,const number<B3,ET3> & b)2032 inline number<B1, ET1>& multiply(number<B1, ET1>& result, const number<B2, ET2>& a, const number<B3, ET3>& b)
2033 {
2034 BOOST_STATIC_ASSERT_MSG((is_convertible<B2, B1>::value), "No conversion to the target of a mixed precision addition exists");
2035 BOOST_STATIC_ASSERT_MSG((is_convertible<B3, B1>::value), "No conversion to the target of a mixed precision addition exists");
2036 using default_ops::eval_multiply;
2037 eval_multiply(result.backend(), a.backend(), b.backend());
2038 return result;
2039 }
2040
2041 template <class B, expression_template_option ET, class I>
2042 inline typename enable_if_c<is_integral<I>::value, number<B, ET>&>::type
add(number<B,ET> & result,const I & a,const I & b)2043 add(number<B, ET>& result, const I& a, const I& b)
2044 {
2045 using default_ops::eval_add;
2046 typedef typename detail::canonical<I, B>::type canonical_type;
2047 eval_add(result.backend(), static_cast<canonical_type>(a), static_cast<canonical_type>(b));
2048 return result;
2049 }
2050
2051 template <class B, expression_template_option ET, class I>
2052 inline typename enable_if_c<is_integral<I>::value, number<B, ET>&>::type
subtract(number<B,ET> & result,const I & a,const I & b)2053 subtract(number<B, ET>& result, const I& a, const I& b)
2054 {
2055 using default_ops::eval_subtract;
2056 typedef typename detail::canonical<I, B>::type canonical_type;
2057 eval_subtract(result.backend(), static_cast<canonical_type>(a), static_cast<canonical_type>(b));
2058 return result;
2059 }
2060
2061 template <class B, expression_template_option ET, class I>
2062 inline typename enable_if_c<is_integral<I>::value, number<B, ET>&>::type
multiply(number<B,ET> & result,const I & a,const I & b)2063 multiply(number<B, ET>& result, const I& a, const I& b)
2064 {
2065 using default_ops::eval_multiply;
2066 typedef typename detail::canonical<I, B>::type canonical_type;
2067 eval_multiply(result.backend(), static_cast<canonical_type>(a), static_cast<canonical_type>(b));
2068 return result;
2069 }
2070
2071 template <class tag, class A1, class A2, class A3, class A4, class Policy>
trunc(const detail::expression<tag,A1,A2,A3,A4> & v,const Policy & pol)2072 inline typename detail::expression<tag, A1, A2, A3, A4>::result_type trunc(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2073 {
2074 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2075 return BOOST_MP_MOVE(trunc(number_type(v), pol));
2076 }
2077
2078 template <class Backend, expression_template_option ExpressionTemplates, class Policy>
trunc(const number<Backend,ExpressionTemplates> & v,const Policy &)2079 inline number<Backend, ExpressionTemplates> trunc(const number<Backend, ExpressionTemplates>& v, const Policy&)
2080 {
2081 using default_ops::eval_trunc;
2082 number<Backend, ExpressionTemplates> result;
2083 eval_trunc(result.backend(), v.backend());
2084 return BOOST_MP_MOVE(result);
2085 }
2086
2087 template <class tag, class A1, class A2, class A3, class A4, class Policy>
itrunc(const detail::expression<tag,A1,A2,A3,A4> & v,const Policy & pol)2088 inline int itrunc(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2089 {
2090 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2091 number_type r = trunc(v, pol);
2092 if((r > (std::numeric_limits<int>::max)()) || r < (std::numeric_limits<int>::min)() || !(boost::math::isfinite)(v))
2093 return boost::math::policies::raise_rounding_error("boost::multiprecision::itrunc<%1%>(%1%)", 0, number_type(v), 0, pol);
2094 return r.template convert_to<int>();
2095 }
2096 template <class tag, class A1, class A2, class A3, class A4>
itrunc(const detail::expression<tag,A1,A2,A3,A4> & v)2097 inline int itrunc(const detail::expression<tag, A1, A2, A3, A4>& v)
2098 {
2099 return itrunc(v, boost::math::policies::policy<>());
2100 }
2101 template <class Backend, expression_template_option ExpressionTemplates, class Policy>
itrunc(const number<Backend,ExpressionTemplates> & v,const Policy & pol)2102 inline int itrunc(const number<Backend, ExpressionTemplates>& v, const Policy& pol)
2103 {
2104 number<Backend, ExpressionTemplates> r = trunc(v, pol);
2105 if((r > (std::numeric_limits<int>::max)()) || r < (std::numeric_limits<int>::min)() || !(boost::math::isfinite)(v))
2106 return boost::math::policies::raise_rounding_error("boost::multiprecision::itrunc<%1%>(%1%)", 0, v, 0, pol);
2107 return r.template convert_to<int>();
2108 }
2109 template <class Backend, expression_template_option ExpressionTemplates>
itrunc(const number<Backend,ExpressionTemplates> & v)2110 inline int itrunc(const number<Backend, ExpressionTemplates>& v)
2111 {
2112 return itrunc(v, boost::math::policies::policy<>());
2113 }
2114 template <class tag, class A1, class A2, class A3, class A4, class Policy>
ltrunc(const detail::expression<tag,A1,A2,A3,A4> & v,const Policy & pol)2115 inline long ltrunc(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2116 {
2117 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2118 number_type r = trunc(v, pol);
2119 if((r > (std::numeric_limits<long>::max)()) || r < (std::numeric_limits<long>::min)() || !(boost::math::isfinite)(v))
2120 return boost::math::policies::raise_rounding_error("boost::multiprecision::ltrunc<%1%>(%1%)", 0, number_type(v), 0L, pol);
2121 return r.template convert_to<long>();
2122 }
2123 template <class tag, class A1, class A2, class A3, class A4>
ltrunc(const detail::expression<tag,A1,A2,A3,A4> & v)2124 inline long ltrunc(const detail::expression<tag, A1, A2, A3, A4>& v)
2125 {
2126 return ltrunc(v, boost::math::policies::policy<>());
2127 }
2128 template <class T, expression_template_option ExpressionTemplates, class Policy>
ltrunc(const number<T,ExpressionTemplates> & v,const Policy & pol)2129 inline long ltrunc(const number<T, ExpressionTemplates>& v, const Policy& pol)
2130 {
2131 number<T, ExpressionTemplates> r = trunc(v, pol);
2132 if((r > (std::numeric_limits<long>::max)()) || r < (std::numeric_limits<long>::min)() || !(boost::math::isfinite)(v))
2133 return boost::math::policies::raise_rounding_error("boost::multiprecision::ltrunc<%1%>(%1%)", 0, v, 0L, pol);
2134 return r.template convert_to<long>();
2135 }
2136 template <class T, expression_template_option ExpressionTemplates>
ltrunc(const number<T,ExpressionTemplates> & v)2137 inline long ltrunc(const number<T, ExpressionTemplates>& v)
2138 {
2139 return ltrunc(v, boost::math::policies::policy<>());
2140 }
2141 #ifndef BOOST_NO_LONG_LONG
2142 template <class tag, class A1, class A2, class A3, class A4, class Policy>
lltrunc(const detail::expression<tag,A1,A2,A3,A4> & v,const Policy & pol)2143 inline boost::long_long_type lltrunc(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2144 {
2145 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2146 number_type r = trunc(v, pol);
2147 if((r > (std::numeric_limits<boost::long_long_type>::max)()) || r < (std::numeric_limits<boost::long_long_type>::min)() || !(boost::math::isfinite)(v))
2148 return boost::math::policies::raise_rounding_error("boost::multiprecision::lltrunc<%1%>(%1%)", 0, number_type(v), 0LL, pol);
2149 return r.template convert_to<boost::long_long_type>();
2150 }
2151 template <class tag, class A1, class A2, class A3, class A4>
lltrunc(const detail::expression<tag,A1,A2,A3,A4> & v)2152 inline boost::long_long_type lltrunc(const detail::expression<tag, A1, A2, A3, A4>& v)
2153 {
2154 return lltrunc(v, boost::math::policies::policy<>());
2155 }
2156 template <class T, expression_template_option ExpressionTemplates, class Policy>
lltrunc(const number<T,ExpressionTemplates> & v,const Policy & pol)2157 inline boost::long_long_type lltrunc(const number<T, ExpressionTemplates>& v, const Policy& pol)
2158 {
2159 number<T, ExpressionTemplates> r = trunc(v, pol);
2160 if((r > (std::numeric_limits<boost::long_long_type>::max)()) || r < (std::numeric_limits<boost::long_long_type>::min)() || !(boost::math::isfinite)(v))
2161 return boost::math::policies::raise_rounding_error("boost::multiprecision::lltrunc<%1%>(%1%)", 0, v, 0LL, pol);
2162 return r.template convert_to<boost::long_long_type>();
2163 }
2164 template <class T, expression_template_option ExpressionTemplates>
lltrunc(const number<T,ExpressionTemplates> & v)2165 inline boost::long_long_type lltrunc(const number<T, ExpressionTemplates>& v)
2166 {
2167 return lltrunc(v, boost::math::policies::policy<>());
2168 }
2169 #endif
2170 template <class tag, class A1, class A2, class A3, class A4, class Policy>
round(const detail::expression<tag,A1,A2,A3,A4> & v,const Policy & pol)2171 inline typename detail::expression<tag, A1, A2, A3, A4>::result_type round(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2172 {
2173 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2174 return BOOST_MP_MOVE(round(static_cast<number_type>(v), pol));
2175 }
2176 template <class T, expression_template_option ExpressionTemplates, class Policy>
round(const number<T,ExpressionTemplates> & v,const Policy &)2177 inline number<T, ExpressionTemplates> round(const number<T, ExpressionTemplates>& v, const Policy&)
2178 {
2179 using default_ops::eval_round;
2180 number<T, ExpressionTemplates> result;
2181 eval_round(result.backend(), v.backend());
2182 return BOOST_MP_MOVE(result);
2183 }
2184
2185 template <class tag, class A1, class A2, class A3, class A4, class Policy>
iround(const detail::expression<tag,A1,A2,A3,A4> & v,const Policy & pol)2186 inline int iround(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2187 {
2188 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2189 number_type r = round(v, pol);
2190 if((r > (std::numeric_limits<int>::max)()) || r < (std::numeric_limits<int>::min)() || !(boost::math::isfinite)(v))
2191 return boost::math::policies::raise_rounding_error("boost::multiprecision::iround<%1%>(%1%)", 0, number_type(v), 0, pol);
2192 return r.template convert_to<int>();
2193 }
2194 template <class tag, class A1, class A2, class A3, class A4>
iround(const detail::expression<tag,A1,A2,A3,A4> & v)2195 inline int iround(const detail::expression<tag, A1, A2, A3, A4>& v)
2196 {
2197 return iround(v, boost::math::policies::policy<>());
2198 }
2199 template <class T, expression_template_option ExpressionTemplates, class Policy>
iround(const number<T,ExpressionTemplates> & v,const Policy & pol)2200 inline int iround(const number<T, ExpressionTemplates>& v, const Policy& pol)
2201 {
2202 number<T, ExpressionTemplates> r = round(v, pol);
2203 if((r > (std::numeric_limits<int>::max)()) || r < (std::numeric_limits<int>::min)() || !(boost::math::isfinite)(v))
2204 return boost::math::policies::raise_rounding_error("boost::multiprecision::iround<%1%>(%1%)", 0, v, 0, pol);
2205 return r.template convert_to<int>();
2206 }
2207 template <class T, expression_template_option ExpressionTemplates>
iround(const number<T,ExpressionTemplates> & v)2208 inline int iround(const number<T, ExpressionTemplates>& v)
2209 {
2210 return iround(v, boost::math::policies::policy<>());
2211 }
2212 template <class tag, class A1, class A2, class A3, class A4, class Policy>
lround(const detail::expression<tag,A1,A2,A3,A4> & v,const Policy & pol)2213 inline long lround(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2214 {
2215 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2216 number_type r = round(v, pol);
2217 if((r > (std::numeric_limits<long>::max)()) || r < (std::numeric_limits<long>::min)() || !(boost::math::isfinite)(v))
2218 return boost::math::policies::raise_rounding_error("boost::multiprecision::lround<%1%>(%1%)", 0, number_type(v), 0L, pol);
2219 return r.template convert_to<long>();
2220 }
2221 template <class tag, class A1, class A2, class A3, class A4>
lround(const detail::expression<tag,A1,A2,A3,A4> & v)2222 inline long lround(const detail::expression<tag, A1, A2, A3, A4>& v)
2223 {
2224 return lround(v, boost::math::policies::policy<>());
2225 }
2226 template <class T, expression_template_option ExpressionTemplates, class Policy>
lround(const number<T,ExpressionTemplates> & v,const Policy & pol)2227 inline long lround(const number<T, ExpressionTemplates>& v, const Policy& pol)
2228 {
2229 number<T, ExpressionTemplates> r = round(v, pol);
2230 if((r > (std::numeric_limits<long>::max)()) || r < (std::numeric_limits<long>::min)() || !(boost::math::isfinite)(v))
2231 return boost::math::policies::raise_rounding_error("boost::multiprecision::lround<%1%>(%1%)", 0, v, 0L, pol);
2232 return r.template convert_to<long>();
2233 }
2234 template <class T, expression_template_option ExpressionTemplates>
lround(const number<T,ExpressionTemplates> & v)2235 inline long lround(const number<T, ExpressionTemplates>& v)
2236 {
2237 return lround(v, boost::math::policies::policy<>());
2238 }
2239 #ifndef BOOST_NO_LONG_LONG
2240 template <class tag, class A1, class A2, class A3, class A4, class Policy>
llround(const detail::expression<tag,A1,A2,A3,A4> & v,const Policy & pol)2241 inline boost::long_long_type llround(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2242 {
2243 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2244 number_type r = round(v, pol);
2245 if((r > (std::numeric_limits<boost::long_long_type>::max)()) || r < (std::numeric_limits<boost::long_long_type>::min)() || !(boost::math::isfinite)(v))
2246 return boost::math::policies::raise_rounding_error("boost::multiprecision::iround<%1%>(%1%)", 0, number_type(v), 0LL, pol);
2247 return r.template convert_to<boost::long_long_type>();
2248 }
2249 template <class tag, class A1, class A2, class A3, class A4>
llround(const detail::expression<tag,A1,A2,A3,A4> & v)2250 inline boost::long_long_type llround(const detail::expression<tag, A1, A2, A3, A4>& v)
2251 {
2252 return llround(v, boost::math::policies::policy<>());
2253 }
2254 template <class T, expression_template_option ExpressionTemplates, class Policy>
llround(const number<T,ExpressionTemplates> & v,const Policy & pol)2255 inline boost::long_long_type llround(const number<T, ExpressionTemplates>& v, const Policy& pol)
2256 {
2257 number<T, ExpressionTemplates> r = round(v, pol);
2258 if((r > (std::numeric_limits<boost::long_long_type>::max)()) || r < (std::numeric_limits<boost::long_long_type>::min)() || !(boost::math::isfinite)(v))
2259 return boost::math::policies::raise_rounding_error("boost::multiprecision::iround<%1%>(%1%)", 0, v, 0LL, pol);
2260 return r.template convert_to<boost::long_long_type>();
2261 }
2262 template <class T, expression_template_option ExpressionTemplates>
llround(const number<T,ExpressionTemplates> & v)2263 inline boost::long_long_type llround(const number<T, ExpressionTemplates>& v)
2264 {
2265 return llround(v, boost::math::policies::policy<>());
2266 }
2267 #endif
2268 //
2269 // frexp does not return an expression template since we require the
2270 // integer argument to be evaluated even if the returned value is
2271 // not assigned to anything...
2272 //
2273 template <class T, expression_template_option ExpressionTemplates>
frexp(const number<T,ExpressionTemplates> & v,short * pint)2274 inline typename enable_if_c<number_category<T>::value == number_kind_floating_point, number<T, ExpressionTemplates> >::type frexp(const number<T, ExpressionTemplates>& v, short* pint)
2275 {
2276 using default_ops::eval_frexp;
2277 number<T, ExpressionTemplates> result;
2278 eval_frexp(result.backend(), v.backend(), pint);
2279 return BOOST_MP_MOVE(result);
2280 }
2281 template <class tag, class A1, class A2, class A3, class A4>
2282 inline typename enable_if_c<number_category<typename detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_floating_point, typename detail::expression<tag, A1, A2, A3, A4>::result_type>::type
frexp(const detail::expression<tag,A1,A2,A3,A4> & v,short * pint)2283 frexp(const detail::expression<tag, A1, A2, A3, A4>& v, short* pint)
2284 {
2285 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2286 return BOOST_MP_MOVE(frexp(static_cast<number_type>(v), pint));
2287 }
2288 template <class T, expression_template_option ExpressionTemplates>
frexp(const number<T,ExpressionTemplates> & v,int * pint)2289 inline typename enable_if_c<number_category<T>::value == number_kind_floating_point, number<T, ExpressionTemplates> >::type frexp(const number<T, ExpressionTemplates>& v, int* pint)
2290 {
2291 using default_ops::eval_frexp;
2292 number<T, ExpressionTemplates> result;
2293 eval_frexp(result.backend(), v.backend(), pint);
2294 return BOOST_MP_MOVE(result);
2295 }
2296 template <class tag, class A1, class A2, class A3, class A4>
2297 inline typename enable_if_c<number_category<typename detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_floating_point, typename detail::expression<tag, A1, A2, A3, A4>::result_type>::type
frexp(const detail::expression<tag,A1,A2,A3,A4> & v,int * pint)2298 frexp(const detail::expression<tag, A1, A2, A3, A4>& v, int* pint)
2299 {
2300 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2301 return BOOST_MP_MOVE(frexp(static_cast<number_type>(v), pint));
2302 }
2303 template <class T, expression_template_option ExpressionTemplates>
frexp(const number<T,ExpressionTemplates> & v,long * pint)2304 inline typename enable_if_c<number_category<T>::value == number_kind_floating_point, number<T, ExpressionTemplates> >::type frexp(const number<T, ExpressionTemplates>& v, long* pint)
2305 {
2306 using default_ops::eval_frexp;
2307 number<T, ExpressionTemplates> result;
2308 eval_frexp(result.backend(), v.backend(), pint);
2309 return BOOST_MP_MOVE(result);
2310 }
2311 template <class tag, class A1, class A2, class A3, class A4>
2312 inline typename enable_if_c<number_category<typename detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_floating_point, typename detail::expression<tag, A1, A2, A3, A4>::result_type>::type
frexp(const detail::expression<tag,A1,A2,A3,A4> & v,long * pint)2313 frexp(const detail::expression<tag, A1, A2, A3, A4>& v, long* pint)
2314 {
2315 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2316 return BOOST_MP_MOVE(frexp(static_cast<number_type>(v), pint));
2317 }
2318 template <class T, expression_template_option ExpressionTemplates>
frexp(const number<T,ExpressionTemplates> & v,boost::long_long_type * pint)2319 inline typename enable_if_c<number_category<T>::value == number_kind_floating_point, number<T, ExpressionTemplates> >::type frexp(const number<T, ExpressionTemplates>& v, boost::long_long_type* pint)
2320 {
2321 using default_ops::eval_frexp;
2322 number<T, ExpressionTemplates> result;
2323 eval_frexp(result.backend(), v.backend(), pint);
2324 return BOOST_MP_MOVE(result);
2325 }
2326 template <class tag, class A1, class A2, class A3, class A4>
2327 inline typename enable_if_c<number_category<typename detail::expression<tag, A1, A2, A3, A4>::result_type>::value == number_kind_floating_point, typename detail::expression<tag, A1, A2, A3, A4>::result_type>::type
frexp(const detail::expression<tag,A1,A2,A3,A4> & v,boost::long_long_type * pint)2328 frexp(const detail::expression<tag, A1, A2, A3, A4>& v, boost::long_long_type* pint)
2329 {
2330 typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2331 return BOOST_MP_MOVE(frexp(static_cast<number_type>(v), pint));
2332 }
2333 //
2334 // modf does not return an expression template since we require the
2335 // second argument to be evaluated even if the returned value is
2336 // not assigned to anything...
2337 //
2338 template <class T, expression_template_option ExpressionTemplates>
modf(const number<T,ExpressionTemplates> & v,number<T,ExpressionTemplates> * pipart)2339 inline typename enable_if_c<number_category<T>::value == number_kind_floating_point, number<T, ExpressionTemplates> >::type modf(const number<T, ExpressionTemplates>& v, number<T, ExpressionTemplates>* pipart)
2340 {
2341 using default_ops::eval_modf;
2342 number<T, ExpressionTemplates> result;
2343 eval_modf(result.backend(), v.backend(), pipart ? &pipart->backend() : 0);
2344 return BOOST_MP_MOVE(result);
2345 }
2346 template <class T, expression_template_option ExpressionTemplates, class tag, class A1, class A2, class A3, class A4>
modf(const detail::expression<tag,A1,A2,A3,A4> & v,number<T,ExpressionTemplates> * pipart)2347 inline typename enable_if_c<number_category<T>::value == number_kind_floating_point, number<T, ExpressionTemplates> >::type modf(const detail::expression<tag, A1, A2, A3, A4>& v, number<T, ExpressionTemplates>* pipart)
2348 {
2349 using default_ops::eval_modf;
2350 number<T, ExpressionTemplates> result, arg(v);
2351 eval_modf(result.backend(), arg.backend(), pipart ? &pipart->backend() : 0);
2352 return BOOST_MP_MOVE(result);
2353 }
2354
2355 //
2356 // Integer square root:
2357 //
2358 template <class B, expression_template_option ExpressionTemplates>
2359 inline typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, ExpressionTemplates> >::type
sqrt(const number<B,ExpressionTemplates> & x)2360 sqrt(const number<B, ExpressionTemplates>& x)
2361 {
2362 using default_ops::eval_integer_sqrt;
2363 number<B, ExpressionTemplates> s, r;
2364 eval_integer_sqrt(s.backend(), r.backend(), x.backend());
2365 return s;
2366 }
2367 //
2368 // fma:
2369 //
2370
2371 namespace default_ops {
2372
2373 struct fma_func
2374 {
2375 template <class B, class T, class U, class V>
operator ()boost::multiprecision::default_ops::fma_func2376 void operator()(B& result, const T& a, const U& b, const V& c)const
2377 {
2378 eval_multiply_add(result, a, b, c);
2379 }
2380 };
2381
2382
2383 }
2384
2385 template <class Backend, class U, class V>
2386 inline typename enable_if<
2387 mpl::and_<
2388 mpl::bool_<number_category<number<Backend, et_on> >::value == number_kind_floating_point>,
2389 mpl::or_<
2390 is_number<U>,
2391 is_number_expression<U>,
2392 is_arithmetic<U>
2393 >,
2394 mpl::or_<
2395 is_number<V>,
2396 is_number_expression<V>,
2397 is_arithmetic<V>
2398 >
2399 >,
2400 detail::expression<detail::function, default_ops::fma_func, number<Backend, et_on>, U, V>
2401 >::type
fma(const number<Backend,et_on> & a,const U & b,const V & c)2402 fma(const number<Backend, et_on>& a, const U& b, const V& c)
2403 {
2404 return detail::expression<detail::function, default_ops::fma_func, number<Backend, et_on>, U, V>(
2405 default_ops::fma_func(), a, b, c);
2406 }
2407
2408 template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class U, class V>
2409 inline typename enable_if<
2410 mpl::and_<
2411 mpl::bool_<number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point>,
2412 mpl::or_<
2413 is_number<U>,
2414 is_number_expression<U>,
2415 is_arithmetic<U>
2416 >,
2417 mpl::or_<
2418 is_number<V>,
2419 is_number_expression<V>,
2420 is_arithmetic<V>
2421 >
2422 >,
2423 detail::expression<detail::function, default_ops::fma_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, V>
2424 >::type
fma(const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & a,const U & b,const V & c)2425 fma(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& a, const U& b, const V& c)
2426 {
2427 return detail::expression<detail::function, default_ops::fma_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, V>(
2428 default_ops::fma_func(), a, b, c);
2429 }
2430
2431 template <class Backend, class U, class V>
2432 inline typename enable_if<
2433 mpl::and_<
2434 mpl::bool_<number_category<number<Backend, et_off> >::value == number_kind_floating_point>,
2435 mpl::or_<
2436 is_number<U>,
2437 is_number_expression<U>,
2438 is_arithmetic<U>
2439 >,
2440 mpl::or_<
2441 is_number<V>,
2442 is_number_expression<V>,
2443 is_arithmetic<V>
2444 >
2445 >,
2446 number<Backend, et_off>
2447 >::type
fma(const number<Backend,et_off> & a,const U & b,const V & c)2448 fma(const number<Backend, et_off>& a, const U& b, const V& c)
2449 {
2450 using default_ops::eval_multiply_add;
2451 number<Backend, et_off> result;
2452 eval_multiply_add(result.backend(), number<Backend, et_off>::canonical_value(a), number<Backend, et_off>::canonical_value(b), number<Backend, et_off>::canonical_value(c));
2453 return BOOST_MP_MOVE(result);
2454 }
2455
2456 template <class U, class Backend, class V>
2457 inline typename enable_if<
2458 mpl::and_<
2459 mpl::bool_<number_category<number<Backend, et_on> >::value == number_kind_floating_point>,
2460 is_arithmetic<U>,
2461 mpl::or_<
2462 is_number<V>,
2463 is_number_expression<V>,
2464 is_arithmetic<V>
2465 >
2466 >,
2467 detail::expression<detail::function, default_ops::fma_func, U, number<Backend, et_on>, V>
2468 >::type
fma(const U & a,const number<Backend,et_on> & b,const V & c)2469 fma(const U& a, const number<Backend, et_on>& b, const V& c)
2470 {
2471 return detail::expression<detail::function, default_ops::fma_func, U, number<Backend, et_on>, V>(
2472 default_ops::fma_func(), a, b, c);
2473 }
2474
2475 template <class U, class tag, class Arg1, class Arg2, class Arg3, class Arg4, class V>
2476 inline typename enable_if<
2477 mpl::and_<
2478 mpl::bool_<number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point>,
2479 is_arithmetic<U>,
2480 mpl::or_<
2481 is_number<V>,
2482 is_number_expression<V>,
2483 is_arithmetic<V>
2484 >
2485 >,
2486 detail::expression<detail::function, default_ops::fma_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, V>
2487 >::type
fma(const U & a,const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & b,const V & c)2488 fma(const U& a, const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& b, const V& c)
2489 {
2490 return detail::expression<detail::function, default_ops::fma_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, V>(
2491 default_ops::fma_func(), a, b, c);
2492 }
2493
2494 template <class U, class Backend, class V>
2495 inline typename enable_if<
2496 mpl::and_<
2497 mpl::bool_<number_category<number<Backend, et_off> >::value == number_kind_floating_point>,
2498 is_arithmetic<U>,
2499 mpl::or_<
2500 is_number<V>,
2501 is_number_expression<V>,
2502 is_arithmetic<V>
2503 >
2504 >,
2505 number<Backend, et_off>
2506 >::type
fma(const U & a,const number<Backend,et_off> & b,const V & c)2507 fma(const U& a, const number<Backend, et_off>& b, const V& c)
2508 {
2509 using default_ops::eval_multiply_add;
2510 number<Backend, et_off> result;
2511 eval_multiply_add(result.backend(), number<Backend, et_off>::canonical_value(a), number<Backend, et_off>::canonical_value(b), number<Backend, et_off>::canonical_value(c));
2512 return BOOST_MP_MOVE(result);
2513 }
2514
2515 template <class U, class V, class Backend>
2516 inline typename enable_if<
2517 mpl::and_<
2518 mpl::bool_<number_category<number<Backend, et_on> >::value == number_kind_floating_point>,
2519 is_arithmetic<U>,
2520 is_arithmetic<V>
2521 >,
2522 detail::expression<detail::function, default_ops::fma_func, U, V, number<Backend, et_on> >
2523 >::type
fma(const U & a,const V & b,const number<Backend,et_on> & c)2524 fma(const U& a, const V& b, const number<Backend, et_on>& c)
2525 {
2526 return detail::expression<detail::function, default_ops::fma_func, U, V, number<Backend, et_on> >(
2527 default_ops::fma_func(), a, b, c);
2528 }
2529
2530 template <class U, class V, class tag, class Arg1, class Arg2, class Arg3, class Arg4>
2531 inline typename enable_if<
2532 mpl::and_<
2533 mpl::bool_<number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point>,
2534 is_arithmetic<U>,
2535 is_arithmetic<V>
2536 >,
2537 detail::expression<detail::function, default_ops::fma_func, U, V, detail::expression<tag, Arg1, Arg2, Arg3, Arg4> >
2538 >::type
fma(const U & a,const V & b,const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & c)2539 fma(const U& a, const V& b, const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& c)
2540 {
2541 return detail::expression<detail::function, default_ops::fma_func, U, V, detail::expression<tag, Arg1, Arg2, Arg3, Arg4> >(
2542 default_ops::fma_func(), a, b, c);
2543 }
2544
2545 template <class U, class V, class Backend>
2546 inline typename enable_if<
2547 mpl::and_<
2548 mpl::bool_<number_category<number<Backend, et_off> >::value == number_kind_floating_point>,
2549 is_arithmetic<U>,
2550 is_arithmetic<V>
2551 >,
2552 number<Backend, et_off>
2553 >::type
fma(const U & a,const V & b,const number<Backend,et_off> & c)2554 fma(const U& a, const V& b, const number<Backend, et_off>& c)
2555 {
2556 using default_ops::eval_multiply_add;
2557 number<Backend, et_off> result;
2558 eval_multiply_add(result.backend(), number<Backend, et_off>::canonical_value(a), number<Backend, et_off>::canonical_value(b), number<Backend, et_off>::canonical_value(c));
2559 return BOOST_MP_MOVE(result);
2560 }
2561
2562 namespace default_ops {
2563
2564 struct remquo_func
2565 {
2566 template <class B, class T, class U>
operator ()boost::multiprecision::default_ops::remquo_func2567 void operator()(B& result, const T& a, const U& b, int* pi)const
2568 {
2569 eval_remquo(result, a, b, pi);
2570 }
2571 };
2572
2573 }
2574
2575 template <class Backend, class U>
2576 inline typename enable_if_c<
2577 number_category<number<Backend, et_on> >::value == number_kind_floating_point,
2578 detail::expression<detail::function, default_ops::remquo_func, number<Backend, et_on>, U, int*>
2579 >::type
remquo(const number<Backend,et_on> & a,const U & b,int * pi)2580 remquo(const number<Backend, et_on>& a, const U& b, int* pi)
2581 {
2582 return detail::expression<detail::function, default_ops::remquo_func, number<Backend, et_on>, U, int*>(
2583 default_ops::remquo_func(), a, b, pi);
2584 }
2585
2586 template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class U>
2587 inline typename enable_if_c<
2588 number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point,
2589 detail::expression<detail::function, default_ops::remquo_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, int*>
2590 >::type
remquo(const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & a,const U & b,int * pi)2591 remquo(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& a, const U& b, int* pi)
2592 {
2593 return detail::expression<detail::function, default_ops::remquo_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, int*>(
2594 default_ops::remquo_func(), a, b, pi);
2595 }
2596
2597 template <class U, class Backend>
2598 inline typename enable_if_c<
2599 (number_category<number<Backend, et_on> >::value == number_kind_floating_point)
2600 && !is_number<U>::value && !is_number_expression<U>::value,
2601 detail::expression<detail::function, default_ops::remquo_func, U, number<Backend, et_on>, int*>
2602 >::type
remquo(const U & a,const number<Backend,et_on> & b,int * pi)2603 remquo(const U& a, const number<Backend, et_on>& b, int* pi)
2604 {
2605 return detail::expression<detail::function, default_ops::remquo_func, U, number<Backend, et_on>, int*>(
2606 default_ops::remquo_func(), a, b, pi);
2607 }
2608
2609 template <class U, class tag, class Arg1, class Arg2, class Arg3, class Arg4>
2610 inline typename enable_if_c<
2611 (number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point)
2612 && !is_number<U>::value && !is_number_expression<U>::value,
2613 detail::expression<detail::function, default_ops::remquo_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, int*>
2614 >::type
remquo(const U & a,const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & b,int * pi)2615 remquo(const U& a, const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& b, int* pi)
2616 {
2617 return detail::expression<detail::function, default_ops::remquo_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, int*>(
2618 default_ops::remquo_func(), a, b, pi);
2619 }
2620
2621 template <class Backend, class U>
2622 inline typename enable_if_c<
2623 number_category<number<Backend, et_on> >::value == number_kind_floating_point,
2624 number<Backend, et_off>
2625 >::type
remquo(const number<Backend,et_off> & a,const U & b,int * pi)2626 remquo(const number<Backend, et_off>& a, const U& b, int* pi)
2627 {
2628 using default_ops::eval_remquo;
2629 number<Backend, et_off> result;
2630 eval_remquo(result.backend(), a.backend(), number<Backend, et_off>::canonical_value(b), pi);
2631 return BOOST_MP_MOVE(result);
2632 }
2633 template <class U, class Backend>
2634 inline typename enable_if_c<
2635 (number_category<number<Backend, et_on> >::value == number_kind_floating_point)
2636 && !is_number<U>::value && !is_number_expression<U>::value,
2637 number<Backend, et_off>
2638 >::type
remquo(const U & a,const number<Backend,et_off> & b,int * pi)2639 remquo(const U& a, const number<Backend, et_off>& b, int* pi)
2640 {
2641 using default_ops::eval_remquo;
2642 number<Backend, et_off> result;
2643 eval_remquo(result.backend(), number<Backend, et_off>::canonical_value(a), b.backend(), pi);
2644 return BOOST_MP_MOVE(result);
2645 }
2646
2647
2648 template <class B, expression_template_option ExpressionTemplates>
2649 inline typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, ExpressionTemplates> >::type
sqrt(const number<B,ExpressionTemplates> & x,number<B,ExpressionTemplates> & r)2650 sqrt(const number<B, ExpressionTemplates>& x, number<B, ExpressionTemplates>& r)
2651 {
2652 using default_ops::eval_integer_sqrt;
2653 number<B, ExpressionTemplates> s;
2654 eval_integer_sqrt(s.backend(), r.backend(), x.backend());
2655 return s;
2656 }
2657
2658 #define UNARY_OP_FUNCTOR(func, category)\
2659 namespace detail{\
2660 template <class Backend> \
2661 struct BOOST_JOIN(func, _funct)\
2662 {\
2663 void operator()(Backend& result, const Backend& arg)const\
2664 {\
2665 using default_ops::BOOST_JOIN(eval_,func);\
2666 BOOST_JOIN(eval_,func)(result, arg);\
2667 }\
2668 };\
2669 \
2670 }\
2671 \
2672 template <class tag, class A1, class A2, class A3, class A4> \
2673 inline typename enable_if_c<number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category,\
2674 detail::expression<\
2675 detail::function\
2676 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2677 , detail::expression<tag, A1, A2, A3, A4> > \
2678 >::type \
2679 func(const detail::expression<tag, A1, A2, A3, A4>& arg)\
2680 {\
2681 return detail::expression<\
2682 detail::function\
2683 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2684 , detail::expression<tag, A1, A2, A3, A4> \
2685 > (\
2686 detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2687 , arg \
2688 );\
2689 }\
2690 template <class Backend> \
2691 inline typename enable_if_c<number_category<Backend>::value == category,\
2692 detail::expression<\
2693 detail::function\
2694 , detail::BOOST_JOIN(func, _funct)<Backend> \
2695 , number<Backend, et_on> > \
2696 >::type \
2697 func(const number<Backend, et_on>& arg)\
2698 {\
2699 return detail::expression<\
2700 detail::function\
2701 , detail::BOOST_JOIN(func, _funct)<Backend> \
2702 , number<Backend, et_on> \
2703 >(\
2704 detail::BOOST_JOIN(func, _funct)<Backend>() \
2705 , arg \
2706 );\
2707 }\
2708 template <class Backend> \
2709 inline typename boost::enable_if_c<\
2710 boost::multiprecision::number_category<Backend>::value == category,\
2711 number<Backend, et_off> >::type \
2712 func(const number<Backend, et_off>& arg)\
2713 {\
2714 number<Backend, et_off> result;\
2715 using default_ops::BOOST_JOIN(eval_,func);\
2716 BOOST_JOIN(eval_,func)(result.backend(), arg.backend());\
2717 return BOOST_MP_MOVE(result);\
2718 }
2719
2720 #define BINARY_OP_FUNCTOR(func, category)\
2721 namespace detail{\
2722 template <class Backend> \
2723 struct BOOST_JOIN(func, _funct)\
2724 {\
2725 void operator()(Backend& result, const Backend& arg, const Backend& a)const\
2726 {\
2727 using default_ops:: BOOST_JOIN(eval_,func);\
2728 BOOST_JOIN(eval_,func)(result, arg, a);\
2729 }\
2730 template <class Arithmetic> \
2731 void operator()(Backend& result, const Backend& arg, const Arithmetic& a)const\
2732 {\
2733 using default_ops:: BOOST_JOIN(eval_,func);\
2734 BOOST_JOIN(eval_,func)(result, arg, a);\
2735 }\
2736 template <class Arithmetic> \
2737 void operator()(Backend& result, const Arithmetic& arg, const Backend& a)const\
2738 {\
2739 using default_ops:: BOOST_JOIN(eval_,func);\
2740 BOOST_JOIN(eval_,func)(result, arg, a);\
2741 }\
2742 };\
2743 \
2744 }\
2745 template <class Backend> \
2746 inline typename enable_if_c<number_category<Backend>::value == category,\
2747 detail::expression<\
2748 detail::function\
2749 , detail::BOOST_JOIN(func, _funct)<Backend> \
2750 , number<Backend, et_on> \
2751 , number<Backend, et_on> > \
2752 >::type \
2753 func(const number<Backend, et_on>& arg, const number<Backend, et_on>& a)\
2754 {\
2755 return detail::expression<\
2756 detail::function\
2757 , detail::BOOST_JOIN(func, _funct)<Backend> \
2758 , number<Backend, et_on> \
2759 , number<Backend, et_on> \
2760 >(\
2761 detail::BOOST_JOIN(func, _funct)<Backend>() \
2762 , arg,\
2763 a\
2764 );\
2765 }\
2766 template <class Backend, class tag, class A1, class A2, class A3, class A4> \
2767 inline typename enable_if_c<\
2768 (number_category<Backend>::value == category) && (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2769 detail::expression<\
2770 detail::function\
2771 , detail::BOOST_JOIN(func, _funct)<Backend> \
2772 , number<Backend, et_on> \
2773 , detail::expression<tag, A1, A2, A3, A4> > \
2774 >::type \
2775 func(const number<Backend, et_on>& arg, const detail::expression<tag, A1, A2, A3, A4>& a)\
2776 {\
2777 return detail::expression<\
2778 detail::function\
2779 , detail::BOOST_JOIN(func, _funct)<Backend> \
2780 , number<Backend, et_on> \
2781 , detail::expression<tag, A1, A2, A3, A4> \
2782 >(\
2783 detail::BOOST_JOIN(func, _funct)<Backend>() \
2784 , arg,\
2785 a\
2786 );\
2787 }\
2788 template <class tag, class A1, class A2, class A3, class A4, class Backend> \
2789 inline typename enable_if_c<\
2790 (number_category<Backend>::value == category) && (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2791 detail::expression<\
2792 detail::function\
2793 , detail::BOOST_JOIN(func, _funct)<Backend> \
2794 , detail::expression<tag, A1, A2, A3, A4> \
2795 , number<Backend, et_on> > \
2796 >::type \
2797 func(const detail::expression<tag, A1, A2, A3, A4>& arg, const number<Backend, et_on>& a)\
2798 {\
2799 return detail::expression<\
2800 detail::function\
2801 , detail::BOOST_JOIN(func, _funct)<Backend> \
2802 , detail::expression<tag, A1, A2, A3, A4> \
2803 , number<Backend, et_on> \
2804 >(\
2805 detail::BOOST_JOIN(func, _funct)<Backend>() \
2806 , arg,\
2807 a\
2808 );\
2809 }\
2810 template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b> \
2811 inline typename enable_if_c<\
2812 (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category) && (number_category<detail::expression<tagb, A1b, A2b, A3b, A4b> >::value == category),\
2813 detail::expression<\
2814 detail::function\
2815 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2816 , detail::expression<tag, A1, A2, A3, A4> \
2817 , detail::expression<tagb, A1b, A2b, A3b, A4b> > \
2818 >::type \
2819 func(const detail::expression<tag, A1, A2, A3, A4>& arg, const detail::expression<tagb, A1b, A2b, A3b, A4b>& a)\
2820 {\
2821 return detail::expression<\
2822 detail::function\
2823 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2824 , detail::expression<tag, A1, A2, A3, A4> \
2825 , detail::expression<tagb, A1b, A2b, A3b, A4b> \
2826 >(\
2827 detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2828 , arg,\
2829 a\
2830 );\
2831 }\
2832 template <class Backend, class Arithmetic> \
2833 inline typename enable_if_c<\
2834 is_arithmetic<Arithmetic>::value && (number_category<Backend>::value == category),\
2835 detail::expression<\
2836 detail::function\
2837 , detail::BOOST_JOIN(func, _funct)<Backend> \
2838 , number<Backend, et_on> \
2839 , Arithmetic\
2840 > \
2841 >::type \
2842 func(const number<Backend, et_on>& arg, const Arithmetic& a)\
2843 {\
2844 return detail::expression<\
2845 detail::function\
2846 , detail::BOOST_JOIN(func, _funct)<Backend> \
2847 , number<Backend, et_on> \
2848 , Arithmetic\
2849 >(\
2850 detail::BOOST_JOIN(func, _funct)<Backend>() \
2851 , arg,\
2852 a\
2853 );\
2854 }\
2855 template <class tag, class A1, class A2, class A3, class A4, class Arithmetic> \
2856 inline typename enable_if_c<\
2857 is_arithmetic<Arithmetic>::value && (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2858 detail::expression<\
2859 detail::function\
2860 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2861 , detail::expression<tag, A1, A2, A3, A4> \
2862 , Arithmetic\
2863 > \
2864 >::type \
2865 func(const detail::expression<tag, A1, A2, A3, A4>& arg, const Arithmetic& a)\
2866 {\
2867 return detail::expression<\
2868 detail::function\
2869 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2870 , detail::expression<tag, A1, A2, A3, A4> \
2871 , Arithmetic\
2872 >(\
2873 detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2874 , arg,\
2875 a\
2876 );\
2877 }\
2878 template <class Backend, class Arithmetic> \
2879 inline typename enable_if_c<\
2880 is_arithmetic<Arithmetic>::value && (number_category<Backend>::value == category),\
2881 detail::expression<\
2882 detail::function\
2883 , detail::BOOST_JOIN(func, _funct)<Backend> \
2884 , Arithmetic \
2885 , number<Backend, et_on> \
2886 > \
2887 >::type \
2888 func(const Arithmetic& arg, const number<Backend, et_on>& a)\
2889 {\
2890 return detail::expression<\
2891 detail::function\
2892 , detail::BOOST_JOIN(func, _funct)<Backend> \
2893 , Arithmetic \
2894 , number<Backend, et_on> \
2895 >(\
2896 detail::BOOST_JOIN(func, _funct)<Backend>() \
2897 , arg,\
2898 a\
2899 );\
2900 }\
2901 template <class tag, class A1, class A2, class A3, class A4, class Arithmetic> \
2902 inline typename enable_if_c<\
2903 is_arithmetic<Arithmetic>::value && (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2904 detail::expression<\
2905 detail::function\
2906 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2907 , Arithmetic \
2908 , detail::expression<tag, A1, A2, A3, A4> \
2909 > \
2910 >::type \
2911 func(const Arithmetic& arg, const detail::expression<tag, A1, A2, A3, A4>& a)\
2912 {\
2913 return detail::expression<\
2914 detail::function\
2915 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2916 , Arithmetic \
2917 , detail::expression<tag, A1, A2, A3, A4> \
2918 >(\
2919 detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2920 , arg,\
2921 a\
2922 );\
2923 }\
2924 template <class Backend> \
2925 inline typename enable_if_c<(number_category<Backend>::value == category),\
2926 number<Backend, et_off> >::type \
2927 func(const number<Backend, et_off>& arg, const number<Backend, et_off>& a)\
2928 {\
2929 number<Backend, et_off> result;\
2930 using default_ops:: BOOST_JOIN(eval_,func);\
2931 BOOST_JOIN(eval_,func)(result.backend(), arg.backend(), a.backend());\
2932 return BOOST_MP_MOVE(result);\
2933 }\
2934 template <class Backend, class Arithmetic> \
2935 inline typename enable_if_c<\
2936 is_arithmetic<Arithmetic>::value && (number_category<Backend>::value == category),\
2937 number<Backend, et_off> \
2938 >::type \
2939 func(const number<Backend, et_off>& arg, const Arithmetic& a)\
2940 {\
2941 typedef typename detail::canonical<Arithmetic, Backend>::type canonical_type;\
2942 number<Backend, et_off> result;\
2943 using default_ops:: BOOST_JOIN(eval_,func);\
2944 BOOST_JOIN(eval_,func)(result.backend(), arg.backend(), static_cast<canonical_type>(a));\
2945 return BOOST_MP_MOVE(result);\
2946 }\
2947 template <class Backend, class Arithmetic> \
2948 inline typename enable_if_c<\
2949 is_arithmetic<Arithmetic>::value && (number_category<Backend>::value == category),\
2950 number<Backend, et_off> \
2951 >::type \
2952 func(const Arithmetic& a, const number<Backend, et_off>& arg)\
2953 {\
2954 typedef typename detail::canonical<Arithmetic, Backend>::type canonical_type;\
2955 number<Backend, et_off> result;\
2956 using default_ops:: BOOST_JOIN(eval_,func);\
2957 BOOST_JOIN(eval_,func)(result.backend(), static_cast<canonical_type>(a), arg.backend());\
2958 return BOOST_MP_MOVE(result);\
2959 }\
2960
2961
2962 #define HETERO_BINARY_OP_FUNCTOR_B(func, Arg2, category)\
2963 template <class tag, class A1, class A2, class A3, class A4> \
2964 inline typename enable_if_c<\
2965 (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2966 detail::expression<\
2967 detail::function\
2968 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2969 , detail::expression<tag, A1, A2, A3, A4> \
2970 , Arg2> \
2971 >::type \
2972 func(const detail::expression<tag, A1, A2, A3, A4>& arg, Arg2 const& a)\
2973 {\
2974 return detail::expression<\
2975 detail::function\
2976 , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2977 , detail::expression<tag, A1, A2, A3, A4> \
2978 , Arg2\
2979 >(\
2980 detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2981 , arg, a \
2982 );\
2983 }\
2984 template <class Backend> \
2985 inline typename enable_if_c<\
2986 (number_category<Backend>::value == category),\
2987 detail::expression<\
2988 detail::function\
2989 , detail::BOOST_JOIN(func, _funct)<Backend> \
2990 , number<Backend, et_on> \
2991 , Arg2> \
2992 >::type \
2993 func(const number<Backend, et_on>& arg, Arg2 const& a)\
2994 {\
2995 return detail::expression<\
2996 detail::function\
2997 , detail::BOOST_JOIN(func, _funct)<Backend> \
2998 , number<Backend, et_on> \
2999 , Arg2\
3000 >(\
3001 detail::BOOST_JOIN(func, _funct)<Backend>() \
3002 , arg,\
3003 a\
3004 );\
3005 }\
3006 template <class Backend> \
3007 inline typename enable_if_c<\
3008 (number_category<Backend>::value == category),\
3009 number<Backend, et_off> >::type \
3010 func(const number<Backend, et_off>& arg, Arg2 const& a)\
3011 {\
3012 number<Backend, et_off> result;\
3013 using default_ops:: BOOST_JOIN(eval_,func);\
3014 BOOST_JOIN(eval_,func)(result.backend(), arg.backend(), a);\
3015 return BOOST_MP_MOVE(result);\
3016 }\
3017
3018 #define HETERO_BINARY_OP_FUNCTOR(func, Arg2, category)\
3019 namespace detail{\
3020 template <class Backend> \
3021 struct BOOST_JOIN(func, _funct)\
3022 {\
3023 template <class Arg>\
3024 void operator()(Backend& result, Backend const& arg, Arg a)const\
3025 {\
3026 using default_ops:: BOOST_JOIN(eval_,func);\
3027 BOOST_JOIN(eval_,func)(result, arg, a);\
3028 }\
3029 };\
3030 \
3031 }\
3032 \
3033 HETERO_BINARY_OP_FUNCTOR_B(func, Arg2, category)
3034
3035 namespace detail{
3036 template <class Backend>
3037 struct abs_funct
3038 {
operator ()boost::multiprecision::detail::abs_funct3039 void operator()(Backend& result, const Backend& arg)const
3040 {
3041 using default_ops::eval_abs;
3042 eval_abs(result, arg);
3043 }
3044 };
3045
3046 }
3047
3048 template <class tag, class A1, class A2, class A3, class A4>
3049 inline detail::expression<
3050 detail::function
3051 , detail::abs_funct<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>
3052 , detail::expression<tag, A1, A2, A3, A4> >
abs(const detail::expression<tag,A1,A2,A3,A4> & arg)3053 abs(const detail::expression<tag, A1, A2, A3, A4>& arg)
3054 {
3055 return detail::expression<
3056 detail::function
3057 , detail::abs_funct<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>
3058 , detail::expression<tag, A1, A2, A3, A4>
3059 > (
3060 detail::abs_funct<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>()
3061 , arg
3062 );
3063 }
3064 template <class Backend>
3065 inline detail::expression<
3066 detail::function
3067 , detail::abs_funct<Backend>
3068 , number<Backend, et_on> >
abs(const number<Backend,et_on> & arg)3069 abs(const number<Backend, et_on>& arg)
3070 {
3071 return detail::expression<
3072 detail::function
3073 , detail::abs_funct<Backend>
3074 , number<Backend, et_on>
3075 >(
3076 detail::abs_funct<Backend>()
3077 , arg
3078 );
3079 }
3080 template <class Backend>
3081 inline number<Backend, et_off>
abs(const number<Backend,et_off> & arg)3082 abs(const number<Backend, et_off>& arg)
3083 {
3084 number<Backend, et_off> result;
3085 using default_ops::eval_abs;
3086 eval_abs(result.backend(), arg.backend());
3087 return BOOST_MP_MOVE(result);
3088 }
3089
UNARY_OP_FUNCTOR(fabs,number_kind_floating_point)3090 UNARY_OP_FUNCTOR(fabs, number_kind_floating_point)
3091 UNARY_OP_FUNCTOR(sqrt, number_kind_floating_point)
3092 UNARY_OP_FUNCTOR(floor, number_kind_floating_point)
3093 UNARY_OP_FUNCTOR(ceil, number_kind_floating_point)
3094 UNARY_OP_FUNCTOR(trunc, number_kind_floating_point)
3095 UNARY_OP_FUNCTOR(round, number_kind_floating_point)
3096 UNARY_OP_FUNCTOR(exp, number_kind_floating_point)
3097 UNARY_OP_FUNCTOR(exp2, number_kind_floating_point)
3098 UNARY_OP_FUNCTOR(log, number_kind_floating_point)
3099 UNARY_OP_FUNCTOR(log10, number_kind_floating_point)
3100 UNARY_OP_FUNCTOR(cos, number_kind_floating_point)
3101 UNARY_OP_FUNCTOR(sin, number_kind_floating_point)
3102 UNARY_OP_FUNCTOR(tan, number_kind_floating_point)
3103 UNARY_OP_FUNCTOR(asin, number_kind_floating_point)
3104 UNARY_OP_FUNCTOR(acos, number_kind_floating_point)
3105 UNARY_OP_FUNCTOR(atan, number_kind_floating_point)
3106 UNARY_OP_FUNCTOR(cosh, number_kind_floating_point)
3107 UNARY_OP_FUNCTOR(sinh, number_kind_floating_point)
3108 UNARY_OP_FUNCTOR(tanh, number_kind_floating_point)
3109 UNARY_OP_FUNCTOR(log2, number_kind_floating_point)
3110 UNARY_OP_FUNCTOR(nearbyint, number_kind_floating_point)
3111 UNARY_OP_FUNCTOR(rint, number_kind_floating_point)
3112
3113 HETERO_BINARY_OP_FUNCTOR(ldexp, short, number_kind_floating_point)
3114 //HETERO_BINARY_OP_FUNCTOR(frexp, short*, number_kind_floating_point)
3115 HETERO_BINARY_OP_FUNCTOR_B(ldexp, int, number_kind_floating_point)
3116 //HETERO_BINARY_OP_FUNCTOR_B(frexp, int*, number_kind_floating_point)
3117 HETERO_BINARY_OP_FUNCTOR_B(ldexp, long, number_kind_floating_point)
3118 //HETERO_BINARY_OP_FUNCTOR_B(frexp, long*, number_kind_floating_point)
3119 HETERO_BINARY_OP_FUNCTOR_B(ldexp, boost::long_long_type, number_kind_floating_point)
3120 //HETERO_BINARY_OP_FUNCTOR_B(frexp, boost::long_long_type*, number_kind_floating_point)
3121 BINARY_OP_FUNCTOR(pow, number_kind_floating_point)
3122 BINARY_OP_FUNCTOR(fmod, number_kind_floating_point)
3123 BINARY_OP_FUNCTOR(fmax, number_kind_floating_point)
3124 BINARY_OP_FUNCTOR(fmin, number_kind_floating_point)
3125 BINARY_OP_FUNCTOR(atan2, number_kind_floating_point)
3126 BINARY_OP_FUNCTOR(fdim, number_kind_floating_point)
3127 BINARY_OP_FUNCTOR(hypot, number_kind_floating_point)
3128 BINARY_OP_FUNCTOR(remainder, number_kind_floating_point)
3129
3130 UNARY_OP_FUNCTOR(logb, number_kind_floating_point)
3131 HETERO_BINARY_OP_FUNCTOR(scalbn, short, number_kind_floating_point)
3132 HETERO_BINARY_OP_FUNCTOR(scalbln, short, number_kind_floating_point)
3133 HETERO_BINARY_OP_FUNCTOR_B(scalbn, int, number_kind_floating_point)
3134 HETERO_BINARY_OP_FUNCTOR_B(scalbln, int, number_kind_floating_point)
3135 HETERO_BINARY_OP_FUNCTOR_B(scalbn, long, number_kind_floating_point)
3136 HETERO_BINARY_OP_FUNCTOR_B(scalbln, long, number_kind_floating_point)
3137 HETERO_BINARY_OP_FUNCTOR_B(scalbn, boost::long_long_type, number_kind_floating_point)
3138 HETERO_BINARY_OP_FUNCTOR_B(scalbln, boost::long_long_type, number_kind_floating_point)
3139
3140 //
3141 // Integer functions:
3142 //
3143 BINARY_OP_FUNCTOR(gcd, number_kind_integer)
3144 BINARY_OP_FUNCTOR(lcm, number_kind_integer)
3145 HETERO_BINARY_OP_FUNCTOR_B(pow, unsigned, number_kind_integer)
3146
3147 #undef BINARY_OP_FUNCTOR
3148 #undef UNARY_OP_FUNCTOR
3149
3150 //
3151 // ilogb:
3152 //
3153 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
3154 inline typename enable_if_c<number_category<Backend>::value == number_kind_floating_point, typename Backend::exponent_type>::type
3155 ilogb(const multiprecision::number<Backend, ExpressionTemplates>& val)
3156 {
3157 using default_ops::eval_ilogb;
3158 return eval_ilogb(val.backend());
3159 }
3160
3161 template <class tag, class A1, class A2, class A3, class A4>
3162 inline typename enable_if_c<number_category<detail::expression<tag, A1, A2, A3, A4> >::value == number_kind_floating_point, typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type::backend_type::exponent_type>::type
ilogb(const detail::expression<tag,A1,A2,A3,A4> & val)3163 ilogb(const detail::expression<tag, A1, A2, A3, A4>& val)
3164 {
3165 using default_ops::eval_ilogb;
3166 typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type arg(val);
3167 return eval_ilogb(arg.backend());
3168 }
3169
3170 } //namespace multiprecision
3171
3172 namespace math{
3173 //
3174 // Overload of Boost.Math functions that find the wrong overload when used with number:
3175 //
3176 namespace detail{
3177 template <class T> T sinc_pi_imp(T);
3178 template <class T> T sinhc_pi_imp(T);
3179 }
3180 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
sinc_pi(const multiprecision::number<Backend,ExpressionTemplates> & x)3181 inline multiprecision::number<Backend, ExpressionTemplates> sinc_pi(const multiprecision::number<Backend, ExpressionTemplates>& x)
3182 {
3183 return BOOST_MP_MOVE(detail::sinc_pi_imp(x));
3184 }
3185
3186 template <class Backend, multiprecision::expression_template_option ExpressionTemplates, class Policy>
sinc_pi(const multiprecision::number<Backend,ExpressionTemplates> & x,const Policy &)3187 inline multiprecision::number<Backend, ExpressionTemplates> sinc_pi(const multiprecision::number<Backend, ExpressionTemplates>& x, const Policy&)
3188 {
3189 return BOOST_MP_MOVE(detail::sinc_pi_imp(x));
3190 }
3191
3192 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
sinhc_pi(const multiprecision::number<Backend,ExpressionTemplates> & x)3193 inline multiprecision::number<Backend, ExpressionTemplates> sinhc_pi(const multiprecision::number<Backend, ExpressionTemplates>& x)
3194 {
3195 return BOOST_MP_MOVE(detail::sinhc_pi_imp(x));
3196 }
3197
3198 template <class Backend, multiprecision::expression_template_option ExpressionTemplates, class Policy>
sinhc_pi(const multiprecision::number<Backend,ExpressionTemplates> & x,const Policy &)3199 inline multiprecision::number<Backend, ExpressionTemplates> sinhc_pi(const multiprecision::number<Backend, ExpressionTemplates>& x, const Policy&)
3200 {
3201 return BOOST_MP_MOVE(boost::math::sinhc_pi(x));
3202 }
3203
3204 using boost::multiprecision::gcd;
3205 using boost::multiprecision::lcm;
3206
3207 #ifdef BOOST_MSVC
3208 #pragma warning(pop)
3209 #endif
3210 } // namespace math
3211
3212 namespace integer {
3213
3214 using boost::multiprecision::gcd;
3215 using boost::multiprecision::lcm;
3216
3217 }
3218
3219 } // namespace boost
3220
3221 //
3222 // This has to come last of all:
3223 //
3224 #include <boost/multiprecision/detail/no_et_ops.hpp>
3225 #include <boost/multiprecision/detail/et_ops.hpp>
3226 //
3227 // min/max overloads:
3228 //
3229 #include <boost/multiprecision/detail/min_max.hpp>
3230
3231 #endif
3232
3233