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    T n;
987    eval_divide(result, a, b);
988    if(eval_get_sign(result) < 0)
989       eval_ceil(n, result);
990    else
991       eval_floor(n, result);
992    eval_multiply(n, b);
993    eval_subtract(result, a, n);
994 }
995 template<class T, class A>
eval_fmod(T & result,const T & x,const A & a)996 inline typename enable_if<is_arithmetic<A>, void>::type eval_fmod(T& result, const T& x, const A& a)
997 {
998    typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
999    typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1000    cast_type c;
1001    c = a;
1002    eval_fmod(result, x, c);
1003 }
1004 
1005 template<class T, class A>
eval_fmod(T & result,const A & x,const T & a)1006 inline typename enable_if<is_arithmetic<A>, void>::type eval_fmod(T& result, const A& x, const T& a)
1007 {
1008    typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
1009    typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1010    cast_type c;
1011    c = x;
1012    eval_fmod(result, c, a);
1013 }
1014 
1015 template <class T>
1016 void eval_round(T& result, const T& a);
1017 
1018 template <class T>
eval_remquo(T & result,const T & a,const T & b,int * pi)1019 inline void eval_remquo(T& result, const T& a, const T& b, int* pi)
1020 {
1021    BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The remquo function is only valid for floating point types.");
1022    if((&result == &a) || (&result == &b))
1023    {
1024       T temp;
1025       eval_remquo(temp, a, b, pi);
1026       result = temp;
1027       return;
1028    }
1029    T n;
1030    eval_divide(result, a, b);
1031    eval_round(n, result);
1032    eval_convert_to(pi, n);
1033    eval_multiply(n, b);
1034    eval_subtract(result, a, n);
1035 }
1036 template<class T, class A>
eval_remquo(T & result,const T & x,const A & a,int * pi)1037 inline typename enable_if<is_arithmetic<A>, void>::type eval_remquo(T& result, const T& x, const A& a, int* pi)
1038 {
1039    typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
1040    typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1041    cast_type c;
1042    c = a;
1043    eval_remquo(result, x, c, pi);
1044 }
1045 template<class T, class A>
eval_remquo(T & result,const A & x,const T & a,int * pi)1046 inline typename enable_if<is_arithmetic<A>, void>::type eval_remquo(T& result, const A& x, const T& a, int* pi)
1047 {
1048    typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
1049    typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
1050    cast_type c;
1051    c = x;
1052    eval_remquo(result, c, a, pi);
1053 }
1054 template <class T, class U, class V>
eval_remainder(T & result,const U & a,const V & b)1055 inline void eval_remainder(T& result, const U& a, const V& b)
1056 {
1057    int i;
1058    eval_remquo(result, a, b, &i);
1059 }
1060 
1061 template <class B>
1062 bool eval_gt(const B& a, const B& b);
1063 template <class T, class U>
1064 bool eval_gt(const T& a, const U& b);
1065 template <class B>
1066 bool eval_lt(const B& a, const B& b);
1067 template <class T, class U>
1068 bool eval_lt(const T& a, const U& b);
1069 
1070 template<class T>
eval_fdim(T & result,const T & a,const T & b)1071 inline void eval_fdim(T& result, const T& a, const T& b)
1072 {
1073    typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1074    static const ui_type zero = 0u;
1075    switch(eval_fpclassify(b))
1076    {
1077    case FP_NAN:
1078    case FP_INFINITE:
1079       result = zero;
1080       return;
1081    }
1082    switch(eval_fpclassify(a))
1083    {
1084    case FP_NAN:
1085       result = zero;
1086       return;
1087    case FP_INFINITE:
1088       result = a;
1089       return;
1090    }
1091    if(eval_gt(a, b))
1092    {
1093       eval_subtract(result, a, b);
1094    }
1095    else
1096       result = zero;
1097 }
1098 
1099 template<class T, class A>
eval_fdim(T & result,const T & a,const A & b)1100 inline typename boost::enable_if_c<boost::is_arithmetic<A>::value>::type eval_fdim(T& result, const T& a, const A& b)
1101 {
1102    typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1103    typedef typename boost::multiprecision::detail::canonical<A, T>::type arithmetic_type;
1104    static const ui_type zero = 0u;
1105    arithmetic_type canonical_b = b;
1106    switch((::boost::math::fpclassify)(b))
1107    {
1108    case FP_NAN:
1109    case FP_INFINITE:
1110       result = zero;
1111       return;
1112    }
1113    switch(eval_fpclassify(a))
1114    {
1115    case FP_NAN:
1116       result = zero;
1117       return;
1118    case FP_INFINITE:
1119       result = a;
1120       return;
1121    }
1122    if(eval_gt(a, canonical_b))
1123    {
1124       eval_subtract(result, a, canonical_b);
1125    }
1126    else
1127       result = zero;
1128 }
1129 
1130 template<class T, class A>
eval_fdim(T & result,const A & a,const T & b)1131 inline typename boost::enable_if_c<boost::is_arithmetic<A>::value>::type eval_fdim(T& result, const A& a, const T& b)
1132 {
1133    typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1134    typedef typename boost::multiprecision::detail::canonical<A, T>::type arithmetic_type;
1135    static const ui_type zero = 0u;
1136    arithmetic_type canonical_a = a;
1137    switch(eval_fpclassify(b))
1138    {
1139    case FP_NAN:
1140    case FP_INFINITE:
1141       result = zero;
1142       return;
1143    }
1144    switch((::boost::math::fpclassify)(a))
1145    {
1146    case FP_NAN:
1147       result = zero;
1148       return;
1149    case FP_INFINITE:
1150       result = std::numeric_limits<number<T> >::infinity().backend();
1151       return;
1152    }
1153    if(eval_gt(canonical_a, b))
1154    {
1155       eval_subtract(result, canonical_a, b);
1156    }
1157    else
1158       result = zero;
1159 }
1160 
1161 template <class T>
eval_trunc(T & result,const T & a)1162 inline void eval_trunc(T& result, const T& a)
1163 {
1164    BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The trunc function is only valid for floating point types.");
1165    int c = eval_fpclassify(a);
1166    if(c == (int)FP_NAN || c == (int)FP_INFINITE)
1167    {
1168       result = boost::math::policies::raise_rounding_error("boost::multiprecision::trunc<%1%>(%1%)", 0, number<T>(a), number<T>(a), boost::math::policies::policy<>()).backend();
1169       return;
1170    }
1171    if(eval_get_sign(a) < 0)
1172       eval_ceil(result, a);
1173    else
1174       eval_floor(result, a);
1175 }
1176 
1177 template <class T>
eval_modf(T & result,T const & arg,T * pipart)1178 inline void eval_modf(T& result, T const& arg, T* pipart)
1179 {
1180    typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1181    int c = eval_fpclassify(arg);
1182    if(c == (int)FP_NAN)
1183    {
1184       if(pipart)
1185          *pipart = arg;
1186       result = arg;
1187       return;
1188    }
1189    else if(c == (int)FP_INFINITE)
1190    {
1191       if(pipart)
1192          *pipart = arg;
1193       result = ui_type(0u);
1194       return;
1195    }
1196    if(pipart)
1197    {
1198       eval_trunc(*pipart, arg);
1199       eval_subtract(result, arg, *pipart);
1200    }
1201    else
1202    {
1203       T ipart;
1204       eval_trunc(ipart, arg);
1205       eval_subtract(result, arg, ipart);
1206    }
1207 }
1208 
1209 template <class T>
eval_round(T & result,const T & a)1210 inline void eval_round(T& result, const T& a)
1211 {
1212    BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The round function is only valid for floating point types.");
1213    typedef typename boost::multiprecision::detail::canonical<float, T>::type fp_type;
1214    int c = eval_fpclassify(a);
1215    if((c == (int)FP_NAN) || (c == (int)FP_INFINITE))
1216    {
1217       result = boost::math::policies::raise_rounding_error("boost::multiprecision::round<%1%>(%1%)", 0, number<T>(a), number<T>(a), boost::math::policies::policy<>()).backend();
1218       return;
1219    }
1220    if(eval_get_sign(a) < 0)
1221    {
1222       eval_subtract(result, a, fp_type(0.5f));
1223       eval_ceil(result, result);
1224    }
1225    else
1226    {
1227       eval_add(result, a, fp_type(0.5f));
1228       eval_floor(result, result);
1229    }
1230 }
1231 
1232 template <class B>
1233 void eval_lcm(B& result, const B& a, const B& b);
1234 template <class B>
1235 void eval_gcd(B& result, const B& a, const B& b);
1236 
1237 template <class T, class Arithmetic>
eval_gcd(T & result,const T & a,const Arithmetic & b)1238 inline typename enable_if<is_integral<Arithmetic> >::type eval_gcd(T& result, const T& a, const Arithmetic& b)
1239 {
1240    typedef typename boost::multiprecision::detail::canonical<Arithmetic, T>::type si_type;
1241    using default_ops::eval_gcd;
1242    T t;
1243    t = static_cast<si_type>(b);
1244    eval_gcd(result, a, t);
1245 }
1246 template <class T, class Arithmetic>
eval_gcd(T & result,const Arithmetic & a,const T & b)1247 inline typename enable_if<is_integral<Arithmetic> >::type eval_gcd(T& result, const Arithmetic& a, const T& b)
1248 {
1249    eval_gcd(result, b, a);
1250 }
1251 template <class T, class Arithmetic>
eval_lcm(T & result,const T & a,const Arithmetic & b)1252 inline typename enable_if<is_integral<Arithmetic> >::type eval_lcm(T& result, const T& a, const Arithmetic& b)
1253 {
1254    typedef typename boost::multiprecision::detail::canonical<Arithmetic, T>::type si_type;
1255    using default_ops::eval_lcm;
1256    T t;
1257    t = static_cast<si_type>(b);
1258    eval_lcm(result, a, t);
1259 }
1260 template <class T, class Arithmetic>
eval_lcm(T & result,const Arithmetic & a,const T & b)1261 inline typename enable_if<is_integral<Arithmetic> >::type eval_lcm(T& result, const Arithmetic& a, const T& b)
1262 {
1263    eval_lcm(result, b, a);
1264 }
1265 
1266 template <class T>
eval_lsb(const T & val)1267 inline unsigned eval_lsb(const T& val)
1268 {
1269    typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1270    int c = eval_get_sign(val);
1271    if(c == 0)
1272    {
1273       BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
1274    }
1275    if(c < 0)
1276    {
1277       BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
1278    }
1279    unsigned result = 0;
1280    T mask, t;
1281    mask = ui_type(1);
1282    do
1283    {
1284       eval_bitwise_and(t, mask, val);
1285       ++result;
1286       eval_left_shift(mask, 1);
1287    }
1288    while(eval_is_zero(t));
1289 
1290    return --result;
1291 }
1292 
1293 template <class T>
eval_msb(const T & val)1294 inline int eval_msb(const T& val)
1295 {
1296    int c = eval_get_sign(val);
1297    if(c == 0)
1298    {
1299       BOOST_THROW_EXCEPTION(std::range_error("No bits were set in the operand."));
1300    }
1301    if(c < 0)
1302    {
1303       BOOST_THROW_EXCEPTION(std::range_error("Testing individual bits in negative values is not supported - results are undefined."));
1304    }
1305    //
1306    // This implementation is really really rubbish - it does
1307    // a linear scan for the most-significant-bit.  We should really
1308    // do a binary search, but as none of our backends actually needs
1309    // this implementation, we'll leave it for now.  In fact for most
1310    // backends it's likely that there will always be a more efficient
1311    // native implementation possible.
1312    //
1313    unsigned result = 0;
1314    T t(val);
1315    while(!eval_is_zero(t))
1316    {
1317       eval_right_shift(t, 1);
1318       ++result;
1319    }
1320    return --result;
1321 }
1322 
1323 template <class T>
eval_bit_test(const T & val,unsigned index)1324 inline bool eval_bit_test(const T& val, unsigned index)
1325 {
1326    typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1327    T mask, t;
1328    mask = ui_type(1);
1329    eval_left_shift(mask, index);
1330    eval_bitwise_and(t, mask, val);
1331    return !eval_is_zero(t);
1332 }
1333 
1334 template <class T>
eval_bit_set(T & val,unsigned index)1335 inline void eval_bit_set(T& val, unsigned index)
1336 {
1337    typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1338    T mask;
1339    mask = ui_type(1);
1340    eval_left_shift(mask, index);
1341    eval_bitwise_or(val, mask);
1342 }
1343 
1344 template <class T>
eval_bit_flip(T & val,unsigned index)1345 inline void eval_bit_flip(T& val, unsigned index)
1346 {
1347    typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1348    T mask;
1349    mask = ui_type(1);
1350    eval_left_shift(mask, index);
1351    eval_bitwise_xor(val, mask);
1352 }
1353 
1354 template <class T>
eval_bit_unset(T & val,unsigned index)1355 inline void eval_bit_unset(T& val, unsigned index)
1356 {
1357    typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
1358    T mask, t;
1359    mask = ui_type(1);
1360    eval_left_shift(mask, index);
1361    eval_bitwise_and(t, mask, val);
1362    if(!eval_is_zero(t))
1363       eval_bitwise_xor(val, mask);
1364 }
1365 
1366 template <class B>
eval_integer_sqrt(B & s,B & r,const B & x)1367 void eval_integer_sqrt(B& s, B& r, const B& x)
1368 {
1369    //
1370    // This is slow bit-by-bit integer square root, see for example
1371    // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
1372    // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
1373    // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
1374    // at some point.
1375    //
1376    typedef typename boost::multiprecision::detail::canonical<unsigned char, B>::type ui_type;
1377 
1378    s = ui_type(0u);
1379    if(eval_get_sign(x) == 0)
1380    {
1381       r = ui_type(0u);
1382       return;
1383    }
1384    int g = eval_msb(x);
1385    if(g == 0)
1386    {
1387       r = ui_type(1);
1388       return;
1389    }
1390 
1391    B t;
1392    r = x;
1393    g /= 2;
1394    int org_g = g;
1395    eval_bit_set(s, g);
1396    eval_bit_set(t, 2 * g);
1397    eval_subtract(r, x, t);
1398    --g;
1399    if(eval_get_sign(r) == 0)
1400       return;
1401    int msbr = eval_msb(r);
1402    do
1403    {
1404       if(msbr >= org_g + g + 1)
1405       {
1406          t = s;
1407          eval_left_shift(t, g + 1);
1408          eval_bit_set(t, 2 * g);
1409          if(t.compare(r) <= 0)
1410          {
1411             eval_bit_set(s, g);
1412             eval_subtract(r, t);
1413             if(eval_get_sign(r) == 0)
1414                return;
1415             msbr = eval_msb(r);
1416          }
1417       }
1418       --g;
1419    }
1420    while(g >= 0);
1421 }
1422 
1423 //
1424 // These have to implemented by the backend, declared here so that our macro generated code compiles OK.
1425 //
1426 template <class T>
1427 typename enable_if_c<sizeof(T) == 0>::type eval_floor();
1428 template <class T>
1429 typename enable_if_c<sizeof(T) == 0>::type eval_ceil();
1430 template <class T>
1431 typename enable_if_c<sizeof(T) == 0>::type eval_trunc();
1432 template <class T>
1433 typename enable_if_c<sizeof(T) == 0>::type eval_sqrt();
1434 template <class T>
1435 typename enable_if_c<sizeof(T) == 0>::type eval_ldexp();
1436 template <class T>
1437 typename enable_if_c<sizeof(T) == 0>::type eval_frexp();
1438 
1439 //
1440 // eval_logb and eval_scalbn simply assume base 2 and forward to
1441 // eval_ldexp and eval_frexp:
1442 //
1443 template <class B>
eval_ilogb(const B & val)1444 inline typename B::exponent_type eval_ilogb(const B& val)
1445 {
1446    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");
1447    typename B::exponent_type e;
1448    switch(eval_fpclassify(val))
1449    {
1450    case FP_NAN:
1451       return (std::numeric_limits<typename B::exponent_type>::min)();
1452    case FP_INFINITE:
1453       return (std::numeric_limits<typename B::exponent_type>::max)();
1454    case FP_ZERO:
1455       return (std::numeric_limits<typename B::exponent_type>::min)();
1456    }
1457    B result;
1458    eval_frexp(result, val, &e);
1459    return e - 1;
1460 }
1461 template <class B>
eval_logb(B & result,const B & val)1462 inline void eval_logb(B& result, const B& val)
1463 {
1464    typedef typename boost::mpl::if_c<boost::is_same<boost::intmax_t, long>::value, boost::long_long_type, boost::intmax_t>::type max_t;
1465    result = static_cast<max_t>(eval_ilogb(val));
1466 }
1467 template <class B, class A>
eval_scalbn(B & result,const B & val,A e)1468 inline void eval_scalbn(B& result, const B& val, A e)
1469 {
1470    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");
1471    eval_ldexp(result, val, static_cast<typename B::exponent_type>(e));
1472 }
1473 template <class B, class A>
eval_scalbln(B & result,const B & val,A e)1474 inline void eval_scalbln(B& result, const B& val, A e)
1475 {
1476    eval_scalbn(result, val, e);
1477 }
1478 
1479 template <class T>
is_arg_nan(const T & val,mpl::true_ const &,const mpl::false_ &)1480 inline bool is_arg_nan(const T& val, mpl::true_ const&, const mpl::false_&)
1481 {
1482    return eval_fpclassify(val) == FP_NAN;
1483 }
1484 template <class T>
is_arg_nan(const T & val,mpl::false_ const &,const mpl::true_ &)1485 inline bool is_arg_nan(const T& val, mpl::false_ const&, const mpl::true_&)
1486 {
1487    return (boost::math::isnan)(val);
1488 }
1489 template <class T>
is_arg_nan(const T &,mpl::false_ const &,const mpl::false_ &)1490 inline bool is_arg_nan(const T&, mpl::false_ const&, const mpl::false_&)
1491 {
1492    return false;
1493 }
1494 
1495 template <class T>
is_arg_nan(const T & val)1496 inline bool is_arg_nan(const T& val)
1497 {
1498    return is_arg_nan(val, mpl::bool_<boost::multiprecision::detail::is_backend<T>::value>(), is_floating_point<T>());
1499 }
1500 
1501 template <class T, class U, class V>
eval_fmax(T & result,const U & a,const V & b)1502 inline void eval_fmax(T& result, const U& a, const V& b)
1503 {
1504    if(is_arg_nan(a))
1505       result = number<T>::canonical_value(b);
1506    else if(is_arg_nan(b))
1507       result = number<T>::canonical_value(a);
1508    else if(eval_lt(number<T>::canonical_value(a), number<T>::canonical_value(b)))
1509       result = number<T>::canonical_value(b);
1510    else
1511       result = number<T>::canonical_value(a);
1512 }
1513 template <class T, class U, class V>
eval_fmin(T & result,const U & a,const V & b)1514 inline void eval_fmin(T& result, const U& a, const V& b)
1515 {
1516    if(is_arg_nan(a))
1517       result = number<T>::canonical_value(b);
1518    else if(is_arg_nan(b))
1519       result = number<T>::canonical_value(a);
1520    else if(eval_lt(number<T>::canonical_value(a), number<T>::canonical_value(b)))
1521       result = number<T>::canonical_value(a);
1522    else
1523       result = number<T>::canonical_value(b);
1524 }
1525 
1526 template <class R, class T, class U>
eval_hypot(R & result,const T & a,const U & b)1527 inline void eval_hypot(R& result, const T& a, const U& b)
1528 {
1529    //
1530    // Normalize x and y, so that both are positive and x >= y:
1531    //
1532    R x, y;
1533    x = number<R>::canonical_value(a);
1534    y = number<R>::canonical_value(b);
1535    if(eval_get_sign(x) < 0)
1536       x.negate();
1537    if(eval_get_sign(y) < 0)
1538       y.negate();
1539 
1540    // Special case, see C99 Annex F.
1541    // The order of the if's is important: do not change!
1542    int c1 = eval_fpclassify(x);
1543    int c2 = eval_fpclassify(y);
1544 
1545    if(c1 == FP_ZERO)
1546    {
1547       result = y;
1548       return;
1549    }
1550    if(c2 == FP_ZERO)
1551    {
1552       result = x;
1553       return;
1554    }
1555    if(c1 == FP_INFINITE)
1556    {
1557       result = x;
1558       return;
1559    }
1560    if((c2 == FP_INFINITE) || (c2 == FP_NAN))
1561    {
1562       result = y;
1563       return;
1564    }
1565    if(c1 == FP_NAN)
1566    {
1567       result = x;
1568       return;
1569    }
1570 
1571    if(eval_gt(y, x))
1572       x.swap(y);
1573 
1574    eval_multiply(result, x, std::numeric_limits<number<R> >::epsilon().backend());
1575 
1576    if(eval_gt(result, y))
1577    {
1578       result = x;
1579       return;
1580    }
1581 
1582    R rat;
1583    eval_divide(rat, y, x);
1584    eval_multiply(result, rat, rat);
1585    eval_increment(result);
1586    eval_sqrt(rat, result);
1587    eval_multiply(result, rat, x);
1588 }
1589 
1590 template <class R, class T>
eval_nearbyint(R & result,const T & a)1591 inline void eval_nearbyint(R& result, const T& a)
1592 {
1593    eval_round(result, a);
1594 }
1595 template <class R, class T>
eval_rint(R & result,const T & a)1596 inline void eval_rint(R& result, const T& a)
1597 {
1598    eval_nearbyint(result, a);
1599 }
1600 
1601 //
1602 // These functions are implemented in separate files, but expanded inline here,
1603 // DO NOT CHANGE THE ORDER OF THESE INCLUDES:
1604 //
1605 #include <boost/multiprecision/detail/functions/constants.hpp>
1606 #include <boost/multiprecision/detail/functions/pow.hpp>
1607 #include <boost/multiprecision/detail/functions/trig.hpp>
1608 
1609 }
1610 
1611 //
1612 // Default versions of floating point classification routines:
1613 //
1614 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1615 inline int fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1616 {
1617    using multiprecision::default_ops::eval_fpclassify;
1618    return eval_fpclassify(arg.backend());
1619 }
1620 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)1621 inline int fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1622 {
1623    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1624    return fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1625 }
1626 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1627 inline bool isfinite BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1628 {
1629    int v = fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(arg);
1630    return (v != (int)FP_INFINITE) && (v != (int)FP_NAN);
1631 }
1632 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)1633 inline bool isfinite BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1634 {
1635    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1636    return isfinite BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1637 }
1638 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1639 inline bool isnan BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1640 {
1641    return fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(arg) == (int)FP_NAN;
1642 }
1643 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)1644 inline bool isnan BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1645 {
1646    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1647    return isnan BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1648 }
1649 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1650 inline bool isinf BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1651 {
1652    return fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(arg) == (int)FP_INFINITE;
1653 }
1654 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)1655 inline bool isinf BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1656 {
1657    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1658    return isinf BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1659 }
1660 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1661 inline bool isnormal BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1662 {
1663    return fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(arg) == (int)FP_NORMAL;
1664 }
1665 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)1666 inline bool isnormal BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1667 {
1668    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1669    return isnormal BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1670 }
1671 
1672 // Default versions of sign manipulation functions, if individual backends can do better than this
1673 // (for example with signed zero), then they should overload these functions further:
1674 
1675 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1676 inline int sign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1677 {
1678    return arg.sign();
1679 }
1680 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)1681 inline int sign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1682 {
1683    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1684    return sign BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1685 }
1686 
1687 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1688 inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1689 {
1690    return arg.sign() < 0;
1691 }
1692 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)1693 inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1694 {
1695    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1696    return signbit BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1697 }
1698 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1699 inline multiprecision::number<Backend, ExpressionTemplates> changesign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1700 {
1701    return -arg;
1702 }
1703 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)1704 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)
1705 {
1706    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1707    return changesign BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(arg));
1708 }
1709 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & a,const multiprecision::number<Backend,ExpressionTemplates> & b)1710 inline multiprecision::number<Backend, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& a, const multiprecision::number<Backend, ExpressionTemplates>& b)
1711 {
1712    return (boost::multiprecision::signbit)(a) != (boost::multiprecision::signbit)(b) ? (boost::multiprecision::changesign)(a) : a;
1713 }
1714 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)1715 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)
1716 {
1717    return copysign BOOST_PREVENT_MACRO_SUBSTITUTION(a, multiprecision::number<Backend, ExpressionTemplates>(b));
1718 }
1719 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)1720 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)
1721 {
1722    return copysign BOOST_PREVENT_MACRO_SUBSTITUTION(multiprecision::number<Backend, ExpressionTemplates>(a), b);
1723 }
1724 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)1725 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)
1726 {
1727    typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1728    return copysign BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(a), value_type(b));
1729 }
1730 
1731 } // namespace multiprecision
1732 
1733 namespace math {
1734 
1735    //
1736    // Import Math functions here, so they can be found by Boost.Math:
1737    //
1738    using boost::multiprecision::signbit;
1739    using boost::multiprecision::sign;
1740    using boost::multiprecision::copysign;
1741    using boost::multiprecision::changesign;
1742    using boost::multiprecision::fpclassify;
1743    using boost::multiprecision::isinf;
1744    using boost::multiprecision::isnan;
1745    using boost::multiprecision::isnormal;
1746    using boost::multiprecision::isfinite;
1747 
1748 }
1749 
1750 namespace multiprecision{
1751 
1752    typedef ::boost::math::policies::policy<
1753       ::boost::math::policies::domain_error< ::boost::math::policies::errno_on_error>,
1754       ::boost::math::policies::pole_error< ::boost::math::policies::errno_on_error>,
1755       ::boost::math::policies::overflow_error< ::boost::math::policies::errno_on_error>,
1756       ::boost::math::policies::evaluation_error< ::boost::math::policies::errno_on_error>,
1757       ::boost::math::policies::rounding_error< ::boost::math::policies::errno_on_error>
1758    > c99_error_policy;
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> asinh BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1762    {
1763       return boost::math::asinh(arg, c99_error_policy());
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 asinh 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 asinh(value_type(arg));
1770    }
1771    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1772    inline multiprecision::number<Backend, ExpressionTemplates> acosh BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1773    {
1774       return boost::math::acosh(arg, c99_error_policy());
1775    }
1776    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)1777    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)
1778    {
1779       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1780       return acosh(value_type(arg));
1781    }
1782    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1783    inline multiprecision::number<Backend, ExpressionTemplates> atanh BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1784    {
1785       return boost::math::atanh(arg, c99_error_policy());
1786    }
1787    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)1788    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)
1789    {
1790       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1791       return atanh(value_type(arg));
1792    }
1793    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1794    inline multiprecision::number<Backend, ExpressionTemplates> cbrt BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1795    {
1796       return boost::math::cbrt(arg, c99_error_policy());
1797    }
1798    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)1799    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)
1800    {
1801       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1802       return cbrt(value_type(arg));
1803    }
1804    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1805    inline multiprecision::number<Backend, ExpressionTemplates> erf BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1806    {
1807       return boost::math::erf(arg, c99_error_policy());
1808    }
1809    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)1810    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)
1811    {
1812       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1813       return erf(value_type(arg));
1814    }
1815    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1816    inline multiprecision::number<Backend, ExpressionTemplates> erfc BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1817    {
1818       return boost::math::erfc(arg, c99_error_policy());
1819    }
1820    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)1821    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)
1822    {
1823       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1824       return erfc(value_type(arg));
1825    }
1826    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1827    inline multiprecision::number<Backend, ExpressionTemplates> expm1 BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1828    {
1829       return boost::math::expm1(arg, c99_error_policy());
1830    }
1831    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)1832    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)
1833    {
1834       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1835       return expm1(value_type(arg));
1836    }
1837    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1838    inline multiprecision::number<Backend, ExpressionTemplates> lgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1839    {
1840       return boost::math::lgamma(arg, c99_error_policy());
1841    }
1842    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)1843    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)
1844    {
1845       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1846       return lgamma(value_type(arg));
1847    }
1848    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1849    inline multiprecision::number<Backend, ExpressionTemplates> tgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1850    {
1851       return boost::math::tgamma(arg, c99_error_policy());
1852    }
1853    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)1854    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)
1855    {
1856       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1857       return tgamma(value_type(arg));
1858    }
1859 
1860    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1861    inline long lrint BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1862    {
1863       return lround(arg);
1864    }
1865    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)1866    inline long lrint BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1867    {
1868       return lround(arg);
1869    }
1870 #ifndef BOOST_NO_LONG_LONG
1871    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1872    inline boost::long_long_type llrint BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1873    {
1874       return llround(arg);
1875    }
1876    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)1877    inline boost::long_long_type llrint BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::detail::expression<tag, A1, A2, A3, A4>& arg)
1878    {
1879       return llround(arg);
1880    }
1881 #endif
1882    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & arg)1883    inline multiprecision::number<Backend, ExpressionTemplates> log1p BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& arg)
1884    {
1885       return boost::math::log1p(arg, c99_error_policy());
1886    }
1887    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)1888    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)
1889    {
1890       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1891       return log1p(value_type(arg));
1892    }
1893 
1894    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & a,const multiprecision::number<Backend,ExpressionTemplates> & b)1895    inline multiprecision::number<Backend, ExpressionTemplates> nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& a, const multiprecision::number<Backend, ExpressionTemplates>& b)
1896    {
1897       return boost::math::nextafter(a, b, c99_error_policy());
1898    }
1899    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)1900    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)
1901    {
1902       return nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(a, multiprecision::number<Backend, ExpressionTemplates>(b));
1903    }
1904    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)1905    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)
1906    {
1907       return nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(multiprecision::number<Backend, ExpressionTemplates>(a), b);
1908    }
1909    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)1910    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)
1911    {
1912       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1913       return nextafter BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(a), value_type(b));
1914    }
1915    template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend,ExpressionTemplates> & a,const multiprecision::number<Backend,ExpressionTemplates> & b)1916    inline multiprecision::number<Backend, ExpressionTemplates> nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(const multiprecision::number<Backend, ExpressionTemplates>& a, const multiprecision::number<Backend, ExpressionTemplates>& b)
1917    {
1918       return boost::math::nextafter(a, b, c99_error_policy());
1919    }
1920    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)1921    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)
1922    {
1923       return nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(a, multiprecision::number<Backend, ExpressionTemplates>(b));
1924    }
1925    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)1926    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)
1927    {
1928       return nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(multiprecision::number<Backend, ExpressionTemplates>(a), b);
1929    }
1930    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)1931    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)
1932    {
1933       typedef typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type value_type;
1934       return nexttoward BOOST_PREVENT_MACRO_SUBSTITUTION(value_type(a), value_type(b));
1935    }
1936 
1937 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)1938 inline number<B1, ET1>& add(number<B1, ET1>& result, const number<B2, ET2>& a, const number<B3, ET3>& b)
1939 {
1940    BOOST_STATIC_ASSERT_MSG((is_convertible<B2, B1>::value), "No conversion to the target of a mixed precision addition exists");
1941    BOOST_STATIC_ASSERT_MSG((is_convertible<B3, B1>::value), "No conversion to the target of a mixed precision addition exists");
1942    using default_ops::eval_add;
1943    eval_add(result.backend(), a.backend(), b.backend());
1944    return result;
1945 }
1946 
1947 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)1948 inline number<B1, ET1>& subtract(number<B1, ET1>& result, const number<B2, ET2>& a, const number<B3, ET3>& b)
1949 {
1950    BOOST_STATIC_ASSERT_MSG((is_convertible<B2, B1>::value), "No conversion to the target of a mixed precision addition exists");
1951    BOOST_STATIC_ASSERT_MSG((is_convertible<B3, B1>::value), "No conversion to the target of a mixed precision addition exists");
1952    using default_ops::eval_subtract;
1953    eval_subtract(result.backend(), a.backend(), b.backend());
1954    return result;
1955 }
1956 
1957 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)1958 inline number<B1, ET1>& multiply(number<B1, ET1>& result, const number<B2, ET2>& a, const number<B3, ET3>& b)
1959 {
1960    BOOST_STATIC_ASSERT_MSG((is_convertible<B2, B1>::value), "No conversion to the target of a mixed precision addition exists");
1961    BOOST_STATIC_ASSERT_MSG((is_convertible<B3, B1>::value), "No conversion to the target of a mixed precision addition exists");
1962    using default_ops::eval_multiply;
1963    eval_multiply(result.backend(), a.backend(), b.backend());
1964    return result;
1965 }
1966 
1967 template <class B, expression_template_option ET, class I>
1968 inline typename enable_if_c<is_integral<I>::value, number<B, ET>&>::type
add(number<B,ET> & result,const I & a,const I & b)1969    add(number<B, ET>& result, const I& a, const I& b)
1970 {
1971    using default_ops::eval_add;
1972    typedef typename detail::canonical<I, B>::type canonical_type;
1973    eval_add(result.backend(), static_cast<canonical_type>(a), static_cast<canonical_type>(b));
1974    return result;
1975 }
1976 
1977 template <class B, expression_template_option ET, class I>
1978 inline typename enable_if_c<is_integral<I>::value, number<B, ET>&>::type
subtract(number<B,ET> & result,const I & a,const I & b)1979    subtract(number<B, ET>& result, const I& a, const I& b)
1980 {
1981    using default_ops::eval_subtract;
1982    typedef typename detail::canonical<I, B>::type canonical_type;
1983    eval_subtract(result.backend(), static_cast<canonical_type>(a), static_cast<canonical_type>(b));
1984    return result;
1985 }
1986 
1987 template <class B, expression_template_option ET, class I>
1988 inline typename enable_if_c<is_integral<I>::value, number<B, ET>&>::type
multiply(number<B,ET> & result,const I & a,const I & b)1989    multiply(number<B, ET>& result, const I& a, const I& b)
1990 {
1991    using default_ops::eval_multiply;
1992    typedef typename detail::canonical<I, B>::type canonical_type;
1993    eval_multiply(result.backend(), static_cast<canonical_type>(a), static_cast<canonical_type>(b));
1994    return result;
1995 }
1996 
1997 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)1998 inline typename detail::expression<tag, A1, A2, A3, A4>::result_type trunc(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
1999 {
2000    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2001    return BOOST_MP_MOVE(trunc(number_type(v), pol));
2002 }
2003 
2004 template <class Backend, expression_template_option ExpressionTemplates, class Policy>
trunc(const number<Backend,ExpressionTemplates> & v,const Policy &)2005 inline number<Backend, ExpressionTemplates> trunc(const number<Backend, ExpressionTemplates>& v, const Policy&)
2006 {
2007    using default_ops::eval_trunc;
2008    number<Backend, ExpressionTemplates> result;
2009    eval_trunc(result.backend(), v.backend());
2010    return BOOST_MP_MOVE(result);
2011 }
2012 
2013 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)2014 inline int itrunc(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2015 {
2016    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2017    number_type r = trunc(v, pol);
2018    if((r > (std::numeric_limits<int>::max)()) || r < (std::numeric_limits<int>::min)() || !(boost::math::isfinite)(v))
2019       return boost::math::policies::raise_rounding_error("boost::multiprecision::itrunc<%1%>(%1%)", 0, number_type(v), 0, pol);
2020    return r.template convert_to<int>();
2021 }
2022 template <class tag, class A1, class A2, class A3, class A4>
itrunc(const detail::expression<tag,A1,A2,A3,A4> & v)2023 inline int itrunc(const detail::expression<tag, A1, A2, A3, A4>& v)
2024 {
2025    return itrunc(v, boost::math::policies::policy<>());
2026 }
2027 template <class Backend, expression_template_option ExpressionTemplates, class Policy>
itrunc(const number<Backend,ExpressionTemplates> & v,const Policy & pol)2028 inline int itrunc(const number<Backend, ExpressionTemplates>& v, const Policy& pol)
2029 {
2030    number<Backend, ExpressionTemplates> r = trunc(v, pol);
2031    if((r > (std::numeric_limits<int>::max)()) || r < (std::numeric_limits<int>::min)() || !(boost::math::isfinite)(v))
2032       return boost::math::policies::raise_rounding_error("boost::multiprecision::itrunc<%1%>(%1%)", 0, v, 0, pol);
2033    return r.template convert_to<int>();
2034 }
2035 template <class Backend, expression_template_option ExpressionTemplates>
itrunc(const number<Backend,ExpressionTemplates> & v)2036 inline int itrunc(const number<Backend, ExpressionTemplates>& v)
2037 {
2038    return itrunc(v, boost::math::policies::policy<>());
2039 }
2040 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)2041 inline long ltrunc(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2042 {
2043    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2044    number_type r = trunc(v, pol);
2045    if((r > (std::numeric_limits<long>::max)()) || r < (std::numeric_limits<long>::min)() || !(boost::math::isfinite)(v))
2046       return boost::math::policies::raise_rounding_error("boost::multiprecision::ltrunc<%1%>(%1%)", 0, number_type(v), 0L, pol);
2047    return r.template convert_to<long>();
2048 }
2049 template <class tag, class A1, class A2, class A3, class A4>
ltrunc(const detail::expression<tag,A1,A2,A3,A4> & v)2050 inline long ltrunc(const detail::expression<tag, A1, A2, A3, A4>& v)
2051 {
2052    return ltrunc(v, boost::math::policies::policy<>());
2053 }
2054 template <class T, expression_template_option ExpressionTemplates, class Policy>
ltrunc(const number<T,ExpressionTemplates> & v,const Policy & pol)2055 inline long ltrunc(const number<T, ExpressionTemplates>& v, const Policy& pol)
2056 {
2057    number<T, ExpressionTemplates> r = trunc(v, pol);
2058    if((r > (std::numeric_limits<long>::max)()) || r < (std::numeric_limits<long>::min)() || !(boost::math::isfinite)(v))
2059       return boost::math::policies::raise_rounding_error("boost::multiprecision::ltrunc<%1%>(%1%)", 0, v, 0L, pol);
2060    return r.template convert_to<long>();
2061 }
2062 template <class T, expression_template_option ExpressionTemplates>
ltrunc(const number<T,ExpressionTemplates> & v)2063 inline long ltrunc(const number<T, ExpressionTemplates>& v)
2064 {
2065    return ltrunc(v, boost::math::policies::policy<>());
2066 }
2067 #ifndef BOOST_NO_LONG_LONG
2068 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)2069 inline boost::long_long_type lltrunc(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2070 {
2071    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2072    number_type r = trunc(v, pol);
2073    if((r > (std::numeric_limits<boost::long_long_type>::max)()) || r < (std::numeric_limits<boost::long_long_type>::min)() || !(boost::math::isfinite)(v))
2074       return boost::math::policies::raise_rounding_error("boost::multiprecision::lltrunc<%1%>(%1%)", 0, number_type(v), 0LL, pol);
2075    return r.template convert_to<boost::long_long_type>();
2076 }
2077 template <class tag, class A1, class A2, class A3, class A4>
lltrunc(const detail::expression<tag,A1,A2,A3,A4> & v)2078 inline boost::long_long_type lltrunc(const detail::expression<tag, A1, A2, A3, A4>& v)
2079 {
2080    return lltrunc(v, boost::math::policies::policy<>());
2081 }
2082 template <class T, expression_template_option ExpressionTemplates, class Policy>
lltrunc(const number<T,ExpressionTemplates> & v,const Policy & pol)2083 inline boost::long_long_type lltrunc(const number<T, ExpressionTemplates>& v, const Policy& pol)
2084 {
2085    number<T, ExpressionTemplates> r = trunc(v, pol);
2086    if((r > (std::numeric_limits<boost::long_long_type>::max)()) || r < (std::numeric_limits<boost::long_long_type>::min)() || !(boost::math::isfinite)(v))
2087       return boost::math::policies::raise_rounding_error("boost::multiprecision::lltrunc<%1%>(%1%)", 0, v, 0LL, pol);
2088    return r.template convert_to<boost::long_long_type>();
2089 }
2090 template <class T, expression_template_option ExpressionTemplates>
lltrunc(const number<T,ExpressionTemplates> & v)2091 inline boost::long_long_type lltrunc(const number<T, ExpressionTemplates>& v)
2092 {
2093    return lltrunc(v, boost::math::policies::policy<>());
2094 }
2095 #endif
2096 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)2097 inline typename detail::expression<tag, A1, A2, A3, A4>::result_type round(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2098 {
2099    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2100    return BOOST_MP_MOVE(round(static_cast<number_type>(v), pol));
2101 }
2102 template <class T, expression_template_option ExpressionTemplates, class Policy>
round(const number<T,ExpressionTemplates> & v,const Policy &)2103 inline number<T, ExpressionTemplates> round(const number<T, ExpressionTemplates>& v, const Policy&)
2104 {
2105    using default_ops::eval_round;
2106    number<T, ExpressionTemplates> result;
2107    eval_round(result.backend(), v.backend());
2108    return BOOST_MP_MOVE(result);
2109 }
2110 
2111 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)2112 inline int iround(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2113 {
2114    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2115    number_type r = round(v, pol);
2116    if((r > (std::numeric_limits<int>::max)()) || r < (std::numeric_limits<int>::min)() || !(boost::math::isfinite)(v))
2117       return boost::math::policies::raise_rounding_error("boost::multiprecision::iround<%1%>(%1%)", 0, number_type(v), 0, pol);
2118    return r.template convert_to<int>();
2119 }
2120 template <class tag, class A1, class A2, class A3, class A4>
iround(const detail::expression<tag,A1,A2,A3,A4> & v)2121 inline int iround(const detail::expression<tag, A1, A2, A3, A4>& v)
2122 {
2123    return iround(v, boost::math::policies::policy<>());
2124 }
2125 template <class T, expression_template_option ExpressionTemplates, class Policy>
iround(const number<T,ExpressionTemplates> & v,const Policy & pol)2126 inline int iround(const number<T, ExpressionTemplates>& v, const Policy& pol)
2127 {
2128    number<T, ExpressionTemplates> r = round(v, pol);
2129    if((r > (std::numeric_limits<int>::max)()) || r < (std::numeric_limits<int>::min)() || !(boost::math::isfinite)(v))
2130       return boost::math::policies::raise_rounding_error("boost::multiprecision::iround<%1%>(%1%)", 0, v, 0, pol);
2131    return r.template convert_to<int>();
2132 }
2133 template <class T, expression_template_option ExpressionTemplates>
iround(const number<T,ExpressionTemplates> & v)2134 inline int iround(const number<T, ExpressionTemplates>& v)
2135 {
2136    return iround(v, boost::math::policies::policy<>());
2137 }
2138 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)2139 inline long lround(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2140 {
2141    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2142    number_type r = round(v, pol);
2143    if((r > (std::numeric_limits<long>::max)()) || r < (std::numeric_limits<long>::min)() || !(boost::math::isfinite)(v))
2144       return boost::math::policies::raise_rounding_error("boost::multiprecision::lround<%1%>(%1%)", 0, number_type(v), 0L, pol);
2145    return r.template convert_to<long>();
2146 }
2147 template <class tag, class A1, class A2, class A3, class A4>
lround(const detail::expression<tag,A1,A2,A3,A4> & v)2148 inline long lround(const detail::expression<tag, A1, A2, A3, A4>& v)
2149 {
2150    return lround(v, boost::math::policies::policy<>());
2151 }
2152 template <class T, expression_template_option ExpressionTemplates, class Policy>
lround(const number<T,ExpressionTemplates> & v,const Policy & pol)2153 inline long lround(const number<T, ExpressionTemplates>& v, const Policy& pol)
2154 {
2155    number<T, ExpressionTemplates> r = round(v, pol);
2156    if((r > (std::numeric_limits<long>::max)()) || r < (std::numeric_limits<long>::min)() || !(boost::math::isfinite)(v))
2157       return boost::math::policies::raise_rounding_error("boost::multiprecision::lround<%1%>(%1%)", 0, v, 0L, pol);
2158    return r.template convert_to<long>();
2159 }
2160 template <class T, expression_template_option ExpressionTemplates>
lround(const number<T,ExpressionTemplates> & v)2161 inline long lround(const number<T, ExpressionTemplates>& v)
2162 {
2163    return lround(v, boost::math::policies::policy<>());
2164 }
2165 #ifndef BOOST_NO_LONG_LONG
2166 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)2167 inline boost::long_long_type llround(const detail::expression<tag, A1, A2, A3, A4>& v, const Policy& pol)
2168 {
2169    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2170    number_type r = round(v, pol);
2171    if((r > (std::numeric_limits<boost::long_long_type>::max)()) || r < (std::numeric_limits<boost::long_long_type>::min)() || !(boost::math::isfinite)(v))
2172       return boost::math::policies::raise_rounding_error("boost::multiprecision::iround<%1%>(%1%)", 0, number_type(v), 0LL, pol);
2173    return r.template convert_to<boost::long_long_type>();
2174 }
2175 template <class tag, class A1, class A2, class A3, class A4>
llround(const detail::expression<tag,A1,A2,A3,A4> & v)2176 inline boost::long_long_type llround(const detail::expression<tag, A1, A2, A3, A4>& v)
2177 {
2178    return llround(v, boost::math::policies::policy<>());
2179 }
2180 template <class T, expression_template_option ExpressionTemplates, class Policy>
llround(const number<T,ExpressionTemplates> & v,const Policy & pol)2181 inline boost::long_long_type llround(const number<T, ExpressionTemplates>& v, const Policy& pol)
2182 {
2183    number<T, ExpressionTemplates> r = round(v, pol);
2184    if((r > (std::numeric_limits<boost::long_long_type>::max)()) || r < (std::numeric_limits<boost::long_long_type>::min)() || !(boost::math::isfinite)(v))
2185       return boost::math::policies::raise_rounding_error("boost::multiprecision::iround<%1%>(%1%)", 0, v, 0LL, pol);
2186    return r.template convert_to<boost::long_long_type>();
2187 }
2188 template <class T, expression_template_option ExpressionTemplates>
llround(const number<T,ExpressionTemplates> & v)2189 inline boost::long_long_type llround(const number<T, ExpressionTemplates>& v)
2190 {
2191    return llround(v, boost::math::policies::policy<>());
2192 }
2193 #endif
2194 //
2195 // frexp does not return an expression template since we require the
2196 // integer argument to be evaluated even if the returned value is
2197 // not assigned to anything...
2198 //
2199 template <class T, expression_template_option ExpressionTemplates>
frexp(const number<T,ExpressionTemplates> & v,short * pint)2200 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)
2201 {
2202    using default_ops::eval_frexp;
2203    number<T, ExpressionTemplates> result;
2204    eval_frexp(result.backend(), v.backend(), pint);
2205    return BOOST_MP_MOVE(result);
2206 }
2207 template <class tag, class A1, class A2, class A3, class A4>
2208 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)2209    frexp(const detail::expression<tag, A1, A2, A3, A4>& v, short* pint)
2210 {
2211    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2212    return BOOST_MP_MOVE(frexp(static_cast<number_type>(v), pint));
2213 }
2214 template <class T, expression_template_option ExpressionTemplates>
frexp(const number<T,ExpressionTemplates> & v,int * pint)2215 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)
2216 {
2217    using default_ops::eval_frexp;
2218    number<T, ExpressionTemplates> result;
2219    eval_frexp(result.backend(), v.backend(), pint);
2220    return BOOST_MP_MOVE(result);
2221 }
2222 template <class tag, class A1, class A2, class A3, class A4>
2223 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)2224 frexp(const detail::expression<tag, A1, A2, A3, A4>& v, int* pint)
2225 {
2226    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2227    return BOOST_MP_MOVE(frexp(static_cast<number_type>(v), pint));
2228 }
2229 template <class T, expression_template_option ExpressionTemplates>
frexp(const number<T,ExpressionTemplates> & v,long * pint)2230 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)
2231 {
2232    using default_ops::eval_frexp;
2233    number<T, ExpressionTemplates> result;
2234    eval_frexp(result.backend(), v.backend(), pint);
2235    return BOOST_MP_MOVE(result);
2236 }
2237 template <class tag, class A1, class A2, class A3, class A4>
2238 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)2239 frexp(const detail::expression<tag, A1, A2, A3, A4>& v, long* pint)
2240 {
2241    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2242    return BOOST_MP_MOVE(frexp(static_cast<number_type>(v), pint));
2243 }
2244 template <class T, expression_template_option ExpressionTemplates>
frexp(const number<T,ExpressionTemplates> & v,boost::long_long_type * pint)2245 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)
2246 {
2247    using default_ops::eval_frexp;
2248    number<T, ExpressionTemplates> result;
2249    eval_frexp(result.backend(), v.backend(), pint);
2250    return BOOST_MP_MOVE(result);
2251 }
2252 template <class tag, class A1, class A2, class A3, class A4>
2253 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)2254 frexp(const detail::expression<tag, A1, A2, A3, A4>& v, boost::long_long_type* pint)
2255 {
2256    typedef typename detail::expression<tag, A1, A2, A3, A4>::result_type number_type;
2257    return BOOST_MP_MOVE(frexp(static_cast<number_type>(v), pint));
2258 }
2259 //
2260 // modf does not return an expression template since we require the
2261 // second argument to be evaluated even if the returned value is
2262 // not assigned to anything...
2263 //
2264 template <class T, expression_template_option ExpressionTemplates>
modf(const number<T,ExpressionTemplates> & v,number<T,ExpressionTemplates> * pipart)2265 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)
2266 {
2267    using default_ops::eval_modf;
2268    number<T, ExpressionTemplates> result;
2269    eval_modf(result.backend(), v.backend(), pipart ? &pipart->backend() : 0);
2270    return BOOST_MP_MOVE(result);
2271 }
2272 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)2273 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)
2274 {
2275    using default_ops::eval_modf;
2276    number<T, ExpressionTemplates> result, arg(v);
2277    eval_modf(result.backend(), arg.backend(), pipart ? &pipart->backend() : 0);
2278    return BOOST_MP_MOVE(result);
2279 }
2280 
2281 //
2282 // Integer square root:
2283 //
2284 template <class B, expression_template_option ExpressionTemplates>
2285 inline typename enable_if_c<number_category<B>::value == number_kind_integer, number<B, ExpressionTemplates> >::type
sqrt(const number<B,ExpressionTemplates> & x)2286    sqrt(const number<B, ExpressionTemplates>& x)
2287 {
2288    using default_ops::eval_integer_sqrt;
2289    number<B, ExpressionTemplates> s, r;
2290    eval_integer_sqrt(s.backend(), r.backend(), x.backend());
2291    return s;
2292 }
2293 //
2294 // fma:
2295 //
2296 
2297 namespace default_ops {
2298 
2299    struct fma_func
2300    {
2301       template <class B, class T, class U, class V>
operator ()boost::multiprecision::default_ops::fma_func2302       void operator()(B& result, const T& a, const U& b, const V& c)const
2303       {
2304          eval_multiply_add(result, a, b, c);
2305       }
2306    };
2307 
2308 
2309 }
2310 
2311 template <class Backend, class U, class V>
2312 inline typename enable_if<
2313    mpl::and_<
2314       mpl::bool_<number_category<number<Backend, et_on> >::value == number_kind_floating_point>,
2315       mpl::or_<
2316          is_number<U>,
2317          is_number_expression<U>,
2318          is_arithmetic<U>
2319       >,
2320       mpl::or_<
2321          is_number<V>,
2322          is_number_expression<V>,
2323          is_arithmetic<V>
2324       >
2325    >,
2326    detail::expression<detail::function, default_ops::fma_func, number<Backend, et_on>, U, V>
2327 >::type
fma(const number<Backend,et_on> & a,const U & b,const V & c)2328 fma(const number<Backend, et_on>& a, const U& b, const V& c)
2329 {
2330    return detail::expression<detail::function, default_ops::fma_func, number<Backend, et_on>, U, V>(
2331       default_ops::fma_func(), a, b, c);
2332 }
2333 
2334 template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class U, class V>
2335 inline typename enable_if<
2336    mpl::and_<
2337    mpl::bool_<number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point>,
2338    mpl::or_<
2339    is_number<U>,
2340    is_number_expression<U>,
2341    is_arithmetic<U>
2342    >,
2343    mpl::or_<
2344    is_number<V>,
2345    is_number_expression<V>,
2346    is_arithmetic<V>
2347    >
2348    >,
2349    detail::expression<detail::function, default_ops::fma_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, V>
2350 >::type
fma(const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & a,const U & b,const V & c)2351 fma(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& a, const U& b, const V& c)
2352 {
2353    return detail::expression<detail::function, default_ops::fma_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, V>(
2354       default_ops::fma_func(), a, b, c);
2355 }
2356 
2357 template <class Backend, class U, class V>
2358 inline typename enable_if<
2359    mpl::and_<
2360    mpl::bool_<number_category<number<Backend, et_off> >::value == number_kind_floating_point>,
2361    mpl::or_<
2362    is_number<U>,
2363    is_number_expression<U>,
2364    is_arithmetic<U>
2365    >,
2366    mpl::or_<
2367    is_number<V>,
2368    is_number_expression<V>,
2369    is_arithmetic<V>
2370    >
2371    >,
2372    number<Backend, et_off>
2373 >::type
fma(const number<Backend,et_off> & a,const U & b,const V & c)2374 fma(const number<Backend, et_off>& a, const U& b, const V& c)
2375 {
2376    using default_ops::eval_multiply_add;
2377    number<Backend, et_off> result;
2378    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));
2379    return BOOST_MP_MOVE(result);
2380 }
2381 
2382 template <class U, class Backend, class V>
2383 inline typename enable_if<
2384    mpl::and_<
2385       mpl::bool_<number_category<number<Backend, et_on> >::value == number_kind_floating_point>,
2386       is_arithmetic<U>,
2387       mpl::or_<
2388          is_number<V>,
2389          is_number_expression<V>,
2390          is_arithmetic<V>
2391       >
2392    >,
2393    detail::expression<detail::function, default_ops::fma_func, U, number<Backend, et_on>, V>
2394 >::type
fma(const U & a,const number<Backend,et_on> & b,const V & c)2395 fma(const U& a, const number<Backend, et_on>& b, const V& c)
2396 {
2397    return detail::expression<detail::function, default_ops::fma_func, U, number<Backend, et_on>, V>(
2398       default_ops::fma_func(), a, b, c);
2399 }
2400 
2401 template <class U, class tag, class Arg1, class Arg2, class Arg3, class Arg4, class V>
2402 inline typename enable_if<
2403    mpl::and_<
2404       mpl::bool_<number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point>,
2405       is_arithmetic<U>,
2406       mpl::or_<
2407          is_number<V>,
2408          is_number_expression<V>,
2409          is_arithmetic<V>
2410       >
2411    >,
2412    detail::expression<detail::function, default_ops::fma_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, V>
2413 >::type
fma(const U & a,const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & b,const V & c)2414 fma(const U& a, const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& b, const V& c)
2415 {
2416    return detail::expression<detail::function, default_ops::fma_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, V>(
2417       default_ops::fma_func(), a, b, c);
2418 }
2419 
2420 template <class U, class Backend, class V>
2421 inline typename enable_if<
2422    mpl::and_<
2423       mpl::bool_<number_category<number<Backend, et_off> >::value == number_kind_floating_point>,
2424       is_arithmetic<U>,
2425       mpl::or_<
2426          is_number<V>,
2427          is_number_expression<V>,
2428          is_arithmetic<V>
2429       >
2430    >,
2431    number<Backend, et_off>
2432 >::type
fma(const U & a,const number<Backend,et_off> & b,const V & c)2433 fma(const U& a, const number<Backend, et_off>& b, const V& c)
2434 {
2435    using default_ops::eval_multiply_add;
2436    number<Backend, et_off> result;
2437    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));
2438    return BOOST_MP_MOVE(result);
2439 }
2440 
2441 template <class U, class V, class Backend>
2442 inline typename enable_if<
2443    mpl::and_<
2444    mpl::bool_<number_category<number<Backend, et_on> >::value == number_kind_floating_point>,
2445       is_arithmetic<U>,
2446       is_arithmetic<V>
2447    >,
2448    detail::expression<detail::function, default_ops::fma_func, U, V, number<Backend, et_on> >
2449 >::type
fma(const U & a,const V & b,const number<Backend,et_on> & c)2450 fma(const U& a, const V& b, const number<Backend, et_on>& c)
2451 {
2452    return detail::expression<detail::function, default_ops::fma_func, U, V, number<Backend, et_on> >(
2453       default_ops::fma_func(), a, b, c);
2454 }
2455 
2456 template <class U, class V, class tag, class Arg1, class Arg2, class Arg3, class Arg4>
2457 inline typename enable_if<
2458    mpl::and_<
2459    mpl::bool_<number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point>,
2460       is_arithmetic<U>,
2461       is_arithmetic<V>
2462    >,
2463    detail::expression<detail::function, default_ops::fma_func, U, V, detail::expression<tag, Arg1, Arg2, Arg3, Arg4> >
2464 >::type
fma(const U & a,const V & b,const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & c)2465 fma(const U& a, const V& b, const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& c)
2466 {
2467    return detail::expression<detail::function, default_ops::fma_func, U, V, detail::expression<tag, Arg1, Arg2, Arg3, Arg4> >(
2468       default_ops::fma_func(), a, b, c);
2469 }
2470 
2471 template <class U, class V, class Backend>
2472 inline typename enable_if<
2473    mpl::and_<
2474    mpl::bool_<number_category<number<Backend, et_off> >::value == number_kind_floating_point>,
2475       is_arithmetic<U>,
2476       is_arithmetic<V>
2477    >,
2478    number<Backend, et_off>
2479 >::type
fma(const U & a,const V & b,const number<Backend,et_off> & c)2480 fma(const U& a, const V& b, const number<Backend, et_off>& c)
2481 {
2482    using default_ops::eval_multiply_add;
2483    number<Backend, et_off> result;
2484    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));
2485    return BOOST_MP_MOVE(result);
2486 }
2487 
2488 namespace default_ops {
2489 
2490    struct remquo_func
2491    {
2492       template <class B, class T, class U>
operator ()boost::multiprecision::default_ops::remquo_func2493       void operator()(B& result, const T& a, const U& b, int* pi)const
2494       {
2495          eval_remquo(result, a, b, pi);
2496       }
2497    };
2498 
2499 }
2500 
2501 template <class Backend, class U>
2502 inline typename enable_if_c<
2503    number_category<number<Backend, et_on> >::value == number_kind_floating_point,
2504    detail::expression<detail::function, default_ops::remquo_func, number<Backend, et_on>, U, int*>
2505 >::type
remquo(const number<Backend,et_on> & a,const U & b,int * pi)2506 remquo(const number<Backend, et_on>& a, const U& b, int* pi)
2507 {
2508    return detail::expression<detail::function, default_ops::remquo_func, number<Backend, et_on>, U, int*>(
2509       default_ops::remquo_func(), a, b, pi);
2510 }
2511 
2512 template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class U>
2513 inline typename enable_if_c<
2514    number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point,
2515    detail::expression<detail::function, default_ops::remquo_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, int*>
2516 >::type
remquo(const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & a,const U & b,int * pi)2517 remquo(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& a, const U& b, int* pi)
2518 {
2519    return detail::expression<detail::function, default_ops::remquo_func, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, U, int*>(
2520       default_ops::remquo_func(), a, b, pi);
2521 }
2522 
2523 template <class U, class Backend>
2524 inline typename enable_if_c<
2525    (number_category<number<Backend, et_on> >::value == number_kind_floating_point)
2526    && !is_number<U>::value && !is_number_expression<U>::value,
2527    detail::expression<detail::function, default_ops::remquo_func, U, number<Backend, et_on>, int*>
2528 >::type
remquo(const U & a,const number<Backend,et_on> & b,int * pi)2529 remquo(const U& a, const number<Backend, et_on>& b, int* pi)
2530 {
2531    return detail::expression<detail::function, default_ops::remquo_func, U, number<Backend, et_on>, int*>(
2532       default_ops::remquo_func(), a, b, pi);
2533 }
2534 
2535 template <class U, class tag, class Arg1, class Arg2, class Arg3, class Arg4>
2536 inline typename enable_if_c<
2537    (number_category<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type >::value == number_kind_floating_point)
2538    && !is_number<U>::value && !is_number_expression<U>::value,
2539    detail::expression<detail::function, default_ops::remquo_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, int*>
2540 >::type
remquo(const U & a,const detail::expression<tag,Arg1,Arg2,Arg3,Arg4> & b,int * pi)2541 remquo(const U& a, const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& b, int* pi)
2542 {
2543    return detail::expression<detail::function, default_ops::remquo_func, U, detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, int*>(
2544       default_ops::remquo_func(), a, b, pi);
2545 }
2546 
2547 template <class Backend, class U>
2548 inline typename enable_if_c<
2549    number_category<number<Backend, et_on> >::value == number_kind_floating_point,
2550    number<Backend, et_off>
2551 >::type
remquo(const number<Backend,et_off> & a,const U & b,int * pi)2552 remquo(const number<Backend, et_off>& a, const U& b, int* pi)
2553 {
2554    using default_ops::eval_remquo;
2555    number<Backend, et_off> result;
2556    eval_remquo(result.backend(), a.backend(), number<Backend, et_off>::canonical_value(b), pi);
2557    return BOOST_MP_MOVE(result);
2558 }
2559 template <class U, class Backend>
2560 inline typename enable_if_c<
2561 (number_category<number<Backend, et_on> >::value == number_kind_floating_point)
2562 && !is_number<U>::value && !is_number_expression<U>::value,
2563 number<Backend, et_off>
2564 >::type
remquo(const U & a,const number<Backend,et_off> & b,int * pi)2565 remquo(const U& a, const number<Backend, et_off>& b, int* pi)
2566 {
2567    using default_ops::eval_remquo;
2568    number<Backend, et_off> result;
2569    eval_remquo(result.backend(), number<Backend, et_off>::canonical_value(a), b.backend(), pi);
2570    return BOOST_MP_MOVE(result);
2571 }
2572 
2573 
2574 template <class B, expression_template_option ExpressionTemplates>
2575 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)2576    sqrt(const number<B, ExpressionTemplates>& x, number<B, ExpressionTemplates>& r)
2577 {
2578    using default_ops::eval_integer_sqrt;
2579    number<B, ExpressionTemplates> s;
2580    eval_integer_sqrt(s.backend(), r.backend(), x.backend());
2581    return s;
2582 }
2583 
2584 #define UNARY_OP_FUNCTOR(func, category)\
2585 namespace detail{\
2586 template <class Backend> \
2587 struct BOOST_JOIN(func, _funct)\
2588 {\
2589    void operator()(Backend& result, const Backend& arg)const\
2590    {\
2591       using default_ops::BOOST_JOIN(eval_,func);\
2592       BOOST_JOIN(eval_,func)(result, arg);\
2593    }\
2594 };\
2595 \
2596 }\
2597 \
2598 template <class tag, class A1, class A2, class A3, class A4> \
2599 inline typename enable_if_c<number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category,\
2600    detail::expression<\
2601     detail::function\
2602   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2603   , detail::expression<tag, A1, A2, A3, A4> > \
2604 >::type \
2605 func(const detail::expression<tag, A1, A2, A3, A4>& arg)\
2606 {\
2607     return detail::expression<\
2608     detail::function\
2609   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2610   , detail::expression<tag, A1, A2, A3, A4> \
2611 > (\
2612         detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2613       , arg   \
2614     );\
2615 }\
2616 template <class Backend> \
2617 inline typename enable_if_c<number_category<Backend>::value == category,\
2618    detail::expression<\
2619     detail::function\
2620   , detail::BOOST_JOIN(func, _funct)<Backend> \
2621   , number<Backend, et_on> > \
2622 >::type \
2623 func(const number<Backend, et_on>& arg)\
2624 {\
2625     return detail::expression<\
2626     detail::function\
2627   , detail::BOOST_JOIN(func, _funct)<Backend> \
2628   , number<Backend, et_on> \
2629   >(\
2630         detail::BOOST_JOIN(func, _funct)<Backend>() \
2631       , arg   \
2632     );\
2633 }\
2634 template <class Backend> \
2635 inline typename boost::enable_if_c<\
2636    boost::multiprecision::number_category<Backend>::value == category,\
2637    number<Backend, et_off> >::type \
2638 func(const number<Backend, et_off>& arg)\
2639 {\
2640    number<Backend, et_off> result;\
2641    using default_ops::BOOST_JOIN(eval_,func);\
2642    BOOST_JOIN(eval_,func)(result.backend(), arg.backend());\
2643    return BOOST_MP_MOVE(result);\
2644 }
2645 
2646 #define BINARY_OP_FUNCTOR(func, category)\
2647 namespace detail{\
2648 template <class Backend> \
2649 struct BOOST_JOIN(func, _funct)\
2650 {\
2651    void operator()(Backend& result, const Backend& arg, const Backend& a)const\
2652    {\
2653       using default_ops:: BOOST_JOIN(eval_,func);\
2654       BOOST_JOIN(eval_,func)(result, arg, a);\
2655    }\
2656    template <class Arithmetic> \
2657    void operator()(Backend& result, const Backend& arg, const Arithmetic& a)const\
2658    {\
2659       using default_ops:: BOOST_JOIN(eval_,func);\
2660       BOOST_JOIN(eval_,func)(result, arg, a);\
2661    }\
2662    template <class Arithmetic> \
2663    void operator()(Backend& result, const Arithmetic& arg, const Backend& a)const\
2664    {\
2665       using default_ops:: BOOST_JOIN(eval_,func);\
2666       BOOST_JOIN(eval_,func)(result, arg, a);\
2667    }\
2668 };\
2669 \
2670 }\
2671 template <class Backend> \
2672 inline typename enable_if_c<number_category<Backend>::value == category,\
2673    detail::expression<\
2674     detail::function\
2675   , detail::BOOST_JOIN(func, _funct)<Backend> \
2676   , number<Backend, et_on> \
2677   , number<Backend, et_on> > \
2678 >::type \
2679 func(const number<Backend, et_on>& arg, const number<Backend, et_on>& a)\
2680 {\
2681     return detail::expression<\
2682     detail::function\
2683   , detail::BOOST_JOIN(func, _funct)<Backend> \
2684   , number<Backend, et_on> \
2685   , number<Backend, et_on> \
2686   >(\
2687         detail::BOOST_JOIN(func, _funct)<Backend>() \
2688       , arg,\
2689       a\
2690     );\
2691 }\
2692 template <class Backend, class tag, class A1, class A2, class A3, class A4> \
2693 inline typename enable_if_c<\
2694    (number_category<Backend>::value == category) && (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2695    detail::expression<\
2696     detail::function\
2697   , detail::BOOST_JOIN(func, _funct)<Backend> \
2698   , number<Backend, et_on> \
2699   , detail::expression<tag, A1, A2, A3, A4> > \
2700 >::type \
2701 func(const number<Backend, et_on>& arg, const detail::expression<tag, A1, A2, A3, A4>& a)\
2702 {\
2703     return detail::expression<\
2704     detail::function\
2705   , detail::BOOST_JOIN(func, _funct)<Backend> \
2706   , number<Backend, et_on> \
2707   , detail::expression<tag, A1, A2, A3, A4> \
2708   >(\
2709         detail::BOOST_JOIN(func, _funct)<Backend>() \
2710       , arg,\
2711       a\
2712     );\
2713 }\
2714 template <class tag, class A1, class A2, class A3, class A4, class Backend> \
2715 inline typename enable_if_c<\
2716    (number_category<Backend>::value == category) && (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2717    detail::expression<\
2718     detail::function\
2719   , detail::BOOST_JOIN(func, _funct)<Backend> \
2720   , detail::expression<tag, A1, A2, A3, A4> \
2721   , number<Backend, et_on> > \
2722 >::type \
2723 func(const detail::expression<tag, A1, A2, A3, A4>& arg, const number<Backend, et_on>& a)\
2724 {\
2725     return detail::expression<\
2726     detail::function\
2727   , detail::BOOST_JOIN(func, _funct)<Backend> \
2728   , detail::expression<tag, A1, A2, A3, A4> \
2729   , number<Backend, et_on> \
2730   >(\
2731         detail::BOOST_JOIN(func, _funct)<Backend>() \
2732       , arg,\
2733       a\
2734     );\
2735 }\
2736 template <class tag, class A1, class A2, class A3, class A4, class tagb, class A1b, class A2b, class A3b, class A4b> \
2737 inline typename enable_if_c<\
2738       (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category) && (number_category<detail::expression<tagb, A1b, A2b, A3b, A4b> >::value == category),\
2739    detail::expression<\
2740     detail::function\
2741   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2742   , detail::expression<tag, A1, A2, A3, A4> \
2743   , detail::expression<tagb, A1b, A2b, A3b, A4b> > \
2744 >::type \
2745 func(const detail::expression<tag, A1, A2, A3, A4>& arg, const detail::expression<tagb, A1b, A2b, A3b, A4b>& a)\
2746 {\
2747     return detail::expression<\
2748     detail::function\
2749   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2750   , detail::expression<tag, A1, A2, A3, A4> \
2751   , detail::expression<tagb, A1b, A2b, A3b, A4b> \
2752   >(\
2753         detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2754       , arg,\
2755       a\
2756     );\
2757 }\
2758 template <class Backend, class Arithmetic> \
2759 inline typename enable_if_c<\
2760    is_arithmetic<Arithmetic>::value && (number_category<Backend>::value == category),\
2761    detail::expression<\
2762     detail::function\
2763   , detail::BOOST_JOIN(func, _funct)<Backend> \
2764   , number<Backend, et_on> \
2765   , Arithmetic\
2766   > \
2767 >::type \
2768 func(const number<Backend, et_on>& arg, const Arithmetic& a)\
2769 {\
2770     return detail::expression<\
2771     detail::function\
2772   , detail::BOOST_JOIN(func, _funct)<Backend> \
2773   , number<Backend, et_on> \
2774   , Arithmetic\
2775   >(\
2776         detail::BOOST_JOIN(func, _funct)<Backend>() \
2777       , arg,\
2778       a\
2779     );\
2780 }\
2781 template <class tag, class A1, class A2, class A3, class A4, class Arithmetic> \
2782 inline typename enable_if_c<\
2783    is_arithmetic<Arithmetic>::value && (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2784    detail::expression<\
2785     detail::function\
2786   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2787   , detail::expression<tag, A1, A2, A3, A4> \
2788   , Arithmetic\
2789   > \
2790 >::type \
2791 func(const detail::expression<tag, A1, A2, A3, A4>& arg, const Arithmetic& a)\
2792 {\
2793     return detail::expression<\
2794     detail::function\
2795   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2796   , detail::expression<tag, A1, A2, A3, A4> \
2797   , Arithmetic\
2798    >(\
2799         detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2800       , arg,\
2801       a\
2802     );\
2803 }\
2804 template <class Backend, class Arithmetic> \
2805 inline typename enable_if_c<\
2806    is_arithmetic<Arithmetic>::value && (number_category<Backend>::value == category),\
2807    detail::expression<\
2808     detail::function\
2809   , detail::BOOST_JOIN(func, _funct)<Backend> \
2810   , Arithmetic \
2811   , number<Backend, et_on> \
2812   > \
2813 >::type \
2814 func(const Arithmetic& arg, const number<Backend, et_on>& a)\
2815 {\
2816     return detail::expression<\
2817     detail::function\
2818   , detail::BOOST_JOIN(func, _funct)<Backend> \
2819   , Arithmetic \
2820   , number<Backend, et_on> \
2821   >(\
2822         detail::BOOST_JOIN(func, _funct)<Backend>() \
2823       , arg,\
2824       a\
2825     );\
2826 }\
2827 template <class tag, class A1, class A2, class A3, class A4, class Arithmetic> \
2828 inline typename enable_if_c<\
2829    is_arithmetic<Arithmetic>::value && (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2830    detail::expression<\
2831     detail::function\
2832   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2833   , Arithmetic \
2834   , detail::expression<tag, A1, A2, A3, A4> \
2835   > \
2836 >::type \
2837 func(const Arithmetic& arg, const detail::expression<tag, A1, A2, A3, A4>& a)\
2838 {\
2839     return detail::expression<\
2840     detail::function\
2841   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2842   , Arithmetic \
2843   , detail::expression<tag, A1, A2, A3, A4> \
2844    >(\
2845         detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2846       , arg,\
2847       a\
2848     );\
2849 }\
2850 template <class Backend> \
2851 inline typename enable_if_c<(number_category<Backend>::value == category),\
2852    number<Backend, et_off> >::type \
2853 func(const number<Backend, et_off>& arg, const number<Backend, et_off>& a)\
2854 {\
2855    number<Backend, et_off> result;\
2856    using default_ops:: BOOST_JOIN(eval_,func);\
2857    BOOST_JOIN(eval_,func)(result.backend(), arg.backend(), a.backend());\
2858    return BOOST_MP_MOVE(result);\
2859 }\
2860 template <class Backend, class Arithmetic> \
2861 inline typename enable_if_c<\
2862    is_arithmetic<Arithmetic>::value && (number_category<Backend>::value == category),\
2863    number<Backend, et_off> \
2864 >::type \
2865 func(const number<Backend, et_off>& arg, const Arithmetic& a)\
2866 {\
2867    typedef typename detail::canonical<Arithmetic, Backend>::type canonical_type;\
2868    number<Backend, et_off> result;\
2869    using default_ops:: BOOST_JOIN(eval_,func);\
2870    BOOST_JOIN(eval_,func)(result.backend(), arg.backend(), static_cast<canonical_type>(a));\
2871    return BOOST_MP_MOVE(result);\
2872 }\
2873 template <class Backend, class Arithmetic> \
2874 inline typename enable_if_c<\
2875    is_arithmetic<Arithmetic>::value && (number_category<Backend>::value == category),\
2876    number<Backend, et_off> \
2877 >::type \
2878 func(const Arithmetic& a, const number<Backend, et_off>& arg)\
2879 {\
2880    typedef typename detail::canonical<Arithmetic, Backend>::type canonical_type;\
2881    number<Backend, et_off> result;\
2882    using default_ops:: BOOST_JOIN(eval_,func);\
2883    BOOST_JOIN(eval_,func)(result.backend(), static_cast<canonical_type>(a), arg.backend());\
2884    return BOOST_MP_MOVE(result);\
2885 }\
2886 
2887 
2888 #define HETERO_BINARY_OP_FUNCTOR_B(func, Arg2, category)\
2889 template <class tag, class A1, class A2, class A3, class A4> \
2890 inline typename enable_if_c<\
2891    (number_category<detail::expression<tag, A1, A2, A3, A4> >::value == category),\
2892    detail::expression<\
2893     detail::function\
2894   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2895   , detail::expression<tag, A1, A2, A3, A4> \
2896   , Arg2> \
2897 >::type \
2898 func(const detail::expression<tag, A1, A2, A3, A4>& arg, Arg2 const& a)\
2899 {\
2900     return detail::expression<\
2901     detail::function\
2902   , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type> \
2903   , detail::expression<tag, A1, A2, A3, A4> \
2904   , Arg2\
2905    >(\
2906         detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>() \
2907       , arg, a   \
2908     );\
2909 }\
2910 template <class Backend> \
2911 inline typename enable_if_c<\
2912    (number_category<Backend>::value == category),\
2913   detail::expression<\
2914     detail::function\
2915   , detail::BOOST_JOIN(func, _funct)<Backend> \
2916   , number<Backend, et_on> \
2917   , Arg2> \
2918 >::type \
2919 func(const number<Backend, et_on>& arg, Arg2 const& a)\
2920 {\
2921     return detail::expression<\
2922     detail::function\
2923   , detail::BOOST_JOIN(func, _funct)<Backend> \
2924   , number<Backend, et_on> \
2925   , Arg2\
2926   >(\
2927         detail::BOOST_JOIN(func, _funct)<Backend>() \
2928       , arg,\
2929       a\
2930     );\
2931 }\
2932 template <class Backend> \
2933 inline typename enable_if_c<\
2934   (number_category<Backend>::value == category),\
2935   number<Backend, et_off> >::type \
2936 func(const number<Backend, et_off>& arg, Arg2 const& a)\
2937 {\
2938    number<Backend, et_off> result;\
2939    using default_ops:: BOOST_JOIN(eval_,func);\
2940    BOOST_JOIN(eval_,func)(result.backend(), arg.backend(), a);\
2941    return BOOST_MP_MOVE(result);\
2942 }\
2943 
2944 #define HETERO_BINARY_OP_FUNCTOR(func, Arg2, category)\
2945 namespace detail{\
2946 template <class Backend> \
2947 struct BOOST_JOIN(func, _funct)\
2948 {\
2949    template <class Arg>\
2950    void operator()(Backend& result, Backend const& arg, Arg a)const\
2951    {\
2952       using default_ops:: BOOST_JOIN(eval_,func);\
2953       BOOST_JOIN(eval_,func)(result, arg, a);\
2954    }\
2955 };\
2956 \
2957 }\
2958 \
2959 HETERO_BINARY_OP_FUNCTOR_B(func, Arg2, category)
2960 
2961 namespace detail{
2962 template <class Backend>
2963 struct abs_funct
2964 {
operator ()boost::multiprecision::detail::abs_funct2965    void operator()(Backend& result, const Backend& arg)const
2966    {
2967       using default_ops::eval_abs;
2968       eval_abs(result, arg);
2969    }
2970 };
2971 
2972 }
2973 
2974 template <class tag, class A1, class A2, class A3, class A4>
2975 inline detail::expression<
2976     detail::function
2977   , detail::abs_funct<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>
2978   , detail::expression<tag, A1, A2, A3, A4> >
abs(const detail::expression<tag,A1,A2,A3,A4> & arg)2979 abs(const detail::expression<tag, A1, A2, A3, A4>& arg)
2980 {
2981     return detail::expression<
2982     detail::function
2983   , detail::abs_funct<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>
2984   , detail::expression<tag, A1, A2, A3, A4>
2985 > (
2986         detail::abs_funct<typename detail::backend_type<detail::expression<tag, A1, A2, A3, A4> >::type>()
2987       , arg
2988     );
2989 }
2990 template <class Backend>
2991 inline detail::expression<
2992     detail::function
2993   , detail::abs_funct<Backend>
2994   , number<Backend, et_on> >
abs(const number<Backend,et_on> & arg)2995 abs(const number<Backend, et_on>& arg)
2996 {
2997     return detail::expression<
2998     detail::function
2999   , detail::abs_funct<Backend>
3000   , number<Backend, et_on>
3001   >(
3002         detail::abs_funct<Backend>()
3003       , arg
3004     );
3005 }
3006 template <class Backend>
3007 inline number<Backend, et_off>
abs(const number<Backend,et_off> & arg)3008 abs(const number<Backend, et_off>& arg)
3009 {
3010    number<Backend, et_off> result;
3011    using default_ops::eval_abs;
3012    eval_abs(result.backend(), arg.backend());
3013    return BOOST_MP_MOVE(result);
3014 }
3015 
UNARY_OP_FUNCTOR(fabs,number_kind_floating_point)3016 UNARY_OP_FUNCTOR(fabs, number_kind_floating_point)
3017 UNARY_OP_FUNCTOR(sqrt, number_kind_floating_point)
3018 UNARY_OP_FUNCTOR(floor, number_kind_floating_point)
3019 UNARY_OP_FUNCTOR(ceil, number_kind_floating_point)
3020 UNARY_OP_FUNCTOR(trunc, number_kind_floating_point)
3021 UNARY_OP_FUNCTOR(round, number_kind_floating_point)
3022 UNARY_OP_FUNCTOR(exp, number_kind_floating_point)
3023 UNARY_OP_FUNCTOR(exp2, number_kind_floating_point)
3024 UNARY_OP_FUNCTOR(log, number_kind_floating_point)
3025 UNARY_OP_FUNCTOR(log10, number_kind_floating_point)
3026 UNARY_OP_FUNCTOR(cos, number_kind_floating_point)
3027 UNARY_OP_FUNCTOR(sin, number_kind_floating_point)
3028 UNARY_OP_FUNCTOR(tan, number_kind_floating_point)
3029 UNARY_OP_FUNCTOR(asin, number_kind_floating_point)
3030 UNARY_OP_FUNCTOR(acos, number_kind_floating_point)
3031 UNARY_OP_FUNCTOR(atan, number_kind_floating_point)
3032 UNARY_OP_FUNCTOR(cosh, number_kind_floating_point)
3033 UNARY_OP_FUNCTOR(sinh, number_kind_floating_point)
3034 UNARY_OP_FUNCTOR(tanh, number_kind_floating_point)
3035 UNARY_OP_FUNCTOR(log2, number_kind_floating_point)
3036 UNARY_OP_FUNCTOR(nearbyint, number_kind_floating_point)
3037 UNARY_OP_FUNCTOR(rint, number_kind_floating_point)
3038 
3039 HETERO_BINARY_OP_FUNCTOR(ldexp, short, number_kind_floating_point)
3040 //HETERO_BINARY_OP_FUNCTOR(frexp, short*, number_kind_floating_point)
3041 HETERO_BINARY_OP_FUNCTOR_B(ldexp, int, number_kind_floating_point)
3042 //HETERO_BINARY_OP_FUNCTOR_B(frexp, int*, number_kind_floating_point)
3043 HETERO_BINARY_OP_FUNCTOR_B(ldexp, long, number_kind_floating_point)
3044 //HETERO_BINARY_OP_FUNCTOR_B(frexp, long*, number_kind_floating_point)
3045 HETERO_BINARY_OP_FUNCTOR_B(ldexp, boost::long_long_type, number_kind_floating_point)
3046 //HETERO_BINARY_OP_FUNCTOR_B(frexp, boost::long_long_type*, number_kind_floating_point)
3047 BINARY_OP_FUNCTOR(pow, number_kind_floating_point)
3048 BINARY_OP_FUNCTOR(fmod, number_kind_floating_point)
3049 BINARY_OP_FUNCTOR(fmax, number_kind_floating_point)
3050 BINARY_OP_FUNCTOR(fmin, number_kind_floating_point)
3051 BINARY_OP_FUNCTOR(atan2, number_kind_floating_point)
3052 BINARY_OP_FUNCTOR(fdim, number_kind_floating_point)
3053 BINARY_OP_FUNCTOR(hypot, number_kind_floating_point)
3054 BINARY_OP_FUNCTOR(remainder, number_kind_floating_point)
3055 
3056 UNARY_OP_FUNCTOR(logb, number_kind_floating_point)
3057 HETERO_BINARY_OP_FUNCTOR(scalbn, short, number_kind_floating_point)
3058 HETERO_BINARY_OP_FUNCTOR(scalbln, short, number_kind_floating_point)
3059 HETERO_BINARY_OP_FUNCTOR_B(scalbn, int, number_kind_floating_point)
3060 HETERO_BINARY_OP_FUNCTOR_B(scalbln, int, number_kind_floating_point)
3061 HETERO_BINARY_OP_FUNCTOR_B(scalbn, long, number_kind_floating_point)
3062 HETERO_BINARY_OP_FUNCTOR_B(scalbln, long, number_kind_floating_point)
3063 HETERO_BINARY_OP_FUNCTOR_B(scalbn, boost::long_long_type, number_kind_floating_point)
3064 HETERO_BINARY_OP_FUNCTOR_B(scalbln, boost::long_long_type, number_kind_floating_point)
3065 
3066 //
3067 // Integer functions:
3068 //
3069 BINARY_OP_FUNCTOR(gcd, number_kind_integer)
3070 BINARY_OP_FUNCTOR(lcm, number_kind_integer)
3071 HETERO_BINARY_OP_FUNCTOR_B(pow, unsigned, number_kind_integer)
3072 
3073 #undef BINARY_OP_FUNCTOR
3074 #undef UNARY_OP_FUNCTOR
3075 
3076 //
3077 // ilogb:
3078 //
3079 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
3080 inline typename enable_if_c<number_category<Backend>::value == number_kind_floating_point, typename Backend::exponent_type>::type
3081    ilogb(const multiprecision::number<Backend, ExpressionTemplates>& val)
3082 {
3083    using default_ops::eval_ilogb;
3084    return eval_ilogb(val.backend());
3085 }
3086 
3087 template <class tag, class A1, class A2, class A3, class A4>
3088 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)3089 ilogb(const detail::expression<tag, A1, A2, A3, A4>& val)
3090 {
3091    using default_ops::eval_ilogb;
3092    typename multiprecision::detail::expression<tag, A1, A2, A3, A4>::result_type arg(val);
3093    return eval_ilogb(arg.backend());
3094 }
3095 
3096 } //namespace multiprecision
3097 
3098 namespace math{
3099 //
3100 // Overload of Boost.Math functions that find the wrong overload when used with number:
3101 //
3102 namespace detail{
3103    template <class T> T sinc_pi_imp(T);
3104    template <class T> T sinhc_pi_imp(T);
3105 }
3106 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
sinc_pi(const multiprecision::number<Backend,ExpressionTemplates> & x)3107 inline multiprecision::number<Backend, ExpressionTemplates> sinc_pi(const multiprecision::number<Backend, ExpressionTemplates>& x)
3108 {
3109    return BOOST_MP_MOVE(detail::sinc_pi_imp(x));
3110 }
3111 
3112 template <class Backend, multiprecision::expression_template_option ExpressionTemplates, class Policy>
sinc_pi(const multiprecision::number<Backend,ExpressionTemplates> & x,const Policy &)3113 inline multiprecision::number<Backend, ExpressionTemplates> sinc_pi(const multiprecision::number<Backend, ExpressionTemplates>& x, const Policy&)
3114 {
3115    return BOOST_MP_MOVE(detail::sinc_pi_imp(x));
3116 }
3117 
3118 template <class Backend, multiprecision::expression_template_option ExpressionTemplates>
sinhc_pi(const multiprecision::number<Backend,ExpressionTemplates> & x)3119 inline multiprecision::number<Backend, ExpressionTemplates> sinhc_pi(const multiprecision::number<Backend, ExpressionTemplates>& x)
3120 {
3121    return BOOST_MP_MOVE(detail::sinhc_pi_imp(x));
3122 }
3123 
3124 template <class Backend, multiprecision::expression_template_option ExpressionTemplates, class Policy>
sinhc_pi(const multiprecision::number<Backend,ExpressionTemplates> & x,const Policy &)3125 inline multiprecision::number<Backend, ExpressionTemplates> sinhc_pi(const multiprecision::number<Backend, ExpressionTemplates>& x, const Policy&)
3126 {
3127    return BOOST_MP_MOVE(boost::math::sinhc_pi(x));
3128 }
3129 
3130 #ifdef BOOST_MSVC
3131 #pragma warning(pop)
3132 #endif
3133 } // namespace math
3134 } // namespace boost
3135 
3136 //
3137 // This has to come last of all:
3138 //
3139 #include <boost/multiprecision/detail/no_et_ops.hpp>
3140 #include <boost/multiprecision/detail/et_ops.hpp>
3141 //
3142 // min/max overloads:
3143 //
3144 #include <boost/multiprecision/detail/min_max.hpp>
3145 
3146 #endif
3147 
3148