update mpreal.h

This commit is contained in:
Gael Guennebaud 2015-10-13 09:58:54 +02:00
parent ea9749fd6c
commit 3e32f6b554

View File

@ -5,14 +5,15 @@
Project homepage: http://www.holoborodko.com/pavel/mpfr
Contact e-mail: pavel@holoborodko.com
Copyright (c) 2008-2014 Pavel Holoborodko
Copyright (c) 2008-2015 Pavel Holoborodko
Contributors:
Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman,
Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen,
Pere Constans, Peter van Hoof, Gael Guennebaud, Tsai Chia Cheng,
Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood,
Petr Aleksandrov, Orion Poplawski, Charles Karney.
Petr Aleksandrov, Orion Poplawski, Charles Karney, Arash Partow,
Rodney James, Jorge Leitao.
Licensing:
(A) MPFR C++ is under GNU General Public License ("GPL").
@ -55,10 +56,11 @@
#include <cmath>
#include <cstring>
#include <limits>
#include <cstdint>
#include <complex>
#include <algorithm>
// Options
// FIXME HAVE_INT64_SUPPORT leads to clashes with long int and int64_t on some systems.
//#define MPREAL_HAVE_INT64_SUPPORT // Enable int64_t support if possible. Available only for MSVC 2010 & GCC.
#define MPREAL_HAVE_MSVC_DEBUGVIEW // Enable Debugger Visualizer for "Debug" builds in MSVC.
#define MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS // Enable extended std::numeric_limits<mpfr::mpreal> specialization.
// Meaning that "digits", "round_style" and similar members are defined as functions, not constants.
@ -66,19 +68,19 @@
// Library version
#define MPREAL_VERSION_MAJOR 3
#define MPREAL_VERSION_MINOR 5
#define MPREAL_VERSION_PATCHLEVEL 9
#define MPREAL_VERSION_STRING "3.5.9"
#define MPREAL_VERSION_MINOR 6
#define MPREAL_VERSION_PATCHLEVEL 2
#define MPREAL_VERSION_STRING "3.6.2"
// Detect compiler using signatures from http://predef.sourceforge.net/
#if defined(__GNUC__) && defined(__INTEL_COMPILER)
#define IsInf(x) (isinf)(x) // Intel ICC compiler on Linux
#define IsInf(x) isinf(x) // Intel ICC compiler on Linux
#elif defined(_MSC_VER) // Microsoft Visual C++
#define IsInf(x) (!_finite(x))
#else
#define IsInf(x) (std::isinf)(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance
#define IsInf(x) std::isinf(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance
#endif
// A Clang feature extension to determine compiler features.
@ -100,40 +102,13 @@
// Detect support for explicit converters.
#if (__has_feature(cxx_explicit_conversions) || \
defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \
(defined(__GXX_EXPERIMENTAL_CXX0X__) && __GNUC_MINOR >= 5) || __cplusplus >= 201103L || \
(defined(_MSC_VER) && _MSC_VER >= 1800))
#define MPREAL_HAVE_EXPLICIT_CONVERTERS
#endif
// Detect available 64-bit capabilities
#if defined(MPREAL_HAVE_INT64_SUPPORT)
#define MPFR_USE_INTMAX_T // Should be defined before mpfr.h
#if defined(_MSC_VER) // MSVC + Windows
#if (_MSC_VER >= 1600)
#include <stdint.h> // <stdint.h> is available only in msvc2010!
#else // MPFR relies on intmax_t which is available only in msvc2010
#undef MPREAL_HAVE_INT64_SUPPORT // Besides, MPFR & MPIR have to be compiled with msvc2010
#undef MPFR_USE_INTMAX_T // Since we cannot detect this, disable x64 by default
// Someone should change this manually if needed.
#endif
#elif defined (__GNUC__) && defined(__linux__)
#if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64) || defined(__ia64) || defined(__itanium__) || defined(_M_IA64) || defined (__PPC64__)
#undef MPREAL_HAVE_INT64_SUPPORT // Remove all shaman dances for x64 builds since
#undef MPFR_USE_INTMAX_T // GCC already supports x64 as of "long int" is 64-bit integer, nothing left to do
#else
#include <stdint.h> // use int64_t, uint64_t otherwise
#endif
#else
#include <stdint.h> // rely on int64_t, uint64_t in all other cases, Mac OSX, etc.
#endif
#endif
#define MPFR_USE_INTMAX_T // Enable 64-bit integer types - should be defined before mpfr.h
#if defined(MPREAL_HAVE_MSVC_DEBUGVIEW) && defined(_MSC_VER) && defined(_DEBUG)
#define MPREAL_MSVC_DEBUGVIEW_CODE DebugView = toString();
@ -153,6 +128,12 @@
#define MPREAL_DOUBLE_BITS_OVERFLOW -1 // Triggers overflow exception during conversion to double if mpreal
// cannot fit in MPREAL_DOUBLE_BITS_OVERFLOW bits
// = -1 disables overflow checks (default)
// Fast replacement for mpfr_set_zero(x, +1):
// (a) uses low-level data members, might not be compatible with new versions of MPFR
// (b) sign is not set, add (x)->_mpfr_sign = 1;
#define mpfr_set_zero_fast(x) ((x)->_mpfr_exp = __MPFR_EXP_ZERO)
#if defined(__GNUC__)
#define MPREAL_PERMISSIVE_EXPR __extension__
#else
@ -179,6 +160,8 @@ public:
mpreal(const mpq_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const double u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const long double u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const unsigned long long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const long long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const unsigned long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const unsigned int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
@ -188,11 +171,6 @@ public:
// shared = true allows to avoid deep copy, so that mpreal and 'u' share the same data & pointers.
mpreal(const mpfr_t u, bool shared = false);
#if defined (MPREAL_HAVE_INT64_SUPPORT)
mpreal(const uint64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const int64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd());
#endif
mpreal(const char* s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, mp_rnd_t mode = mpreal::get_default_rnd());
mpreal(const std::string& s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, mp_rnd_t mode = mpreal::get_default_rnd());
@ -217,11 +195,14 @@ public:
mpreal& operator=(const long double v);
mpreal& operator=(const double v);
mpreal& operator=(const unsigned long int v);
mpreal& operator=(const unsigned long long int v);
mpreal& operator=(const long long int v);
mpreal& operator=(const unsigned int v);
mpreal& operator=(const long int v);
mpreal& operator=(const int v);
mpreal& operator=(const char* s);
mpreal& operator=(const std::string& s);
template <typename real_t> mpreal& operator= (const std::complex<real_t>& z);
// +
mpreal& operator+=(const mpreal& v);
@ -235,16 +216,14 @@ public:
mpreal& operator+=(const long int u);
mpreal& operator+=(const int u);
#if defined (MPREAL_HAVE_INT64_SUPPORT)
mpreal& operator+=(const int64_t u);
mpreal& operator+=(const uint64_t u);
mpreal& operator-=(const int64_t u);
mpreal& operator-=(const uint64_t u);
mpreal& operator*=(const int64_t u);
mpreal& operator*=(const uint64_t u);
mpreal& operator/=(const int64_t u);
mpreal& operator/=(const uint64_t u);
#endif
mpreal& operator+=(const long long int u);
mpreal& operator+=(const unsigned long long int u);
mpreal& operator-=(const long long int u);
mpreal& operator-=(const unsigned long long int u);
mpreal& operator*=(const long long int u);
mpreal& operator*=(const unsigned long long int u);
mpreal& operator/=(const long long int u);
mpreal& operator/=(const unsigned long long int u);
const mpreal operator+() const;
mpreal& operator++ ();
@ -308,26 +287,12 @@ public:
mpreal& operator>>=(const long int u);
mpreal& operator>>=(const int u);
// Boolean Operators
friend bool operator > (const mpreal& a, const mpreal& b);
friend bool operator >= (const mpreal& a, const mpreal& b);
friend bool operator < (const mpreal& a, const mpreal& b);
friend bool operator <= (const mpreal& a, const mpreal& b);
friend bool operator == (const mpreal& a, const mpreal& b);
friend bool operator != (const mpreal& a, const mpreal& b);
// Optimized specializations for boolean operators
friend bool operator == (const mpreal& a, const unsigned long int b);
friend bool operator == (const mpreal& a, const unsigned int b);
friend bool operator == (const mpreal& a, const long int b);
friend bool operator == (const mpreal& a, const int b);
friend bool operator == (const mpreal& a, const long double b);
friend bool operator == (const mpreal& a, const double b);
// Type Conversion operators
bool toBool (mp_rnd_t mode = GMP_RNDZ) const;
bool toBool ( ) const;
long toLong (mp_rnd_t mode = GMP_RNDZ) const;
unsigned long toULong (mp_rnd_t mode = GMP_RNDZ) const;
long long toLLong (mp_rnd_t mode = GMP_RNDZ) const;
unsigned long long toULLong (mp_rnd_t mode = GMP_RNDZ) const;
float toFloat (mp_rnd_t mode = GMP_RNDN) const;
double toDouble (mp_rnd_t mode = GMP_RNDN) const;
long double toLDouble (mp_rnd_t mode = GMP_RNDN) const;
@ -336,25 +301,15 @@ public:
explicit operator bool () const { return toBool(); }
explicit operator int () const { return int(toLong()); }
explicit operator long () const { return toLong(); }
explicit operator long long () const { return toLong(); }
explicit operator long long () const { return toLLong(); }
explicit operator unsigned () const { return unsigned(toULong()); }
explicit operator unsigned long () const { return toULong(); }
explicit operator unsigned long long () const { return toULong(); }
explicit operator unsigned long long () const { return toULLong(); }
explicit operator float () const { return toFloat(); }
explicit operator double () const { return toDouble(); }
explicit operator long double () const { return toLDouble(); }
#endif
#if defined (MPREAL_HAVE_INT64_SUPPORT)
int64_t toInt64 (mp_rnd_t mode = GMP_RNDZ) const;
uint64_t toUInt64 (mp_rnd_t mode = GMP_RNDZ) const;
#if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS)
explicit operator int64_t () const { return toInt64(); }
explicit operator uint64_t () const { return toUInt64(); }
#endif
#endif
// Get raw pointers so that mpreal can be directly used in raw mpfr_* functions
::mpfr_ptr mpfr_ptr();
::mpfr_srcptr mpfr_ptr() const;
@ -394,6 +349,7 @@ public:
friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal logb (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode);
@ -436,6 +392,7 @@ public:
friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal tgamma (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode);
friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode);
@ -450,7 +407,7 @@ public:
friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode);
friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode);
friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode);
friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode);
friend const mpreal sum (const mpreal tab[], const unsigned long int n, int& status, mp_rnd_t rnd_mode);
friend int sgn(const mpreal& v); // returns -1 or +1
// MPFR 2.4.0 Specifics
@ -465,11 +422,13 @@ public:
friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Modulus after division
#endif
// MPFR 3.0.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
friend const mpreal digamma (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal ai (const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear
#endif
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0))
friend const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear
friend const mpreal grandom (unsigned int seed);
#endif
@ -480,10 +439,6 @@ public:
// Check urandom() for more precise control.
friend const mpreal random(unsigned int seed);
// Exponent and mantissa manipulation
friend const mpreal frexp(const mpreal& v, mp_exp_t* exp);
friend const mpreal ldexp(const mpreal& v, mp_exp_t exp);
// Splits mpreal value into fractional and integer parts.
// Returns fractional part and stores integer part in n.
friend const mpreal modf(const mpreal& v, mpreal& n);
@ -530,9 +485,9 @@ public:
#endif
// Instance Checkers
friend bool (isnan) (const mpreal& v);
friend bool (isinf) (const mpreal& v);
friend bool (isfinite) (const mpreal& v);
friend bool isnan (const mpreal& v);
friend bool isinf (const mpreal& v);
friend bool isfinite (const mpreal& v);
friend bool isnum (const mpreal& v);
friend bool iszero (const mpreal& v);
@ -612,7 +567,7 @@ public:
inline mpreal::mpreal()
{
mpfr_init2(mpfr_ptr(), mpreal::get_default_prec());
mpfr_set_ui(mpfr_ptr(), 0, mpreal::get_default_rnd());
mpfr_set_zero_fast(mpfr_ptr());
MPREAL_MSVC_DEBUGVIEW_CODE;
}
@ -628,7 +583,7 @@ inline mpreal::mpreal(const mpreal& u)
#ifdef MPREAL_HAVE_MOVE_SUPPORT
inline mpreal::mpreal(mpreal&& other)
{
mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds no pinter to actual data
mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds no pointer to actual data
mpfr_swap(mpfr_ptr(), other.mpfr_ptr());
MPREAL_MSVC_DEBUGVIEW_CODE;
@ -707,6 +662,22 @@ inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const unsigned long long int u, mp_prec_t prec, mp_rnd_t mode)
{
mpfr_init2 (mpfr_ptr(), prec);
mpfr_set_uj(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const long long int u, mp_prec_t prec, mp_rnd_t mode)
{
mpfr_init2 (mpfr_ptr(), prec);
mpfr_set_sj(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
{
mpfr_init2 (mpfr_ptr(), prec);
@ -739,24 +710,6 @@ inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
MPREAL_MSVC_DEBUGVIEW_CODE;
}
#if defined (MPREAL_HAVE_INT64_SUPPORT)
inline mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode)
{
mpfr_init2 (mpfr_ptr(), prec);
mpfr_set_uj(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode)
{
mpfr_init2 (mpfr_ptr(), prec);
mpfr_set_sj(mpfr_ptr(), u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
#endif
inline mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
{
mpfr_init2 (mpfr_ptr(), prec);
@ -803,11 +756,8 @@ namespace internal{
template <> struct result_type<unsigned int> {typedef mpreal type;};
template <> struct result_type<long int> {typedef mpreal type;};
template <> struct result_type<int> {typedef mpreal type;};
#if defined (MPREAL_HAVE_INT64_SUPPORT)
template <> struct result_type<int64_t > {typedef mpreal type;};
template <> struct result_type<uint64_t > {typedef mpreal type;};
#endif
template <> struct result_type<long long> {typedef mpreal type;};
template <> struct result_type<unsigned long long> {typedef mpreal type;};
}
// + Addition
@ -1040,6 +990,22 @@ inline mpreal& mpreal::operator=(const unsigned int v)
return *this;
}
inline mpreal& mpreal::operator=(const unsigned long long int v)
{
mpfr_set_uj(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator=(const long long int v)
{
mpfr_set_sj(mpfr_ptr(), v, mpreal::get_default_rnd());
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
inline mpreal& mpreal::operator=(const long int v)
{
mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());
@ -1102,6 +1068,11 @@ inline mpreal& mpreal::operator=(const std::string& s)
return *this;
}
template <typename real_t>
inline mpreal& mpreal::operator= (const std::complex<real_t>& z)
{
return *this = z.real();
}
//////////////////////////////////////////////////////////////////////////
// + Addition
@ -1180,16 +1151,14 @@ inline mpreal& mpreal::operator+=(const int u)
return *this;
}
#if defined (MPREAL_HAVE_INT64_SUPPORT)
inline mpreal& mpreal::operator+=(const int64_t u){ *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator+=(const uint64_t u){ *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator-=(const int64_t u){ *this -= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator-=(const uint64_t u){ *this -= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator*=(const int64_t u){ *this *= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator*=(const uint64_t u){ *this *= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator/=(const int64_t u){ *this /= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator/=(const uint64_t u){ *this /= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
#endif
inline mpreal& mpreal::operator+=(const long long int u) { *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator+=(const unsigned long long int u){ *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator-=(const long long int u) { *this -= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator-=(const unsigned long long int u){ *this -= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator*=(const long long int u) { *this *= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator*=(const unsigned long long int u){ *this *= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator/=(const long long int u) { *this /= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline mpreal& mpreal::operator/=(const unsigned long long int u){ *this /= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; }
inline const mpreal mpreal::operator+()const { return mpreal(*this); }
@ -1671,25 +1640,69 @@ inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)
}
//////////////////////////////////////////////////////////////////////////
//Boolean operators
//Relational operators
// WARNING:
//
// Please note that following checks for double-NaN are guaranteed to work only in IEEE math mode:
//
// isnan(b) = (b != b)
// isnan(b) = !(b == b) (we use in code below)
//
// Be cautions if you use compiler options which break strict IEEE compliance (e.g. -ffast-math in GCC).
// Use std::isnan instead (C++11).
inline bool operator > (const mpreal& a, const mpreal& b ){ return (mpfr_greater_p(a.mpfr_srcptr(),b.mpfr_srcptr()) != 0 ); }
inline bool operator > (const mpreal& a, const unsigned long int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) > 0 ); }
inline bool operator > (const mpreal& a, const unsigned int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) > 0 ); }
inline bool operator > (const mpreal& a, const long int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) > 0 ); }
inline bool operator > (const mpreal& a, const int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) > 0 ); }
inline bool operator > (const mpreal& a, const long double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(),b) > 0 ); }
inline bool operator > (const mpreal& a, const double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_d (a.mpfr_srcptr(),b) > 0 ); }
inline bool operator >= (const mpreal& a, const mpreal& b ){ return (mpfr_greaterequal_p(a.mpfr_srcptr(),b.mpfr_srcptr()) != 0 ); }
inline bool operator >= (const mpreal& a, const unsigned long int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator >= (const mpreal& a, const unsigned int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator >= (const mpreal& a, const long int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator >= (const mpreal& a, const int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator >= (const mpreal& a, const long double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator >= (const mpreal& a, const double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_d (a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator < (const mpreal& a, const mpreal& b ){ return (mpfr_less_p(a.mpfr_srcptr(),b.mpfr_srcptr()) != 0 ); }
inline bool operator < (const mpreal& a, const unsigned long int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) < 0 ); }
inline bool operator < (const mpreal& a, const unsigned int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) < 0 ); }
inline bool operator < (const mpreal& a, const long int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) < 0 ); }
inline bool operator < (const mpreal& a, const int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) < 0 ); }
inline bool operator < (const mpreal& a, const long double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(),b) < 0 ); }
inline bool operator < (const mpreal& a, const double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_d (a.mpfr_srcptr(),b) < 0 ); }
inline bool operator <= (const mpreal& a, const mpreal& b ){ return (mpfr_lessequal_p(a.mpfr_srcptr(),b.mpfr_srcptr()) != 0 ); }
inline bool operator <= (const mpreal& a, const unsigned long int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) <= 0 ); }
inline bool operator <= (const mpreal& a, const unsigned int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) <= 0 ); }
inline bool operator <= (const mpreal& a, const long int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) <= 0 ); }
inline bool operator <= (const mpreal& a, const int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) <= 0 ); }
inline bool operator <= (const mpreal& a, const long double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(),b) <= 0 ); }
inline bool operator <= (const mpreal& a, const double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_d (a.mpfr_srcptr(),b) <= 0 ); }
inline bool operator == (const mpreal& a, const mpreal& b ){ return (mpfr_equal_p(a.mpfr_srcptr(),b.mpfr_srcptr()) != 0 ); }
inline bool operator != (const mpreal& a, const mpreal& b){ return (mpfr_lessgreater_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); }
inline bool operator == (const mpreal& a, const unsigned long int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const unsigned int b ){ return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const long int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const int b ){ return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const long double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const double b ){ return !isnan(a) && (b == b) && (mpfr_cmp_d (a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const unsigned long int b ){ return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const unsigned int b ){ return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const long int b ){ return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const int b ){ return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const long double b ){ return (mpfr_cmp_ld(a.mpfr_srcptr(),b) == 0 ); }
inline bool operator == (const mpreal& a, const double b ){ return (mpfr_cmp_d (a.mpfr_srcptr(),b) == 0 ); }
inline bool operator != (const mpreal& a, const mpreal& b ){ return !(a == b); }
inline bool operator != (const mpreal& a, const unsigned long int b ){ return !(a == b); }
inline bool operator != (const mpreal& a, const unsigned int b ){ return !(a == b); }
inline bool operator != (const mpreal& a, const long int b ){ return !(a == b); }
inline bool operator != (const mpreal& a, const int b ){ return !(a == b); }
inline bool operator != (const mpreal& a, const long double b ){ return !(a == b); }
inline bool operator != (const mpreal& a, const double b ){ return !(a == b); }
inline bool (isnan) (const mpreal& op){ return (mpfr_nan_p (op.mpfr_srcptr()) != 0 ); }
inline bool (isinf) (const mpreal& op){ return (mpfr_inf_p (op.mpfr_srcptr()) != 0 ); }
inline bool (isfinite) (const mpreal& op){ return (mpfr_number_p (op.mpfr_srcptr()) != 0 ); }
inline bool isnan (const mpreal& op){ return (mpfr_nan_p (op.mpfr_srcptr()) != 0 ); }
inline bool isinf (const mpreal& op){ return (mpfr_inf_p (op.mpfr_srcptr()) != 0 ); }
inline bool isfinite (const mpreal& op){ return (mpfr_number_p (op.mpfr_srcptr()) != 0 ); }
inline bool iszero (const mpreal& op){ return (mpfr_zero_p (op.mpfr_srcptr()) != 0 ); }
inline bool isint (const mpreal& op){ return (mpfr_integer_p(op.mpfr_srcptr()) != 0 ); }
@ -1699,17 +1712,14 @@ inline bool isregular(const mpreal& op){ return (mpfr_regular_p(op.mpfr_srcpt
//////////////////////////////////////////////////////////////////////////
// Type Converters
inline bool mpreal::toBool (mp_rnd_t /*mode*/) const { return mpfr_zero_p (mpfr_srcptr()) == 0; }
inline bool mpreal::toBool ( ) const { return mpfr_zero_p (mpfr_srcptr()) == 0; }
inline long mpreal::toLong (mp_rnd_t mode) const { return mpfr_get_si (mpfr_srcptr(), mode); }
inline unsigned long mpreal::toULong (mp_rnd_t mode) const { return mpfr_get_ui (mpfr_srcptr(), mode); }
inline float mpreal::toFloat (mp_rnd_t mode) const { return mpfr_get_flt(mpfr_srcptr(), mode); }
inline double mpreal::toDouble (mp_rnd_t mode) const { return mpfr_get_d (mpfr_srcptr(), mode); }
inline long double mpreal::toLDouble(mp_rnd_t mode) const { return mpfr_get_ld (mpfr_srcptr(), mode); }
#if defined (MPREAL_HAVE_INT64_SUPPORT)
inline int64_t mpreal::toInt64 (mp_rnd_t mode) const{ return mpfr_get_sj(mpfr_srcptr(), mode); }
inline uint64_t mpreal::toUInt64(mp_rnd_t mode) const{ return mpfr_get_uj(mpfr_srcptr(), mode); }
#endif
inline long long mpreal::toLLong (mp_rnd_t mode) const { return mpfr_get_sj (mpfr_srcptr(), mode); }
inline unsigned long long mpreal::toULLong (mp_rnd_t mode) const { return mpfr_get_uj (mpfr_srcptr(), mode); }
inline ::mpfr_ptr mpreal::mpfr_ptr() { return mp; }
inline ::mpfr_srcptr mpreal::mpfr_ptr() const { return mp; }
@ -1755,7 +1765,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
std::ostringstream format;
int digits = (n >= 0) ? n : bits2digits(mpfr_get_prec(mpfr_srcptr()));
int digits = (n >= 0) ? n : 1 + bits2digits(mpfr_get_prec(mpfr_srcptr()));
format << "%." << digits << "RNg";
@ -1926,8 +1936,7 @@ inline int bits2digits(mp_prec_t b)
// Set/Get number properties
inline int sgn(const mpreal& op)
{
int r = mpfr_signbit(op.mpfr_srcptr());
return (r > 0? -1 : 1);
return mpfr_sgn(op.mpfr_srcptr());
}
inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode)
@ -1965,7 +1974,6 @@ inline mpreal& mpreal::setNan()
inline mpreal& mpreal::setZero(int sign)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
mpfr_set_zero(mpfr_ptr(), sign);
#else
@ -2000,23 +2008,32 @@ inline int mpreal::set_exp (mp_exp_t e)
return x;
}
inline const mpreal frexp(const mpreal& v, mp_exp_t* exp)
inline const mpreal frexp(const mpreal& x, mp_exp_t* exp, mp_rnd_t mode = mpreal::get_default_rnd())
{
mpreal x(v);
*exp = x.get_exp();
x.set_exp(0);
return x;
mpreal y(x);
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0))
mpfr_frexp(exp,y.mpfr_ptr(),x.mpfr_srcptr(),mode);
#else
*exp = mpfr_get_exp(y.mpfr_srcptr());
mpfr_set_exp(y.mpfr_ptr(),0);
#endif
return y;
}
inline const mpreal ldexp(const mpreal& v, mp_exp_t exp)
{
mpreal x(v);
// rounding is not important since we just increasing the exponent
// rounding is not important since we are just increasing the exponent (= exact operation)
mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd());
return x;
}
inline const mpreal scalbn(const mpreal& v, mp_exp_t exp)
{
return ldexp(v, exp);
}
inline mpreal machine_epsilon(mp_prec_t prec)
{
/* the smallest eps such that 1 + eps != 1 */
@ -2063,6 +2080,20 @@ inline bool isEqualFuzzy(const mpreal& a, const mpreal& b)
return isEqualFuzzy(a, b, machine_epsilon((max)(1, (min)(abs(a), abs(b)))));
}
//////////////////////////////////////////////////////////////////////////
// C++11 sign functions.
inline mpreal copysign(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal rop(0, mpfr_get_prec(x.mpfr_ptr()));
mpfr_setsign(rop.mpfr_ptr(), x.mpfr_srcptr(), mpfr_signbit(y.mpfr_srcptr()), rnd_mode);
return rop;
}
inline bool signbit(const mpreal& x)
{
return mpfr_signbit(x.mpfr_srcptr());
}
inline const mpreal modf(const mpreal& v, mpreal& n)
{
mpreal f(v);
@ -2209,6 +2240,8 @@ inline const mpreal acos (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd
inline const mpreal asin (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(asin ); }
inline const mpreal atan (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(atan ); }
inline const mpreal logb (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { return log2 (abs(x),r); }
inline const mpreal acot (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return atan (1/v, r); }
inline const mpreal asec (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return acos (1/v, r); }
inline const mpreal acsc (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return asin (1/v, r); }
@ -2230,6 +2263,7 @@ inline const mpreal log1p (const mpreal& x, mp_rnd_t r = mpreal::get_default_r
inline const mpreal expm1 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(expm1 ); }
inline const mpreal eint (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(eint ); }
inline const mpreal gamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(gamma ); }
inline const mpreal tgamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(gamma ); }
inline const mpreal lngamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma); }
inline const mpreal zeta (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(zeta ); }
inline const mpreal erf (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(erf ); }
@ -2347,16 +2381,17 @@ inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode =
return a;
}
inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
inline const mpreal sum (const mpreal tab[], const unsigned long int n, int& status, mp_rnd_t mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_ptr* t;
unsigned long int i;
mpfr_srcptr *p = new mpfr_srcptr[n];
t = new mpfr_ptr[n];
for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp;
mpfr_sum(x.mp,t,n,rnd_mode);
delete[] t;
for (unsigned long int i = 0; i < n; i++)
p[i] = tab[i].mpfr_srcptr();
mpreal x;
status = mpfr_sum(x.mpfr_ptr(), (mpfr_ptr*)p, n, mode);
delete [] p;
return x;
}
@ -2553,33 +2588,24 @@ inline const mpreal nextbelow (const mpreal& x)
inline const mpreal urandomb (gmp_randstate_t& state)
{
mpreal x;
mpfr_urandomb(x.mp,state);
mpfr_urandomb(x.mpfr_ptr(),state);
return x;
}
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0))
// use gmp_randinit_default() to init state, gmp_randclear() to clear
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_urandom(x.mp,state,rnd_mode);
mpfr_urandom(x.mpfr_ptr(), state, rnd_mode);
return x;
}
inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_grandom(x.mp, NULL, state, rnd_mode);
return x;
}
#endif
#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2))
inline const mpreal random2 (mp_size_t size, mp_exp_t exp)
{
mpreal x;
mpfr_random2(x.mp,size,exp);
mpfr_random2(x.mpfr_ptr(),size,exp);
return x;
}
#endif
@ -2590,16 +2616,15 @@ inline const mpreal random2 (mp_size_t size, mp_exp_t exp)
// seed != 0
inline const mpreal random(unsigned int seed = 0)
{
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
static gmp_randstate_t state;
static bool isFirstTime = true;
static bool initialize = true;
if(isFirstTime)
if(initialize)
{
gmp_randinit_default(state);
gmp_randseed_ui(state,0);
isFirstTime = false;
initialize = false;
}
if(seed != 0) gmp_randseed_ui(state,seed);
@ -2612,17 +2637,25 @@ inline const mpreal random(unsigned int seed = 0)
}
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0))
inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_grandom(x.mpfr_ptr(), NULL, state, rnd_mode);
return x;
}
inline const mpreal grandom(unsigned int seed = 0)
{
static gmp_randstate_t state;
static bool isFirstTime = true;
static bool initialize = true;
if(isFirstTime)
if(initialize)
{
gmp_randinit_default(state);
gmp_randseed_ui(state,0);
isFirstTime = false;
initialize = false;
}
if(seed != 0) gmp_randseed_ui(state,seed);