Make file formatting comply with POSIX and Unix standards

UTF-8, LF, no BOM, and newlines at the end of files

4 files changed, 3187 insertions, 3187 deletions

diff --git a/unsupported/doc/examples/SYCL/CwiseMul.cpp b/unsupported/doc/examples/SYCL/CwiseMul.cpp
index 31eb104c6..a7c33140e 100644
--- a/unsupported/doc/examples/SYCL/CwiseMul.cpp
+++ b/unsupported/doc/examples/SYCL/CwiseMul.cpp
@@ -60,4 +60,4 @@ int main()
     }
   }
   printf("c=a*b Done\n");
-}
\ No newline at end of file
+}
diff --git a/unsupported/test/cxx11_tensor_block_io.cpp b/unsupported/test/cxx11_tensor_block_io.cpp
index b8600eaea..a6f7a44cd 100644
--- a/unsupported/test/cxx11_tensor_block_io.cpp
+++ b/unsupported/test/cxx11_tensor_block_io.cpp
@@ -434,4 +434,4 @@ EIGEN_DECLARE_TEST(cxx11_tensor_block_io) {
   CALL_SUBTEST(test_block_io_squeeze_ones<RowMajor>());
   CALL_SUBTEST(test_block_io_squeeze_ones<ColMajor>());
   // clang-format on
-}
\ No newline at end of file
+}
diff --git a/unsupported/test/cxx11_tensor_trace.cpp b/unsupported/test/cxx11_tensor_trace.cpp
index 77e5dca20..009722895 100644
--- a/unsupported/test/cxx11_tensor_trace.cpp
+++ b/unsupported/test/cxx11_tensor_trace.cpp
@@ -169,4 +169,4 @@ EIGEN_DECLARE_TEST(cxx11_tensor_trace) {
   CALL_SUBTEST(test_simple_trace<RowMajor>());
   CALL_SUBTEST(test_trace_in_expr<ColMajor>());
   CALL_SUBTEST(test_trace_in_expr<RowMajor>());
-}
\ No newline at end of file
+}
diff --git a/unsupported/test/mpreal/mpreal.h b/unsupported/test/mpreal/mpreal.h
index 2b4f8e80f..5cfd66a10 100644
--- a/unsupported/test/mpreal/mpreal.h
+++ b/unsupported/test/mpreal/mpreal.h
@@ -1,3184 +1,3184 @@
-/*

-    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-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, 

-    Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood,

-    Petr Aleksandrov, Orion Poplawski, Charles Karney, Arash Partow,

-    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 

-        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 

-        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

-    it under the terms of the GNU General Public License as published by

-    the Free Software Foundation, either version 3 of the License, or

-    (at your option) any later version.

-

-    This program is distributed in the hope that it will be useful,

-    but WITHOUT ANY WARRANTY; without even the implied warranty of

-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

-    GNU General Public License for more details.

-

-    You should have received a copy of the GNU General Public License

-    along with this program.  If not, see <http://www.gnu.org/licenses/>.

-*/

-

-#ifndef __MPREAL_H__

-#define __MPREAL_H__

-

-#include <string>

-#include <iostream>

-#include <sstream>

-#include <stdexcept>

-#include <cfloat>

-#include <cmath>

-#include <cstring>

-#include <limits>

-#include <complex>

-#include <algorithm>

-#include <stdint.h>

-

-// Options

-#define MPREAL_HAVE_MSVC_DEBUGVIEW              // Enable Debugger Visualizer for "Debug" builds in MSVC.

-#define MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS  // Enable extended std::numeric_limits<mpfr::mpreal> specialization.

-                                                // Meaning that "digits", "round_style" and similar members are defined as functions, not constants.

-                                                // See std::numeric_limits<mpfr::mpreal> at the end of the file for more information.

-

-// Library version

-#define MPREAL_VERSION_MAJOR 3

-#define MPREAL_VERSION_MINOR 6

-#define MPREAL_VERSION_PATCHLEVEL 5

-#define MPREAL_VERSION_STRING "3.6.5"

-

-// Detect compiler using signatures from http://predef.sourceforge.net/

-#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 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.

-#ifndef __has_feature

-    #define __has_feature(x) 0

-#endif

-

-// 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(MPREAL_DISABLE_MOVE_SEMANTIC))

-

-    #define MPREAL_HAVE_MOVE_SUPPORT

-

-    // 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. 

-#if (__has_feature(cxx_explicit_conversions) || \

-       (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

-

-#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 

-#endif

-

-#include <mpfr.h>

-

-#if (MPFR_VERSION < MPFR_VERSION_NUM(3,0,0))

-    #include <cstdlib>                          // Needed for random()

-#endif

-

-// Less important options

-#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 forward compatible

-// (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

-  #define MPREAL_PERMISSIVE_EXPR

-#endif

-

-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)();                          }

-

-    // 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 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());

-    

-    // 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();                           

-

-#ifdef MPREAL_HAVE_MOVE_SUPPORT

-    mpreal& operator=(mpreal&& v);

-    mpreal(mpreal&& u);

-#endif

-

-    // Operations

-    // =

-    // +, -, *, /, ++, --, <<, >> 

-    // *=, +=, -=, /=,

-    // <, >, ==, <=, >=

-

-    // =

-    mpreal& operator=(const mpreal& v);

-    mpreal& operator=(const mpf_t v);

-    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 unsigned long int v);

-    mpreal& operator=(const unsigned long long int v);

-    mpreal& operator=(const long long int v);

-    mpreal& operator=(const unsigned int v);

-    mpreal& operator=(const long int v);

-    mpreal& operator=(const int v);

-    mpreal& operator=(const char* s);

-    mpreal& operator=(const std::string& s);

-    template <typename real_t> mpreal& operator= (const std::complex<real_t>& z);

-

-    // +

-    mpreal& operator+=(const mpreal& v);

-    mpreal& operator+=(const mpf_t v);

-    mpreal& operator+=(const mpz_t v);

-    mpreal& operator+=(const mpq_t v);

-    mpreal& operator+=(const long double u);

-    mpreal& operator+=(const double u);

-    mpreal& operator+=(const unsigned long int u);

-    mpreal& operator+=(const unsigned int u);

-    mpreal& operator+=(const long int u);

-    mpreal& operator+=(const 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);

-    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); 

-

-    // -

-    mpreal& operator-=(const mpreal& v);

-    mpreal& operator-=(const mpz_t v);

-    mpreal& operator-=(const mpq_t v);

-    mpreal& operator-=(const long double u);

-    mpreal& operator-=(const double u);

-    mpreal& operator-=(const unsigned long int u);

-    mpreal& operator-=(const unsigned int u);

-    mpreal& operator-=(const long int u);

-    mpreal& operator-=(const int u);

-    const mpreal operator-() const;

-    friend const mpreal operator-(const unsigned long int b, const mpreal& a);

-    friend const mpreal operator-(const unsigned int b,      const mpreal& a);

-    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-- ();    

-    const mpreal  operator-- (int);

-

-    // *

-    mpreal& operator*=(const mpreal& v);

-    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 unsigned long int v);

-    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);

-    mpreal& operator/=(const mpq_t v);

-    mpreal& operator/=(const long double v);

-    mpreal& operator/=(const double v);

-    mpreal& operator/=(const unsigned long int v);

-    mpreal& operator/=(const unsigned int v);

-    mpreal& operator/=(const long int v);

-    mpreal& operator/=(const int v);

-    friend const mpreal operator/(const unsigned long int b, const mpreal& a);

-    friend const mpreal operator/(const unsigned int b,      const mpreal& a);

-    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);

-

-    //<<= Fast Multiplication by 2^u

-    mpreal& operator<<=(const unsigned long int u);

-    mpreal& operator<<=(const unsigned int u);

-    mpreal& operator<<=(const long int u);

-    mpreal& operator<<=(const int u);

-

-    //>>= Fast Division by 2^u

-    mpreal& operator>>=(const unsigned long int u);

-    mpreal& operator>>=(const unsigned int u);

-    mpreal& operator>>=(const long int u);

-    mpreal& operator>>=(const int u);

-

-    // Type Conversion operators

-    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 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

-    ::mpfr_ptr    mpfr_ptr();

-    ::mpfr_srcptr mpfr_ptr()    const;

-    ::mpfr_srcptr mpfr_srcptr() const;

-

-    // Convert mpreal to string with n significant digits in base b

-    // n = -1 -> convert with the maximum available digits

-    std::string toString(int n = -1, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const;

-

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    std::string toString(const std::string& format) const;

-#endif

-

-    std::ostream& output(std::ostream& os) const;

-

-    // Math Functions

-    friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode);

-    friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);

-    friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);

-    friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode);

-    friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode);

-    friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode);

-    friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode);

-    friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode);

-    friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode);

-

-    friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);

-    friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);

-    friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode);

-    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 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);

-    friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode);

-    friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);

-

-    friend const mpreal acos  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal asin  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal atan  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode);

-    friend const mpreal acot  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal asec  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal acsc  (const mpreal& v, mp_rnd_t rnd_mode);

-

-    friend const mpreal cosh  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal sinh  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal tanh  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal sech  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal csch  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal coth  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode);

-

-    friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);

-

-    friend const mpreal fac_ui (unsigned long int v,  mp_prec_t prec, mp_rnd_t rnd_mode);

-    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 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 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);

-

-// MPFR 2.4.0 Specifics

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    friend int          sinh_cosh   (mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal li2         (const mpreal& v,                       mp_rnd_t rnd_mode);

-    friend const mpreal fmod        (const mpreal& x, const mpreal& y,      mp_rnd_t rnd_mode);

-    friend const mpreal rec_sqrt    (const mpreal& v,                       mp_rnd_t rnd_mode);

-

-    // MATLAB's semantic equivalents

-    friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Remainder after division

-    friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Modulus after division

-#endif

-

-#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);

-

-    // 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);    

-

-    // Constants

-    // don't forget to call mpfr_free_cache() for every thread where you are using const-functions

-    friend const mpreal const_log2      (mp_prec_t prec, mp_rnd_t rnd_mode);

-    friend const mpreal const_pi        (mp_prec_t prec, mp_rnd_t rnd_mode);

-    friend const mpreal const_euler     (mp_prec_t prec, mp_rnd_t rnd_mode);

-    friend const mpreal const_catalan   (mp_prec_t prec, mp_rnd_t rnd_mode);

-

-    // returns +inf iff sign>=0 otherwise -inf

-    friend const mpreal const_infinity(int sign, mp_prec_t prec);

-

-    // Output/ Input

-    friend std::ostream& operator<<(std::ostream& os, const mpreal& v);

-    friend std::istream& operator>>(std::istream& is, mpreal& v);

-

-    // Integer Related Functions

-    friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal ceil (const mpreal& v);

-    friend const mpreal floor(const mpreal& v);

-    friend const mpreal round(const mpreal& v);

-    friend const mpreal trunc(const mpreal& v);

-    friend const mpreal rint_ceil   (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal rint_floor  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal rint_round  (const mpreal& v, mp_rnd_t rnd_mode);

-    friend const mpreal rint_trunc  (const mpreal& v, mp_rnd_t rnd_mode);

-    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); 

-

-// MPFR < 2.4.2 Specifics

-#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2))

-    friend const mpreal random2 (mp_size_t size, mp_exp_t exp);

-#endif

-

-    // Instance Checkers

-    friend bool isnan EIGEN_NOT_A_MACRO     (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);

-    friend bool isint    (const mpreal& v);

-

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))

-    friend bool isregular(const mpreal& v);

-#endif

-

-    // Set/Get instance properties

-    inline mp_prec_t    get_prec() const;

-    inline void         set_prec(mp_prec_t prec, mp_rnd_t rnd_mode = get_default_rnd());    // Change precision with rounding mode

-

-    // Aliases for get_prec(), set_prec() - needed for compatibility with std::complex<mpreal> 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&        setNan  ();

-    mpreal&        setZero (int Sign = +1);

-    mpreal&        setSign (int Sign, mp_rnd_t RoundingMode = get_default_rnd());

-

-    //Exponent

-    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());

-

-    // Inexact conversion from float

-    inline bool fits_in_bits(double x, int n);

-

-    // Set/Get global properties

-    static void            set_default_prec(mp_prec_t prec);

-    static void            set_default_rnd(mp_rnd_t rnd_mode);

-

-    static mp_exp_t  get_emin (void);

-    static mp_exp_t  get_emax (void);

-    static mp_exp_t  get_emin_min (void);

-    static mp_exp_t  get_emin_max (void);

-    static mp_exp_t  get_emax_min (void);

-    static mp_exp_t  get_emax_max (void);

-    static int       set_emin (mp_exp_t exp);

-    static int       set_emax (mp_exp_t exp);

-

-    // 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);

-

-private:

-    // Human friendly Debug Preview in Visual Studio.

-    // Put one of these lines:

-    //

-    // mpfr::mpreal=<DebugView>                              ; Show value only

-    // mpfr::mpreal=<DebugView>, <mp[0]._mpfr_prec,u>bits    ; Show value & precision

-    // 

-    // at the beginning of

-    // [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat

-    MPREAL_MSVC_DEBUGVIEW_DATA

-

-    // "Smart" resources deallocation. Checks if instance initialized before deletion.

-    void clear(::mpfr_ptr);

-};

-

-//////////////////////////////////////////////////////////////////////////

-// Exceptions

-class conversion_overflow : public std::exception {

-public:

-    std::string why() { return "inexact conversion from floating point"; }

-};

-

-//////////////////////////////////////////////////////////////////////////

-// 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_zero_fast(mpfr_ptr());

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-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());

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-#ifdef MPREAL_HAVE_MOVE_SUPPORT

-inline mpreal::mpreal(mpreal&& other)

-{

-    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;

-}

-

-inline mpreal& mpreal::operator=(mpreal&& other)

-{

-    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

-

-inline mpreal::mpreal(const mpfr_t  u, bool shared)

-{

-    if(shared)

-    {

-        std::memcpy(mpfr_ptr(), u, sizeof(mpfr_t));

-    }

-    else

-    {

-        mpfr_init2(mpfr_ptr(), mpfr_get_prec(u));

-        mpfr_set  (mpfr_ptr(), u, mpreal::get_default_rnd());

-    }

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-inline mpreal::mpreal(const mpf_t u)

-{

-    mpfr_init2(mpfr_ptr(),(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)

-    mpfr_set_f(mpfr_ptr(),u,mpreal::get_default_rnd());

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-inline mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)

-{

-    mpfr_init2(mpfr_ptr(), prec);

-    mpfr_set_z(mpfr_ptr(), u, mode);

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-inline mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)

-{

-    mpfr_init2(mpfr_ptr(), prec);

-    mpfr_set_q(mpfr_ptr(), u, mode);

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-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();

-#else

-	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);

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-inline mpreal::mpreal(const unsigned long long int u, mp_prec_t prec, mp_rnd_t mode)

-{ 

-    mpfr_init2 (mpfr_ptr(), prec);

-    mpfr_set_uj(mpfr_ptr(), u, mode);

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-inline mpreal::mpreal(const long long int u, mp_prec_t prec, mp_rnd_t mode)

-{ 

-    mpfr_init2 (mpfr_ptr(), prec);

-    mpfr_set_sj(mpfr_ptr(), u, mode);

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)

-{ 

-    mpfr_init2 (mpfr_ptr(), prec);

-    mpfr_set_ui(mpfr_ptr(), u, mode);

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-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);

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-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);

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-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);

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-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); 

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-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); 

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-inline void mpreal::clear(::mpfr_ptr x)

-{

-#ifdef MPREAL_HAVE_MOVE_SUPPORT

-    if(mpfr_is_initialized(x)) 

-#endif

-    mpfr_clear(x);

-}

-

-inline mpreal::~mpreal() 

-{ 

-    clear(mpfr_ptr());

-}                           

-

-// internal namespace needed for template magic

-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 <typename ArgumentType> struct result_type {};    

-    

-    template <> struct result_type<mpreal>              {typedef mpreal type;};    

-    template <> struct result_type<mpz_t>               {typedef mpreal type;};    

-    template <> struct result_type<mpq_t>               {typedef mpreal type;};    

-    template <> struct result_type<long double>         {typedef mpreal type;};    

-    template <> struct result_type<double>              {typedef mpreal type;};    

-    template <> struct result_type<unsigned long int>   {typedef mpreal type;};    

-    template <> struct result_type<unsigned int>        {typedef mpreal type;};    

-    template <> struct result_type<long int>            {typedef mpreal type;};    

-    template <> struct result_type<int>                 {typedef mpreal type;};    

-    template <> struct result_type<long long>           {typedef mpreal type;};    

-    template <> struct result_type<unsigned long long>  {typedef mpreal type;};    

-}

-

-// + Addition

-template <typename Rhs> 

-inline const typename internal::result_type<Rhs>::type 

-    operator+(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) += rhs;    }

-

-template <typename Lhs> 

-inline const typename internal::result_type<Lhs>::type 

-    operator+(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) += lhs;    } 

-

-// - Subtraction

-template <typename Rhs> 

-inline const typename internal::result_type<Rhs>::type 

-    operator-(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) -= rhs;    }

-

-template <typename Lhs> 

-inline const typename internal::result_type<Lhs>::type 

-    operator-(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) -= rhs;    }

-

-// * Multiplication

-template <typename Rhs> 

-inline const typename internal::result_type<Rhs>::type 

-    operator*(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) *= rhs;    }

-

-template <typename Lhs> 

-inline const typename internal::result_type<Lhs>::type 

-    operator*(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) *= lhs;    } 

-

-// / Division

-template <typename Rhs> 

-inline const typename internal::result_type<Rhs>::type 

-    operator/(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) /= rhs;    }

-

-template <typename Lhs> 

-inline const typename internal::result_type<Lhs>::type 

-    operator/(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) /= rhs;    }

-

-//////////////////////////////////////////////////////////////////////////

-// sqrt

-const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal sqrt(const long int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-// abs

-inline const mpreal abs(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd());

-

-//////////////////////////////////////////////////////////////////////////

-// pow

-const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-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 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 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 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());

-const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());

-

-//////////////////////////////////////////////////////////////////////////

-// Estimate machine epsilon for the given precision

-// Returns smallest eps such that 1.0 + eps != 1.0

-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);        

-

-// 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());

-inline mpreal maxval(mp_prec_t prec = mpreal::get_default_prec());

-

-// 'Dirty' equality check 1: |a-b| < min{|a|,|b|} * eps

-inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps);

-

-// 'Dirty' equality check 2: |a-b| < min{|a|,|b|} * eps( min{|a|,|b|} )

-inline bool isEqualFuzzy(const mpreal& a, const mpreal& b);

-

-// 'Bitwise' equality check

-//  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))

-

-inline mp_prec_t digits2bits(int d);

-inline int       bits2digits(mp_prec_t b);

-

-//////////////////////////////////////////////////////////////////////////

-// min, max

-const mpreal (max)(const mpreal& x, const mpreal& y);

-const mpreal (min)(const mpreal& x, const mpreal& y);

-

-//////////////////////////////////////////////////////////////////////////

-// Implementation

-//////////////////////////////////////////////////////////////////////////

-

-//////////////////////////////////////////////////////////////////////////

-// Operators - Assignment

-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());

-

-		if(tp != vp){

-			clear(mpfr_ptr());

-			mpfr_init2(mpfr_ptr(), vp);

-		}

-

-        mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());

-

-        MPREAL_MSVC_DEBUGVIEW_CODE;

-    }

-    return *this;

-}

-

-inline mpreal& mpreal::operator=(const mpf_t v)

-{

-    mpfr_set_f(mpfr_ptr(), v, mpreal::get_default_rnd());

-    

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator=(const mpz_t v)

-{

-    mpfr_set_z(mpfr_ptr(), v, mpreal::get_default_rnd());

-    

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator=(const mpq_t v)

-{

-    mpfr_set_q(mpfr_ptr(), v, mpreal::get_default_rnd());

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-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)                

-{   

-#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();

-#else

-	mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());

-#endif

-

-	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());    

-

-    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 long long int v)    

-{    

-    mpfr_set_uj(mpfr_ptr(), v, mpreal::get_default_rnd());    

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator=(const long long int v)    

-{    

-    mpfr_set_sj(mpfr_ptr(), v, mpreal::get_default_rnd());    

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator=(const long int v)            

-{    

-    mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());    

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator=(const int v)

-{    

-    mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());    

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator=(const char* s)

-{

-    // Use other converters for more precise control on base & precision & rounding:

-    //

-    //        mpreal(const char* s,        mp_prec_t prec, int base, mp_rnd_t mode)

-    //        mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode)

-    //

-    // Here we assume base = 10 and we use precision of target variable.

-

-    mpfr_t t;

-

-    mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));

-

-    if(0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd()))

-    {

-        mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); 

-        MPREAL_MSVC_DEBUGVIEW_CODE;

-    }

-

-    clear(t);

-    return *this;

-}

-

-inline mpreal& mpreal::operator=(const std::string& s)

-{

-    // Use other converters for more precise control on base & precision & rounding:

-    //

-    //        mpreal(const char* s,        mp_prec_t prec, int base, mp_rnd_t mode)

-    //        mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode)

-    //

-    // Here we assume base = 10 and we use precision of target variable.

-

-    mpfr_t t;

-

-    mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));

-

-    if(0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd()))

-    {

-        mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); 

-        MPREAL_MSVC_DEBUGVIEW_CODE;

-    }

-

-    clear(t);

-    return *this;

-}

-

-template <typename real_t> 

-inline mpreal& mpreal::operator= (const std::complex<real_t>& z)

-{

-    return *this = z.real();

-}

-

-//////////////////////////////////////////////////////////////////////////

-// + Addition

-inline mpreal& mpreal::operator+=(const mpreal& v)

-{

-    mpfr_add(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator+=(const mpf_t u)

-{

-    *this += mpreal(u);

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator+=(const mpz_t u)

-{

-    mpfr_add_z(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator+=(const mpq_t u)

-{

-    mpfr_add_q(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator+= (const long double u)

-{

-    *this += mpreal(u);    

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;    

-}

-

-inline mpreal& mpreal::operator+= (const double u)

-{

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    mpfr_add_d(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-#else

-    *this += mpreal(u);

-#endif

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator+=(const unsigned long int u)

-{

-    mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator+=(const unsigned int u)

-{

-    mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator+=(const long int u)

-{

-    mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator+=(const int u)

-{

-    mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    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 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); }

-

-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;

-}

-

-inline mpreal& mpreal::operator++() 

-{

-    return *this += 1;

-}

-

-inline const mpreal mpreal::operator++ (int)

-{

-    mpreal x(*this);

-    *this += 1;

-    return x;

-}

-

-inline mpreal& mpreal::operator--() 

-{

-    return *this -= 1;

-}

-

-inline const mpreal mpreal::operator-- (int)

-{

-    mpreal x(*this);

-    *this -= 1;

-    return x;

-}

-

-//////////////////////////////////////////////////////////////////////////

-// - Subtraction

-inline mpreal& mpreal::operator-=(const mpreal& v)

-{

-    mpfr_sub(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator-=(const mpz_t v)

-{

-    mpfr_sub_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator-=(const mpq_t v)

-{

-    mpfr_sub_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator-=(const long double v)

-{

-    *this -= mpreal(v);    

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;    

-}

-

-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);    

-#endif

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator-=(const unsigned long int v)

-{

-    mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator-=(const unsigned int v)

-{

-    mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator-=(const long int v)

-{

-    mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator-=(const int v)

-{

-    mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline const mpreal mpreal::operator-()const

-{

-    mpreal u(*this);

-    mpfr_neg(u.mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd());

-    return u;

-}

-

-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;

-}

-

-inline const mpreal operator-(const double  b, const mpreal& a)

-{

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

-    mpfr_d_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-#else

-    mpreal x(b, mpfr_get_prec(a.mpfr_ptr()));

-    x -= a;

-    return x;

-#endif

-}

-

-inline const mpreal operator-(const unsigned long int b, const mpreal& a)

-{

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

-    mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-}

-

-inline const mpreal operator-(const unsigned int b, const mpreal& a)

-{

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

-    mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-}

-

-inline const mpreal operator-(const long int b, const mpreal& a)

-{

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

-    mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-}

-

-inline const mpreal operator-(const int b, const mpreal& a)

-{

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

-    mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-}

-

-//////////////////////////////////////////////////////////////////////////

-// * Multiplication

-inline mpreal& mpreal::operator*= (const mpreal& v)

-{

-    mpfr_mul(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator*=(const mpz_t v)

-{

-    mpfr_mul_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator*=(const mpq_t v)

-{

-    mpfr_mul_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator*=(const long double v)

-{

-    *this *= mpreal(v);    

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;    

-}

-

-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);    

-#endif

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator*=(const unsigned long int v)

-{

-    mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator*=(const unsigned int v)

-{

-    mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator*=(const long int v)

-{

-    mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator*=(const int v)

-{

-    mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *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_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());

-	return c;

-}

-

-//////////////////////////////////////////////////////////////////////////

-// / Division

-inline mpreal& mpreal::operator/=(const mpreal& v)

-{

-    mpfr_div(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator/=(const mpz_t v)

-{

-    mpfr_div_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator/=(const mpq_t v)

-{

-    mpfr_div_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator/=(const long double v)

-{

-    *this /= mpreal(v);

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;    

-}

-

-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);    

-#endif

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator/=(const unsigned long int v)

-{

-    mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator/=(const unsigned int v)

-{

-    mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator/=(const long int v)

-{

-    mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator/=(const int v)

-{

-    mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-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;

-}

-

-inline const mpreal operator/(const unsigned long int b, const mpreal& a)

-{

-    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

-    mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-}

-

-inline const mpreal operator/(const unsigned int b, const mpreal& a)

-{

-    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

-    mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-}

-

-inline const mpreal operator/(const long int b, const mpreal& a)

-{

-    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

-    mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-}

-

-inline const mpreal operator/(const int b, const mpreal& a)

-{

-    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

-    mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-}

-

-inline const mpreal operator/(const double  b, const mpreal& a)

-{

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));

-    mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());

-    return x;

-#else

-    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));

-    x /= a;

-    return x;

-#endif

-}

-

-//////////////////////////////////////////////////////////////////////////

-// Shifts operators - Multiplication/Division by power of 2

-inline mpreal& mpreal::operator<<=(const unsigned long int u)

-{

-    mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator<<=(const unsigned int u)

-{

-    mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator<<=(const long int u)

-{

-    mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator<<=(const int u)

-{

-    mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator>>=(const unsigned long int u)

-{

-    mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator>>=(const unsigned int u)

-{

-    mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator>>=(const long int u)

-{

-    mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::operator>>=(const int u)

-{

-    mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline const mpreal operator<<(const mpreal& v, const unsigned long int k)

-{

-    return mul_2ui(v,k);

-}

-

-inline const mpreal operator<<(const mpreal& v, const unsigned int k)

-{

-    return mul_2ui(v,static_cast<unsigned long int>(k));

-}

-

-inline const mpreal operator<<(const mpreal& v, const long int k)

-{

-    return mul_2si(v,k);

-}

-

-inline const mpreal operator<<(const mpreal& v, const int k)

-{

-    return mul_2si(v,static_cast<long int>(k));

-}

-

-inline const mpreal operator>>(const mpreal& v, const unsigned long int k)

-{

-    return div_2ui(v,k);

-}

-

-inline const mpreal operator>>(const mpreal& v, const long int k)

-{

-    return div_2si(v,k);

-}

-

-inline const mpreal operator>>(const mpreal& v, const unsigned int k)

-{

-    return div_2ui(v,static_cast<unsigned long int>(k));

-}

-

-inline const mpreal operator>>(const mpreal& v, const int k)

-{

-    return div_2si(v,static_cast<long int>(k));

-}

-

-// mul_2ui

-inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)

-{

-    mpreal x(v);

-    mpfr_mul_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);

-    return x;

-}

-

-// mul_2si

-inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)

-{

-    mpreal x(v);

-    mpfr_mul_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);

-    return x;

-}

-

-inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)

-{

-    mpreal x(v);

-    mpfr_div_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);

-    return x;

-}

-

-inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)

-{

-    mpreal x(v);

-    mpfr_div_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);

-    return x;

-}

-

-//////////////////////////////////////////////////////////////////////////

-//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 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 (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 );    }

-inline bool operator >  (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 (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 );   }

-inline bool operator >= (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 (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 );    }

-inline bool operator <  (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 (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 );   }

-inline bool operator <= (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 (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 );   }

-inline bool operator == (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 );    }

-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 

-

-//////////////////////////////////////////////////////////////////////////

-// Type Converters

-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; }

-inline ::mpfr_srcptr  mpreal::mpfr_srcptr() const    { return mp; }

-

-template <class T>

-inline std::string toString(T t, std::ios_base & (*f)(std::ios_base&))

-{

-    std::ostringstream oss;

-    oss << f << t;

-    return oss.str();

-}

-

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-

-inline std::string mpreal::toString(const std::string& format) const

-{

-    char *s = NULL;

-    std::string out;

-

-    if( !format.empty() )

-    {

-        if(!(mpfr_asprintf(&s, format.c_str(), mpfr_srcptr()) < 0))

-        {

-            out = std::string(s);

-

-            mpfr_free_str(s);

-        }

-    }

-

-    return out;

-}

-

-#endif

-

-inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const

-{

-    // TODO: Add extended format specification (f, e, rounding mode) as it done in output operator

-    (void)b;

-    (void)mode;

-

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-

-    std::ostringstream format;

-

-    int digits = (n >= 0) ? n : 2 + bits2digits(mpfr_get_prec(mpfr_srcptr()));

-    

-    format << "%." << digits << "RNg";

-

-    return toString(format.str());

-

-#else

-

-    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";

-    }

-

-    if(mpfr_zero_p(mp)) return "0";

-    if(mpfr_nan_p(mp))  return "NaN";

-

-    s  = mpfr_get_str(NULL, &exp, b, 0, mp, mode);

-    ns = mpfr_get_str(NULL, &exp, b, (std::max)(0,n), mp, mode);

-

-    if(s!=NULL && ns!=NULL)

-    {

-        slen  = strlen(s);

-        nslen = strlen(ns);

-        if(nslen<=slen) 

-        {

-            mpfr_free_str(s);

-            s = ns;

-            slen = nslen;

-        }

-        else {

-            mpfr_free_str(ns);

-        }

-

-        // Make human eye-friendly formatting if possible

-        if (exp>0 && static_cast<size_t>(exp)<slen)

-        {

-            if(s[0]=='-')

-            {

-                // Remove zeros starting from right end

-                char* ptr = s+slen-1;

-                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);

-

-                //out = string(s,exp+1)+'.'+string(s+exp+1);

-            }

-            else

-            {

-                // Remove zeros starting from right end

-                char* ptr = s+slen-1;

-                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);

-

-                //out = string(s,exp)+'.'+string(s+exp);

-            }

-

-        }else{ // exp<0 || exp>slen

-            if(s[0]=='-')

-            {

-                // Remove zeros starting from right end

-                char* ptr = s+slen-1;

-                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);

-

-                //out = string(s,2)+'.'+string(s+2);

-            }

-            else

-            {

-                // Remove zeros starting from right end

-                char* ptr = s+slen-1;

-                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);

-

-                //out = string(s,1)+'.'+string(s+1);

-            }

-

-            // Make final string

-            if(--exp)

-            {

-                if(exp>0) out += "e+"+mpfr::toString<mp_exp_t>(exp,std::dec);

-                else       out += "e"+mpfr::toString<mp_exp_t>(exp,std::dec);

-            }

-        }

-

-        mpfr_free_str(s);

-        return out;

-    }else{

-        return "conversion error!";

-    }

-#endif

-}

-

-

-//////////////////////////////////////////////////////////////////////////

-// I/O

-inline std::ostream& mpreal::output(std::ostream& os) const 

-{

-    std::ostringstream format;

-    const std::ios::fmtflags flags = os.flags();

-

-    format << ((flags & std::ios::showpos) ? "%+" : "%");

-    if (os.precision() >= 0)

-        format << '.' << os.precision() << "R*"

-               << ((flags & std::ios::floatfield) == std::ios::fixed ? 'f' :

-                   (flags & std::ios::floatfield) == std::ios::scientific ? 'e' :

-                   'g');

-    else

-        format << "R*e";

-

-    char *s = NULL;

-    if(!(mpfr_asprintf(&s, format.str().c_str(),

-                        mpfr::mpreal::get_default_rnd(),

-                        mpfr_srcptr())

-        < 0))

-    {

-        os << std::string(s);

-        mpfr_free_str(s);

-    }

-    return os;

-}

-

-inline std::ostream& operator<<(std::ostream& os, const mpreal& v)

-{

-    return v.output(os);

-}

-

-inline std::istream& operator>>(std::istream &is, mpreal& v)

-{

-    // TODO: use cout::hexfloat and other flags to setup base

-    std::string tmp;

-    is >> tmp;

-    mpfr_set_str(v.mpfr_ptr(), tmp.c_str(), 10, mpreal::get_default_rnd());

-    return is;

-}

-

-//////////////////////////////////////////////////////////////////////////

-//     Bits - decimal digits relation

-//        bits   = ceil(digits*log[2](10))

-//        digits = floor(bits*log[10](2))

-

-inline mp_prec_t digits2bits(int d)

-{

-    const double LOG2_10 = 3.3219280948873624;

-

-    return mp_prec_t(std::ceil( d * LOG2_10 ));

-}

-

-inline int bits2digits(mp_prec_t b)

-{

-    const double LOG10_2 = 0.30102999566398119;

-

-    return int(std::floor( b * LOG10_2 ));

-}

-

-//////////////////////////////////////////////////////////////////////////

-// Set/Get number properties

-inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode)

-{

-    mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), sign < 0, RoundingMode);

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline int mpreal::getPrecision() const

-{

-    return int(mpfr_get_prec(mpfr_srcptr()));

-}

-

-inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode)

-{

-    mpfr_prec_round(mpfr_ptr(), Precision, RoundingMode);

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mpreal& mpreal::setInf(int sign) 

-{ 

-    mpfr_set_inf(mpfr_ptr(), sign);

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}    

-

-inline mpreal& mpreal::setNan() 

-{

-    mpfr_set_nan(mpfr_ptr());

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-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 

-

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return *this;

-}

-

-inline mp_prec_t mpreal::get_prec() const

-{

-    return mpfr_get_prec(mpfr_srcptr());

-}

-

-inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode)

-{

-    mpfr_prec_round(mpfr_ptr(),prec,rnd_mode);

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-}

-

-inline mp_exp_t mpreal::get_exp () const

-{

-    return mpfr_get_exp(mpfr_srcptr());

-}

-

-inline int mpreal::set_exp (mp_exp_t e)

-{

-    int x = mpfr_set_exp(mpfr_ptr(), e);

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return x;

-}

-

-inline const mpreal frexp(const mpreal& x, mp_exp_t* exp, mp_rnd_t mode = mpreal::get_default_rnd())

-{

-    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 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 */

-    return machine_epsilon(mpreal(1, prec));

-}

-

-inline mpreal machine_epsilon(const mpreal& x)

-{    

-    /* the smallest eps such that x + eps != x */

-    if( x < 0)

-    {

-        return nextabove(-x) + x;

-    }else{

-        return nextabove( x) - x;

-    }

-}

-

-// minval is 'safe' meaning 1 / minval does not overflow

-inline mpreal minval(mp_prec_t prec)

-{

-    /* min = 1/2 * 2^emin = 2^(emin - 1) */

-    return mpreal(1, prec) << mpreal::get_emin()-1;

-}

-

-// maxval is 'safe' meaning 1 / maxval does not underflow

-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(); 

-}

-

-inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps)

-{

-    return abs(a - b) <= machine_epsilon((max)(abs(a), abs(b))) * maxUlps;

-}

-

-inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps)

-{

-    return abs(a - b) <= eps;

-}

-

-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 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_trunc(n.mpfr_ptr(),v.mpfr_srcptr());

-    return f;

-}

-

-inline int mpreal::check_range (int t, mp_rnd_t rnd_mode)

-{

-    return mpfr_check_range(mpfr_ptr(),t,rnd_mode);

-}

-

-inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode)

-{

-    int r = mpfr_subnormalize(mpfr_ptr(),t,rnd_mode);

-    MPREAL_MSVC_DEBUGVIEW_CODE;

-    return r;

-}

-

-inline mp_exp_t mpreal::get_emin (void)

-{

-    return mpfr_get_emin();

-}

-

-inline int mpreal::set_emin (mp_exp_t exp)

-{

-    return mpfr_set_emin(exp);

-}

-

-inline mp_exp_t mpreal::get_emax (void)

-{

-    return mpfr_get_emax();

-}

-

-inline int mpreal::set_emax (mp_exp_t exp)

-{

-    return mpfr_set_emax(exp);

-}

-

-inline mp_exp_t mpreal::get_emin_min (void)

-{

-    return mpfr_get_emin_min();

-}

-

-inline mp_exp_t mpreal::get_emin_max (void)

-{

-    return mpfr_get_emin_max();

-}

-

-inline mp_exp_t mpreal::get_emax_min (void)

-{

-    return mpfr_get_emax_min();

-}

-

-inline mp_exp_t mpreal::get_emax_max (void)

-{

-    return mpfr_get_emax_max();

-}

-

-//////////////////////////////////////////////////////////////////////////

-// Mathematical Functions

-//////////////////////////////////////////////////////////////////////////

-#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; 

-

-inline const mpreal sqr  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())

-{   MPREAL_UNARY_MATH_FUNCTION_BODY(sqr );    }

-

-inline const mpreal sqrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())

-{   MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt);    }

-

-inline const mpreal sqrt(const unsigned long int x, mp_rnd_t r)

-{

-    mpreal y;

-    mpfr_sqrt_ui(y.mpfr_ptr(), x, r);

-    return y;

-}

-

-inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode)

-{

-    return sqrt(static_cast<unsigned long int>(v),rnd_mode);

-}

-

-inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode)

-{

-    if (v>=0)   return sqrt(static_cast<unsigned long int>(v),rnd_mode);

-    else        return mpreal().setNan(); // NaN  

-}

-

-inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode)

-{

-    if (v>=0)   return sqrt(static_cast<unsigned long int>(v),rnd_mode);

-    else        return mpreal().setNan(); // NaN

-}

-

-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())); 

-    #if (MPFR_VERSION >= MPFR_VERSION_NUM(4,0,0))

-    mpfr_rootn_ui(y.mpfr_ptr(), x.mpfr_srcptr(), k, r);

-    #else

-    mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r);

-    #endif

-    return y; 

-}

-

-inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r = mpreal::get_default_rnd())

-{

-    mpreal y(0, mpfr_get_prec(a.mpfr_srcptr()));

-    mpfr_dim(y.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), r);

-    return y;

-}

-

-inline int cmpabs(const mpreal& a,const mpreal& b)

-{

-    return mpfr_cmpabs(a.mpfr_ptr(), b.mpfr_srcptr());

-}

-

-inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    return mpfr_sin_cos(s.mpfr_ptr(), c.mpfr_ptr(), v.mpfr_srcptr(), rnd_mode);

-}

-

-inline const mpreal sqrt  (const long double v, mp_rnd_t rnd_mode)    {   return sqrt(mpreal(v),rnd_mode);    }

-inline const mpreal sqrt  (const double v, mp_rnd_t rnd_mode)         {   return sqrt(mpreal(v),rnd_mode);    }

-

-inline const mpreal cbrt  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt );    }

-inline const mpreal fabs  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(abs  );    }

-inline const mpreal abs   (const mpreal& x, mp_rnd_t r)                             {   MPREAL_UNARY_MATH_FUNCTION_BODY(abs  );    }

-inline const mpreal log   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log  );    }

-inline const mpreal log2  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log2 );    }

-inline const mpreal log10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log10);    }

-inline const mpreal exp   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp  );    }

-inline const mpreal exp2  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 );    }

-inline const mpreal exp10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp10);    }

-inline const mpreal cos   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cos  );    }

-inline const mpreal sin   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sin  );    }

-inline const mpreal tan   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(tan  );    }

-inline const mpreal sec   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sec  );    }

-inline const mpreal csc   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(csc  );    }

-inline const mpreal cot   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cot  );    }

-inline const mpreal acos  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(acos );    }

-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);                      }

-inline const mpreal acoth (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return atanh(1/v, r);                      }

-inline const mpreal asech (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return acosh(1/v, r);                      }

-inline const mpreal acsch (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return asinh(1/v, r);                      }

-

-inline const mpreal cosh  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cosh );    }

-inline const mpreal sinh  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sinh );    }

-inline const mpreal tanh  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(tanh );    }

-inline const mpreal sech  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sech );    }

-inline const mpreal csch  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(csch );    }

-inline const mpreal coth  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(coth );    }

-inline const mpreal acosh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(acosh);    }

-inline const mpreal asinh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(asinh);    }

-inline const mpreal atanh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(atanh);    }

-

-inline const mpreal log1p   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log1p  );    }

-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    );    }

-inline const mpreal erfc    (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(erfc   );    }

-inline const mpreal besselj0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(j0     );    }

-inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(j1     );    }

-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()));

-    mpfr_atan2(a.mpfr_ptr(), y.mpfr_srcptr(), x.mpfr_srcptr(), rnd_mode);

-    return a;

-}

-

-inline const mpreal hypot (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_hypot(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);

-    return a;

-}

-

-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;

-}

-

-inline const mpreal remquo (long* q, 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_remquo(a.mpfr_ptr(),q, x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);

-    return a;

-}

-

-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())

-{

-    mpreal x(0, prec);

-    mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode);

-    return x;

-}

-

-

-inline const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal x(v);

-    int tsignp;

-

-    if(signp)   mpfr_lgamma(x.mpfr_ptr(),  signp,v.mpfr_srcptr(),rnd_mode);

-    else        mpfr_lgamma(x.mpfr_ptr(),&tsignp,v.mpfr_srcptr(),rnd_mode);

-

-    return x;

-}

-

-

-inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())

-{

-    mpreal  y(0, x.getPrecision());

-    mpfr_jn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);

-    return y;

-}

-

-inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())

-{

-    mpreal  y(0, x.getPrecision());

-    mpfr_yn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);

-    return y;

-}

-

-inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal a;

-    mp_prec_t p1, p2, p3;

-

-    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));

-

-    mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);

-    return a;

-}

-

-inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal a;

-    mp_prec_t p1, p2, p3;

-

-    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));

-

-    mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);

-    return a;

-}

-

-inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal a;

-    mp_prec_t p1, p2;

-

-    p1 = v1.get_prec(); 

-    p2 = v2.get_prec(); 

-

-    a.set_prec(p1>p2?p1:p2);

-

-    mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode);

-

-    return a;

-}

-

-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;

-    status = mpfr_sum(x.mpfr_ptr(), (mpfr_ptr*)p, n, mode);

-    

-    delete [] p;

-    return x;

-}

-

-//////////////////////////////////////////////////////////////////////////

-// MPFR 2.4.0 Specifics

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))

-

-inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    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 rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    /*  R = rem(X,Y) if Y != 0, returns X - n * Y where n = trunc(X/Y). */

-    return fmod(x, y, rnd_mode);

-}

-

-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.    

-    

-    */

-

-    if(iszero(y)) return x;

-    if(x == y) return 0;

-

-    mpreal m = x - floor(x / y) * y;

-

-    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())

-{

-    mpreal a;

-    mp_prec_t yp, xp;

-

-    yp = y.get_prec(); 

-    xp = x.get_prec(); 

-

-    a.set_prec(yp>xp?yp:xp);

-

-    mpfr_fmod(a.mp, x.mp, y.mp, rnd_mode);

-

-    return a;

-}

-

-inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal x(v);

-    mpfr_rec_sqrt(x.mp,v.mp,rnd_mode);

-    return x;

-}

-#endif //  MPFR 2.4.0 Specifics

-

-//////////////////////////////////////////////////////////////////////////

-// MPFR 3.0.0 Specifics

-#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))

-inline const mpreal digamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(digamma);     }

-inline const mpreal ai      (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(ai);          }

-#endif // MPFR 3.0.0 Specifics

-

-//////////////////////////////////////////////////////////////////////////

-// Constants

-inline const mpreal const_log2 (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())

-{

-    mpreal x(0, p);

-    mpfr_const_log2(x.mpfr_ptr(), r);

-    return x;

-}

-

-inline const mpreal const_pi (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())

-{

-    mpreal x(0, p);

-    mpfr_const_pi(x.mpfr_ptr(), r);

-    return x;

-}

-

-inline const mpreal const_euler (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())

-{

-    mpreal x(0, p);

-    mpfr_const_euler(x.mpfr_ptr(), r);

-    return x;

-}

-

-inline const mpreal const_catalan (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())

-{

-    mpreal x(0, p);

-    mpfr_const_catalan(x.mpfr_ptr(), r);

-    return x;

-}

-

-inline const mpreal const_infinity (int sign = 1, mp_prec_t p = mpreal::get_default_prec())

-{

-    mpreal x(0, p);

-    mpfr_set_inf(x.mpfr_ptr(), sign);

-    return x;

-}

-

-//////////////////////////////////////////////////////////////////////////

-// Integer Related Functions

-inline const mpreal ceil(const mpreal& v)

-{

-    mpreal x(v);

-    mpfr_ceil(x.mp,v.mp);

-    return x;

-}

-

-inline const mpreal floor(const mpreal& v)

-{

-    mpreal x(v);

-    mpfr_floor(x.mp,v.mp);

-    return x;

-}

-

-inline const mpreal round(const mpreal& v)

-{

-    mpreal x(v);

-    mpfr_round(x.mp,v.mp);

-    return x;

-}

-

-inline const mpreal trunc(const mpreal& v)

-{

-    mpreal x(v);

-    mpfr_trunc(x.mp,v.mp);

-    return x;

-}

-

-inline const mpreal rint       (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint      );     }

-inline const mpreal rint_ceil  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil );     }

-inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor);     }

-inline const mpreal rint_round (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round);     }

-inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc);     }

-inline const mpreal frac       (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(frac      );     }

-

-//////////////////////////////////////////////////////////////////////////

-// Miscellaneous Functions

-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<y?x:y);       }

-

-inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal a;

-    mpfr_max(a.mp,x.mp,y.mp,rnd_mode);

-    return a;

-}

-

-inline const mpreal fmin(const mpreal& x, const mpreal& y,  mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal a;

-    mpfr_min(a.mp,x.mp,y.mp,rnd_mode);

-    return a;

-}

-

-inline const mpreal nexttoward (const mpreal& x, const mpreal& y)

-{

-    mpreal a(x);

-    mpfr_nexttoward(a.mp,y.mp);

-    return a;

-}

-

-inline const mpreal nextabove  (const mpreal& x)

-{

-    mpreal a(x);

-    mpfr_nextabove(a.mp);

-    return a;

-}

-

-inline const mpreal nextbelow  (const mpreal& x)

-{

-    mpreal a(x);

-    mpfr_nextbelow(a.mp);

-    return a;

-}

-

-inline const mpreal urandomb (gmp_randstate_t& state)

-{

-    mpreal x;

-    mpfr_urandomb(x.mpfr_ptr(),state);

-    return x;

-}

-

-#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.mpfr_ptr(), state, rnd_mode);

-    return x;

-}

-#endif

-

-#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2))

-inline const mpreal random2 (mp_size_t size, mp_exp_t exp)

-{

-    mpreal x;

-    mpfr_random2(x.mpfr_ptr(),size,exp);

-    return x;

-}

-#endif

-

-// Uniformly distributed random number generation

-// a = random(seed); <- initialization & first random number generation

-// a = random();     <- next random numbers generation

-// 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 initialize = true;

-

-    if(initialize)

-    {

-        gmp_randinit_default(state);

-        gmp_randseed_ui(state,0);

-        initialize = false;

-    }

-

-    if(seed != 0)    gmp_randseed_ui(state,seed);

-

-    return mpfr::urandom(state);

-#else

-    if(seed != 0)    std::srand(seed);

-    return mpfr::mpreal(std::rand()/(double)RAND_MAX);

-#endif

-

-}

-

-#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)

-{

-    static gmp_randstate_t state;

-    static bool initialize = true;

-

-    if(initialize)

-    {

-        gmp_randinit_default(state);

-        gmp_randseed_ui(state,0);

-        initialize = false;

-    }

-

-    if(seed != 0) gmp_randseed_ui(state,seed);

-

-    return mpfr::grandom(state);

-}

-#endif

-

-//////////////////////////////////////////////////////////////////////////

-// Set/Get global properties

-inline void mpreal::set_default_prec(mp_prec_t prec)

-{ 

-    mpfr_set_default_prec(prec); 

-}

-

-inline void mpreal::set_default_rnd(mp_rnd_t 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);

-}

-

-inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal x(a);

-    mpfr_pow(x.mp,x.mp,b.mp,rnd_mode);

-    return x;

-}

-

-inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal x(a);

-    mpfr_pow_z(x.mp,x.mp,b,rnd_mode);

-    return x;

-}

-

-inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal x(a);

-    mpfr_pow_ui(x.mp,x.mp,b,rnd_mode);

-    return x;

-}

-

-inline const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode)

-{

-    return pow(a,static_cast<unsigned long int>(b),rnd_mode);

-}

-

-inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal x(a);

-    mpfr_pow_si(x.mp,x.mp,b,rnd_mode);

-    return x;

-}

-

-inline const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode)

-{

-    return pow(a,static_cast<long int>(b),rnd_mode);

-}

-

-inline const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode)

-{

-    return pow(a,mpreal(b),rnd_mode);

-}

-

-inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode)

-{

-    return pow(a,mpreal(b),rnd_mode);

-}

-

-inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())

-{

-    mpreal x(a);

-    mpfr_ui_pow(x.mp,a,b.mp,rnd_mode);

-    return x;

-}

-

-inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode)

-{

-    return pow(static_cast<unsigned long int>(a),b,rnd_mode);

-}

-

-inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode)

-{

-    if (a>=0)     return pow(static_cast<unsigned long int>(a),b,rnd_mode);

-    else          return pow(mpreal(a),b,rnd_mode);

-}

-

-inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode)

-{

-    if (a>=0)     return pow(static_cast<unsigned long int>(a),b,rnd_mode);

-    else          return pow(mpreal(a),b,rnd_mode);

-}

-

-inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),b,rnd_mode);

-}

-

-inline const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),b,rnd_mode);

-}

-

-// pow unsigned long int

-inline const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode)

-{

-    mpreal x(a);

-    mpfr_ui_pow_ui(x.mp,a,b,rnd_mode);

-    return x;

-}

-

-inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode)

-{

-    return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-}

-

-inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode)

-{

-    if(b>0)    return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-    else       return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow

-}

-

-inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode)

-{

-    if(b>0)    return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-    else       return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow

-}

-

-inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode)

-{

-    return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow

-}

-

-inline const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode)

-{

-    return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow

-}

-

-// pow unsigned int

-inline const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode)

-{

-    return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui

-}

-

-inline const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode)

-{

-    return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-}

-

-inline const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode)

-{

-    if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-    else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-}

-

-inline const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode)

-{

-    if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-    else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-}

-

-inline const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode)

-{

-    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-}

-

-inline const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode)

-{

-    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-}

-

-// pow long int

-inline const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode)

-{

-    if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui

-    else     return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui

-}

-

-inline const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode)

-{

-    if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode);  //mpfr_ui_pow_ui

-    else     return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui

-}

-

-inline const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode)

-{

-    if (a>0)

-    {

-        if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-        else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-    }else{

-        return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si

-    }

-}

-

-inline const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode)

-{

-    if (a>0)

-    {

-        if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-        else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-    }else{

-        return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si

-    }

-}

-

-inline const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode)

-{

-    if (a>=0)   return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-    else        return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow

-}

-

-inline const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode)

-{

-    if (a>=0)   return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-    else        return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow

-}

-

-// pow int

-inline const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode)

-{

-    if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui

-    else     return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui

-}

-

-inline const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode)

-{

-    if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode);  //mpfr_ui_pow_ui

-    else     return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui

-}

-

-inline const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode)

-{

-    if (a>0)

-    {

-        if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-        else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-    }else{

-        return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si

-    }

-}

-

-inline const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode)

-{

-    if (a>0)

-    {

-        if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui

-        else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-    }else{

-        return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si

-    }

-}

-

-inline const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode)

-{

-    if (a>=0)   return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-    else        return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow

-}

-

-inline const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode)

-{

-    if (a>=0)   return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow

-    else        return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow

-}

-

-// 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);

-}

-

-inline const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui

-}

-

-inline const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui

-}

-

-inline const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si

-}

-

-inline const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si

-}

-

-inline const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),mpreal(b),rnd_mode);

-}

-

-inline const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),b,rnd_mode); // mpfr_pow_ui

-}

-

-inline const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); // mpfr_pow_ui

-}

-

-inline const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si

-}

-

-inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode)

-{

-    return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si

-}

-} // End of mpfr namespace

-

-// Explicit specialization of std::swap for mpreal numbers

-// Thus standard algorithms will use efficient version of swap (due to Koenig lookup)

-// 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

-    template <>

-    inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) 

-    { 

-        return mpfr::swap(x, y); 

-    }

-

-    template<>

-    class numeric_limits<mpfr::mpreal>

-    {

-    public:

-        static const bool is_specialized    = true;

-        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 bool has_infinity      = true;

-        static const bool has_quiet_NaN     = true;

-        static const bool has_signaling_NaN = true;

-

-        static const bool is_iec559         = true;        // = IEEE 754

-        static const bool is_bounded        = true;

-        static const bool is_modulo         = false;

-        static const bool traps             = true;

-        static const bool tinyness_before   = true;

-

-        static const float_denorm_style has_denorm  = denorm_absent;

-

-        inline static mpfr::mpreal (min)    (mp_prec_t precision = mpfr::mpreal::get_default_prec()) {  return  mpfr::minval(precision);  }

-        inline static mpfr::mpreal (max)    (mp_prec_t precision = mpfr::mpreal::get_default_prec()) {  return  mpfr::maxval(precision);  }

-        inline static mpfr::mpreal lowest   (mp_prec_t precision = mpfr::mpreal::get_default_prec()) {  return -mpfr::maxval(precision);  }

-

-        // 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);  }

-

-        inline static mpfr::mpreal round_error(mp_prec_t precision = mpfr::mpreal::get_default_prec())

-        {

-            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);    

-        }

-

-        inline static const mpfr::mpreal infinity()         { return mpfr::const_infinity();     }

-        inline static const mpfr::mpreal quiet_NaN()        { return mpfr::mpreal().setNan();    }

-        inline static const mpfr::mpreal signaling_NaN()    { return mpfr::mpreal().setNan();    }

-        inline static const mpfr::mpreal denorm_min()       { return (min)();                    }

-

-        // 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); 

-

-#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. 

-        //

-        // This is preferable way for std::numeric_limits<mpfr::mpreal> 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. 

-        inline static float_round_style round_style()

-        {

-            mp_rnd_t r = mpfr::mpreal::get_default_rnd();

-

-            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; 

-            default: return round_indeterminate;

-            }

-        }

-

-        inline static int digits()                        {    return int(mpfr::mpreal::get_default_prec());    }

-        inline static int digits(const mpfr::mpreal& x)   {    return x.getPrecision();                         }

-

-        inline static int digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec())

-        {

-            return mpfr::bits2digits(precision);

-        }

-

-        inline static int digits10(const mpfr::mpreal& x)

-        {

-            return mpfr::bits2digits(x.getPrecision());

-        }

-

-        inline static int max_digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec())

-        {

-            return digits10(precision);

-        }

-#else

-        // Digits and round_style are NOT constants when it comes to mpreal.

-        // 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. 

-        //

-        // For example, if you use 256 bits of precision uniformly in your program, then:

-        // digits       = 256

-        // 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;

-        static const int digits       = 53;

-        static const int digits10     = 15;

-        static const int max_digits10 = 16;

-#endif

-    };

-

-}

-

-#endif /* __MPREAL_H__ */

+/*
+    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-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, 
+    Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood,
+    Petr Aleksandrov, Orion Poplawski, Charles Karney, Arash Partow,
+    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 
+        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 
+        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
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __MPREAL_H__
+#define __MPREAL_H__
+
+#include <string>
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+#include <cfloat>
+#include <cmath>
+#include <cstring>
+#include <limits>
+#include <complex>
+#include <algorithm>
+#include <stdint.h>
+
+// Options
+#define MPREAL_HAVE_MSVC_DEBUGVIEW              // Enable Debugger Visualizer for "Debug" builds in MSVC.
+#define MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS  // Enable extended std::numeric_limits<mpfr::mpreal> specialization.
+                                                // Meaning that "digits", "round_style" and similar members are defined as functions, not constants.
+                                                // See std::numeric_limits<mpfr::mpreal> at the end of the file for more information.
+
+// Library version
+#define MPREAL_VERSION_MAJOR 3
+#define MPREAL_VERSION_MINOR 6
+#define MPREAL_VERSION_PATCHLEVEL 5
+#define MPREAL_VERSION_STRING "3.6.5"
+
+// Detect compiler using signatures from http://predef.sourceforge.net/
+#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 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.
+#ifndef __has_feature
+    #define __has_feature(x) 0
+#endif
+
+// 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(MPREAL_DISABLE_MOVE_SEMANTIC))
+
+    #define MPREAL_HAVE_MOVE_SUPPORT
+
+    // 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. 
+#if (__has_feature(cxx_explicit_conversions) || \
+       (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
+
+#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 
+#endif
+
+#include <mpfr.h>
+
+#if (MPFR_VERSION < MPFR_VERSION_NUM(3,0,0))
+    #include <cstdlib>                          // Needed for random()
+#endif
+
+// Less important options
+#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 forward compatible
+// (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
+  #define MPREAL_PERMISSIVE_EXPR
+#endif
+
+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)();                          }
+
+    // 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 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());
+    
+    // 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();                           
+
+#ifdef MPREAL_HAVE_MOVE_SUPPORT
+    mpreal& operator=(mpreal&& v);
+    mpreal(mpreal&& u);
+#endif
+
+    // Operations
+    // =
+    // +, -, *, /, ++, --, <<, >> 
+    // *=, +=, -=, /=,
+    // <, >, ==, <=, >=
+
+    // =
+    mpreal& operator=(const mpreal& v);
+    mpreal& operator=(const mpf_t v);
+    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 unsigned long int v);
+    mpreal& operator=(const unsigned long long int v);
+    mpreal& operator=(const long long int v);
+    mpreal& operator=(const unsigned int v);
+    mpreal& operator=(const long int v);
+    mpreal& operator=(const int v);
+    mpreal& operator=(const char* s);
+    mpreal& operator=(const std::string& s);
+    template <typename real_t> mpreal& operator= (const std::complex<real_t>& z);
+
+    // +
+    mpreal& operator+=(const mpreal& v);
+    mpreal& operator+=(const mpf_t v);
+    mpreal& operator+=(const mpz_t v);
+    mpreal& operator+=(const mpq_t v);
+    mpreal& operator+=(const long double u);
+    mpreal& operator+=(const double u);
+    mpreal& operator+=(const unsigned long int u);
+    mpreal& operator+=(const unsigned int u);
+    mpreal& operator+=(const long int u);
+    mpreal& operator+=(const 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);
+    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); 
+
+    // -
+    mpreal& operator-=(const mpreal& v);
+    mpreal& operator-=(const mpz_t v);
+    mpreal& operator-=(const mpq_t v);
+    mpreal& operator-=(const long double u);
+    mpreal& operator-=(const double u);
+    mpreal& operator-=(const unsigned long int u);
+    mpreal& operator-=(const unsigned int u);
+    mpreal& operator-=(const long int u);
+    mpreal& operator-=(const int u);
+    const mpreal operator-() const;
+    friend const mpreal operator-(const unsigned long int b, const mpreal& a);
+    friend const mpreal operator-(const unsigned int b,      const mpreal& a);
+    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-- ();    
+    const mpreal  operator-- (int);
+
+    // *
+    mpreal& operator*=(const mpreal& v);
+    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 unsigned long int v);
+    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);
+    mpreal& operator/=(const mpq_t v);
+    mpreal& operator/=(const long double v);
+    mpreal& operator/=(const double v);
+    mpreal& operator/=(const unsigned long int v);
+    mpreal& operator/=(const unsigned int v);
+    mpreal& operator/=(const long int v);
+    mpreal& operator/=(const int v);
+    friend const mpreal operator/(const unsigned long int b, const mpreal& a);
+    friend const mpreal operator/(const unsigned int b,      const mpreal& a);
+    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);
+
+    //<<= Fast Multiplication by 2^u
+    mpreal& operator<<=(const unsigned long int u);
+    mpreal& operator<<=(const unsigned int u);
+    mpreal& operator<<=(const long int u);
+    mpreal& operator<<=(const int u);
+
+    //>>= Fast Division by 2^u
+    mpreal& operator>>=(const unsigned long int u);
+    mpreal& operator>>=(const unsigned int u);
+    mpreal& operator>>=(const long int u);
+    mpreal& operator>>=(const int u);
+
+    // Type Conversion operators
+    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 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
+    ::mpfr_ptr    mpfr_ptr();
+    ::mpfr_srcptr mpfr_ptr()    const;
+    ::mpfr_srcptr mpfr_srcptr() const;
+
+    // Convert mpreal to string with n significant digits in base b
+    // n = -1 -> convert with the maximum available digits
+    std::string toString(int n = -1, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const;
+
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+    std::string toString(const std::string& format) const;
+#endif
+
+    std::ostream& output(std::ostream& os) const;
+
+    // Math Functions
+    friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode);
+    friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);
+    friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);
+    friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode);
+    friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode);
+    friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode);
+    friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode);
+    friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode);
+    friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode);
+
+    friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode);
+    friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode);
+    friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode);
+    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 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);
+    friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode);
+    friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);
+
+    friend const mpreal acos  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal asin  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal atan  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode);
+    friend const mpreal acot  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal asec  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal acsc  (const mpreal& v, mp_rnd_t rnd_mode);
+
+    friend const mpreal cosh  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal sinh  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal tanh  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal sech  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal csch  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal coth  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode);
+
+    friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode);
+
+    friend const mpreal fac_ui (unsigned long int v,  mp_prec_t prec, mp_rnd_t rnd_mode);
+    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 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 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);
+
+// MPFR 2.4.0 Specifics
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+    friend int          sinh_cosh   (mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal li2         (const mpreal& v,                       mp_rnd_t rnd_mode);
+    friend const mpreal fmod        (const mpreal& x, const mpreal& y,      mp_rnd_t rnd_mode);
+    friend const mpreal rec_sqrt    (const mpreal& v,                       mp_rnd_t rnd_mode);
+
+    // MATLAB's semantic equivalents
+    friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Remainder after division
+    friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Modulus after division
+#endif
+
+#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);
+
+    // 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);    
+
+    // Constants
+    // don't forget to call mpfr_free_cache() for every thread where you are using const-functions
+    friend const mpreal const_log2      (mp_prec_t prec, mp_rnd_t rnd_mode);
+    friend const mpreal const_pi        (mp_prec_t prec, mp_rnd_t rnd_mode);
+    friend const mpreal const_euler     (mp_prec_t prec, mp_rnd_t rnd_mode);
+    friend const mpreal const_catalan   (mp_prec_t prec, mp_rnd_t rnd_mode);
+
+    // returns +inf iff sign>=0 otherwise -inf
+    friend const mpreal const_infinity(int sign, mp_prec_t prec);
+
+    // Output/ Input
+    friend std::ostream& operator<<(std::ostream& os, const mpreal& v);
+    friend std::istream& operator>>(std::istream& is, mpreal& v);
+
+    // Integer Related Functions
+    friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal ceil (const mpreal& v);
+    friend const mpreal floor(const mpreal& v);
+    friend const mpreal round(const mpreal& v);
+    friend const mpreal trunc(const mpreal& v);
+    friend const mpreal rint_ceil   (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal rint_floor  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal rint_round  (const mpreal& v, mp_rnd_t rnd_mode);
+    friend const mpreal rint_trunc  (const mpreal& v, mp_rnd_t rnd_mode);
+    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); 
+
+// MPFR < 2.4.2 Specifics
+#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2))
+    friend const mpreal random2 (mp_size_t size, mp_exp_t exp);
+#endif
+
+    // Instance Checkers
+    friend bool isnan EIGEN_NOT_A_MACRO     (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);
+    friend bool isint    (const mpreal& v);
+
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
+    friend bool isregular(const mpreal& v);
+#endif
+
+    // Set/Get instance properties
+    inline mp_prec_t    get_prec() const;
+    inline void         set_prec(mp_prec_t prec, mp_rnd_t rnd_mode = get_default_rnd());    // Change precision with rounding mode
+
+    // Aliases for get_prec(), set_prec() - needed for compatibility with std::complex<mpreal> 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&        setNan  ();
+    mpreal&        setZero (int Sign = +1);
+    mpreal&        setSign (int Sign, mp_rnd_t RoundingMode = get_default_rnd());
+
+    //Exponent
+    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());
+
+    // Inexact conversion from float
+    inline bool fits_in_bits(double x, int n);
+
+    // Set/Get global properties
+    static void            set_default_prec(mp_prec_t prec);
+    static void            set_default_rnd(mp_rnd_t rnd_mode);
+
+    static mp_exp_t  get_emin (void);
+    static mp_exp_t  get_emax (void);
+    static mp_exp_t  get_emin_min (void);
+    static mp_exp_t  get_emin_max (void);
+    static mp_exp_t  get_emax_min (void);
+    static mp_exp_t  get_emax_max (void);
+    static int       set_emin (mp_exp_t exp);
+    static int       set_emax (mp_exp_t exp);
+
+    // 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);
+
+private:
+    // Human friendly Debug Preview in Visual Studio.
+    // Put one of these lines:
+    //
+    // mpfr::mpreal=<DebugView>                              ; Show value only
+    // mpfr::mpreal=<DebugView>, <mp[0]._mpfr_prec,u>bits    ; Show value & precision
+    // 
+    // at the beginning of
+    // [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat
+    MPREAL_MSVC_DEBUGVIEW_DATA
+
+    // "Smart" resources deallocation. Checks if instance initialized before deletion.
+    void clear(::mpfr_ptr);
+};
+
+//////////////////////////////////////////////////////////////////////////
+// Exceptions
+class conversion_overflow : public std::exception {
+public:
+    std::string why() { return "inexact conversion from floating point"; }
+};
+
+//////////////////////////////////////////////////////////////////////////
+// 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_zero_fast(mpfr_ptr());
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+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());
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+#ifdef MPREAL_HAVE_MOVE_SUPPORT
+inline mpreal::mpreal(mpreal&& other)
+{
+    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;
+}
+
+inline mpreal& mpreal::operator=(mpreal&& other)
+{
+    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
+
+inline mpreal::mpreal(const mpfr_t  u, bool shared)
+{
+    if(shared)
+    {
+        std::memcpy(mpfr_ptr(), u, sizeof(mpfr_t));
+    }
+    else
+    {
+        mpfr_init2(mpfr_ptr(), mpfr_get_prec(u));
+        mpfr_set  (mpfr_ptr(), u, mpreal::get_default_rnd());
+    }
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+inline mpreal::mpreal(const mpf_t u)
+{
+    mpfr_init2(mpfr_ptr(),(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)
+    mpfr_set_f(mpfr_ptr(),u,mpreal::get_default_rnd());
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+inline mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
+{
+    mpfr_init2(mpfr_ptr(), prec);
+    mpfr_set_z(mpfr_ptr(), u, mode);
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+inline mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
+{
+    mpfr_init2(mpfr_ptr(), prec);
+    mpfr_set_q(mpfr_ptr(), u, mode);
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+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();
+#else
+	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);
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+inline mpreal::mpreal(const unsigned long long int u, mp_prec_t prec, mp_rnd_t mode)
+{ 
+    mpfr_init2 (mpfr_ptr(), prec);
+    mpfr_set_uj(mpfr_ptr(), u, mode);
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+inline mpreal::mpreal(const long long int u, mp_prec_t prec, mp_rnd_t mode)
+{ 
+    mpfr_init2 (mpfr_ptr(), prec);
+    mpfr_set_sj(mpfr_ptr(), u, mode);
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
+{ 
+    mpfr_init2 (mpfr_ptr(), prec);
+    mpfr_set_ui(mpfr_ptr(), u, mode);
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+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);
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+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);
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+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);
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+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); 
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+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); 
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+inline void mpreal::clear(::mpfr_ptr x)
+{
+#ifdef MPREAL_HAVE_MOVE_SUPPORT
+    if(mpfr_is_initialized(x)) 
+#endif
+    mpfr_clear(x);
+}
+
+inline mpreal::~mpreal() 
+{ 
+    clear(mpfr_ptr());
+}                           
+
+// internal namespace needed for template magic
+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 <typename ArgumentType> struct result_type {};    
+    
+    template <> struct result_type<mpreal>              {typedef mpreal type;};    
+    template <> struct result_type<mpz_t>               {typedef mpreal type;};    
+    template <> struct result_type<mpq_t>               {typedef mpreal type;};    
+    template <> struct result_type<long double>         {typedef mpreal type;};    
+    template <> struct result_type<double>              {typedef mpreal type;};    
+    template <> struct result_type<unsigned long int>   {typedef mpreal type;};    
+    template <> struct result_type<unsigned int>        {typedef mpreal type;};    
+    template <> struct result_type<long int>            {typedef mpreal type;};    
+    template <> struct result_type<int>                 {typedef mpreal type;};    
+    template <> struct result_type<long long>           {typedef mpreal type;};    
+    template <> struct result_type<unsigned long long>  {typedef mpreal type;};    
+}
+
+// + Addition
+template <typename Rhs> 
+inline const typename internal::result_type<Rhs>::type 
+    operator+(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) += rhs;    }
+
+template <typename Lhs> 
+inline const typename internal::result_type<Lhs>::type 
+    operator+(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) += lhs;    } 
+
+// - Subtraction
+template <typename Rhs> 
+inline const typename internal::result_type<Rhs>::type 
+    operator-(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) -= rhs;    }
+
+template <typename Lhs> 
+inline const typename internal::result_type<Lhs>::type 
+    operator-(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) -= rhs;    }
+
+// * Multiplication
+template <typename Rhs> 
+inline const typename internal::result_type<Rhs>::type 
+    operator*(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) *= rhs;    }
+
+template <typename Lhs> 
+inline const typename internal::result_type<Lhs>::type 
+    operator*(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) *= lhs;    } 
+
+// / Division
+template <typename Rhs> 
+inline const typename internal::result_type<Rhs>::type 
+    operator/(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) /= rhs;    }
+
+template <typename Lhs> 
+inline const typename internal::result_type<Lhs>::type 
+    operator/(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) /= rhs;    }
+
+//////////////////////////////////////////////////////////////////////////
+// sqrt
+const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal sqrt(const long int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
+// abs
+inline const mpreal abs(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd());
+
+//////////////////////////////////////////////////////////////////////////
+// pow
+const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
+const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
+const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
+const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
+const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
+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 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 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 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());
+const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
+inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd());
+
+//////////////////////////////////////////////////////////////////////////
+// Estimate machine epsilon for the given precision
+// Returns smallest eps such that 1.0 + eps != 1.0
+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);        
+
+// 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());
+inline mpreal maxval(mp_prec_t prec = mpreal::get_default_prec());
+
+// 'Dirty' equality check 1: |a-b| < min{|a|,|b|} * eps
+inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps);
+
+// 'Dirty' equality check 2: |a-b| < min{|a|,|b|} * eps( min{|a|,|b|} )
+inline bool isEqualFuzzy(const mpreal& a, const mpreal& b);
+
+// 'Bitwise' equality check
+//  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))
+
+inline mp_prec_t digits2bits(int d);
+inline int       bits2digits(mp_prec_t b);
+
+//////////////////////////////////////////////////////////////////////////
+// min, max
+const mpreal (max)(const mpreal& x, const mpreal& y);
+const mpreal (min)(const mpreal& x, const mpreal& y);
+
+//////////////////////////////////////////////////////////////////////////
+// Implementation
+//////////////////////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////////////////////
+// Operators - Assignment
+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());
+
+		if(tp != vp){
+			clear(mpfr_ptr());
+			mpfr_init2(mpfr_ptr(), vp);
+		}
+
+        mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());
+
+        MPREAL_MSVC_DEBUGVIEW_CODE;
+    }
+    return *this;
+}
+
+inline mpreal& mpreal::operator=(const mpf_t v)
+{
+    mpfr_set_f(mpfr_ptr(), v, mpreal::get_default_rnd());
+    
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator=(const mpz_t v)
+{
+    mpfr_set_z(mpfr_ptr(), v, mpreal::get_default_rnd());
+    
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator=(const mpq_t v)
+{
+    mpfr_set_q(mpfr_ptr(), v, mpreal::get_default_rnd());
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+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)                
+{   
+#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();
+#else
+	mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());
+#endif
+
+	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());    
+
+    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 long long int v)    
+{    
+    mpfr_set_uj(mpfr_ptr(), v, mpreal::get_default_rnd());    
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator=(const long long int v)    
+{    
+    mpfr_set_sj(mpfr_ptr(), v, mpreal::get_default_rnd());    
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator=(const long int v)            
+{    
+    mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());    
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator=(const int v)
+{    
+    mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd());    
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator=(const char* s)
+{
+    // Use other converters for more precise control on base & precision & rounding:
+    //
+    //        mpreal(const char* s,        mp_prec_t prec, int base, mp_rnd_t mode)
+    //        mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode)
+    //
+    // Here we assume base = 10 and we use precision of target variable.
+
+    mpfr_t t;
+
+    mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));
+
+    if(0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd()))
+    {
+        mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); 
+        MPREAL_MSVC_DEBUGVIEW_CODE;
+    }
+
+    clear(t);
+    return *this;
+}
+
+inline mpreal& mpreal::operator=(const std::string& s)
+{
+    // Use other converters for more precise control on base & precision & rounding:
+    //
+    //        mpreal(const char* s,        mp_prec_t prec, int base, mp_rnd_t mode)
+    //        mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode)
+    //
+    // Here we assume base = 10 and we use precision of target variable.
+
+    mpfr_t t;
+
+    mpfr_init2(t, mpfr_get_prec(mpfr_srcptr()));
+
+    if(0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd()))
+    {
+        mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); 
+        MPREAL_MSVC_DEBUGVIEW_CODE;
+    }
+
+    clear(t);
+    return *this;
+}
+
+template <typename real_t> 
+inline mpreal& mpreal::operator= (const std::complex<real_t>& z)
+{
+    return *this = z.real();
+}
+
+//////////////////////////////////////////////////////////////////////////
+// + Addition
+inline mpreal& mpreal::operator+=(const mpreal& v)
+{
+    mpfr_add(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator+=(const mpf_t u)
+{
+    *this += mpreal(u);
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator+=(const mpz_t u)
+{
+    mpfr_add_z(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator+=(const mpq_t u)
+{
+    mpfr_add_q(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator+= (const long double u)
+{
+    *this += mpreal(u);    
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;    
+}
+
+inline mpreal& mpreal::operator+= (const double u)
+{
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+    mpfr_add_d(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+#else
+    *this += mpreal(u);
+#endif
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator+=(const unsigned long int u)
+{
+    mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator+=(const unsigned int u)
+{
+    mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator+=(const long int u)
+{
+    mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator+=(const int u)
+{
+    mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    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 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); }
+
+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;
+}
+
+inline mpreal& mpreal::operator++() 
+{
+    return *this += 1;
+}
+
+inline const mpreal mpreal::operator++ (int)
+{
+    mpreal x(*this);
+    *this += 1;
+    return x;
+}
+
+inline mpreal& mpreal::operator--() 
+{
+    return *this -= 1;
+}
+
+inline const mpreal mpreal::operator-- (int)
+{
+    mpreal x(*this);
+    *this -= 1;
+    return x;
+}
+
+//////////////////////////////////////////////////////////////////////////
+// - Subtraction
+inline mpreal& mpreal::operator-=(const mpreal& v)
+{
+    mpfr_sub(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator-=(const mpz_t v)
+{
+    mpfr_sub_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator-=(const mpq_t v)
+{
+    mpfr_sub_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator-=(const long double v)
+{
+    *this -= mpreal(v);    
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;    
+}
+
+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);    
+#endif
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator-=(const unsigned long int v)
+{
+    mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator-=(const unsigned int v)
+{
+    mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator-=(const long int v)
+{
+    mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator-=(const int v)
+{
+    mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline const mpreal mpreal::operator-()const
+{
+    mpreal u(*this);
+    mpfr_neg(u.mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd());
+    return u;
+}
+
+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;
+}
+
+inline const mpreal operator-(const double  b, const mpreal& a)
+{
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+    mpfr_d_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+#else
+    mpreal x(b, mpfr_get_prec(a.mpfr_ptr()));
+    x -= a;
+    return x;
+#endif
+}
+
+inline const mpreal operator-(const unsigned long int b, const mpreal& a)
+{
+    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+    mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+}
+
+inline const mpreal operator-(const unsigned int b, const mpreal& a)
+{
+    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+    mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+}
+
+inline const mpreal operator-(const long int b, const mpreal& a)
+{
+    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+    mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+}
+
+inline const mpreal operator-(const int b, const mpreal& a)
+{
+    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+    mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+}
+
+//////////////////////////////////////////////////////////////////////////
+// * Multiplication
+inline mpreal& mpreal::operator*= (const mpreal& v)
+{
+    mpfr_mul(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator*=(const mpz_t v)
+{
+    mpfr_mul_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator*=(const mpq_t v)
+{
+    mpfr_mul_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator*=(const long double v)
+{
+    *this *= mpreal(v);    
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;    
+}
+
+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);    
+#endif
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator*=(const unsigned long int v)
+{
+    mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator*=(const unsigned int v)
+{
+    mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator*=(const long int v)
+{
+    mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator*=(const int v)
+{
+    mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *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_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
+	return c;
+}
+
+//////////////////////////////////////////////////////////////////////////
+// / Division
+inline mpreal& mpreal::operator/=(const mpreal& v)
+{
+    mpfr_div(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator/=(const mpz_t v)
+{
+    mpfr_div_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator/=(const mpq_t v)
+{
+    mpfr_div_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator/=(const long double v)
+{
+    *this /= mpreal(v);
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;    
+}
+
+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);    
+#endif
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator/=(const unsigned long int v)
+{
+    mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator/=(const unsigned int v)
+{
+    mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator/=(const long int v)
+{
+    mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator/=(const int v)
+{
+    mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+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;
+}
+
+inline const mpreal operator/(const unsigned long int b, const mpreal& a)
+{
+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+    mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+}
+
+inline const mpreal operator/(const unsigned int b, const mpreal& a)
+{
+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+    mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+}
+
+inline const mpreal operator/(const long int b, const mpreal& a)
+{
+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+    mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+}
+
+inline const mpreal operator/(const int b, const mpreal& a)
+{
+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+    mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+}
+
+inline const mpreal operator/(const double  b, const mpreal& a)
+{
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+    mpreal x(0, mpfr_get_prec(a.mpfr_srcptr()));
+    mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd());
+    return x;
+#else
+    mpreal x(0, mpfr_get_prec(a.mpfr_ptr()));
+    x /= a;
+    return x;
+#endif
+}
+
+//////////////////////////////////////////////////////////////////////////
+// Shifts operators - Multiplication/Division by power of 2
+inline mpreal& mpreal::operator<<=(const unsigned long int u)
+{
+    mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator<<=(const unsigned int u)
+{
+    mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator<<=(const long int u)
+{
+    mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator<<=(const int u)
+{
+    mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator>>=(const unsigned long int u)
+{
+    mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator>>=(const unsigned int u)
+{
+    mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator>>=(const long int u)
+{
+    mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::operator>>=(const int u)
+{
+    mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline const mpreal operator<<(const mpreal& v, const unsigned long int k)
+{
+    return mul_2ui(v,k);
+}
+
+inline const mpreal operator<<(const mpreal& v, const unsigned int k)
+{
+    return mul_2ui(v,static_cast<unsigned long int>(k));
+}
+
+inline const mpreal operator<<(const mpreal& v, const long int k)
+{
+    return mul_2si(v,k);
+}
+
+inline const mpreal operator<<(const mpreal& v, const int k)
+{
+    return mul_2si(v,static_cast<long int>(k));
+}
+
+inline const mpreal operator>>(const mpreal& v, const unsigned long int k)
+{
+    return div_2ui(v,k);
+}
+
+inline const mpreal operator>>(const mpreal& v, const long int k)
+{
+    return div_2si(v,k);
+}
+
+inline const mpreal operator>>(const mpreal& v, const unsigned int k)
+{
+    return div_2ui(v,static_cast<unsigned long int>(k));
+}
+
+inline const mpreal operator>>(const mpreal& v, const int k)
+{
+    return div_2si(v,static_cast<long int>(k));
+}
+
+// mul_2ui
+inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)
+{
+    mpreal x(v);
+    mpfr_mul_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
+    return x;
+}
+
+// mul_2si
+inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)
+{
+    mpreal x(v);
+    mpfr_mul_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
+    return x;
+}
+
+inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)
+{
+    mpreal x(v);
+    mpfr_div_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
+    return x;
+}
+
+inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)
+{
+    mpreal x(v);
+    mpfr_div_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode);
+    return x;
+}
+
+//////////////////////////////////////////////////////////////////////////
+//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 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 (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 );    }
+inline bool operator >  (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 (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 );   }
+inline bool operator >= (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 (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 );    }
+inline bool operator <  (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 (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 );   }
+inline bool operator <= (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 (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 );   }
+inline bool operator == (const mpreal& a, const double b            ){  return !isnan EIGEN_NOT_A_MACRO (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 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 );    }
+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 
+
+//////////////////////////////////////////////////////////////////////////
+// Type Converters
+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; }
+inline ::mpfr_srcptr  mpreal::mpfr_srcptr() const    { return mp; }
+
+template <class T>
+inline std::string toString(T t, std::ios_base & (*f)(std::ios_base&))
+{
+    std::ostringstream oss;
+    oss << f << t;
+    return oss.str();
+}
+
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+
+inline std::string mpreal::toString(const std::string& format) const
+{
+    char *s = NULL;
+    std::string out;
+
+    if( !format.empty() )
+    {
+        if(!(mpfr_asprintf(&s, format.c_str(), mpfr_srcptr()) < 0))
+        {
+            out = std::string(s);
+
+            mpfr_free_str(s);
+        }
+    }
+
+    return out;
+}
+
+#endif
+
+inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
+{
+    // TODO: Add extended format specification (f, e, rounding mode) as it done in output operator
+    (void)b;
+    (void)mode;
+
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+
+    std::ostringstream format;
+
+    int digits = (n >= 0) ? n : 2 + bits2digits(mpfr_get_prec(mpfr_srcptr()));
+    
+    format << "%." << digits << "RNg";
+
+    return toString(format.str());
+
+#else
+
+    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";
+    }
+
+    if(mpfr_zero_p(mp)) return "0";
+    if(mpfr_nan_p(mp))  return "NaN";
+
+    s  = mpfr_get_str(NULL, &exp, b, 0, mp, mode);
+    ns = mpfr_get_str(NULL, &exp, b, (std::max)(0,n), mp, mode);
+
+    if(s!=NULL && ns!=NULL)
+    {
+        slen  = strlen(s);
+        nslen = strlen(ns);
+        if(nslen<=slen) 
+        {
+            mpfr_free_str(s);
+            s = ns;
+            slen = nslen;
+        }
+        else {
+            mpfr_free_str(ns);
+        }
+
+        // Make human eye-friendly formatting if possible
+        if (exp>0 && static_cast<size_t>(exp)<slen)
+        {
+            if(s[0]=='-')
+            {
+                // Remove zeros starting from right end
+                char* ptr = s+slen-1;
+                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);
+
+                //out = string(s,exp+1)+'.'+string(s+exp+1);
+            }
+            else
+            {
+                // Remove zeros starting from right end
+                char* ptr = s+slen-1;
+                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);
+
+                //out = string(s,exp)+'.'+string(s+exp);
+            }
+
+        }else{ // exp<0 || exp>slen
+            if(s[0]=='-')
+            {
+                // Remove zeros starting from right end
+                char* ptr = s+slen-1;
+                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);
+
+                //out = string(s,2)+'.'+string(s+2);
+            }
+            else
+            {
+                // Remove zeros starting from right end
+                char* ptr = s+slen-1;
+                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);
+
+                //out = string(s,1)+'.'+string(s+1);
+            }
+
+            // Make final string
+            if(--exp)
+            {
+                if(exp>0) out += "e+"+mpfr::toString<mp_exp_t>(exp,std::dec);
+                else       out += "e"+mpfr::toString<mp_exp_t>(exp,std::dec);
+            }
+        }
+
+        mpfr_free_str(s);
+        return out;
+    }else{
+        return "conversion error!";
+    }
+#endif
+}
+
+
+//////////////////////////////////////////////////////////////////////////
+// I/O
+inline std::ostream& mpreal::output(std::ostream& os) const 
+{
+    std::ostringstream format;
+    const std::ios::fmtflags flags = os.flags();
+
+    format << ((flags & std::ios::showpos) ? "%+" : "%");
+    if (os.precision() >= 0)
+        format << '.' << os.precision() << "R*"
+               << ((flags & std::ios::floatfield) == std::ios::fixed ? 'f' :
+                   (flags & std::ios::floatfield) == std::ios::scientific ? 'e' :
+                   'g');
+    else
+        format << "R*e";
+
+    char *s = NULL;
+    if(!(mpfr_asprintf(&s, format.str().c_str(),
+                        mpfr::mpreal::get_default_rnd(),
+                        mpfr_srcptr())
+        < 0))
+    {
+        os << std::string(s);
+        mpfr_free_str(s);
+    }
+    return os;
+}
+
+inline std::ostream& operator<<(std::ostream& os, const mpreal& v)
+{
+    return v.output(os);
+}
+
+inline std::istream& operator>>(std::istream &is, mpreal& v)
+{
+    // TODO: use cout::hexfloat and other flags to setup base
+    std::string tmp;
+    is >> tmp;
+    mpfr_set_str(v.mpfr_ptr(), tmp.c_str(), 10, mpreal::get_default_rnd());
+    return is;
+}
+
+//////////////////////////////////////////////////////////////////////////
+//     Bits - decimal digits relation
+//        bits   = ceil(digits*log[2](10))
+//        digits = floor(bits*log[10](2))
+
+inline mp_prec_t digits2bits(int d)
+{
+    const double LOG2_10 = 3.3219280948873624;
+
+    return mp_prec_t(std::ceil( d * LOG2_10 ));
+}
+
+inline int bits2digits(mp_prec_t b)
+{
+    const double LOG10_2 = 0.30102999566398119;
+
+    return int(std::floor( b * LOG10_2 ));
+}
+
+//////////////////////////////////////////////////////////////////////////
+// Set/Get number properties
+inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode)
+{
+    mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), sign < 0, RoundingMode);
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline int mpreal::getPrecision() const
+{
+    return int(mpfr_get_prec(mpfr_srcptr()));
+}
+
+inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode)
+{
+    mpfr_prec_round(mpfr_ptr(), Precision, RoundingMode);
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mpreal& mpreal::setInf(int sign) 
+{ 
+    mpfr_set_inf(mpfr_ptr(), sign);
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}    
+
+inline mpreal& mpreal::setNan() 
+{
+    mpfr_set_nan(mpfr_ptr());
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+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 
+
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return *this;
+}
+
+inline mp_prec_t mpreal::get_prec() const
+{
+    return mpfr_get_prec(mpfr_srcptr());
+}
+
+inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode)
+{
+    mpfr_prec_round(mpfr_ptr(),prec,rnd_mode);
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+}
+
+inline mp_exp_t mpreal::get_exp () const
+{
+    return mpfr_get_exp(mpfr_srcptr());
+}
+
+inline int mpreal::set_exp (mp_exp_t e)
+{
+    int x = mpfr_set_exp(mpfr_ptr(), e);
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return x;
+}
+
+inline const mpreal frexp(const mpreal& x, mp_exp_t* exp, mp_rnd_t mode = mpreal::get_default_rnd())
+{
+    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 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 */
+    return machine_epsilon(mpreal(1, prec));
+}
+
+inline mpreal machine_epsilon(const mpreal& x)
+{    
+    /* the smallest eps such that x + eps != x */
+    if( x < 0)
+    {
+        return nextabove(-x) + x;
+    }else{
+        return nextabove( x) - x;
+    }
+}
+
+// minval is 'safe' meaning 1 / minval does not overflow
+inline mpreal minval(mp_prec_t prec)
+{
+    /* min = 1/2 * 2^emin = 2^(emin - 1) */
+    return mpreal(1, prec) << mpreal::get_emin()-1;
+}
+
+// maxval is 'safe' meaning 1 / maxval does not underflow
+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(); 
+}
+
+inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps)
+{
+    return abs(a - b) <= machine_epsilon((max)(abs(a), abs(b))) * maxUlps;
+}
+
+inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps)
+{
+    return abs(a - b) <= eps;
+}
+
+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 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_trunc(n.mpfr_ptr(),v.mpfr_srcptr());
+    return f;
+}
+
+inline int mpreal::check_range (int t, mp_rnd_t rnd_mode)
+{
+    return mpfr_check_range(mpfr_ptr(),t,rnd_mode);
+}
+
+inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode)
+{
+    int r = mpfr_subnormalize(mpfr_ptr(),t,rnd_mode);
+    MPREAL_MSVC_DEBUGVIEW_CODE;
+    return r;
+}
+
+inline mp_exp_t mpreal::get_emin (void)
+{
+    return mpfr_get_emin();
+}
+
+inline int mpreal::set_emin (mp_exp_t exp)
+{
+    return mpfr_set_emin(exp);
+}
+
+inline mp_exp_t mpreal::get_emax (void)
+{
+    return mpfr_get_emax();
+}
+
+inline int mpreal::set_emax (mp_exp_t exp)
+{
+    return mpfr_set_emax(exp);
+}
+
+inline mp_exp_t mpreal::get_emin_min (void)
+{
+    return mpfr_get_emin_min();
+}
+
+inline mp_exp_t mpreal::get_emin_max (void)
+{
+    return mpfr_get_emin_max();
+}
+
+inline mp_exp_t mpreal::get_emax_min (void)
+{
+    return mpfr_get_emax_min();
+}
+
+inline mp_exp_t mpreal::get_emax_max (void)
+{
+    return mpfr_get_emax_max();
+}
+
+//////////////////////////////////////////////////////////////////////////
+// Mathematical Functions
+//////////////////////////////////////////////////////////////////////////
+#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; 
+
+inline const mpreal sqr  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
+{   MPREAL_UNARY_MATH_FUNCTION_BODY(sqr );    }
+
+inline const mpreal sqrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
+{   MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt);    }
+
+inline const mpreal sqrt(const unsigned long int x, mp_rnd_t r)
+{
+    mpreal y;
+    mpfr_sqrt_ui(y.mpfr_ptr(), x, r);
+    return y;
+}
+
+inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode)
+{
+    return sqrt(static_cast<unsigned long int>(v),rnd_mode);
+}
+
+inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode)
+{
+    if (v>=0)   return sqrt(static_cast<unsigned long int>(v),rnd_mode);
+    else        return mpreal().setNan(); // NaN  
+}
+
+inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode)
+{
+    if (v>=0)   return sqrt(static_cast<unsigned long int>(v),rnd_mode);
+    else        return mpreal().setNan(); // NaN
+}
+
+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())); 
+    #if (MPFR_VERSION >= MPFR_VERSION_NUM(4,0,0))
+    mpfr_rootn_ui(y.mpfr_ptr(), x.mpfr_srcptr(), k, r);
+    #else
+    mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r);
+    #endif
+    return y; 
+}
+
+inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r = mpreal::get_default_rnd())
+{
+    mpreal y(0, mpfr_get_prec(a.mpfr_srcptr()));
+    mpfr_dim(y.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), r);
+    return y;
+}
+
+inline int cmpabs(const mpreal& a,const mpreal& b)
+{
+    return mpfr_cmpabs(a.mpfr_ptr(), b.mpfr_srcptr());
+}
+
+inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    return mpfr_sin_cos(s.mpfr_ptr(), c.mpfr_ptr(), v.mpfr_srcptr(), rnd_mode);
+}
+
+inline const mpreal sqrt  (const long double v, mp_rnd_t rnd_mode)    {   return sqrt(mpreal(v),rnd_mode);    }
+inline const mpreal sqrt  (const double v, mp_rnd_t rnd_mode)         {   return sqrt(mpreal(v),rnd_mode);    }
+
+inline const mpreal cbrt  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt );    }
+inline const mpreal fabs  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(abs  );    }
+inline const mpreal abs   (const mpreal& x, mp_rnd_t r)                             {   MPREAL_UNARY_MATH_FUNCTION_BODY(abs  );    }
+inline const mpreal log   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log  );    }
+inline const mpreal log2  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log2 );    }
+inline const mpreal log10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log10);    }
+inline const mpreal exp   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp  );    }
+inline const mpreal exp2  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 );    }
+inline const mpreal exp10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(exp10);    }
+inline const mpreal cos   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cos  );    }
+inline const mpreal sin   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sin  );    }
+inline const mpreal tan   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(tan  );    }
+inline const mpreal sec   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sec  );    }
+inline const mpreal csc   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(csc  );    }
+inline const mpreal cot   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cot  );    }
+inline const mpreal acos  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(acos );    }
+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);                      }
+inline const mpreal acoth (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return atanh(1/v, r);                      }
+inline const mpreal asech (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return acosh(1/v, r);                      }
+inline const mpreal acsch (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) {   return asinh(1/v, r);                      }
+
+inline const mpreal cosh  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(cosh );    }
+inline const mpreal sinh  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sinh );    }
+inline const mpreal tanh  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(tanh );    }
+inline const mpreal sech  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(sech );    }
+inline const mpreal csch  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(csch );    }
+inline const mpreal coth  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(coth );    }
+inline const mpreal acosh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(acosh);    }
+inline const mpreal asinh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(asinh);    }
+inline const mpreal atanh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(atanh);    }
+
+inline const mpreal log1p   (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(log1p  );    }
+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    );    }
+inline const mpreal erfc    (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(erfc   );    }
+inline const mpreal besselj0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(j0     );    }
+inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(j1     );    }
+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()));
+    mpfr_atan2(a.mpfr_ptr(), y.mpfr_srcptr(), x.mpfr_srcptr(), rnd_mode);
+    return a;
+}
+
+inline const mpreal hypot (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_hypot(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
+    return a;
+}
+
+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;
+}
+
+inline const mpreal remquo (long* q, 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_remquo(a.mpfr_ptr(),q, x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode);
+    return a;
+}
+
+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())
+{
+    mpreal x(0, prec);
+    mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode);
+    return x;
+}
+
+
+inline const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal x(v);
+    int tsignp;
+
+    if(signp)   mpfr_lgamma(x.mpfr_ptr(),  signp,v.mpfr_srcptr(),rnd_mode);
+    else        mpfr_lgamma(x.mpfr_ptr(),&tsignp,v.mpfr_srcptr(),rnd_mode);
+
+    return x;
+}
+
+
+inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
+{
+    mpreal  y(0, x.getPrecision());
+    mpfr_jn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);
+    return y;
+}
+
+inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
+{
+    mpreal  y(0, x.getPrecision());
+    mpfr_yn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r);
+    return y;
+}
+
+inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal a;
+    mp_prec_t p1, p2, p3;
+
+    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));
+
+    mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
+    return a;
+}
+
+inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal a;
+    mp_prec_t p1, p2, p3;
+
+    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));
+
+    mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
+    return a;
+}
+
+inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal a;
+    mp_prec_t p1, p2;
+
+    p1 = v1.get_prec(); 
+    p2 = v2.get_prec(); 
+
+    a.set_prec(p1>p2?p1:p2);
+
+    mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode);
+
+    return a;
+}
+
+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;
+    status = mpfr_sum(x.mpfr_ptr(), (mpfr_ptr*)p, n, mode);
+    
+    delete [] p;
+    return x;
+}
+
+//////////////////////////////////////////////////////////////////////////
+// MPFR 2.4.0 Specifics
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
+
+inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    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 rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    /*  R = rem(X,Y) if Y != 0, returns X - n * Y where n = trunc(X/Y). */
+    return fmod(x, y, rnd_mode);
+}
+
+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.    
+    
+    */
+
+    if(iszero(y)) return x;
+    if(x == y) return 0;
+
+    mpreal m = x - floor(x / y) * y;
+
+    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())
+{
+    mpreal a;
+    mp_prec_t yp, xp;
+
+    yp = y.get_prec(); 
+    xp = x.get_prec(); 
+
+    a.set_prec(yp>xp?yp:xp);
+
+    mpfr_fmod(a.mp, x.mp, y.mp, rnd_mode);
+
+    return a;
+}
+
+inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal x(v);
+    mpfr_rec_sqrt(x.mp,v.mp,rnd_mode);
+    return x;
+}
+#endif //  MPFR 2.4.0 Specifics
+
+//////////////////////////////////////////////////////////////////////////
+// MPFR 3.0.0 Specifics
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
+inline const mpreal digamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(digamma);     }
+inline const mpreal ai      (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(ai);          }
+#endif // MPFR 3.0.0 Specifics
+
+//////////////////////////////////////////////////////////////////////////
+// Constants
+inline const mpreal const_log2 (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
+{
+    mpreal x(0, p);
+    mpfr_const_log2(x.mpfr_ptr(), r);
+    return x;
+}
+
+inline const mpreal const_pi (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
+{
+    mpreal x(0, p);
+    mpfr_const_pi(x.mpfr_ptr(), r);
+    return x;
+}
+
+inline const mpreal const_euler (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
+{
+    mpreal x(0, p);
+    mpfr_const_euler(x.mpfr_ptr(), r);
+    return x;
+}
+
+inline const mpreal const_catalan (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd())
+{
+    mpreal x(0, p);
+    mpfr_const_catalan(x.mpfr_ptr(), r);
+    return x;
+}
+
+inline const mpreal const_infinity (int sign = 1, mp_prec_t p = mpreal::get_default_prec())
+{
+    mpreal x(0, p);
+    mpfr_set_inf(x.mpfr_ptr(), sign);
+    return x;
+}
+
+//////////////////////////////////////////////////////////////////////////
+// Integer Related Functions
+inline const mpreal ceil(const mpreal& v)
+{
+    mpreal x(v);
+    mpfr_ceil(x.mp,v.mp);
+    return x;
+}
+
+inline const mpreal floor(const mpreal& v)
+{
+    mpreal x(v);
+    mpfr_floor(x.mp,v.mp);
+    return x;
+}
+
+inline const mpreal round(const mpreal& v)
+{
+    mpreal x(v);
+    mpfr_round(x.mp,v.mp);
+    return x;
+}
+
+inline const mpreal trunc(const mpreal& v)
+{
+    mpreal x(v);
+    mpfr_trunc(x.mp,v.mp);
+    return x;
+}
+
+inline const mpreal rint       (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint      );     }
+inline const mpreal rint_ceil  (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil );     }
+inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor);     }
+inline const mpreal rint_round (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round);     }
+inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc);     }
+inline const mpreal frac       (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) {   MPREAL_UNARY_MATH_FUNCTION_BODY(frac      );     }
+
+//////////////////////////////////////////////////////////////////////////
+// Miscellaneous Functions
+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<y?x:y);       }
+
+inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal a;
+    mpfr_max(a.mp,x.mp,y.mp,rnd_mode);
+    return a;
+}
+
+inline const mpreal fmin(const mpreal& x, const mpreal& y,  mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal a;
+    mpfr_min(a.mp,x.mp,y.mp,rnd_mode);
+    return a;
+}
+
+inline const mpreal nexttoward (const mpreal& x, const mpreal& y)
+{
+    mpreal a(x);
+    mpfr_nexttoward(a.mp,y.mp);
+    return a;
+}
+
+inline const mpreal nextabove  (const mpreal& x)
+{
+    mpreal a(x);
+    mpfr_nextabove(a.mp);
+    return a;
+}
+
+inline const mpreal nextbelow  (const mpreal& x)
+{
+    mpreal a(x);
+    mpfr_nextbelow(a.mp);
+    return a;
+}
+
+inline const mpreal urandomb (gmp_randstate_t& state)
+{
+    mpreal x;
+    mpfr_urandomb(x.mpfr_ptr(),state);
+    return x;
+}
+
+#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.mpfr_ptr(), state, rnd_mode);
+    return x;
+}
+#endif
+
+#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2))
+inline const mpreal random2 (mp_size_t size, mp_exp_t exp)
+{
+    mpreal x;
+    mpfr_random2(x.mpfr_ptr(),size,exp);
+    return x;
+}
+#endif
+
+// Uniformly distributed random number generation
+// a = random(seed); <- initialization & first random number generation
+// a = random();     <- next random numbers generation
+// 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 initialize = true;
+
+    if(initialize)
+    {
+        gmp_randinit_default(state);
+        gmp_randseed_ui(state,0);
+        initialize = false;
+    }
+
+    if(seed != 0)    gmp_randseed_ui(state,seed);
+
+    return mpfr::urandom(state);
+#else
+    if(seed != 0)    std::srand(seed);
+    return mpfr::mpreal(std::rand()/(double)RAND_MAX);
+#endif
+
+}
+
+#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)
+{
+    static gmp_randstate_t state;
+    static bool initialize = true;
+
+    if(initialize)
+    {
+        gmp_randinit_default(state);
+        gmp_randseed_ui(state,0);
+        initialize = false;
+    }
+
+    if(seed != 0) gmp_randseed_ui(state,seed);
+
+    return mpfr::grandom(state);
+}
+#endif
+
+//////////////////////////////////////////////////////////////////////////
+// Set/Get global properties
+inline void mpreal::set_default_prec(mp_prec_t prec)
+{ 
+    mpfr_set_default_prec(prec); 
+}
+
+inline void mpreal::set_default_rnd(mp_rnd_t 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);
+}
+
+inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal x(a);
+    mpfr_pow(x.mp,x.mp,b.mp,rnd_mode);
+    return x;
+}
+
+inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal x(a);
+    mpfr_pow_z(x.mp,x.mp,b,rnd_mode);
+    return x;
+}
+
+inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal x(a);
+    mpfr_pow_ui(x.mp,x.mp,b,rnd_mode);
+    return x;
+}
+
+inline const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode)
+{
+    return pow(a,static_cast<unsigned long int>(b),rnd_mode);
+}
+
+inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal x(a);
+    mpfr_pow_si(x.mp,x.mp,b,rnd_mode);
+    return x;
+}
+
+inline const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode)
+{
+    return pow(a,static_cast<long int>(b),rnd_mode);
+}
+
+inline const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode)
+{
+    return pow(a,mpreal(b),rnd_mode);
+}
+
+inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode)
+{
+    return pow(a,mpreal(b),rnd_mode);
+}
+
+inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
+{
+    mpreal x(a);
+    mpfr_ui_pow(x.mp,a,b.mp,rnd_mode);
+    return x;
+}
+
+inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode)
+{
+    return pow(static_cast<unsigned long int>(a),b,rnd_mode);
+}
+
+inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode)
+{
+    if (a>=0)     return pow(static_cast<unsigned long int>(a),b,rnd_mode);
+    else          return pow(mpreal(a),b,rnd_mode);
+}
+
+inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode)
+{
+    if (a>=0)     return pow(static_cast<unsigned long int>(a),b,rnd_mode);
+    else          return pow(mpreal(a),b,rnd_mode);
+}
+
+inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),b,rnd_mode);
+}
+
+inline const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),b,rnd_mode);
+}
+
+// pow unsigned long int
+inline const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode)
+{
+    mpreal x(a);
+    mpfr_ui_pow_ui(x.mp,a,b,rnd_mode);
+    return x;
+}
+
+inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode)
+{
+    return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+}
+
+inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode)
+{
+    if(b>0)    return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+    else       return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
+}
+
+inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode)
+{
+    if(b>0)    return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+    else       return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
+}
+
+inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode)
+{
+    return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
+}
+
+inline const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode)
+{
+    return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
+}
+
+// pow unsigned int
+inline const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode)
+{
+    return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui
+}
+
+inline const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode)
+{
+    return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+}
+
+inline const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode)
+{
+    if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+    else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+}
+
+inline const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode)
+{
+    if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+    else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+}
+
+inline const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode)
+{
+    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+}
+
+inline const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode)
+{
+    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+}
+
+// pow long int
+inline const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode)
+{
+    if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui
+    else     return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui
+}
+
+inline const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode)
+{
+    if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode);  //mpfr_ui_pow_ui
+    else     return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui
+}
+
+inline const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode)
+{
+    if (a>0)
+    {
+        if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+        else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+    }else{
+        return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si
+    }
+}
+
+inline const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode)
+{
+    if (a>0)
+    {
+        if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+        else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+    }else{
+        return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si
+    }
+}
+
+inline const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode)
+{
+    if (a>=0)   return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+    else        return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow
+}
+
+inline const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode)
+{
+    if (a>=0)   return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+    else        return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow
+}
+
+// pow int
+inline const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode)
+{
+    if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui
+    else     return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui
+}
+
+inline const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode)
+{
+    if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode);  //mpfr_ui_pow_ui
+    else     return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui
+}
+
+inline const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode)
+{
+    if (a>0)
+    {
+        if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+        else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+    }else{
+        return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si
+    }
+}
+
+inline const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode)
+{
+    if (a>0)
+    {
+        if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
+        else    return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+    }else{
+        return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si
+    }
+}
+
+inline const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode)
+{
+    if (a>=0)   return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+    else        return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow
+}
+
+inline const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode)
+{
+    if (a>=0)   return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
+    else        return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow
+}
+
+// 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);
+}
+
+inline const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui
+}
+
+inline const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui
+}
+
+inline const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si
+}
+
+inline const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si
+}
+
+inline const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),mpreal(b),rnd_mode);
+}
+
+inline const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),b,rnd_mode); // mpfr_pow_ui
+}
+
+inline const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); // mpfr_pow_ui
+}
+
+inline const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si
+}
+
+inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode)
+{
+    return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si
+}
+} // End of mpfr namespace
+
+// Explicit specialization of std::swap for mpreal numbers
+// Thus standard algorithms will use efficient version of swap (due to Koenig lookup)
+// 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
+    template <>
+    inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) 
+    { 
+        return mpfr::swap(x, y); 
+    }
+
+    template<>
+    class numeric_limits<mpfr::mpreal>
+    {
+    public:
+        static const bool is_specialized    = true;
+        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 bool has_infinity      = true;
+        static const bool has_quiet_NaN     = true;
+        static const bool has_signaling_NaN = true;
+
+        static const bool is_iec559         = true;        // = IEEE 754
+        static const bool is_bounded        = true;
+        static const bool is_modulo         = false;
+        static const bool traps             = true;
+        static const bool tinyness_before   = true;
+
+        static const float_denorm_style has_denorm  = denorm_absent;
+
+        inline static mpfr::mpreal (min)    (mp_prec_t precision = mpfr::mpreal::get_default_prec()) {  return  mpfr::minval(precision);  }
+        inline static mpfr::mpreal (max)    (mp_prec_t precision = mpfr::mpreal::get_default_prec()) {  return  mpfr::maxval(precision);  }
+        inline static mpfr::mpreal lowest   (mp_prec_t precision = mpfr::mpreal::get_default_prec()) {  return -mpfr::maxval(precision);  }
+
+        // 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);  }
+
+        inline static mpfr::mpreal round_error(mp_prec_t precision = mpfr::mpreal::get_default_prec())
+        {
+            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);    
+        }
+
+        inline static const mpfr::mpreal infinity()         { return mpfr::const_infinity();     }
+        inline static const mpfr::mpreal quiet_NaN()        { return mpfr::mpreal().setNan();    }
+        inline static const mpfr::mpreal signaling_NaN()    { return mpfr::mpreal().setNan();    }
+        inline static const mpfr::mpreal denorm_min()       { return (min)();                    }
+
+        // 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); 
+
+#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. 
+        //
+        // This is preferable way for std::numeric_limits<mpfr::mpreal> 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. 
+        inline static float_round_style round_style()
+        {
+            mp_rnd_t r = mpfr::mpreal::get_default_rnd();
+
+            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; 
+            default: return round_indeterminate;
+            }
+        }
+
+        inline static int digits()                        {    return int(mpfr::mpreal::get_default_prec());    }
+        inline static int digits(const mpfr::mpreal& x)   {    return x.getPrecision();                         }
+
+        inline static int digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec())
+        {
+            return mpfr::bits2digits(precision);
+        }
+
+        inline static int digits10(const mpfr::mpreal& x)
+        {
+            return mpfr::bits2digits(x.getPrecision());
+        }
+
+        inline static int max_digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec())
+        {
+            return digits10(precision);
+        }
+#else
+        // Digits and round_style are NOT constants when it comes to mpreal.
+        // 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. 
+        //
+        // For example, if you use 256 bits of precision uniformly in your program, then:
+        // digits       = 256
+        // 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;
+        static const int digits       = 53;
+        static const int digits10     = 15;
+        static const int max_digits10 = 16;
+#endif
+    };
+
+}
+
+#endif /* __MPREAL_H__ */