ceph/src/boost/libs/math/include/boost/math/bindings/rr.hpp

   1 //  Copyright John Maddock 2007.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost 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_NTL_RR_HPP
   7 #define BOOST_MATH_NTL_RR_HPP
   8
   9 #include <boost/config.hpp>
  10 #include <boost/limits.hpp>
  11 #include <boost/math/tools/real_cast.hpp>
  12 #include <boost/math/tools/precision.hpp>
  13 #include <boost/math/constants/constants.hpp>
  14 #include <boost/math/tools/roots.hpp>
  15 #include <boost/math/special_functions/fpclassify.hpp>
  16 #include <boost/math/bindings/detail/big_digamma.hpp>
  17 #include <boost/math/bindings/detail/big_lanczos.hpp>
  18
  19 #include <ostream>
  20 #include <istream>
  21 #include <boost/config/no_tr1/cmath.hpp>
  22 #include <NTL/RR.h>
  23
  24 namespace boost{ namespace math{
  25
  26 namespace ntl
  27 {
  28
  29 class RR;
  30
  31 RR ldexp(RR r, int exp);
  32 RR frexp(RR r, int* exp);
  33
  34 class RR
  35 {
  36 public:
  37    // Constructors:
  38    RR() {}
  39    RR(const ::NTL::RR& c) : m_value(c){}
  40    RR(char c)
  41    {
  42       m_value = c;
  43    }
  44 #ifndef BOOST_NO_INTRINSIC_WCHAR_T
  45    RR(wchar_t c)
  46    {
  47       m_value = c;
  48    }
  49 #endif
  50    RR(unsigned char c)
  51    {
  52       m_value = c;
  53    }
  54    RR(signed char c)
  55    {
  56       m_value = c;
  57    }
  58    RR(unsigned short c)
  59    {
  60       m_value = c;
  61    }
  62    RR(short c)
  63    {
  64       m_value = c;
  65    }
  66    RR(unsigned int c)
  67    {
  68       assign_large_int(c);
  69    }
  70    RR(int c)
  71    {
  72       assign_large_int(c);
  73    }
  74    RR(unsigned long c)
  75    {
  76       assign_large_int(c);
  77    }
  78    RR(long c)
  79    {
  80       assign_large_int(c);
  81    }
  82 #ifdef BOOST_HAS_LONG_LONG
  83    RR(boost::ulong_long_type c)
  84    {
  85       assign_large_int(c);
  86    }
  87    RR(boost::long_long_type c)
  88    {
  89       assign_large_int(c);
  90    }
  91 #endif
  92    RR(float c)
  93    {
  94       m_value = c;
  95    }
  96    RR(double c)
  97    {
  98       m_value = c;
  99    }
 100    RR(long double c)
 101    {
 102       assign_large_real(c);
 103    }
 104
 105    // Assignment:
 106    RR& operator=(char c) { m_value = c; return *this; }
 107    RR& operator=(unsigned char c) { m_value = c; return *this; }
 108    RR& operator=(signed char c) { m_value = c; return *this; }
 109 #ifndef BOOST_NO_INTRINSIC_WCHAR_T
 110    RR& operator=(wchar_t c) { m_value = c; return *this; }
 111 #endif
 112    RR& operator=(short c) { m_value = c; return *this; }
 113    RR& operator=(unsigned short c) { m_value = c; return *this; }
 114    RR& operator=(int c) { assign_large_int(c); return *this; }
 115    RR& operator=(unsigned int c) { assign_large_int(c); return *this; }
 116    RR& operator=(long c) { assign_large_int(c); return *this; }
 117    RR& operator=(unsigned long c) { assign_large_int(c); return *this; }
 118 #ifdef BOOST_HAS_LONG_LONG
 119    RR& operator=(boost::long_long_type c) { assign_large_int(c); return *this; }
 120    RR& operator=(boost::ulong_long_type c) { assign_large_int(c); return *this; }
 121 #endif
 122    RR& operator=(float c) { m_value = c; return *this; }
 123    RR& operator=(double c) { m_value = c; return *this; }
 124    RR& operator=(long double c) { assign_large_real(c); return *this; }
 125
 126    // Access:
 127    NTL::RR& value(){ return m_value; }
 128    NTL::RR const& value()const{ return m_value; }
 129
 130    // Member arithmetic:
 131    RR& operator+=(const RR& other)
 132    { m_value += other.value(); return *this; }
 133    RR& operator-=(const RR& other)
 134    { m_value -= other.value(); return *this; }
 135    RR& operator*=(const RR& other)
 136    { m_value *= other.value(); return *this; }
 137    RR& operator/=(const RR& other)
 138    { m_value /= other.value(); return *this; }
 139    RR operator-()const
 140    { return -m_value; }
 141    RR const& operator+()const
 142    { return *this; }
 143
 144    // RR compatibity:
 145    const ::NTL::ZZ& mantissa() const
 146    { return m_value.mantissa(); }
 147    long exponent() const
 148    { return m_value.exponent(); }
 149
 150    static void SetPrecision(long p)
 151    { ::NTL::RR::SetPrecision(p); }
 152
 153    static long precision()
 154    { return ::NTL::RR::precision(); }
 155
 156    static void SetOutputPrecision(long p)
 157    { ::NTL::RR::SetOutputPrecision(p); }
 158    static long OutputPrecision()
 159    { return ::NTL::RR::OutputPrecision(); }
 160
 161
 162 private:
 163    ::NTL::RR m_value;
 164
 165    template <class V>
 166    void assign_large_real(const V& a)
 167    {
 168       using std::frexp;
 169       using std::ldexp;
 170       using std::floor;
 171       if (a == 0) {
 172          clear(m_value);
 173          return;
 174       }
 175
 176       if (a == 1) {
 177          NTL::set(m_value);
 178          return;
 179       }
 180
 181       if (!(boost::math::isfinite)(a))
 182       {
 183          throw std::overflow_error("Cannot construct an instance of NTL::RR with an infinite value.");
 184       }
 185
 186       int e;
 187       long double f, term;
 188       ::NTL::RR t;
 189       clear(m_value);
 190
 191       f = frexp(a, &e);
 192
 193       while(f)
 194       {
 195          // extract 30 bits from f:
 196          f = ldexp(f, 30);
 197          term = floor(f);
 198          e -= 30;
 199          conv(t.x, (int)term);
 200          t.e = e;
 201          m_value += t;
 202          f -= term;
 203       }
 204    }
 205
 206    template <class V>
 207    void assign_large_int(V a)
 208    {
 209 #ifdef BOOST_MSVC
 210 #pragma warning(push)
 211 #pragma warning(disable:4146)
 212 #endif
 213       clear(m_value);
 214       int exp = 0;
 215       NTL::RR t;
 216       bool neg = a < V(0) ? true : false;
 217       if(neg)
 218          a = -a;
 219       while(a)
 220       {
 221          t = static_cast<double>(a & 0xffff);
 222          m_value += ldexp(RR(t), exp).value();
 223          a >>= 16;
 224          exp += 16;
 225       }
 226       if(neg)
 227          m_value = -m_value;
 228 #ifdef BOOST_MSVC
 229 #pragma warning(pop)
 230 #endif
 231    }
 232 };
 233
 234 // Non-member arithmetic:
 235 inline RR operator+(const RR& a, const RR& b)
 236 {
 237    RR result(a);
 238    result += b;
 239    return result;
 240 }
 241 inline RR operator-(const RR& a, const RR& b)
 242 {
 243    RR result(a);
 244    result -= b;
 245    return result;
 246 }
 247 inline RR operator*(const RR& a, const RR& b)
 248 {
 249    RR result(a);
 250    result *= b;
 251    return result;
 252 }
 253 inline RR operator/(const RR& a, const RR& b)
 254 {
 255    RR result(a);
 256    result /= b;
 257    return result;
 258 }
 259
 260 // Comparison:
 261 inline bool operator == (const RR& a, const RR& b)
 262 { return a.value() == b.value() ? true : false; }
 263 inline bool operator != (const RR& a, const RR& b)
 264 { return a.value() != b.value() ? true : false;}
 265 inline bool operator < (const RR& a, const RR& b)
 266 { return a.value() < b.value() ? true : false; }
 267 inline bool operator <= (const RR& a, const RR& b)
 268 { return a.value() <= b.value() ? true : false; }
 269 inline bool operator > (const RR& a, const RR& b)
 270 { return a.value() > b.value() ? true : false; }
 271 inline bool operator >= (const RR& a, const RR& b)
 272 { return a.value() >= b.value() ? true : false; }
 273
 274 #if 0
 275 // Non-member mixed compare:
 276 template <class T>
 277 inline bool operator == (const T& a, const RR& b)
 278 {
 279    return a == b.value();
 280 }
 281 template <class T>
 282 inline bool operator != (const T& a, const RR& b)
 283 {
 284    return a != b.value();
 285 }
 286 template <class T>
 287 inline bool operator < (const T& a, const RR& b)
 288 {
 289    return a < b.value();
 290 }
 291 template <class T>
 292 inline bool operator > (const T& a, const RR& b)
 293 {
 294    return a > b.value();
 295 }
 296 template <class T>
 297 inline bool operator <= (const T& a, const RR& b)
 298 {
 299    return a <= b.value();
 300 }
 301 template <class T>
 302 inline bool operator >= (const T& a, const RR& b)
 303 {
 304    return a >= b.value();
 305 }
 306 #endif  // Non-member mixed compare:
 307
 308 // Non-member functions:
 309 /*
 310 inline RR acos(RR a)
 311 { return ::NTL::acos(a.value()); }
 312 */
 313 inline RR cos(RR a)
 314 { return ::NTL::cos(a.value()); }
 315 /*
 316 inline RR asin(RR a)
 317 { return ::NTL::asin(a.value()); }
 318 inline RR atan(RR a)
 319 { return ::NTL::atan(a.value()); }
 320 inline RR atan2(RR a, RR b)
 321 { return ::NTL::atan2(a.value(), b.value()); }
 322 */
 323 inline RR ceil(RR a)
 324 { return ::NTL::ceil(a.value()); }
 325 /*
 326 inline RR fmod(RR a, RR b)
 327 { return ::NTL::fmod(a.value(), b.value()); }
 328 inline RR cosh(RR a)
 329 { return ::NTL::cosh(a.value()); }
 330 */
 331 inline RR exp(RR a)
 332 { return ::NTL::exp(a.value()); }
 333 inline RR fabs(RR a)
 334 { return ::NTL::fabs(a.value()); }
 335 inline RR abs(RR a)
 336 { return ::NTL::abs(a.value()); }
 337 inline RR floor(RR a)
 338 { return ::NTL::floor(a.value()); }
 339 /*
 340 inline RR modf(RR a, RR* ipart)
 341 {
 342    ::NTL::RR ip;
 343    RR result = modf(a.value(), &ip);
 344    *ipart = ip;
 345    return result;
 346 }
 347 inline RR frexp(RR a, int* expon)
 348 { return ::NTL::frexp(a.value(), expon); }
 349 inline RR ldexp(RR a, int expon)
 350 { return ::NTL::ldexp(a.value(), expon); }
 351 */
 352 inline RR log(RR a)
 353 { return ::NTL::log(a.value()); }
 354 inline RR log10(RR a)
 355 { return ::NTL::log10(a.value()); }
 356 /*
 357 inline RR tan(RR a)
 358 { return ::NTL::tan(a.value()); }
 359 */
 360 inline RR pow(RR a, RR b)
 361 { return ::NTL::pow(a.value(), b.value()); }
 362 inline RR pow(RR a, int b)
 363 { return ::NTL::power(a.value(), b); }
 364 inline RR sin(RR a)
 365 { return ::NTL::sin(a.value()); }
 366 /*
 367 inline RR sinh(RR a)
 368 { return ::NTL::sinh(a.value()); }
 369 */
 370 inline RR sqrt(RR a)
 371 { return ::NTL::sqrt(a.value()); }
 372 /*
 373 inline RR tanh(RR a)
 374 { return ::NTL::tanh(a.value()); }
 375 */
 376    inline RR pow(const RR& r, long l)
 377    {
 378       return ::NTL::power(r.value(), l);
 379    }
 380    inline RR tan(const RR& a)
 381    {
 382       return sin(a)/cos(a);
 383    }
 384    inline RR frexp(RR r, int* exp)
 385    {
 386       *exp = r.value().e;
 387       r.value().e = 0;
 388       while(r >= 1)
 389       {
 390          *exp += 1;
 391          r.value().e -= 1;
 392       }
 393       while(r < 0.5)
 394       {
 395          *exp -= 1;
 396          r.value().e += 1;
 397       }
 398       BOOST_ASSERT(r < 1);
 399       BOOST_ASSERT(r >= 0.5);
 400       return r;
 401    }
 402    inline RR ldexp(RR r, int exp)
 403    {
 404       r.value().e += exp;
 405       return r;
 406    }
 407
 408 // Streaming:
 409 template <class charT, class traits>
 410 inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& os, const RR& a)
 411 {
 412    return os << a.value();
 413 }
 414 template <class charT, class traits>
 415 inline std::basic_istream<charT, traits>& operator>>(std::basic_istream<charT, traits>& is, RR& a)
 416 {
 417    ::NTL::RR v;
 418    is >> v;
 419    a = v;
 420    return is;
 421 }
 422
 423 } // namespace ntl
 424
 425 namespace lanczos{
 426
 427 struct ntl_lanczos
 428 {
 429    static ntl::RR lanczos_sum(const ntl::RR& z)
 430    {
 431       unsigned long p = ntl::RR::precision();
 432       if(p <= 72)
 433          return lanczos13UDT::lanczos_sum(z);
 434       else if(p <= 120)
 435          return lanczos22UDT::lanczos_sum(z);
 436       else if(p <= 170)
 437          return lanczos31UDT::lanczos_sum(z);
 438       else //if(p <= 370) approx 100 digit precision:
 439          return lanczos61UDT::lanczos_sum(z);
 440    }
 441    static ntl::RR lanczos_sum_expG_scaled(const ntl::RR& z)
 442    {
 443       unsigned long p = ntl::RR::precision();
 444       if(p <= 72)
 445          return lanczos13UDT::lanczos_sum_expG_scaled(z);
 446       else if(p <= 120)
 447          return lanczos22UDT::lanczos_sum_expG_scaled(z);
 448       else if(p <= 170)
 449          return lanczos31UDT::lanczos_sum_expG_scaled(z);
 450       else //if(p <= 370) approx 100 digit precision:
 451          return lanczos61UDT::lanczos_sum_expG_scaled(z);
 452    }
 453    static ntl::RR lanczos_sum_near_1(const ntl::RR& z)
 454    {
 455       unsigned long p = ntl::RR::precision();
 456       if(p <= 72)
 457          return lanczos13UDT::lanczos_sum_near_1(z);
 458       else if(p <= 120)
 459          return lanczos22UDT::lanczos_sum_near_1(z);
 460       else if(p <= 170)
 461          return lanczos31UDT::lanczos_sum_near_1(z);
 462       else //if(p <= 370) approx 100 digit precision:
 463          return lanczos61UDT::lanczos_sum_near_1(z);
 464    }
 465    static ntl::RR lanczos_sum_near_2(const ntl::RR& z)
 466    {
 467       unsigned long p = ntl::RR::precision();
 468       if(p <= 72)
 469          return lanczos13UDT::lanczos_sum_near_2(z);
 470       else if(p <= 120)
 471          return lanczos22UDT::lanczos_sum_near_2(z);
 472       else if(p <= 170)
 473          return lanczos31UDT::lanczos_sum_near_2(z);
 474       else //if(p <= 370) approx 100 digit precision:
 475          return lanczos61UDT::lanczos_sum_near_2(z);
 476    }
 477    static ntl::RR g()
 478    {
 479       unsigned long p = ntl::RR::precision();
 480       if(p <= 72)
 481          return lanczos13UDT::g();
 482       else if(p <= 120)
 483          return lanczos22UDT::g();
 484       else if(p <= 170)
 485          return lanczos31UDT::g();
 486       else //if(p <= 370) approx 100 digit precision:
 487          return lanczos61UDT::g();
 488    }
 489 };
 490
 491 template<class Policy>
 492 struct lanczos<ntl::RR, Policy>
 493 {
 494    typedef ntl_lanczos type;
 495 };
 496
 497 } // namespace lanczos
 498
 499 namespace tools
 500 {
 501
 502 template<>
 503 inline int digits<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
 504 {
 505    return ::NTL::RR::precision();
 506 }
 507
 508 template <>
 509 inline float real_cast<float, boost::math::ntl::RR>(boost::math::ntl::RR t)
 510 {
 511    double r;
 512    conv(r, t.value());
 513    return static_cast<float>(r);
 514 }
 515 template <>
 516 inline double real_cast<double, boost::math::ntl::RR>(boost::math::ntl::RR t)
 517 {
 518    double r;
 519    conv(r, t.value());
 520    return r;
 521 }
 522
 523 namespace detail{
 524
 525 template<class I>
 526 void convert_to_long_result(NTL::RR const& r, I& result)
 527 {
 528    result = 0;
 529    I last_result(0);
 530    NTL::RR t(r);
 531    double term;
 532    do
 533    {
 534       conv(term, t);
 535       last_result = result;
 536       result += static_cast<I>(term);
 537       t -= term;
 538    }while(result != last_result);
 539 }
 540
 541 }
 542
 543 template <>
 544 inline long double real_cast<long double, boost::math::ntl::RR>(boost::math::ntl::RR t)
 545 {
 546    long double result(0);
 547    detail::convert_to_long_result(t.value(), result);
 548    return result;
 549 }
 550 template <>
 551 inline boost::math::ntl::RR real_cast<boost::math::ntl::RR, boost::math::ntl::RR>(boost::math::ntl::RR t)
 552 {
 553    return t;
 554 }
 555 template <>
 556 inline unsigned real_cast<unsigned, boost::math::ntl::RR>(boost::math::ntl::RR t)
 557 {
 558    unsigned result;
 559    detail::convert_to_long_result(t.value(), result);
 560    return result;
 561 }
 562 template <>
 563 inline int real_cast<int, boost::math::ntl::RR>(boost::math::ntl::RR t)
 564 {
 565    int result;
 566    detail::convert_to_long_result(t.value(), result);
 567    return result;
 568 }
 569 template <>
 570 inline long real_cast<long, boost::math::ntl::RR>(boost::math::ntl::RR t)
 571 {
 572    long result;
 573    detail::convert_to_long_result(t.value(), result);
 574    return result;
 575 }
 576 template <>
 577 inline long long real_cast<long long, boost::math::ntl::RR>(boost::math::ntl::RR t)
 578 {
 579    long long result;
 580    detail::convert_to_long_result(t.value(), result);
 581    return result;
 582 }
 583
 584 template <>
 585 inline boost::math::ntl::RR max_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
 586 {
 587    static bool has_init = false;
 588    static NTL::RR val;
 589    if(!has_init)
 590    {
 591       val = 1;
 592       val.e = NTL_OVFBND-20;
 593       has_init = true;
 594    }
 595    return val;
 596 }
 597
 598 template <>
 599 inline boost::math::ntl::RR min_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
 600 {
 601    static bool has_init = false;
 602    static NTL::RR val;
 603    if(!has_init)
 604    {
 605       val = 1;
 606       val.e = -NTL_OVFBND+20;
 607       has_init = true;
 608    }
 609    return val;
 610 }
 611
 612 template <>
 613 inline boost::math::ntl::RR log_max_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
 614 {
 615    static bool has_init = false;
 616    static NTL::RR val;
 617    if(!has_init)
 618    {
 619       val = 1;
 620       val.e = NTL_OVFBND-20;
 621       val = log(val);
 622       has_init = true;
 623    }
 624    return val;
 625 }
 626
 627 template <>
 628 inline boost::math::ntl::RR log_min_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
 629 {
 630    static bool has_init = false;
 631    static NTL::RR val;
 632    if(!has_init)
 633    {
 634       val = 1;
 635       val.e = -NTL_OVFBND+20;
 636       val = log(val);
 637       has_init = true;
 638    }
 639    return val;
 640 }
 641
 642 template <>
 643 inline boost::math::ntl::RR epsilon<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
 644 {
 645    return ldexp(boost::math::ntl::RR(1), 1-boost::math::policies::digits<boost::math::ntl::RR, boost::math::policies::policy<> >());
 646 }
 647
 648 } // namespace tools
 649
 650 //
 651 // The number of digits precision in RR can vary with each call
 652 // so we need to recalculate these with each call:
 653 //
 654 namespace constants{
 655
 656 template<> inline boost::math::ntl::RR pi<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
 657 {
 658     NTL::RR result;
 659     ComputePi(result);
 660     return result;
 661 }
 662 template<> inline boost::math::ntl::RR e<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
 663 {
 664     NTL::RR result;
 665     result = 1;
 666     return exp(result);
 667 }
 668
 669 } // namespace constants
 670
 671 namespace ntl{
 672    //
 673    // These are some fairly brain-dead versions of the math
 674    // functions that NTL fails to provide.
 675    //
 676
 677
 678    //
 679    // Inverse trig functions:
 680    //
 681    struct asin_root
 682    {
 683       asin_root(RR const& target) : t(target){}
 684
 685       boost::math::tuple<RR, RR, RR> operator()(RR const& p)
 686       {
 687          RR f0 = sin(p);
 688          RR f1 = cos(p);
 689          RR f2 = -f0;
 690          f0 -= t;
 691          return boost::math::make_tuple(f0, f1, f2);
 692       }
 693    private:
 694       RR t;
 695    };
 696
 697    inline RR asin(RR z)
 698    {
 699       double r;
 700       conv(r, z.value());
 701       return boost::math::tools::halley_iterate(
 702          asin_root(z),
 703          RR(std::asin(r)),
 704          RR(-boost::math::constants::pi<RR>()/2),
 705          RR(boost::math::constants::pi<RR>()/2),
 706          NTL::RR::precision());
 707    }
 708
 709    struct acos_root
 710    {
 711       acos_root(RR const& target) : t(target){}
 712
 713       boost::math::tuple<RR, RR, RR> operator()(RR const& p)
 714       {
 715          RR f0 = cos(p);
 716          RR f1 = -sin(p);
 717          RR f2 = -f0;
 718          f0 -= t;
 719          return boost::math::make_tuple(f0, f1, f2);
 720       }
 721    private:
 722       RR t;
 723    };
 724
 725    inline RR acos(RR z)
 726    {
 727       double r;
 728       conv(r, z.value());
 729       return boost::math::tools::halley_iterate(
 730          acos_root(z),
 731          RR(std::acos(r)),
 732          RR(-boost::math::constants::pi<RR>()/2),
 733          RR(boost::math::constants::pi<RR>()/2),
 734          NTL::RR::precision());
 735    }
 736
 737    struct atan_root
 738    {
 739       atan_root(RR const& target) : t(target){}
 740
 741       boost::math::tuple<RR, RR, RR> operator()(RR const& p)
 742       {
 743          RR c = cos(p);
 744          RR ta = tan(p);
 745          RR f0 = ta - t;
 746          RR f1 = 1 / (c * c);
 747          RR f2 = 2 * ta / (c * c);
 748          return boost::math::make_tuple(f0, f1, f2);
 749       }
 750    private:
 751       RR t;
 752    };
 753
 754    inline RR atan(RR z)
 755    {
 756       double r;
 757       conv(r, z.value());
 758       return boost::math::tools::halley_iterate(
 759          atan_root(z),
 760          RR(std::atan(r)),
 761          -boost::math::constants::pi<RR>()/2,
 762          boost::math::constants::pi<RR>()/2,
 763          NTL::RR::precision());
 764    }
 765
 766    inline RR atan2(RR y, RR x)
 767    {
 768       if(x > 0)
 769          return atan(y / x);
 770       if(x < 0)
 771       {
 772          return y < 0 ? atan(y / x) - boost::math::constants::pi<RR>() : atan(y / x) + boost::math::constants::pi<RR>();
 773       }
 774       return y < 0 ? -boost::math::constants::half_pi<RR>() : boost::math::constants::half_pi<RR>() ;
 775    }
 776
 777    inline RR sinh(RR z)
 778    {
 779       return (expm1(z.value()) - expm1(-z.value())) / 2;
 780    }
 781
 782    inline RR cosh(RR z)
 783    {
 784       return (exp(z) + exp(-z)) / 2;
 785    }
 786
 787    inline RR tanh(RR z)
 788    {
 789       return sinh(z) / cosh(z);
 790    }
 791
 792    inline RR fmod(RR x, RR y)
 793    {
 794       // This is a really crummy version of fmod, we rely on lots
 795       // of digits to get us out of trouble...
 796       RR factor = floor(x/y);
 797       return x - factor * y;
 798    }
 799
 800    template <class Policy>
 801    inline int iround(RR const& x, const Policy& pol)
 802    {
 803       return tools::real_cast<int>(round(x, pol));
 804    }
 805
 806    template <class Policy>
 807    inline long lround(RR const& x, const Policy& pol)
 808    {
 809       return tools::real_cast<long>(round(x, pol));
 810    }
 811
 812    template <class Policy>
 813    inline long long llround(RR const& x, const Policy& pol)
 814    {
 815       return tools::real_cast<long long>(round(x, pol));
 816    }
 817
 818    template <class Policy>
 819    inline int itrunc(RR const& x, const Policy& pol)
 820    {
 821       return tools::real_cast<int>(trunc(x, pol));
 822    }
 823
 824    template <class Policy>
 825    inline long ltrunc(RR const& x, const Policy& pol)
 826    {
 827       return tools::real_cast<long>(trunc(x, pol));
 828    }
 829
 830    template <class Policy>
 831    inline long long lltrunc(RR const& x, const Policy& pol)
 832    {
 833       return tools::real_cast<long long>(trunc(x, pol));
 834    }
 835
 836 } // namespace ntl
 837
 838 namespace detail{
 839
 840 template <class Policy>
 841 ntl::RR digamma_imp(ntl::RR x, const mpl::int_<0>* , const Policy& pol)
 842 {
 843    //
 844    // This handles reflection of negative arguments, and all our
 845    // error handling, then forwards to the T-specific approximation.
 846    //
 847    BOOST_MATH_STD_USING // ADL of std functions.
 848
 849    ntl::RR result = 0;
 850    //
 851    // Check for negative arguments and use reflection:
 852    //
 853    if(x < 0)
 854    {
 855       // Reflect:
 856       x = 1 - x;
 857       // Argument reduction for tan:
 858       ntl::RR remainder = x - floor(x);
 859       // Shift to negative if > 0.5:
 860       if(remainder > 0.5)
 861       {
 862          remainder -= 1;
 863       }
 864       //
 865       // check for evaluation at a negative pole:
 866       //
 867       if(remainder == 0)
 868       {
 869          return policies::raise_pole_error<ntl::RR>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol);
 870       }
 871       result = constants::pi<ntl::RR>() / tan(constants::pi<ntl::RR>() * remainder);
 872    }
 873    result += big_digamma(x);
 874    return result;
 875 }
 876
 877 } // namespace detail
 878
 879 } // namespace math
 880 } // namespace boost
 881
 882 #endif // BOOST_MATH_REAL_CONCEPT_HPP
 883
 884