From 3e32f6b554d3213ef2af0eed9223b0e579af3056 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 13 Oct 2015 09:58:54 +0200 Subject: update mpreal.h --- unsupported/test/mpreal/mpreal.h | 853 ++++++++++++++++++++------------------- 1 file changed, 443 insertions(+), 410 deletions(-) (limited to 'unsupported/test/mpreal') diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h index 104cb686f..f896515aa 100644 --- a/unsupported/test/mpreal/mpreal.h +++ b/unsupported/test/mpreal/mpreal.h @@ -1,33 +1,34 @@ /* - MPFR C++: Multi-precision floating point number class for C++. + MPFR C++: Multi-precision floating point number class for C++. Based on MPFR library: http://mpfr.org 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, + 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"). - - (B) Non-free licenses may also be purchased from the author, for users who + + (B) Non-free licenses may also be purchased from the author, for users who do not want their programs protected by the GPL. - The non-free licenses are for users that wish to use MPFR C++ in - their products but are unwilling to release their software - under the GPL (which would require them to release source code + The non-free licenses are for users that wish to use MPFR C++ in + their products but are unwilling to release their software + under the GPL (which would require them to release source code and allow free redistribution). Such users can purchase an unlimited-use license from the author. Contact us for more details. - + GNU General Public License ("GPL") copyright permissions statement: ************************************************************************** This program is free software: you can redistribute it and/or modify @@ -55,10 +56,11 @@ #include #include #include +#include +#include +#include // 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 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)) +#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. @@ -93,54 +95,27 @@ #define MPREAL_HAVE_MOVE_SUPPORT - // Use fields in mpfr_t structure to check if it was initialized / set dummy initialization + // Use fields in mpfr_t structure to check if it was initialized / set dummy initialization #define mpfr_is_initialized(x) (0 != (x)->_mpfr_d) #define mpfr_set_uninitialized(x) ((x)->_mpfr_d = 0 ) #endif -// Detect support for explicit converters. +// Detect support for explicit converters. #if (__has_feature(cxx_explicit_conversions) || \ - defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \ - (defined(_MSC_VER) && _MSC_VER >= 1800)) + (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 // 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 // use int64_t, uint64_t otherwise - #endif - - #else - #include // 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(); #define MPREAL_MSVC_DEBUGVIEW_DATA std::string DebugView; #else - #define MPREAL_MSVC_DEBUGVIEW_CODE - #define MPREAL_MSVC_DEBUGVIEW_DATA + #define MPREAL_MSVC_DEBUGVIEW_CODE + #define MPREAL_MSVC_DEBUGVIEW_DATA #endif #include @@ -150,9 +125,15 @@ #endif // Less important options -#define MPREAL_DOUBLE_BITS_OVERFLOW -1 // Triggers overflow exception during conversion to double if mpreal +#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 @@ -164,9 +145,9 @@ namespace mpfr { class mpreal { private: mpfr_t mp; - + public: - + // Get default rounding mode & precision inline static mp_rnd_t get_default_rnd() { return (mp_rnd_t)(mpfr_get_default_rounding_mode()); } inline static mp_prec_t get_default_prec() { return mpfr_get_default_prec(); } @@ -174,29 +155,26 @@ public: // Constructors && type conversions mpreal(); mpreal(const mpreal& u); - mpreal(const mpf_t u); - mpreal(const mpz_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); - 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 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()); - mpreal(const int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); - - // Construct mpreal from mpfr_t structure. - // 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); + mpreal(const mpf_t u); + mpreal(const mpz_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); + 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()); + mpreal(const int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); -#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 + // Construct mpreal from mpfr_t structure. + // 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); 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()); - ~mpreal(); + ~mpreal(); #ifdef MPREAL_HAVE_MOVE_SUPPORT mpreal& operator=(mpreal&& v); @@ -205,7 +183,7 @@ public: // Operations // = - // +, -, *, /, ++, --, <<, >> + // +, -, *, /, ++, --, <<, >> // *=, +=, -=, /=, // <, >, ==, <=, >= @@ -215,13 +193,16 @@ public: mpreal& operator=(const mpz_t v); mpreal& operator=(const mpq_t v); mpreal& operator=(const long double v); - mpreal& operator=(const 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 mpreal& operator= (const std::complex& z); // + mpreal& operator+=(const mpreal& v); @@ -235,20 +216,18 @@ 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++ (); - const mpreal operator++ (int); + const mpreal operator++ (int); // - mpreal& operator-=(const mpreal& v); @@ -266,7 +245,7 @@ public: friend const mpreal operator-(const long int b, const mpreal& a); friend const mpreal operator-(const int b, const mpreal& a); friend const mpreal operator-(const double b, const mpreal& a); - mpreal& operator-- (); + mpreal& operator-- (); const mpreal operator-- (int); // * @@ -279,7 +258,7 @@ public: mpreal& operator*=(const unsigned int v); mpreal& operator*=(const long int v); mpreal& operator*=(const int v); - + // / mpreal& operator/=(const mpreal& v); mpreal& operator/=(const mpz_t v); @@ -308,51 +287,27 @@ 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; - long toLong (mp_rnd_t mode = GMP_RNDZ) const; - unsigned long toULong (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; + 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; #if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS) - 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 unsigned () const { return unsigned(toULong()); } - explicit operator unsigned long () const { return toULong(); } - explicit operator unsigned long long () const { return toULong(); } - 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 + 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 toLLong(); } + explicit operator unsigned () const { return unsigned(toULong()); } + explicit operator unsigned 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 // Get raw pointers so that mpreal can be directly used in raw mpfr_* functions @@ -391,11 +346,12 @@ public: friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode); friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode); friend int cmpabs(const mpreal& a,const mpreal& b); - + 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 exp (const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode); @@ -436,21 +392,22 @@ 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); friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode); + friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode); + friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal besseljn (long n, const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal bessely0 (const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal bessely1 (const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode); + friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode); 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,28 +422,26 @@ 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 - + // Uniformly distributed random number generation in [0,1] using // Mersenne-Twister algorithm by default. // Use parameter to setup seed, e.g.: random((unsigned)time(NULL)) // 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); + friend const mpreal modf(const mpreal& v, mpreal& n); // Constants // don't forget to call mpfr_free_cache() for every thread where you are using const-functions @@ -515,14 +470,14 @@ public: friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal remainder ( const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); - + // Miscellaneous Functions friend const mpreal nexttoward (const mpreal& x, const mpreal& y); friend const mpreal nextabove (const mpreal& x); friend const mpreal nextbelow (const mpreal& x); // use gmp_randinit_default() to init state, gmp_randclear() to clear - friend const mpreal urandomb (gmp_randstate_t& state); + friend const mpreal urandomb (gmp_randstate_t& state); // MPFR < 2.4.2 Specifics #if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) @@ -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); @@ -549,9 +504,9 @@ public: // Aliases for get_prec(), set_prec() - needed for compatibility with std::complex interface inline mpreal& setPrecision(int Precision, mp_rnd_t RoundingMode = get_default_rnd()); inline int getPrecision() const; - + // Set mpreal to +/- inf, NaN, +/-0 - mpreal& setInf (int Sign = +1); + mpreal& setInf (int Sign = +1); mpreal& setNan (); mpreal& setZero (int Sign = +1); mpreal& setSign (int Sign, mp_rnd_t RoundingMode = get_default_rnd()); @@ -560,7 +515,7 @@ public: mp_exp_t get_exp(); int set_exp(mp_exp_t e); int check_range (int t, mp_rnd_t rnd_mode = get_default_rnd()); - int subnormalize (int t,mp_rnd_t rnd_mode = get_default_rnd()); + int subnormalize (int t, mp_rnd_t rnd_mode = get_default_rnd()); // Inexact conversion from float inline bool fits_in_bits(double x, int n); @@ -580,7 +535,7 @@ public: // Efficient swapping of two mpreal values - needed for std algorithms friend void swap(mpreal& x, mpreal& y); - + friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); @@ -590,7 +545,7 @@ private: // // mpfr::mpreal= ; Show value only // mpfr::mpreal=, bits ; Show value & precision - // + // // at the beginning of // [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat MPREAL_MSVC_DEBUGVIEW_DATA @@ -609,15 +564,15 @@ public: ////////////////////////////////////////////////////////////////////////// // Constructors & converters // Default constructor: creates mp number and initializes it to 0. -inline mpreal::mpreal() -{ - mpfr_init2 (mpfr_ptr(), mpreal::get_default_prec()); - mpfr_set_ui(mpfr_ptr(), 0, mpreal::get_default_rnd()); +inline mpreal::mpreal() +{ + mpfr_init2(mpfr_ptr(), mpreal::get_default_prec()); + mpfr_set_zero_fast(mpfr_ptr()); MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mpreal::mpreal(const mpreal& u) +inline mpreal::mpreal(const mpreal& u) { mpfr_init2(mpfr_ptr(),mpfr_get_prec(u.mpfr_srcptr())); mpfr_set (mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd()); @@ -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; @@ -700,67 +655,65 @@ inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) } inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) -{ +{ mpfr_init2 (mpfr_ptr(), prec); mpfr_set_ld(mpfr_ptr(), u, mode); MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) -{ +inline mpreal::mpreal(const unsigned long long int u, mp_prec_t prec, mp_rnd_t mode) +{ mpfr_init2 (mpfr_ptr(), prec); - mpfr_set_ui(mpfr_ptr(), u, mode); + mpfr_set_uj(mpfr_ptr(), u, mode); MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) -{ +inline mpreal::mpreal(const long long int u, mp_prec_t prec, mp_rnd_t mode) +{ mpfr_init2 (mpfr_ptr(), prec); - mpfr_set_ui(mpfr_ptr(), u, mode); + mpfr_set_sj(mpfr_ptr(), u, mode); MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) -{ +inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) +{ mpfr_init2 (mpfr_ptr(), prec); - mpfr_set_si(mpfr_ptr(), u, mode); + mpfr_set_ui(mpfr_ptr(), u, mode); MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) -{ +inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) +{ mpfr_init2 (mpfr_ptr(), prec); - mpfr_set_si(mpfr_ptr(), u, mode); + mpfr_set_ui(mpfr_ptr(), u, 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) +inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) { mpfr_init2 (mpfr_ptr(), prec); - mpfr_set_uj(mpfr_ptr(), u, mode); + mpfr_set_si(mpfr_ptr(), u, mode); MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode) +inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) { mpfr_init2 (mpfr_ptr(), prec); - mpfr_set_sj(mpfr_ptr(), u, mode); + mpfr_set_si(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); - mpfr_set_str(mpfr_ptr(), s, base, mode); + mpfr_set_str(mpfr_ptr(), s, base, mode); MPREAL_MSVC_DEBUGVIEW_CODE; } @@ -768,7 +721,7 @@ inline mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode) { mpfr_init2 (mpfr_ptr(), prec); - mpfr_set_str(mpfr_ptr(), s.c_str(), base, mode); + mpfr_set_str(mpfr_ptr(), s.c_str(), base, mode); MPREAL_MSVC_DEBUGVIEW_CODE; } @@ -776,15 +729,15 @@ inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t m inline void mpreal::clear(::mpfr_ptr x) { #ifdef MPREAL_HAVE_MOVE_SUPPORT - if(mpfr_is_initialized(x)) + if(mpfr_is_initialized(x)) #endif mpfr_clear(x); } -inline mpreal::~mpreal() -{ +inline mpreal::~mpreal() +{ clear(mpfr_ptr()); -} +} // internal namespace needed for template magic namespace internal{ @@ -792,58 +745,55 @@ namespace internal{ // Use SFINAE to restrict arithmetic operations instantiation only for numeric types // This is needed for smooth integration with libraries based on expression templates, like Eigen. // TODO: Do the same for boolean operators. - template struct result_type {}; - - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; - -#if defined (MPREAL_HAVE_INT64_SUPPORT) - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; -#endif + template struct result_type {}; + + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; + template <> struct result_type {typedef mpreal type;}; } // + Addition -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator+(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) += rhs; } -template -inline const typename internal::result_type::type - operator+(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) += lhs; } +template +inline const typename internal::result_type::type + operator+(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) += lhs; } // - Subtraction -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator-(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) -= rhs; } -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator-(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) -= rhs; } // * Multiplication -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator*(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) *= rhs; } -template -inline const typename internal::result_type::type - operator*(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) *= lhs; } +template +inline const typename internal::result_type::type + operator*(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) *= lhs; } // / Division -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator/(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) /= rhs; } -template -inline const typename internal::result_type::type +template +inline const typename internal::result_type::type operator/(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) /= rhs; } ////////////////////////////////////////////////////////////////////////// @@ -893,17 +843,17 @@ const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::g const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); +const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); +const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); +const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); +const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); @@ -920,9 +870,9 @@ inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpr inline mpreal machine_epsilon(mp_prec_t prec = mpreal::get_default_prec()); // Returns smallest eps such that x + eps != x (relative machine epsilon) -inline mpreal machine_epsilon(const mpreal& x); +inline mpreal machine_epsilon(const mpreal& x); -// Gives max & min values for the required precision, +// Gives max & min values for the required precision, // minval is 'safe' meaning 1 / minval does not overflow // maxval is 'safe' meaning 1 / maxval does not underflow inline mpreal minval(mp_prec_t prec = mpreal::get_default_prec()); @@ -935,13 +885,13 @@ inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps); inline bool isEqualFuzzy(const mpreal& a, const mpreal& b); // 'Bitwise' equality check -// maxUlps - a and b can be apart by maxUlps binary numbers. +// maxUlps - a and b can be apart by maxUlps binary numbers. inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps); ////////////////////////////////////////////////////////////////////////// -// Convert precision in 'bits' to decimal digits and vice versa. -// bits = ceil(digits*log[2](10)) -// digits = floor(bits*log[10](2)) +// Convert precision in 'bits' to decimal digits and vice versa. +// bits = ceil(digits*log[2](10)) +// digits = floor(bits*log[10](2)) inline mp_prec_t digits2bits(int d); inline int bits2digits(mp_prec_t b); @@ -979,7 +929,7 @@ inline mpreal& mpreal::operator=(const mpreal& v) inline mpreal& mpreal::operator=(const mpf_t v) { mpfr_set_f(mpfr_ptr(), v, mpreal::get_default_rnd()); - + MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } @@ -987,7 +937,7 @@ inline mpreal& mpreal::operator=(const mpf_t v) inline mpreal& mpreal::operator=(const mpz_t v) { mpfr_set_z(mpfr_ptr(), v, mpreal::get_default_rnd()); - + MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } @@ -1000,16 +950,16 @@ inline mpreal& mpreal::operator=(const mpq_t v) return *this; } -inline mpreal& mpreal::operator=(const long double v) -{ +inline mpreal& mpreal::operator=(const long double v) +{ mpfr_set_ld(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } -inline mpreal& mpreal::operator=(const double v) -{ +inline mpreal& mpreal::operator=(const double v) +{ #if (MPREAL_DOUBLE_BITS_OVERFLOW > -1) if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW)) { @@ -1024,33 +974,49 @@ inline mpreal& mpreal::operator=(const double v) return *this; } -inline mpreal& mpreal::operator=(const unsigned long int v) -{ - mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); +inline mpreal& mpreal::operator=(const unsigned long int v) +{ + mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); + + MPREAL_MSVC_DEBUGVIEW_CODE; + return *this; +} + +inline mpreal& mpreal::operator=(const unsigned int v) +{ + mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } -inline mpreal& mpreal::operator=(const unsigned int v) -{ - mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); +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 int v) -{ - mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); +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()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } inline mpreal& mpreal::operator=(const int v) -{ - mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); +{ + mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; @@ -1071,7 +1037,7 @@ inline mpreal& mpreal::operator=(const char* s) if(0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd())) { - mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); + mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; } @@ -1094,7 +1060,7 @@ inline mpreal& mpreal::operator=(const std::string& s) if(0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd())) { - mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); + mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); MPREAL_MSVC_DEBUGVIEW_CODE; } @@ -1102,6 +1068,11 @@ inline mpreal& mpreal::operator=(const std::string& s) return *this; } +template +inline mpreal& mpreal::operator= (const std::complex& z) +{ + return *this = z.real(); +} ////////////////////////////////////////////////////////////////////////// // + Addition @@ -1135,9 +1106,9 @@ inline mpreal& mpreal::operator+=(const mpq_t u) inline mpreal& mpreal::operator+= (const long double u) { - *this += mpreal(u); + *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; + return *this; } inline mpreal& mpreal::operator+= (const double u) @@ -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); } @@ -1200,7 +1169,7 @@ inline const mpreal operator+(const mpreal& a, const mpreal& b) return c; } -inline mpreal& mpreal::operator++() +inline mpreal& mpreal::operator++() { return *this += 1; } @@ -1212,7 +1181,7 @@ inline const mpreal mpreal::operator++ (int) return x; } -inline mpreal& mpreal::operator--() +inline mpreal& mpreal::operator--() { return *this -= 1; } @@ -1249,9 +1218,9 @@ inline mpreal& mpreal::operator-=(const mpq_t v) inline mpreal& mpreal::operator-=(const long double v) { - *this -= mpreal(v); + *this -= mpreal(v); MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; + return *this; } inline mpreal& mpreal::operator-=(const double v) @@ -1259,7 +1228,7 @@ inline mpreal& mpreal::operator-=(const double v) #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) mpfr_sub_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); #else - *this -= mpreal(v); + *this -= mpreal(v); #endif MPREAL_MSVC_DEBUGVIEW_CODE; @@ -1374,9 +1343,9 @@ inline mpreal& mpreal::operator*=(const mpq_t v) inline mpreal& mpreal::operator*=(const long double v) { - *this *= mpreal(v); + *this *= mpreal(v); MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; + return *this; } inline mpreal& mpreal::operator*=(const double v) @@ -1384,7 +1353,7 @@ inline mpreal& mpreal::operator*=(const double v) #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) mpfr_mul_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); #else - *this *= mpreal(v); + *this *= mpreal(v); #endif MPREAL_MSVC_DEBUGVIEW_CODE; return *this; @@ -1452,7 +1421,7 @@ inline mpreal& mpreal::operator/=(const long double v) { *this /= mpreal(v); MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; + return *this; } inline mpreal& mpreal::operator/=(const double v) @@ -1460,7 +1429,7 @@ inline mpreal& mpreal::operator/=(const double v) #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) mpfr_div_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); #else - *this /= mpreal(v); + *this /= mpreal(v); #endif MPREAL_MSVC_DEBUGVIEW_CODE; return *this; @@ -1671,45 +1640,86 @@ inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode) } ////////////////////////////////////////////////////////////////////////// -//Boolean operators -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 mpreal& b){ return (mpfr_greaterequal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=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 mpreal& b){ return (mpfr_lessequal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=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 (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 (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 ); } +//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 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 !(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 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 ); } #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) inline bool isregular(const mpreal& op){ return (mpfr_regular_p(op.mpfr_srcptr()));} -#endif +#endif ////////////////////////////////////////////////////////////////////////// // Type Converters -inline bool mpreal::toBool (mp_rnd_t /*mode*/) 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 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); } +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,21 +1765,21 @@ 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"; return toString(format.str()); #else - char *s, *ns = NULL; + char *s, *ns = NULL; size_t slen, nslen; mp_exp_t exp; std::string out; if(mpfr_inf_p(mp)) - { + { if(mpfr_sgn(mp)>0) return "+Inf"; else return "-Inf"; } @@ -1784,7 +1794,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { slen = strlen(s); nslen = strlen(ns); - if(nslen<=slen) + if(nslen<=slen) { mpfr_free_str(s); s = ns; @@ -1801,7 +1811,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { // Remove zeros starting from right end char* ptr = s+slen-1; - while (*ptr=='0' && ptr>s+exp) ptr--; + while (*ptr=='0' && ptr>s+exp) ptr--; if(ptr==s+exp) out = std::string(s,exp+1); else out = std::string(s,exp+1)+'.'+std::string(s+exp+1,ptr-(s+exp+1)+1); @@ -1812,7 +1822,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { // Remove zeros starting from right end char* ptr = s+slen-1; - while (*ptr=='0' && ptr>s+exp-1) ptr--; + while (*ptr=='0' && ptr>s+exp-1) ptr--; if(ptr==s+exp-1) out = std::string(s,exp); else out = std::string(s,exp)+'.'+std::string(s+exp,ptr-(s+exp)+1); @@ -1825,7 +1835,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { // Remove zeros starting from right end char* ptr = s+slen-1; - while (*ptr=='0' && ptr>s+1) ptr--; + while (*ptr=='0' && ptr>s+1) ptr--; if(ptr==s+1) out = std::string(s,2); else out = std::string(s,2)+'.'+std::string(s+2,ptr-(s+2)+1); @@ -1836,7 +1846,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { // Remove zeros starting from right end char* ptr = s+slen-1; - while (*ptr=='0' && ptr>s) ptr--; + while (*ptr=='0' && ptr>s) ptr--; if(ptr==s) out = std::string(s,1); else out = std::string(s,1)+'.'+std::string(s+1,ptr-(s+1)+1); @@ -1863,7 +1873,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const ////////////////////////////////////////////////////////////////////////// // I/O -inline std::ostream& mpreal::output(std::ostream& os) const +inline std::ostream& mpreal::output(std::ostream& os) const { std::ostringstream format; const std::ios::fmtflags flags = os.flags(); @@ -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) @@ -1949,29 +1958,28 @@ inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode) return *this; } -inline mpreal& mpreal::setInf(int sign) -{ +inline mpreal& mpreal::setInf(int sign) +{ mpfr_set_inf(mpfr_ptr(), sign); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; -} +} -inline mpreal& mpreal::setNan() +inline mpreal& mpreal::setNan() { mpfr_set_nan(mpfr_ptr()); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } -inline mpreal& mpreal::setZero(int sign) +inline mpreal& mpreal::setZero(int sign) { - #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) mpfr_set_zero(mpfr_ptr(), sign); #else mpfr_set_si(mpfr_ptr(), 0, (mpfr_get_default_rounding_mode)()); setSign(sign); -#endif +#endif MPREAL_MSVC_DEBUGVIEW_CODE; return *this; @@ -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 - mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd()); + // 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 */ @@ -2024,7 +2041,7 @@ inline mpreal machine_epsilon(mp_prec_t prec) } inline mpreal machine_epsilon(const mpreal& x) -{ +{ /* the smallest eps such that x + eps != x */ if( x < 0) { @@ -2045,7 +2062,7 @@ inline mpreal minval(mp_prec_t prec) inline mpreal maxval(mp_prec_t prec) { /* max = (1 - eps) * 2^emax, eps is machine epsilon */ - return (mpreal(1, prec) - machine_epsilon(prec)) << mpreal::get_emax(); + return (mpreal(1, prec) - machine_epsilon(prec)) << mpreal::get_emax(); } inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps) @@ -2063,12 +2080,26 @@ 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); // rounding is not important since we are using the same number - mpfr_frac (f.mpfr_ptr(),f.mpfr_srcptr(),mpreal::get_default_rnd()); + mpfr_frac (f.mpfr_ptr(),f.mpfr_srcptr(),mpreal::get_default_rnd()); mpfr_trunc(n.mpfr_ptr(),v.mpfr_srcptr()); return f; } @@ -2131,7 +2162,7 @@ inline mp_exp_t mpreal::get_emax_max (void) #define MPREAL_UNARY_MATH_FUNCTION_BODY(f) \ mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); \ mpfr_##f(y.mpfr_ptr(), x.mpfr_srcptr(), r); \ - return y; + return y; inline const mpreal sqr (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sqr ); } @@ -2154,7 +2185,7 @@ inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode) inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode) { if (v>=0) return sqrt(static_cast(v),rnd_mode); - else return mpreal().setNan(); // NaN + else return mpreal().setNan(); // NaN } inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) @@ -2165,9 +2196,9 @@ inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r = mpreal::get_default_rnd()) { - mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); - mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r); - return y; + mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); + mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r); + return y; } inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r = mpreal::get_default_rnd()) @@ -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 ); } @@ -2254,7 +2288,7 @@ inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = } inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) -{ +{ mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision())); mpfr_remainder(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode); return a; @@ -2307,9 +2341,9 @@ inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, m mpreal a; mp_prec_t p1, p2, p3; - p1 = v1.get_prec(); - p2 = v2.get_prec(); - p3 = v3.get_prec(); + p1 = v1.get_prec(); + p2 = v2.get_prec(); + p3 = v3.get_prec(); a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); @@ -2322,9 +2356,9 @@ inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, m mpreal a; mp_prec_t p1, p2, p3; - p1 = v1.get_prec(); - p2 = v2.get_prec(); - p3 = v3.get_prec(); + p1 = v1.get_prec(); + p2 = v2.get_prec(); + p3 = v3.get_prec(); a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); @@ -2337,8 +2371,8 @@ inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal a; mp_prec_t p1, p2; - p1 = v1.get_prec(); - p2 = v2.get_prec(); + p1 = v1.get_prec(); + p2 = v2.get_prec(); a.set_prec(p1>p2?p1:p2); @@ -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()) { + mpfr_srcptr *p = new mpfr_srcptr[n]; + + for (unsigned long int i = 0; i < n; i++) + p[i] = tab[i].mpfr_srcptr(); + mpreal x; - mpfr_ptr* t; - unsigned long int i; + status = mpfr_sum(x.mpfr_ptr(), (mpfr_ptr*)p, n, mode); - t = new mpfr_ptr[n]; - for (i=0;ixp?yp:xp); @@ -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); - 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); + mpfr_urandom(x.mpfr_ptr(), state, rnd_mode); return x; } - -#endif +#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); @@ -2634,17 +2667,17 @@ inline const mpreal grandom(unsigned int seed = 0) ////////////////////////////////////////////////////////////////////////// // Set/Get global properties inline void mpreal::set_default_prec(mp_prec_t prec) -{ - mpfr_set_default_prec(prec); +{ + mpfr_set_default_prec(prec); } inline void mpreal::set_default_rnd(mp_rnd_t rnd_mode) -{ - mpfr_set_default_rounding_mode(rnd_mode); +{ + mpfr_set_default_rounding_mode(rnd_mode); } inline bool mpreal::fits_in_bits(double x, int n) -{ +{ int i; double t; return IsInf(x) || (std::modf ( std::ldexp ( std::frexp ( x, &i ), n ), &t ) == 0.0); @@ -2894,7 +2927,7 @@ inline const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode) else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow } -// pow long double +// pow long double inline const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode) { return pow(mpreal(a),mpreal(b),rnd_mode); @@ -2953,9 +2986,9 @@ namespace std { // we are allowed to extend namespace std with specializations only template <> - inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) - { - return mpfr::swap(x, y); + inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) + { + return mpfr::swap(x, y); } template<> @@ -2966,7 +2999,7 @@ namespace std static const bool is_signed = true; static const bool is_integer = false; static const bool is_exact = false; - static const int radix = 2; + static const int radix = 2; static const bool has_infinity = true; static const bool has_quiet_NaN = true; @@ -2986,7 +3019,7 @@ namespace std // Returns smallest eps such that 1 + eps != 1 (classic machine epsilon) inline static mpfr::mpreal epsilon(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(precision); } - + // Returns smallest eps such that x + eps != x (relative machine epsilon) inline static mpfr::mpreal epsilon(const mpfr::mpreal& x) { return mpfr::machine_epsilon(x); } @@ -2994,8 +3027,8 @@ namespace std { mp_rnd_t r = mpfr::mpreal::get_default_rnd(); - if(r == GMP_RNDN) return mpfr::mpreal(0.5, precision); - else return mpfr::mpreal(1.0, precision); + if(r == GMP_RNDN) return mpfr::mpreal(0.5, precision); + else return mpfr::mpreal(1.0, precision); } inline static const mpfr::mpreal infinity() { return mpfr::const_infinity(); } @@ -3006,17 +3039,17 @@ namespace std // Please note, exponent range is not fixed in MPFR static const int min_exponent = MPFR_EMIN_DEFAULT; static const int max_exponent = MPFR_EMAX_DEFAULT; - MPREAL_PERMISSIVE_EXPR static const int min_exponent10 = (int) (MPFR_EMIN_DEFAULT * 0.3010299956639811); - MPREAL_PERMISSIVE_EXPR static const int max_exponent10 = (int) (MPFR_EMAX_DEFAULT * 0.3010299956639811); + MPREAL_PERMISSIVE_EXPR static const int min_exponent10 = (int) (MPFR_EMIN_DEFAULT * 0.3010299956639811); + MPREAL_PERMISSIVE_EXPR static const int max_exponent10 = (int) (MPFR_EMAX_DEFAULT * 0.3010299956639811); #ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS // Following members should be constant according to standard, but they can be variable in MPFR - // So we define them as functions here. + // So we define them as functions here. // // This is preferable way for std::numeric_limits specialization. - // But it is incompatible with standard std::numeric_limits and might not work with other libraries, e.g. boost. - // See below for compatible implementation. + // But it is incompatible with standard std::numeric_limits and might not work with other libraries, e.g. boost. + // See below for compatible implementation. inline static float_round_style round_style() { mp_rnd_t r = mpfr::mpreal::get_default_rnd(); @@ -3024,9 +3057,9 @@ namespace std switch (r) { case GMP_RNDN: return round_to_nearest; - case GMP_RNDZ: return round_toward_zero; - case GMP_RNDU: return round_toward_infinity; - case GMP_RNDD: return round_toward_neg_infinity; + case GMP_RNDZ: return round_toward_zero; + case GMP_RNDU: return round_toward_infinity; + case GMP_RNDD: return round_toward_neg_infinity; default: return round_indeterminate; } } @@ -3053,13 +3086,13 @@ namespace std // If possible, please use functions digits() and round_style() defined above. // // These (default) values are preserved for compatibility with existing libraries, e.g. boost. - // Change them accordingly to your application. + // Change them accordingly to your application. // // For example, if you use 256 bits of precision uniformly in your program, then: // digits = 256 - // digits10 = 77 + // digits10 = 77 // max_digits10 = 78 - // + // // Approximate formula for decimal digits is: digits10 = floor(log10(2) * digits). See bits2digits() for more details. static const std::float_round_style round_style = round_to_nearest; @@ -3071,4 +3104,4 @@ namespace std } -#endif /* __MPREAL_H__ */ +#endif /* __MPREAL_H__ */ \ No newline at end of file -- cgit v1.2.3