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