Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:39:31 

0001 //  Copyright John Maddock 2007.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005 
0006 #ifndef BOOST_MATH_NTL_RR_HPP
0007 #define BOOST_MATH_NTL_RR_HPP
0008 
0009 #include <boost/math/tools/real_cast.hpp>
0010 #include <boost/math/tools/precision.hpp>
0011 #include <boost/math/tools/config.hpp>
0012 #include <boost/math/constants/constants.hpp>
0013 #include <boost/math/tools/roots.hpp>
0014 #include <boost/math/special_functions/fpclassify.hpp>
0015 #include <boost/math/bindings/detail/big_digamma.hpp>
0016 #include <boost/math/bindings/detail/big_lanczos.hpp>
0017 #include <stdexcept>
0018 #include <ostream>
0019 #include <istream>
0020 #include <cmath>
0021 #include <limits>
0022 #include <NTL/RR.h>
0023 
0024 namespace boost{ namespace math{
0025 
0026 namespace ntl
0027 {
0028 
0029 class RR;
0030 
0031 RR ldexp(RR r, int exp);
0032 RR frexp(RR r, int* exp);
0033 
0034 class RR
0035 {
0036 public:
0037    // Constructors:
0038    RR() {}
0039    RR(const ::NTL::RR& c) : m_value(c){}
0040    RR(char c)
0041    {
0042       m_value = c;
0043    }
0044    RR(wchar_t c)
0045    {
0046       m_value = c;
0047    }
0048    RR(unsigned char c)
0049    {
0050       m_value = c;
0051    }
0052    RR(signed char c)
0053    {
0054       m_value = c;
0055    }
0056    RR(unsigned short c)
0057    {
0058       m_value = c;
0059    }
0060    RR(short c)
0061    {
0062       m_value = c;
0063    }
0064    RR(unsigned int c)
0065    {
0066       assign_large_int(c);
0067    }
0068    RR(int c)
0069    {
0070       assign_large_int(c);
0071    }
0072    RR(unsigned long c)
0073    {
0074       assign_large_int(c);
0075    }
0076    RR(long c)
0077    {
0078       assign_large_int(c);
0079    }
0080    RR(unsigned long long c)
0081    {
0082       assign_large_int(c);
0083    }
0084    RR(long long c)
0085    {
0086       assign_large_int(c);
0087    }
0088    RR(float c)
0089    {
0090       m_value = c;
0091    }
0092    RR(double c)
0093    {
0094       m_value = c;
0095    }
0096    RR(long double c)
0097    {
0098       assign_large_real(c);
0099    }
0100 
0101    // Assignment:
0102    RR& operator=(char c) { m_value = c; return *this; }
0103    RR& operator=(unsigned char c) { m_value = c; return *this; }
0104    RR& operator=(signed char c) { m_value = c; return *this; }
0105    RR& operator=(wchar_t c) { m_value = c; return *this; }
0106    RR& operator=(short c) { m_value = c; return *this; }
0107    RR& operator=(unsigned short c) { m_value = c; return *this; }
0108    RR& operator=(int c) { assign_large_int(c); return *this; }
0109    RR& operator=(unsigned int c) { assign_large_int(c); return *this; }
0110    RR& operator=(long c) { assign_large_int(c); return *this; }
0111    RR& operator=(unsigned long c) { assign_large_int(c); return *this; }
0112    RR& operator=(long long c) { assign_large_int(c); return *this; }
0113    RR& operator=(unsigned long long c) { assign_large_int(c); return *this; }
0114    RR& operator=(float c) { m_value = c; return *this; }
0115    RR& operator=(double c) { m_value = c; return *this; }
0116    RR& operator=(long double c) { assign_large_real(c); return *this; }
0117 
0118    // Access:
0119    NTL::RR& value(){ return m_value; }
0120    NTL::RR const& value()const{ return m_value; }
0121 
0122    // Member arithmetic:
0123    RR& operator+=(const RR& other)
0124    { m_value += other.value(); return *this; }
0125    RR& operator-=(const RR& other)
0126    { m_value -= other.value(); return *this; }
0127    RR& operator*=(const RR& other)
0128    { m_value *= other.value(); return *this; }
0129    RR& operator/=(const RR& other)
0130    { m_value /= other.value(); return *this; }
0131    RR operator-()const
0132    { return -m_value; }
0133    RR const& operator+()const
0134    { return *this; }
0135 
0136    // RR compatibility:
0137    const ::NTL::ZZ& mantissa() const
0138    { return m_value.mantissa(); }
0139    long exponent() const
0140    { return m_value.exponent(); }
0141 
0142    static void SetPrecision(long p)
0143    { ::NTL::RR::SetPrecision(p); }
0144 
0145    static long precision()
0146    { return ::NTL::RR::precision(); }
0147 
0148    static void SetOutputPrecision(long p)
0149    { ::NTL::RR::SetOutputPrecision(p); }
0150    static long OutputPrecision()
0151    { return ::NTL::RR::OutputPrecision(); }
0152 
0153 
0154 private:
0155    ::NTL::RR m_value;
0156 
0157    template <class V>
0158    void assign_large_real(const V& a)
0159    {
0160       using std::frexp;
0161       using std::ldexp;
0162       using std::floor;
0163       if (a == 0) {
0164          clear(m_value);
0165          return;
0166       }
0167 
0168       if (a == 1) {
0169          NTL::set(m_value);
0170          return;
0171       }
0172 
0173       if (!(boost::math::isfinite)(a))
0174       {
0175          throw std::overflow_error("Cannot construct an instance of NTL::RR with an infinite value.");
0176       }
0177 
0178       int e;
0179       long double f, term;
0180       ::NTL::RR t;
0181       clear(m_value);
0182 
0183       f = frexp(a, &e);
0184 
0185       while(f)
0186       {
0187          // extract 30 bits from f:
0188          f = ldexp(f, 30);
0189          term = floor(f);
0190          e -= 30;
0191          conv(t.x, (int)term);
0192          t.e = e;
0193          m_value += t;
0194          f -= term;
0195       }
0196    }
0197 
0198    template <class V>
0199    void assign_large_int(V a)
0200    {
0201 #ifdef _MSC_VER
0202 #pragma warning(push)
0203 #pragma warning(disable:4146)
0204 #endif
0205       clear(m_value);
0206       int exp = 0;
0207       NTL::RR t;
0208       bool neg = a < V(0) ? true : false;
0209       if(neg)
0210          a = -a;
0211       while(a)
0212       {
0213          t = static_cast<double>(a & 0xffff);
0214          m_value += ldexp(RR(t), exp).value();
0215          a >>= 16;
0216          exp += 16;
0217       }
0218       if(neg)
0219          m_value = -m_value;
0220 #ifdef _MSC_VER
0221 #pragma warning(pop)
0222 #endif
0223    }
0224 };
0225 
0226 // Non-member arithmetic:
0227 inline RR operator+(const RR& a, const RR& b)
0228 {
0229    RR result(a);
0230    result += b;
0231    return result;
0232 }
0233 inline RR operator-(const RR& a, const RR& b)
0234 {
0235    RR result(a);
0236    result -= b;
0237    return result;
0238 }
0239 inline RR operator*(const RR& a, const RR& b)
0240 {
0241    RR result(a);
0242    result *= b;
0243    return result;
0244 }
0245 inline RR operator/(const RR& a, const RR& b)
0246 {
0247    RR result(a);
0248    result /= b;
0249    return result;
0250 }
0251 
0252 // Comparison:
0253 inline bool operator == (const RR& a, const RR& b)
0254 { return a.value() == b.value() ? true : false; }
0255 inline bool operator != (const RR& a, const RR& b)
0256 { return a.value() != b.value() ? true : false;}
0257 inline bool operator < (const RR& a, const RR& b)
0258 { return a.value() < b.value() ? true : false; }
0259 inline bool operator <= (const RR& a, const RR& b)
0260 { return a.value() <= b.value() ? true : false; }
0261 inline bool operator > (const RR& a, const RR& b)
0262 { return a.value() > b.value() ? true : false; }
0263 inline bool operator >= (const RR& a, const RR& b)
0264 { return a.value() >= b.value() ? true : false; }
0265 
0266 #if 0
0267 // Non-member mixed compare:
0268 template <class T>
0269 inline bool operator == (const T& a, const RR& b)
0270 {
0271    return a == b.value();
0272 }
0273 template <class T>
0274 inline bool operator != (const T& a, const RR& b)
0275 {
0276    return a != b.value();
0277 }
0278 template <class T>
0279 inline bool operator < (const T& a, const RR& b)
0280 {
0281    return a < b.value();
0282 }
0283 template <class T>
0284 inline bool operator > (const T& a, const RR& b)
0285 {
0286    return a > b.value();
0287 }
0288 template <class T>
0289 inline bool operator <= (const T& a, const RR& b)
0290 {
0291    return a <= b.value();
0292 }
0293 template <class T>
0294 inline bool operator >= (const T& a, const RR& b)
0295 {
0296    return a >= b.value();
0297 }
0298 #endif  // Non-member mixed compare:
0299 
0300 // Non-member functions:
0301 /*
0302 inline RR acos(RR a)
0303 { return ::NTL::acos(a.value()); }
0304 */
0305 inline RR cos(RR a)
0306 { return ::NTL::cos(a.value()); }
0307 /*
0308 inline RR asin(RR a)
0309 { return ::NTL::asin(a.value()); }
0310 inline RR atan(RR a)
0311 { return ::NTL::atan(a.value()); }
0312 inline RR atan2(RR a, RR b)
0313 { return ::NTL::atan2(a.value(), b.value()); }
0314 */
0315 inline RR ceil(RR a)
0316 { return ::NTL::ceil(a.value()); }
0317 /*
0318 inline RR fmod(RR a, RR b)
0319 { return ::NTL::fmod(a.value(), b.value()); }
0320 inline RR cosh(RR a)
0321 { return ::NTL::cosh(a.value()); }
0322 */
0323 inline RR exp(RR a)
0324 { return ::NTL::exp(a.value()); }
0325 inline RR fabs(RR a)
0326 { return ::NTL::fabs(a.value()); }
0327 inline RR abs(RR a)
0328 { return ::NTL::abs(a.value()); }
0329 inline RR floor(RR a)
0330 { return ::NTL::floor(a.value()); }
0331 /*
0332 inline RR modf(RR a, RR* ipart)
0333 {
0334    ::NTL::RR ip;
0335    RR result = modf(a.value(), &ip);
0336    *ipart = ip;
0337    return result;
0338 }
0339 inline RR frexp(RR a, int* expon)
0340 { return ::NTL::frexp(a.value(), expon); }
0341 inline RR ldexp(RR a, int expon)
0342 { return ::NTL::ldexp(a.value(), expon); }
0343 */
0344 inline RR log(RR a)
0345 { return ::NTL::log(a.value()); }
0346 inline RR log10(RR a)
0347 { return ::NTL::log10(a.value()); }
0348 /*
0349 inline RR tan(RR a)
0350 { return ::NTL::tan(a.value()); }
0351 */
0352 inline RR pow(RR a, RR b)
0353 { return ::NTL::pow(a.value(), b.value()); }
0354 inline RR pow(RR a, int b)
0355 { return ::NTL::power(a.value(), b); }
0356 inline RR sin(RR a)
0357 { return ::NTL::sin(a.value()); }
0358 /*
0359 inline RR sinh(RR a)
0360 { return ::NTL::sinh(a.value()); }
0361 */
0362 inline RR sqrt(RR a)
0363 { return ::NTL::sqrt(a.value()); }
0364 /*
0365 inline RR tanh(RR a)
0366 { return ::NTL::tanh(a.value()); }
0367 */
0368    inline RR pow(const RR& r, long l)
0369    {
0370       return ::NTL::power(r.value(), l);
0371    }
0372    inline RR tan(const RR& a)
0373    {
0374       return sin(a)/cos(a);
0375    }
0376    inline RR frexp(RR r, int* exp)
0377    {
0378       *exp = r.value().e;
0379       r.value().e = 0;
0380       while(r >= 1)
0381       {
0382          *exp += 1;
0383          r.value().e -= 1;
0384       }
0385       while(r < 0.5)
0386       {
0387          *exp -= 1;
0388          r.value().e += 1;
0389       }
0390       BOOST_MATH_ASSERT(r < 1);
0391       BOOST_MATH_ASSERT(r >= 0.5);
0392       return r;
0393    }
0394    inline RR ldexp(RR r, int exp)
0395    {
0396       r.value().e += exp;
0397       return r;
0398    }
0399 
0400 // Streaming:
0401 template <class charT, class traits>
0402 inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& os, const RR& a)
0403 {
0404    return os << a.value();
0405 }
0406 template <class charT, class traits>
0407 inline std::basic_istream<charT, traits>& operator>>(std::basic_istream<charT, traits>& is, RR& a)
0408 {
0409    ::NTL::RR v;
0410    is >> v;
0411    a = v;
0412    return is;
0413 }
0414 
0415 } // namespace ntl
0416 
0417 namespace lanczos{
0418 
0419 struct ntl_lanczos
0420 {
0421    static ntl::RR lanczos_sum(const ntl::RR& z)
0422    {
0423       unsigned long p = ntl::RR::precision();
0424       if(p <= 72)
0425          return lanczos13UDT::lanczos_sum(z);
0426       else if(p <= 120)
0427          return lanczos22UDT::lanczos_sum(z);
0428       else if(p <= 170)
0429          return lanczos31UDT::lanczos_sum(z);
0430       else //if(p <= 370) approx 100 digit precision:
0431          return lanczos61UDT::lanczos_sum(z);
0432    }
0433    static ntl::RR lanczos_sum_expG_scaled(const ntl::RR& z)
0434    {
0435       unsigned long p = ntl::RR::precision();
0436       if(p <= 72)
0437          return lanczos13UDT::lanczos_sum_expG_scaled(z);
0438       else if(p <= 120)
0439          return lanczos22UDT::lanczos_sum_expG_scaled(z);
0440       else if(p <= 170)
0441          return lanczos31UDT::lanczos_sum_expG_scaled(z);
0442       else //if(p <= 370) approx 100 digit precision:
0443          return lanczos61UDT::lanczos_sum_expG_scaled(z);
0444    }
0445    static ntl::RR lanczos_sum_near_1(const ntl::RR& z)
0446    {
0447       unsigned long p = ntl::RR::precision();
0448       if(p <= 72)
0449          return lanczos13UDT::lanczos_sum_near_1(z);
0450       else if(p <= 120)
0451          return lanczos22UDT::lanczos_sum_near_1(z);
0452       else if(p <= 170)
0453          return lanczos31UDT::lanczos_sum_near_1(z);
0454       else //if(p <= 370) approx 100 digit precision:
0455          return lanczos61UDT::lanczos_sum_near_1(z);
0456    }
0457    static ntl::RR lanczos_sum_near_2(const ntl::RR& z)
0458    {
0459       unsigned long p = ntl::RR::precision();
0460       if(p <= 72)
0461          return lanczos13UDT::lanczos_sum_near_2(z);
0462       else if(p <= 120)
0463          return lanczos22UDT::lanczos_sum_near_2(z);
0464       else if(p <= 170)
0465          return lanczos31UDT::lanczos_sum_near_2(z);
0466       else //if(p <= 370) approx 100 digit precision:
0467          return lanczos61UDT::lanczos_sum_near_2(z);
0468    }
0469    static ntl::RR g()
0470    {
0471       unsigned long p = ntl::RR::precision();
0472       if(p <= 72)
0473          return lanczos13UDT::g();
0474       else if(p <= 120)
0475          return lanczos22UDT::g();
0476       else if(p <= 170)
0477          return lanczos31UDT::g();
0478       else //if(p <= 370) approx 100 digit precision:
0479          return lanczos61UDT::g();
0480    }
0481 };
0482 
0483 template<class Policy>
0484 struct lanczos<ntl::RR, Policy>
0485 {
0486    typedef ntl_lanczos type;
0487 };
0488 
0489 } // namespace lanczos
0490 
0491 namespace tools
0492 {
0493 
0494 template<>
0495 inline int digits<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
0496 {
0497    return ::NTL::RR::precision();
0498 }
0499 
0500 template <>
0501 inline float real_cast<float, boost::math::ntl::RR>(boost::math::ntl::RR t)
0502 {
0503    double r;
0504    conv(r, t.value());
0505    return static_cast<float>(r);
0506 }
0507 template <>
0508 inline double real_cast<double, boost::math::ntl::RR>(boost::math::ntl::RR t)
0509 {
0510    double r;
0511    conv(r, t.value());
0512    return r;
0513 }
0514 
0515 namespace detail{
0516 
0517 template<class I>
0518 void convert_to_long_result(NTL::RR const& r, I& result)
0519 {
0520    result = 0;
0521    I last_result(0);
0522    NTL::RR t(r);
0523    double term;
0524    do
0525    {
0526       conv(term, t);
0527       last_result = result;
0528       result += static_cast<I>(term);
0529       t -= term;
0530    }while(result != last_result);
0531 }
0532 
0533 }
0534 
0535 template <>
0536 inline long double real_cast<long double, boost::math::ntl::RR>(boost::math::ntl::RR t)
0537 {
0538    long double result(0);
0539    detail::convert_to_long_result(t.value(), result);
0540    return result;
0541 }
0542 template <>
0543 inline boost::math::ntl::RR real_cast<boost::math::ntl::RR, boost::math::ntl::RR>(boost::math::ntl::RR t)
0544 {
0545    return t;
0546 }
0547 template <>
0548 inline unsigned real_cast<unsigned, boost::math::ntl::RR>(boost::math::ntl::RR t)
0549 {
0550    unsigned result;
0551    detail::convert_to_long_result(t.value(), result);
0552    return result;
0553 }
0554 template <>
0555 inline int real_cast<int, boost::math::ntl::RR>(boost::math::ntl::RR t)
0556 {
0557    int result;
0558    detail::convert_to_long_result(t.value(), result);
0559    return result;
0560 }
0561 template <>
0562 inline long real_cast<long, boost::math::ntl::RR>(boost::math::ntl::RR t)
0563 {
0564    long result;
0565    detail::convert_to_long_result(t.value(), result);
0566    return result;
0567 }
0568 template <>
0569 inline long long real_cast<long long, boost::math::ntl::RR>(boost::math::ntl::RR t)
0570 {
0571    long long result;
0572    detail::convert_to_long_result(t.value(), result);
0573    return result;
0574 }
0575 
0576 template <>
0577 inline boost::math::ntl::RR max_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
0578 {
0579    static bool has_init = false;
0580    static NTL::RR val;
0581    if(!has_init)
0582    {
0583       val = 1;
0584       val.e = NTL_OVFBND-20;
0585       has_init = true;
0586    }
0587    return val;
0588 }
0589 
0590 template <>
0591 inline boost::math::ntl::RR min_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
0592 {
0593    static bool has_init = false;
0594    static NTL::RR val;
0595    if(!has_init)
0596    {
0597       val = 1;
0598       val.e = -NTL_OVFBND+20;
0599       has_init = true;
0600    }
0601    return val;
0602 }
0603 
0604 template <>
0605 inline boost::math::ntl::RR log_max_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
0606 {
0607    static bool has_init = false;
0608    static NTL::RR val;
0609    if(!has_init)
0610    {
0611       val = 1;
0612       val.e = NTL_OVFBND-20;
0613       val = log(val);
0614       has_init = true;
0615    }
0616    return val;
0617 }
0618 
0619 template <>
0620 inline boost::math::ntl::RR log_min_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
0621 {
0622    static bool has_init = false;
0623    static NTL::RR val;
0624    if(!has_init)
0625    {
0626       val = 1;
0627       val.e = -NTL_OVFBND+20;
0628       val = log(val);
0629       has_init = true;
0630    }
0631    return val;
0632 }
0633 
0634 template <>
0635 inline boost::math::ntl::RR epsilon<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
0636 {
0637    return ldexp(boost::math::ntl::RR(1), 1-boost::math::policies::digits<boost::math::ntl::RR, boost::math::policies::policy<> >());
0638 }
0639 
0640 } // namespace tools
0641 
0642 //
0643 // The number of digits precision in RR can vary with each call
0644 // so we need to recalculate these with each call:
0645 //
0646 namespace constants{
0647 
0648 template<> inline boost::math::ntl::RR pi<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
0649 {
0650     NTL::RR result;
0651     ComputePi(result);
0652     return result;
0653 }
0654 template<> inline boost::math::ntl::RR e<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR))
0655 {
0656     NTL::RR result;
0657     result = 1;
0658     return exp(result);
0659 }
0660 
0661 } // namespace constants
0662 
0663 namespace ntl{
0664    //
0665    // These are some fairly brain-dead versions of the math
0666    // functions that NTL fails to provide.
0667    //
0668 
0669 
0670    //
0671    // Inverse trig functions:
0672    //
0673    struct asin_root
0674    {
0675       asin_root(RR const& target) : t(target){}
0676 
0677       boost::math::tuple<RR, RR, RR> operator()(RR const& p)
0678       {
0679          RR f0 = sin(p);
0680          RR f1 = cos(p);
0681          RR f2 = -f0;
0682          f0 -= t;
0683          return boost::math::make_tuple(f0, f1, f2);
0684       }
0685    private:
0686       RR t;
0687    };
0688 
0689    inline RR asin(RR z)
0690    {
0691       double r;
0692       conv(r, z.value());
0693       return boost::math::tools::halley_iterate(
0694          asin_root(z),
0695          RR(std::asin(r)),
0696          RR(-boost::math::constants::pi<RR>()/2),
0697          RR(boost::math::constants::pi<RR>()/2),
0698          NTL::RR::precision());
0699    }
0700 
0701    struct acos_root
0702    {
0703       acos_root(RR const& target) : t(target){}
0704 
0705       boost::math::tuple<RR, RR, RR> operator()(RR const& p)
0706       {
0707          RR f0 = cos(p);
0708          RR f1 = -sin(p);
0709          RR f2 = -f0;
0710          f0 -= t;
0711          return boost::math::make_tuple(f0, f1, f2);
0712       }
0713    private:
0714       RR t;
0715    };
0716 
0717    inline RR acos(RR z)
0718    {
0719       double r;
0720       conv(r, z.value());
0721       return boost::math::tools::halley_iterate(
0722          acos_root(z),
0723          RR(std::acos(r)),
0724          RR(-boost::math::constants::pi<RR>()/2),
0725          RR(boost::math::constants::pi<RR>()/2),
0726          NTL::RR::precision());
0727    }
0728 
0729    struct atan_root
0730    {
0731       atan_root(RR const& target) : t(target){}
0732 
0733       boost::math::tuple<RR, RR, RR> operator()(RR const& p)
0734       {
0735          RR c = cos(p);
0736          RR ta = tan(p);
0737          RR f0 = ta - t;
0738          RR f1 = 1 / (c * c);
0739          RR f2 = 2 * ta / (c * c);
0740          return boost::math::make_tuple(f0, f1, f2);
0741       }
0742    private:
0743       RR t;
0744    };
0745 
0746    inline RR atan(RR z)
0747    {
0748       double r;
0749       conv(r, z.value());
0750       return boost::math::tools::halley_iterate(
0751          atan_root(z),
0752          RR(std::atan(r)),
0753          -boost::math::constants::pi<RR>()/2,
0754          boost::math::constants::pi<RR>()/2,
0755          NTL::RR::precision());
0756    }
0757 
0758    inline RR atan2(RR y, RR x)
0759    {
0760       if(x > 0)
0761          return atan(y / x);
0762       if(x < 0)
0763       {
0764          return y < 0 ? atan(y / x) - boost::math::constants::pi<RR>() : atan(y / x) + boost::math::constants::pi<RR>();
0765       }
0766       return y < 0 ? -boost::math::constants::half_pi<RR>() : boost::math::constants::half_pi<RR>() ;
0767    }
0768 
0769    inline RR sinh(RR z)
0770    {
0771       return (expm1(z.value()) - expm1(-z.value())) / 2;
0772    }
0773 
0774    inline RR cosh(RR z)
0775    {
0776       return (exp(z) + exp(-z)) / 2;
0777    }
0778 
0779    inline RR tanh(RR z)
0780    {
0781       return sinh(z) / cosh(z);
0782    }
0783 
0784    inline RR fmod(RR x, RR y)
0785    {
0786       // This is a really crummy version of fmod, we rely on lots
0787       // of digits to get us out of trouble...
0788       RR factor = floor(x/y);
0789       return x - factor * y;
0790    }
0791 
0792    template <class Policy>
0793    inline int iround(RR const& x, const Policy& pol)
0794    {
0795       return tools::real_cast<int>(round(x, pol));
0796    }
0797 
0798    template <class Policy>
0799    inline long lround(RR const& x, const Policy& pol)
0800    {
0801       return tools::real_cast<long>(round(x, pol));
0802    }
0803 
0804    template <class Policy>
0805    inline long long llround(RR const& x, const Policy& pol)
0806    {
0807       return tools::real_cast<long long>(round(x, pol));
0808    }
0809 
0810    template <class Policy>
0811    inline int itrunc(RR const& x, const Policy& pol)
0812    {
0813       return tools::real_cast<int>(trunc(x, pol));
0814    }
0815 
0816    template <class Policy>
0817    inline long ltrunc(RR const& x, const Policy& pol)
0818    {
0819       return tools::real_cast<long>(trunc(x, pol));
0820    }
0821 
0822    template <class Policy>
0823    inline long long lltrunc(RR const& x, const Policy& pol)
0824    {
0825       return tools::real_cast<long long>(trunc(x, pol));
0826    }
0827 
0828 } // namespace ntl
0829 
0830 namespace detail{
0831 
0832 template <class Policy>
0833 ntl::RR digamma_imp(ntl::RR x, const std::integral_constant<int, 0>* , const Policy& pol)
0834 {
0835    //
0836    // This handles reflection of negative arguments, and all our
0837    // error handling, then forwards to the T-specific approximation.
0838    //
0839    BOOST_MATH_STD_USING // ADL of std functions.
0840 
0841    ntl::RR result = 0;
0842    //
0843    // Check for negative arguments and use reflection:
0844    //
0845    if(x < 0)
0846    {
0847       // Reflect:
0848       x = 1 - x;
0849       // Argument reduction for tan:
0850       ntl::RR remainder = x - floor(x);
0851       // Shift to negative if > 0.5:
0852       if(remainder > 0.5)
0853       {
0854          remainder -= 1;
0855       }
0856       //
0857       // check for evaluation at a negative pole:
0858       //
0859       if(remainder == 0)
0860       {
0861          return policies::raise_pole_error<ntl::RR>("boost::math::digamma<%1%>(%1%)", nullptr, (1-x), pol);
0862       }
0863       result = constants::pi<ntl::RR>() / tan(constants::pi<ntl::RR>() * remainder);
0864    }
0865    result += big_digamma(x);
0866    return result;
0867 }
0868 
0869 } // namespace detail
0870 
0871 } // namespace math
0872 } // namespace boost
0873 
0874 #endif // BOOST_MATH_REAL_CONCEPT_HPP
0875 
0876
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!