From e29bfe84796370c6b12c09f013ef86e37222f743 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 8 Oct 2018 17:09:23 +0200 Subject: Update included mpreal header to 3.6.5 and fix deprecated warnings. --- unsupported/test/mpreal/mpreal.h | 622 ++++++++++++++++++++++----------------- 1 file changed, 349 insertions(+), 273 deletions(-) (limited to 'unsupported/test/mpreal') diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h index 8404f1ff8..d68d7a78f 100644 --- a/unsupported/test/mpreal/mpreal.h +++ b/unsupported/test/mpreal/mpreal.h @@ -1,34 +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-2015 Pavel Holoborodko + Copyright (c) 2008-2016 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, Arash Partow, - Rodney James, Jorge Leitao. + Rodney James, Jorge Leitao, Jerome Benoit. 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 @@ -58,6 +58,7 @@ #include #include #include +#include // Options #define MPREAL_HAVE_MSVC_DEBUGVIEW // Enable Debugger Visualizer for "Debug" builds in MSVC. @@ -68,16 +69,18 @@ // Library version #define MPREAL_VERSION_MAJOR 3 #define MPREAL_VERSION_MINOR 6 -#define MPREAL_VERSION_PATCHLEVEL 2 -#define MPREAL_VERSION_STRING "3.6.2" +#define MPREAL_VERSION_PATCHLEVEL 5 +#define MPREAL_VERSION_STRING "3.6.5" // Detect compiler using signatures from http://predef.sourceforge.net/ -#if defined(__GNUC__) - #define IsInf(x) (isinf)(x) // GNU C++/Intel ICC compiler on Linux -#elif defined(_MSC_VER) // Microsoft Visual C++ - #define IsInf(x) (!_finite(x)) +#if defined(__GNUC__) && defined(__INTEL_COMPILER) + #define IsInf(x) isinf EIGEN_NOT_A_MACRO (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 EIGEN_NOT_A_MACRO (x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance #endif // A Clang feature extension to determine compiler features. @@ -85,22 +88,26 @@ #define __has_feature(x) 0 #endif -// Detect support for r-value references (move semantic). Borrowed from Eigen. +// Detect support for r-value references (move semantic). +// Move semantic should be enabled with great care in multi-threading environments, +// especially if MPFR uses custom memory allocators. +// Everything should be thread-safe and support passing ownership over thread boundary. #if (__has_feature(cxx_rvalue_references) || \ defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \ - (defined(_MSC_VER) && _MSC_VER >= 1600)) + (defined(_MSC_VER) && _MSC_VER >= 1600) && !defined(MPREAL_DISABLE_MOVE_SEMANTIC)) #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__) && __GNUC_MINOR__ >= 5) || __cplusplus >= 201103L || \ - (defined(_MSC_VER) && _MSC_VER >= 1800)) + (defined(__GXX_EXPERIMENTAL_CXX0X__) && __GNUC_MINOR >= 5) || __cplusplus >= 201103L || \ + (defined(_MSC_VER) && _MSC_VER >= 1800) || \ + (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1300)) #define MPREAL_HAVE_EXPLICIT_CONVERTERS #endif @@ -111,8 +118,8 @@ #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 @@ -122,12 +129,12 @@ #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 +// (a) uses low-level data members, might not be forward compatible // (b) sign is not set, add (x)->_mpfr_sign = 1; #define mpfr_set_zero_fast(x) ((x)->_mpfr_exp = __MPFR_EXP_ZERO) @@ -142,19 +149,19 @@ 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(); } + inline static mp_prec_t get_default_prec() { return (mpfr_get_default_prec)(); } // 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 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()); @@ -163,15 +170,15 @@ public: 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); + // 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); @@ -180,7 +187,7 @@ public: // Operations // = - // +, -, *, /, ++, --, <<, >> + // +, -, *, /, ++, --, <<, >> // *=, +=, -=, /=, // <, >, ==, <=, >= @@ -190,7 +197,7 @@ 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); @@ -224,7 +231,7 @@ public: const mpreal operator+() const; mpreal& operator++ (); - const mpreal operator++ (int); + const mpreal operator++ (int); // - mpreal& operator-=(const mpreal& v); @@ -242,7 +249,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); // * @@ -255,7 +262,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); @@ -343,17 +350,19 @@ 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); friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode); + friend const mpreal nextpow2(const mpreal& v, mp_rnd_t rnd_mode); + friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode); friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode); @@ -395,17 +404,17 @@ public: 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[], const unsigned long int n, int& status, mp_rnd_t rnd_mode); - friend int sgn(const mpreal& v); // returns -1 or +1 + friend int sgn (const mpreal& v); // MPFR 2.4.0 Specifics #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) @@ -438,7 +447,7 @@ public: // 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 @@ -467,14 +476,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)) @@ -482,7 +491,7 @@ public: #endif // Instance Checkers - friend bool (isnan) (const mpreal& v); + friend bool isnan EIGEN_NOT_A_MACRO (const mpreal& v); friend bool (isinf) (const mpreal& v); friend bool (isfinite) (const mpreal& v); @@ -501,15 +510,15 @@ 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()); //Exponent - mp_exp_t get_exp(); + mp_exp_t get_exp() const; 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()); @@ -532,7 +541,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); @@ -542,7 +551,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 @@ -561,15 +570,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()); +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()); @@ -580,7 +589,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 pointer to actual data + mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds null-pointer (in uninitialized state) mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); MPREAL_MSVC_DEBUGVIEW_CODE; @@ -588,9 +597,11 @@ inline mpreal::mpreal(mpreal&& other) inline mpreal& mpreal::operator=(mpreal&& other) { - mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); - - MPREAL_MSVC_DEBUGVIEW_CODE; + if (this != &other) + { + mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); // destructor for "other" will be called just afterwards + MPREAL_MSVC_DEBUGVIEW_CODE; + } return *this; } #endif @@ -639,20 +650,20 @@ inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) mpfr_init2(mpfr_ptr(), prec); #if (MPREAL_DOUBLE_BITS_OVERFLOW > -1) - if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW)) - { - mpfr_set_d(mpfr_ptr(), u, mode); - }else - throw conversion_overflow(); + if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW)) + { + mpfr_set_d(mpfr_ptr(), u, mode); + }else + throw conversion_overflow(); #else - mpfr_set_d(mpfr_ptr(), u, mode); + mpfr_set_d(mpfr_ptr(), u, mode); #endif MPREAL_MSVC_DEBUGVIEW_CODE; } 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); @@ -660,7 +671,7 @@ inline mpreal::mpreal(const long double 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_uj(mpfr_ptr(), u, mode); @@ -668,7 +679,7 @@ inline mpreal::mpreal(const unsigned long long int u, mp_prec_t prec, mp_rnd_t m } 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); @@ -676,7 +687,7 @@ inline mpreal::mpreal(const long 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_ui(mpfr_ptr(), u, mode); @@ -684,7 +695,7 @@ inline mpreal::mpreal(const unsigned long 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_ui(mpfr_ptr(), u, mode); @@ -692,7 +703,7 @@ inline mpreal::mpreal(const unsigned int 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_si(mpfr_ptr(), u, mode); @@ -700,7 +711,7 @@ inline mpreal::mpreal(const long int 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_si(mpfr_ptr(), u, mode); @@ -710,7 +721,7 @@ inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) 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; } @@ -718,7 +729,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; } @@ -726,15 +737,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{ @@ -742,55 +753,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;}; - template <> struct result_type {typedef mpreal type;}; - template <> struct result_type {typedef mpreal type;}; + 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; } ////////////////////////////////////////////////////////////////////////// @@ -840,17 +851,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()); @@ -867,9 +878,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()); @@ -882,7 +893,7 @@ 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); ////////////////////////////////////////////////////////////////////////// @@ -908,13 +919,13 @@ inline mpreal& mpreal::operator=(const mpreal& v) { if (this != &v) { - mp_prec_t tp = mpfr_get_prec( mpfr_srcptr()); - mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr()); + mp_prec_t tp = mpfr_get_prec( mpfr_srcptr()); + mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr()); - if(tp != vp){ - clear(mpfr_ptr()); - mpfr_init2(mpfr_ptr(), vp); - } + if(tp != vp){ + clear(mpfr_ptr()); + mpfr_init2(mpfr_ptr(), vp); + } mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd()); @@ -926,7 +937,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; } @@ -934,7 +945,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; } @@ -947,73 +958,73 @@ 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)) - { - mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd()); - }else - throw conversion_overflow(); + if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW)) + { + mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd()); + }else + throw conversion_overflow(); #else - mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd()); + mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd()); #endif - MPREAL_MSVC_DEBUGVIEW_CODE; + MPREAL_MSVC_DEBUGVIEW_CODE; 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()); +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 long long int v) -{ - mpfr_set_uj(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 long int v) -{ - mpfr_set_sj(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()); +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; @@ -1034,7 +1045,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; } @@ -1057,7 +1068,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; } @@ -1065,7 +1076,7 @@ inline mpreal& mpreal::operator=(const std::string& s) return *this; } -template +template inline mpreal& mpreal::operator= (const std::complex& z) { return *this = z.real(); @@ -1103,9 +1114,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) @@ -1161,12 +1172,12 @@ inline const mpreal mpreal::operator+()const { return mpreal(*this); } inline const mpreal operator+(const mpreal& a, const mpreal& b) { - mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); - mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); - return c; + mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); + mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); + return c; } -inline mpreal& mpreal::operator++() +inline mpreal& mpreal::operator++() { return *this += 1; } @@ -1178,7 +1189,7 @@ inline const mpreal mpreal::operator++ (int) return x; } -inline mpreal& mpreal::operator--() +inline mpreal& mpreal::operator--() { return *this -= 1; } @@ -1215,9 +1226,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) @@ -1225,7 +1236,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; @@ -1269,9 +1280,9 @@ inline const mpreal mpreal::operator-()const inline const mpreal operator-(const mpreal& a, const mpreal& b) { - mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); - mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); - return c; + mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); + mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); + return c; } inline const mpreal operator-(const double b, const mpreal& a) @@ -1340,9 +1351,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) @@ -1350,7 +1361,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; @@ -1386,9 +1397,9 @@ inline mpreal& mpreal::operator*=(const int v) inline const mpreal operator*(const mpreal& a, const mpreal& b) { - mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); - mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); - return c; + mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); + mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); + return c; } ////////////////////////////////////////////////////////////////////////// @@ -1418,7 +1429,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) @@ -1426,7 +1437,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; @@ -1462,9 +1473,9 @@ inline mpreal& mpreal::operator/=(const int v) inline const mpreal operator/(const mpreal& a, const mpreal& b) { - mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_srcptr()), mpfr_get_prec(b.mpfr_srcptr()))); - mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); - return c; + mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_srcptr()), mpfr_get_prec(b.mpfr_srcptr()))); + mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); + return c; } inline const mpreal operator/(const unsigned long int b, const mpreal& a) @@ -1639,9 +1650,9 @@ inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode) ////////////////////////////////////////////////////////////////////////// //Relational operators -// WARNING: +// WARNING: // -// Please note that following checks for double-NaN are guaranteed to work only in IEEE math mode: +// 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) @@ -1659,7 +1670,7 @@ inline bool operator > (const mpreal& a, const double b ){ return ! 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 EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); } -// inline bool operator >= (const mpreal& a, const unsigned int b ){ return !isnan EIGEN_NOT_A_MACRO (isnan()a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); } +inline bool operator >= (const mpreal& a, const unsigned int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); } inline bool operator >= (const mpreal& a, const long int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) >= 0 ); } inline bool operator >= (const mpreal& a, const int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) >= 0 ); } inline bool operator >= (const mpreal& a, const long double b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(),b) >= 0 ); } @@ -1697,7 +1708,7 @@ inline bool operator != (const mpreal& a, const int b ){ return ! 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 isnan EIGEN_NOT_A_MACRO (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 ); } @@ -1705,7 +1716,7 @@ inline bool isint (const mpreal& op){ return (mpfr_integer_p(op.mpfr_srcpt #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 @@ -1762,21 +1773,21 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const std::ostringstream format; - int digits = (n >= 0) ? n : 1 + bits2digits(mpfr_get_prec(mpfr_srcptr())); - + int digits = (n >= 0) ? n : 2 + 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"; } @@ -1791,7 +1802,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; @@ -1808,7 +1819,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); @@ -1819,7 +1830,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); @@ -1832,7 +1843,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); @@ -1843,7 +1854,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); @@ -1870,7 +1881,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(); @@ -1931,14 +1942,9 @@ inline int bits2digits(mp_prec_t b) ////////////////////////////////////////////////////////////////////////// // Set/Get number properties -inline int sgn(const mpreal& op) -{ - return mpfr_sgn(op.mpfr_srcptr()); -} - inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode) { - mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), (sign < 0 ? 1 : 0), RoundingMode); + mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), sign < 0, RoundingMode); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } @@ -1955,14 +1961,14 @@ 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; @@ -1976,7 +1982,7 @@ inline mpreal& mpreal::setZero(int sign) #else mpfr_set_si(mpfr_ptr(), 0, (mpfr_get_default_rounding_mode)()); setSign(sign); -#endif +#endif MPREAL_MSVC_DEBUGVIEW_CODE; return *this; @@ -1993,7 +1999,7 @@ inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode) MPREAL_MSVC_DEBUGVIEW_CODE; } -inline mp_exp_t mpreal::get_exp () +inline mp_exp_t mpreal::get_exp () const { return mpfr_get_exp(mpfr_srcptr()); } @@ -2022,7 +2028,7 @@ inline const mpreal ldexp(const mpreal& v, mp_exp_t exp) mpreal x(v); // 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()); + mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd()); return x; } @@ -2038,7 +2044,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) { @@ -2059,7 +2065,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) @@ -2091,12 +2097,18 @@ inline bool signbit(const mpreal& x) return mpfr_signbit(x.mpfr_srcptr()); } +inline mpreal& setsignbit(mpreal& x, bool minus, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) +{ + mpfr_setsign(x.mpfr_ptr(), x.mpfr_srcptr(), minus, rnd_mode); + return x; +} + 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; } @@ -2159,7 +2171,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 ); } @@ -2182,7 +2194,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) @@ -2193,9 +2205,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_rootn_ui(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()) @@ -2270,6 +2282,16 @@ inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_r inline const mpreal bessely0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y0 ); } inline const mpreal bessely1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y1 ); } +inline const mpreal nextpow2(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) +{ + mpreal y(0, x.getPrecision()); + + if(!iszero(x)) + y = ceil(log2(abs(x,r),r)); + + return y; +} + inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision())); @@ -2284,8 +2306,46 @@ inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = return a; } -inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) +inline const mpreal hypot(const mpreal& a, const mpreal& b, const mpreal& c) +{ + if(isnan EIGEN_NOT_A_MACRO (a) || isnan EIGEN_NOT_A_MACRO (b) || isnan EIGEN_NOT_A_MACRO(c)) return mpreal().setNan(); + else + { + mpreal absa = abs(a), absb = abs(b), absc = abs(c); + mpreal w = (std::max)(absa, (std::max)(absb, absc)); + mpreal r; + + if (!iszero(w)) + { + mpreal iw = 1/w; + r = w * sqrt(sqr(absa*iw) + sqr(absb*iw) + sqr(absc*iw)); + } + + return r; + } +} + +inline const mpreal hypot(const mpreal& a, const mpreal& b, const mpreal& c, const mpreal& d) { + if(isnan EIGEN_NOT_A_MACRO (a) || isnan EIGEN_NOT_A_MACRO (b) || isnan EIGEN_NOT_A_MACRO (c) || isnan EIGEN_NOT_A_MACRO (d)) return mpreal().setNan(); + else + { + mpreal absa = abs(a), absb = abs(b), absc = abs(c), absd = abs(d); + mpreal w = (std::max)(absa, (std::max)(absb, (std::max)(absc, absd))); + mpreal r; + + if (!iszero(w)) + { + mpreal iw = 1/w; + r = w * sqrt(sqr(absa*iw) + sqr(absb*iw) + sqr(absc*iw) + sqr(absd*iw)); + } + + return r; + } +} + +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; @@ -2299,7 +2359,7 @@ inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t } inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t rnd_mode = mpreal::get_default_rnd()) + mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { mpreal x(0, prec); mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode); @@ -2338,9 +2398,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)); @@ -2353,9 +2413,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)); @@ -2368,8 +2428,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); @@ -2387,7 +2447,7 @@ inline const mpreal sum (const mpreal tab[], const unsigned long int n, int& sta mpreal x; status = mpfr_sum(x.mpfr_ptr(), (mpfr_ptr*)p, n, mode); - + delete [] p; return x; } @@ -2401,9 +2461,9 @@ inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode); } -inline const mpreal li2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) -{ - MPREAL_UNARY_MATH_FUNCTION_BODY(li2); +inline const mpreal li2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) +{ + MPREAL_UNARY_MATH_FUNCTION_BODY(li2); } inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) @@ -2415,16 +2475,16 @@ inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = m inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { (void)rnd_mode; - - /* + + /* m = mod(x,y) if y != 0, returns x - n*y where n = floor(x/y) The following are true by convention: - mod(x,0) is x - mod(x,x) is 0 - - mod(x,y) for x != y and y != 0 has the same sign as y. - + - mod(x,y) for x != y and y != 0 has the same sign as y. + */ if(iszero(y)) return x; @@ -2432,9 +2492,7 @@ inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = m mpreal m = x - floor(x / y) * y; - m.setSign(sgn(y)); // make sure result has the same sign as Y - - return m; + return copysign(abs(m),y); // make sure result has the same sign as Y } inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) @@ -2442,8 +2500,8 @@ inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal a; mp_prec_t yp, xp; - yp = y.get_prec(); - xp = x.get_prec(); + yp = y.get_prec(); + xp = x.get_prec(); a.set_prec(yp>xp?yp:xp); @@ -2543,7 +2601,16 @@ inline const mpreal frac (const mpreal& x, mp_rnd_t r = mpreal::get_defaul ////////////////////////////////////////////////////////////////////////// // Miscellaneous Functions -inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mp,b.mp); } +inline int sgn(const mpreal& op) +{ + // Please note, this is classic signum function which ignores sign of zero. + // Use signbit if you need sign of zero. + return mpfr_sgn(op.mpfr_srcptr()); +} + +////////////////////////////////////////////////////////////////////////// +// Miscellaneous Functions +inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mpfr_ptr(),b.mpfr_ptr()); } inline const mpreal (max)(const mpreal& x, const mpreal& y){ return (x>y?x:y); } inline const mpreal (min)(const mpreal& x, const mpreal& y){ return (x= MPFR_VERSION_NUM(3,1,0)) +#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0) && MPFR_VERSION < MPFR_VERSION_NUM(4,0,0)) +// TODO: +// Use mpfr_nrandom since mpfr_grandom is deprecated +#if defined(_MSC_VER) +#pragma warning( push ) +#pragma warning( disable : 1478) +#endif 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; } +#if defined(_MSC_VER) +#pragma warning( pop ) +#endif inline const mpreal grandom(unsigned int seed = 0) { @@ -2664,17 +2740,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); @@ -2924,7 +3000,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); @@ -2981,11 +3057,11 @@ inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode) // Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap namespace std { - // we are allowed to extend namespace std with specializations only + // 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<> @@ -2996,7 +3072,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; @@ -3016,7 +3092,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); } @@ -3024,8 +3100,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(); } @@ -3036,17 +3112,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(); @@ -3054,9 +3130,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; } } @@ -3083,13 +3159,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; -- cgit v1.2.3