boost/multiprecision/detail/default_ops.hpp

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