From 316dadc8e4add23129c6eec663479ccdd4d4b308 Mon Sep 17 00:00:00 2001 From: Hauke Heibel Date: Fri, 24 Sep 2010 17:32:44 +0200 Subject: [PATCH] Fixed some SVD issues. Make the SVD's output unitary. Improved unit tests. Added an assert to the SVD ctor to check whether rows>=cols. --- Eigen/src/SVD/SVD.h | 13 +- test/svd.cpp | 5 +- unsupported/test/mpreal.h | 6244 ++++++++++++++++++------------------- 3 files changed, 3132 insertions(+), 3130 deletions(-) diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index b8f0c6496..9207fbd6f 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -83,11 +83,14 @@ template class SVD * \sa JacobiSVD() */ SVD(Index rows, Index cols) : m_matU(rows, rows), - m_matV(cols,cols), - m_sigma(std::min(rows, cols)), - m_workMatrix(rows, cols), - m_rv1(cols), - m_isInitialized(false) {} + m_matV(cols,cols), + m_sigma(std::min(rows, cols)), + m_workMatrix(rows, cols), + m_rv1(cols), + m_isInitialized(false) + { + ei_assert(rows >= cols && "SVD is only defined if rows>=cols."); + } SVD(const MatrixType& matrix) : m_matU(matrix.rows(), matrix.rows()), m_matV(matrix.cols(),matrix.cols()), diff --git a/test/svd.cpp b/test/svd.cpp index 5aa766fee..c151a5486 100644 --- a/test/svd.cpp +++ b/test/svd.cpp @@ -95,6 +95,8 @@ template void svd_verify_assert() VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp)) VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp)) VERIFY_RAISES_ASSERT(svd.computeScalingRotation(&tmp,&tmp)) + + VERIFY_RAISES_ASSERT(SVD(10, 20)) } void test_svd() @@ -116,7 +118,4 @@ void test_svd() CALL_SUBTEST_2( svd_verify_assert() ); CALL_SUBTEST_3( svd_verify_assert() ); CALL_SUBTEST_4( svd_verify_assert() ); - - // Test problem size constructors - CALL_SUBTEST_9( SVD(10, 20) ); } diff --git a/unsupported/test/mpreal.h b/unsupported/test/mpreal.h index 29a0cf7ff..d8f315111 100644 --- a/unsupported/test/mpreal.h +++ b/unsupported/test/mpreal.h @@ -1,3122 +1,3122 @@ -/* - Multi-precision real number class. C++ wrapper fo MPFR library. - Project homepage: http://www.holoborodko.com/pavel/ - Contact e-mail: pavel@holoborodko.com - - Copyright (c) 2008-2010 Pavel Holoborodko - - This library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - This library 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 - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with this library; if not, write to the Free Software - Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA - - Contributors: - Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, - Heinz van Saanen, Pere Constans, Dmitriy Gubanov -*/ - -#ifndef __MP_REAL_H__ -#define __MP_REAL_H__ - -#include -#include -#include -#include -#include -#include - -#include - -// Detect compiler using signatures from http://predef.sourceforge.net/ -#if defined(__GNUC__) && defined(__INTEL_COMPILER) - #define IsInf(x) isinf(x) // GNU C/C++ + Intel ICC compiler - -#elif defined(__GNUC__) - #define IsInf(x) std::isinf(x) // GNU C/C++ - -#elif defined(_MSC_VER) - #define IsInf(x) (!_finite(x)) // Microsoft Visual C++ - -#else - #define IsInf(x) std::isinf(x) // C99 conformance -#endif - -namespace mpfr { - -class mpreal { -private: - mpfr_t mp; - -public: - static mp_rnd_t default_rnd; - static mp_prec_t default_prec; - static int default_base; - static int double_bits; - -public: - // Constructors && type conversion - mpreal(); - mpreal(const mpreal& u); - - mpreal(const mpfr_t u); - mpreal(const mpf_t u); - - mpreal(const mpz_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); - mpreal(const mpq_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); - mpreal(const double u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); - mpreal(const long double u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); - mpreal(const unsigned long int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); - mpreal(const unsigned int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); - mpreal(const long int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); - mpreal(const int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); - mpreal(const char* s, mp_prec_t prec = default_prec, int base = default_base, mp_rnd_t mode = default_rnd); - - ~mpreal(); - - // 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 int v); - mpreal& operator=(const long int v); - mpreal& operator=(const int v); - mpreal& operator=(const char* s); - - // + - 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); - 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); - - // Boolean Operators - friend bool operator > (const mpreal& a, const mpreal& b); - friend bool operator >= (const mpreal& a, const mpreal& b); - friend bool operator < (const mpreal& a, const mpreal& b); - friend bool operator <= (const mpreal& a, const mpreal& b); - friend bool operator == (const mpreal& a, const mpreal& b); - friend bool operator != (const mpreal& a, const mpreal& b); - - // Type Conversion operators - inline operator long double() const; - inline operator double() const; - inline operator float() const; - inline operator unsigned long() const; - inline operator unsigned int() const; - inline operator long() const; - operator std::string() const; - inline operator mpfr_ptr(); - - // Math Functions - friend const mpreal sqr(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend int cmpabs(const mpreal& a,const mpreal& b); - - friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - - friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - - friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - - friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal fac_ui (unsigned long int v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal log1p (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal expm1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - - friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal _j0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal _j1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal _jn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal _y0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal _y1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal _yn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend int sgn(const mpreal& v); // -1 or +1 - -// 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 = mpreal::default_rnd); - friend const mpreal li2(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); -#endif - -// MPFR 3.0.0 Specifics -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) - friend const mpreal digamma(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal urandom (gmp_randstate_t& state,mp_rnd_t rnd_mode = mpreal::default_rnd); // use gmp_randinit_default() to init state, gmp_randclear() to clear - friend bool _isregular(const mpreal& v); -#endif - - // Exponent and mantissa manipulation - friend const mpreal frexp(const mpreal& v, mp_exp_t* exp); - friend const mpreal ldexp(const mpreal& v, mp_exp_t exp); - - // Splits mpreal value into fractional and integer parts. - // Returns fractional part and stores integer part in n. - friend const mpreal modf(const mpreal& v, mpreal& n); - - // 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 = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal const_pi (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal const_euler (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal const_catalan (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); - // returns +inf iff sign>=0 otherwise -inf - friend const mpreal const_infinity(int sign = 1, mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); - - // 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 = mpreal::default_rnd); - 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 = mpreal::default_rnd); - friend const mpreal rint_floor(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal rint_round(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal rint_trunc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); - friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); - - // 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(const mpreal& v); - friend bool _isinf(const mpreal& v); - friend bool _isnum(const mpreal& v); - friend bool _iszero(const mpreal& v); - friend bool _isint(const mpreal& v); - - // Set/Get instance properties - inline mp_prec_t get_prec() const; - inline void set_prec(mp_prec_t prec, mp_rnd_t rnd_mode = default_rnd); // Change precision with rounding mode - - // Set mpreal to +-inf, NaN - void set_inf(int sign = +1); - void set_nan(); - - // sign = -1 or +1 - void set_sign(int sign, mp_rnd_t rnd_mode = default_rnd); - - //Exponent - mp_exp_t get_exp(); - int set_exp(mp_exp_t e); - int check_range (int t, mp_rnd_t rnd_mode = default_rnd); - int subnormalize (int t,mp_rnd_t rnd_mode = 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 mp_prec_t get_default_prec(); - static void set_default_base(int base); - static int get_default_base(); - static void set_double_bits(int dbits); - static int get_double_bits(); - static void set_default_rnd(mp_rnd_t rnd_mode); - static mp_rnd_t get_default_rnd(); - 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); - - // Get/Set conversions - // Convert mpreal to string with n significant digits in base b - // n = 0 -> convert with the maximum available digits - std::string to_string(size_t n = 0, int b = default_base, mp_rnd_t mode = default_rnd) const; - - // Efficient swapping of two mpreal values - friend void swap(mpreal& x, mpreal& y); - - //Min Max - macros is evil. Needed for systems which defines max and min globally as macros (e.g. Windows) - //Hope that globally defined macros use > < operations only - #ifndef max - friend const mpreal max(const mpreal& x, const mpreal& y); - #endif - - #ifndef min - friend const mpreal min(const mpreal& x, const mpreal& y); - #endif -}; - -////////////////////////////////////////////////////////////////////////// -// Exceptions -class conversion_overflow : public std::exception { -public: - std::string why() { return "inexact conversion from floating point"; } -}; - -////////////////////////////////////////////////////////////////////////// -// + Addition -const mpreal operator+(const mpreal& a, const mpreal& b); - -// + Fast specialized addition - implemented through fast += operations -const mpreal operator+(const mpreal& a, const mpz_t b); -const mpreal operator+(const mpreal& a, const mpq_t b); -const mpreal operator+(const mpreal& a, const long double b); -const mpreal operator+(const mpreal& a, const double b); -const mpreal operator+(const mpreal& a, const unsigned long int b); -const mpreal operator+(const mpreal& a, const unsigned int b); -const mpreal operator+(const mpreal& a, const long int b); -const mpreal operator+(const mpreal& a, const int b); -const mpreal operator+(const mpreal& a, const char* b); -const mpreal operator+(const char* a, const mpreal& b); -const std::string operator+(const mpreal& a, const std::string b); -const std::string operator+(const std::string a, const mpreal& b); - -const mpreal operator+(const mpz_t b, const mpreal& a); -const mpreal operator+(const mpq_t b, const mpreal& a); -const mpreal operator+(const long double b, const mpreal& a); -const mpreal operator+(const double b, const mpreal& a); -const mpreal operator+(const unsigned long int b, const mpreal& a); -const mpreal operator+(const unsigned int b, const mpreal& a); -const mpreal operator+(const long int b, const mpreal& a); -const mpreal operator+(const int b, const mpreal& a); - -////////////////////////////////////////////////////////////////////////// -// - Subtraction -const mpreal operator-(const mpreal& a, const mpreal& b); - -// - Fast specialized subtraction - implemented through fast -= operations -const mpreal operator-(const mpreal& a, const mpz_t b); -const mpreal operator-(const mpreal& a, const mpq_t b); -const mpreal operator-(const mpreal& a, const long double b); -const mpreal operator-(const mpreal& a, const double b); -const mpreal operator-(const mpreal& a, const unsigned long int b); -const mpreal operator-(const mpreal& a, const unsigned int b); -const mpreal operator-(const mpreal& a, const long int b); -const mpreal operator-(const mpreal& a, const int b); -const mpreal operator-(const mpreal& a, const char* b); -const mpreal operator-(const char* a, const mpreal& b); - -const mpreal operator-(const mpz_t b, const mpreal& a); -const mpreal operator-(const mpq_t b, const mpreal& a); -const mpreal operator-(const long double b, const mpreal& a); -//const mpreal operator-(const double b, const mpreal& a); - -////////////////////////////////////////////////////////////////////////// -// * Multiplication -const mpreal operator*(const mpreal& a, const mpreal& b); - -// * Fast specialized multiplication - implemented through fast *= operations -const mpreal operator*(const mpreal& a, const mpz_t b); -const mpreal operator*(const mpreal& a, const mpq_t b); -const mpreal operator*(const mpreal& a, const long double b); -const mpreal operator*(const mpreal& a, const double b); -const mpreal operator*(const mpreal& a, const unsigned long int b); -const mpreal operator*(const mpreal& a, const unsigned int b); -const mpreal operator*(const mpreal& a, const long int b); -const mpreal operator*(const mpreal& a, const int b); - -const mpreal operator*(const mpz_t b, const mpreal& a); -const mpreal operator*(const mpq_t b, const mpreal& a); -const mpreal operator*(const long double b, const mpreal& a); -const mpreal operator*(const double b, const mpreal& a); -const mpreal operator*(const unsigned long int b, const mpreal& a); -const mpreal operator*(const unsigned int b, const mpreal& a); -const mpreal operator*(const long int b, const mpreal& a); -const mpreal operator*(const int b, const mpreal& a); - -////////////////////////////////////////////////////////////////////////// -// / Division -const mpreal operator/(const mpreal& a, const mpreal& b); - -// / Fast specialized division - implemented through fast /= operations -const mpreal operator/(const mpreal& a, const mpz_t b); -const mpreal operator/(const mpreal& a, const mpq_t b); -const mpreal operator/(const mpreal& a, const long double b); -const mpreal operator/(const mpreal& a, const double b); -const mpreal operator/(const mpreal& a, const unsigned long int b); -const mpreal operator/(const mpreal& a, const unsigned int b); -const mpreal operator/(const mpreal& a, const long int b); -const mpreal operator/(const mpreal& a, const int b); - -const mpreal operator/(const long double b, const mpreal& a); - -////////////////////////////////////////////////////////////////////////// -// Shifts operators - Multiplication/Division by a power of 2 -const mpreal operator<<(const mpreal& v, const unsigned long int k); -const mpreal operator<<(const mpreal& v, const unsigned int k); -const mpreal operator<<(const mpreal& v, const long int k); -const mpreal operator<<(const mpreal& v, const int k); - -const mpreal operator>>(const mpreal& v, const unsigned long int k); -const mpreal operator>>(const mpreal& v, const unsigned int k); -const mpreal operator>>(const mpreal& v, const long int k); -const mpreal operator>>(const mpreal& v, const int k); - -////////////////////////////////////////////////////////////////////////// -// Boolean operators -bool operator < (const mpreal& a, const unsigned long int b); -bool operator < (const mpreal& a, const unsigned int b); -bool operator < (const mpreal& a, const long int b); -bool operator < (const mpreal& a, const int b); -bool operator < (const mpreal& a, const long double b); -bool operator < (const mpreal& a, const double b); - -bool operator < (const unsigned long int a,const mpreal& b); -bool operator < (const unsigned int a, const mpreal& b); -bool operator < (const long int a, const mpreal& b); -bool operator < (const int a, const mpreal& b); -bool operator < (const long double a, const mpreal& b); -bool operator < (const double a, const mpreal& b); - -bool operator > (const mpreal& a, const unsigned long int b); -bool operator > (const mpreal& a, const unsigned int b); -bool operator > (const mpreal& a, const long int b); -bool operator > (const mpreal& a, const int b); -bool operator > (const mpreal& a, const long double b); -bool operator > (const mpreal& a, const double b); - -bool operator > (const unsigned long int a,const mpreal& b); -bool operator > (const unsigned int a, const mpreal& b); -bool operator > (const long int a, const mpreal& b); -bool operator > (const int a, const mpreal& b); -bool operator > (const long double a, const mpreal& b); -bool operator > (const double a, const mpreal& b); - -bool operator >= (const mpreal& a, const unsigned long int b); -bool operator >= (const mpreal& a, const unsigned int b); -bool operator >= (const mpreal& a, const long int b); -bool operator >= (const mpreal& a, const int b); -bool operator >= (const mpreal& a, const long double b); -bool operator >= (const mpreal& a, const double b); - -bool operator >= (const unsigned long int a,const mpreal& b); -bool operator >= (const unsigned int a, const mpreal& b); -bool operator >= (const long int a, const mpreal& b); -bool operator >= (const int a, const mpreal& b); -bool operator >= (const long double a, const mpreal& b); -bool operator >= (const double a, const mpreal& b); - -bool operator <= (const mpreal& a, const unsigned long int b); -bool operator <= (const mpreal& a, const unsigned int b); -bool operator <= (const mpreal& a, const long int b); -bool operator <= (const mpreal& a, const int b); -bool operator <= (const mpreal& a, const long double b); -bool operator <= (const mpreal& a, const double b); - -bool operator <= (const unsigned long int a,const mpreal& b); -bool operator <= (const unsigned int a, const mpreal& b); -bool operator <= (const long int a, const mpreal& b); -bool operator <= (const int a, const mpreal& b); -bool operator <= (const long double a, const mpreal& b); -bool operator <= (const double a, const mpreal& b); - -bool operator == (const mpreal& a, const unsigned long int b); -bool operator == (const mpreal& a, const unsigned int b); -bool operator == (const mpreal& a, const long int b); -bool operator == (const mpreal& a, const int b); -bool operator == (const mpreal& a, const long double b); -bool operator == (const mpreal& a, const double b); - -bool operator == (const unsigned long int a,const mpreal& b); -bool operator == (const unsigned int a, const mpreal& b); -bool operator == (const long int a, const mpreal& b); -bool operator == (const int a, const mpreal& b); -bool operator == (const long double a, const mpreal& b); -bool operator == (const double a, const mpreal& b); - -bool operator != (const mpreal& a, const unsigned long int b); -bool operator != (const mpreal& a, const unsigned int b); -bool operator != (const mpreal& a, const long int b); -bool operator != (const mpreal& a, const int b); -bool operator != (const mpreal& a, const long double b); -bool operator != (const mpreal& a, const double b); - -bool operator != (const unsigned long int a,const mpreal& b); -bool operator != (const unsigned int a, const mpreal& b); -bool operator != (const long int a, const mpreal& b); -bool operator != (const int a, const mpreal& b); -bool operator != (const long double a, const mpreal& b); -bool operator != (const double a, const mpreal& b); - -////////////////////////////////////////////////////////////////////////// -// sqrt -const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal sqrt(const long int v, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::default_rnd); - -////////////////////////////////////////////////////////////////////////// -// pow -const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); - -const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); - -const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); - -const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); - -const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); - -const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); - -const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); - -const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); -const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); - -////////////////////////////////////////////////////////////////////////// -// Estimate machine epsilon for the given precision -inline const mpreal machine_epsilon(mp_prec_t prec); -inline const mpreal mpreal_min(mp_prec_t prec); -inline const mpreal mpreal_max(mp_prec_t prec); - -////////////////////////////////////////////////////////////////////////// -// Implementation of inline functions -////////////////////////////////////////////////////////////////////////// - -////////////////////////////////////////////////////////////////////////// -// Operators - Assignment -inline mpreal& mpreal::operator=(const mpreal& v) -{ - if (this!= &v) mpfr_set(mp,v.mp,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator=(const mpf_t v) -{ - mpfr_set_f(mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator=(const mpz_t v) -{ - mpfr_set_z(mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator=(const mpq_t v) -{ - mpfr_set_q(mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator=(const long double v) -{ - mpfr_set_ld(mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator=(const double v) -{ - if(double_bits == -1 || fits_in_bits(v, double_bits)) - { - mpfr_set_d(mp,v,default_rnd); - } - else - throw conversion_overflow(); - - return *this; -} - -inline mpreal& mpreal::operator=(const unsigned long int v) -{ - mpfr_set_ui(mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator=(const unsigned int v) -{ - mpfr_set_ui(mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator=(const long int v) -{ - mpfr_set_si(mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator=(const int v) -{ - mpfr_set_si(mp,v,default_rnd); - return *this; -} - -////////////////////////////////////////////////////////////////////////// -// + Addition -inline mpreal& mpreal::operator+=(const mpreal& v) -{ - mpfr_add(mp,mp,v.mp,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator+=(const mpf_t u) -{ - *this += mpreal(u); - return *this; -} - -inline mpreal& mpreal::operator+=(const mpz_t u) -{ - mpfr_add_z(mp,mp,u,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator+=(const mpq_t u) -{ - mpfr_add_q(mp,mp,u,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator+= (const long double u) -{ - return *this += mpreal(u); -} - -inline mpreal& mpreal::operator+= (const double u) -{ -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) - mpfr_add_d(mp,mp,u,default_rnd); - return *this; -#else - return *this += mpreal(u); -#endif -} - -inline mpreal& mpreal::operator+=(const unsigned long int u) -{ - mpfr_add_ui(mp,mp,u,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator+=(const unsigned int u) -{ - mpfr_add_ui(mp,mp,u,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator+=(const long int u) -{ - mpfr_add_si(mp,mp,u,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator+=(const int u) -{ - mpfr_add_si(mp,mp,u,default_rnd); - return *this; -} - -inline const mpreal mpreal::operator+()const -{ - return mpreal(*this); -} - -inline const mpreal operator+(const mpreal& a, const mpreal& b) -{ - // prec(a+b) = max(prec(a),prec(b)) - if(a.get_prec()>b.get_prec()) return mpreal(a) += b; - else return mpreal(b) += a; -} - -inline const std::string operator+(const mpreal& a, const std::string b) -{ - return (std::string)a+b; -} - -inline const std::string operator+(const std::string a, const mpreal& b) -{ - return a+(std::string)b; -} - -inline const mpreal operator+(const mpreal& a, const mpz_t b) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const mpreal& a, const char* b) -{ - return a+mpreal(b); -} - -inline const mpreal operator+(const char* a, const mpreal& b) -{ - return mpreal(a)+b; - -} - -inline const mpreal operator+(const mpreal& a, const mpq_t b) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const mpreal& a, const long double b) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const mpreal& a, const double b) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const mpreal& a, const unsigned long int b) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const mpreal& a, const unsigned int b) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const mpreal& a, const long int b) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const mpreal& a, const int b) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const mpz_t b, const mpreal& a) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const mpq_t b, const mpreal& a) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const long double b, const mpreal& a) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const double b, const mpreal& a) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const unsigned long int b, const mpreal& a) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const unsigned int b, const mpreal& a) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const long int b, const mpreal& a) -{ - return mpreal(a) += b; -} - -inline const mpreal operator+(const int b, const mpreal& a) -{ - return mpreal(a) += b; -} - -inline mpreal& mpreal::operator++() -{ - *this += 1; - return *this; -} - -inline const mpreal mpreal::operator++ (int) -{ - mpreal x(*this); - *this += 1; - return x; -} - -inline mpreal& mpreal::operator--() -{ - *this -= 1; - return *this; -} - -inline const mpreal mpreal::operator-- (int) -{ - mpreal x(*this); - *this -= 1; - return x; -} - -////////////////////////////////////////////////////////////////////////// -// - Subtraction -inline mpreal& mpreal::operator-= (const mpreal& v) -{ - mpfr_sub(mp,mp,v.mp,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator-=(const mpz_t v) -{ - mpfr_sub_z(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator-=(const mpq_t v) -{ - mpfr_sub_q(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator-=(const long double v) -{ - return *this -= mpreal(v); -} - -inline mpreal& mpreal::operator-=(const double v) -{ -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) - mpfr_sub_d(mp,mp,v,default_rnd); - return *this; -#else - return *this -= mpreal(v); -#endif -} - -inline mpreal& mpreal::operator-=(const unsigned long int v) -{ - mpfr_sub_ui(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator-=(const unsigned int v) -{ - mpfr_sub_ui(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator-=(const long int v) -{ - mpfr_sub_si(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator-=(const int v) -{ - mpfr_sub_si(mp,mp,v,default_rnd); - return *this; -} - -inline const mpreal mpreal::operator-()const -{ - mpreal u(*this); - mpfr_neg(u.mp,u.mp,default_rnd); - return u; -} - -inline const mpreal operator-(const mpreal& a, const mpreal& b) -{ - // prec(a-b) = max(prec(a),prec(b)) - if(a.get_prec()>b.get_prec()) return mpreal(a) -= b; - else return -(mpreal(b) -= a); -} - -inline const mpreal operator-(const mpreal& a, const mpz_t b) -{ - return mpreal(a) -= b; -} - -inline const mpreal operator-(const mpreal& a, const mpq_t b) -{ - return mpreal(a) -= b; -} - -inline const mpreal operator-(const mpreal& a, const long double b) -{ - return mpreal(a) -= b; -} - -inline const mpreal operator-(const mpreal& a, const double b) -{ - return mpreal(a) -= b; -} - -inline const mpreal operator-(const mpreal& a, const unsigned long int b) -{ - return mpreal(a) -= b; -} - -inline const mpreal operator-(const mpreal& a, const unsigned int b) -{ - return mpreal(a) -= b; -} - -inline const mpreal operator-(const mpreal& a, const long int b) -{ - return mpreal(a) -= b; -} - -inline const mpreal operator-(const mpreal& a, const int b) -{ - return mpreal(a) -= b; -} - -inline const mpreal operator-(const mpz_t b, const mpreal& a) -{ - return -(mpreal(a) -= b); -} - -inline const mpreal operator-(const mpq_t b, const mpreal& a) -{ - return -(mpreal(a) -= b); -} - -inline const mpreal operator-(const long double b, const mpreal& a) -{ - return -(mpreal(a) -= b); -} - -inline const mpreal operator-(const double b, const mpreal& a) -{ -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) - mpreal x(a); - mpfr_d_sub(x.mp,b,a.mp,mpreal::default_rnd); - return x; -#else - return -(mpreal(a) -= b); -#endif -} - -inline const mpreal operator-(const unsigned long int b, const mpreal& a) -{ - mpreal x(a); - mpfr_ui_sub(x.mp,b,a.mp,mpreal::default_rnd); - return x; -} - -inline const mpreal operator-(const unsigned int b, const mpreal& a) -{ - mpreal x(a); - mpfr_ui_sub(x.mp,b,a.mp,mpreal::default_rnd); - return x; -} - -inline const mpreal operator-(const long int b, const mpreal& a) -{ - mpreal x(a); - mpfr_si_sub(x.mp,b,a.mp,mpreal::default_rnd); - return x; -} - -inline const mpreal operator-(const int b, const mpreal& a) -{ - mpreal x(a); - mpfr_si_sub(x.mp,b,a.mp,mpreal::default_rnd); - return x; -} - -inline const mpreal operator-(const mpreal& a, const char* b) -{ - return a-mpreal(b); -} - -inline const mpreal operator-(const char* a, const mpreal& b) -{ - return mpreal(a)-b; -} - -////////////////////////////////////////////////////////////////////////// -// * Multiplication -inline mpreal& mpreal::operator*= (const mpreal& v) -{ - mpfr_mul(mp,mp,v.mp,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator*=(const mpz_t v) -{ - mpfr_mul_z(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator*=(const mpq_t v) -{ - mpfr_mul_q(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator*=(const long double v) -{ - return *this *= mpreal(v); -} - -inline mpreal& mpreal::operator*=(const double v) -{ -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) - mpfr_mul_d(mp,mp,v,default_rnd); - return *this; -#else - return *this *= mpreal(v); -#endif -} - -inline mpreal& mpreal::operator*=(const unsigned long int v) -{ - mpfr_mul_ui(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator*=(const unsigned int v) -{ - mpfr_mul_ui(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator*=(const long int v) -{ - mpfr_mul_si(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator*=(const int v) -{ - mpfr_mul_si(mp,mp,v,default_rnd); - return *this; -} - -inline const mpreal operator*(const mpreal& a, const mpreal& b) -{ - // prec(a*b) = max(prec(a),prec(b)) - if(a.get_prec()>b.get_prec()) return mpreal(a) *= b; - else return mpreal(b) *= a; -} - -inline const mpreal operator*(const mpreal& a, const mpz_t b) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const mpreal& a, const mpq_t b) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const mpreal& a, const long double b) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const mpreal& a, const double b) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const mpreal& a, const unsigned long int b) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const mpreal& a, const unsigned int b) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const mpreal& a, const long int b) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const mpreal& a, const int b) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const mpz_t b, const mpreal& a) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const mpq_t b, const mpreal& a) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const long double b, const mpreal& a) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const double b, const mpreal& a) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const unsigned long int b, const mpreal& a) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const unsigned int b, const mpreal& a) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const long int b, const mpreal& a) -{ - return mpreal(a) *= b; -} - -inline const mpreal operator*(const int b, const mpreal& a) -{ - return mpreal(a) *= b; -} - -////////////////////////////////////////////////////////////////////////// -// / Division -inline mpreal& mpreal::operator/=(const mpreal& v) -{ - mpfr_div(mp,mp,v.mp,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator/=(const mpz_t v) -{ - mpfr_div_z(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator/=(const mpq_t v) -{ - mpfr_div_q(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator/=(const long double v) -{ - return *this /= mpreal(v); -} - -inline mpreal& mpreal::operator/=(const double v) -{ -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) - mpfr_div_d(mp,mp,v,default_rnd); - return *this; -#else - return *this /= mpreal(v); -#endif -} - -inline mpreal& mpreal::operator/=(const unsigned long int v) -{ - mpfr_div_ui(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator/=(const unsigned int v) -{ - mpfr_div_ui(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator/=(const long int v) -{ - mpfr_div_si(mp,mp,v,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator/=(const int v) -{ - mpfr_div_si(mp,mp,v,default_rnd); - return *this; -} - -inline const mpreal operator/(const mpreal& a, const mpreal& b) -{ - mpreal x(a); - mp_prec_t pb; - mp_prec_t pa; - - // prec(a/b) = max(prec(a),prec(b)) - pa = a.get_prec(); - pb = b.get_prec(); - if(pb>pa) x.set_prec(pb); - - return x /= b; -} - -inline const mpreal operator/(const mpreal& a, const mpz_t b) -{ - return mpreal(a) /= b; -} - -inline const mpreal operator/(const mpreal& a, const mpq_t b) -{ - return mpreal(a) /= b; -} - -inline const mpreal operator/(const mpreal& a, const long double b) -{ - return mpreal(a) /= b; -} - -inline const mpreal operator/(const mpreal& a, const double b) -{ - return mpreal(a) /= b; -} - -inline const mpreal operator/(const mpreal& a, const unsigned long int b) -{ - return mpreal(a) /= b; -} - -inline const mpreal operator/(const mpreal& a, const unsigned int b) -{ - return mpreal(a) /= b; -} - -inline const mpreal operator/(const mpreal& a, const long int b) -{ - return mpreal(a) /= b; -} - -inline const mpreal operator/(const mpreal& a, const int b) -{ - return mpreal(a) /= b; -} - -inline const mpreal operator/(const unsigned long int b, const mpreal& a) -{ - mpreal x(a); - mpfr_ui_div(x.mp,b,a.mp,mpreal::default_rnd); - return x; -} - -inline const mpreal operator/(const unsigned int b, const mpreal& a) -{ - mpreal x(a); - mpfr_ui_div(x.mp,b,a.mp,mpreal::default_rnd); - return x; -} - -inline const mpreal operator/(const long int b, const mpreal& a) -{ - mpreal x(a); - mpfr_si_div(x.mp,b,a.mp,mpreal::default_rnd); - return x; -} - -inline const mpreal operator/(const int b, const mpreal& a) -{ - mpreal x(a); - mpfr_si_div(x.mp,b,a.mp,mpreal::default_rnd); - return x; -} - -inline const mpreal operator/(const long double b, const mpreal& a) -{ - mpreal x(b); - return x/a; -} - -inline const mpreal operator/(const double b, const mpreal& a) -{ -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) - mpreal x(a); - mpfr_d_div(x.mp,b,a.mp,mpreal::default_rnd); - return x; -#else - mpreal x(b); - return x/a; -#endif -} - -////////////////////////////////////////////////////////////////////////// -// Shifts operators - Multiplication/Division by power of 2 -inline mpreal& mpreal::operator<<=(const unsigned long int u) -{ - mpfr_mul_2ui(mp,mp,u,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator<<=(const unsigned int u) -{ - mpfr_mul_2ui(mp,mp,static_cast(u),default_rnd); - return *this; -} - -inline mpreal& mpreal::operator<<=(const long int u) -{ - mpfr_mul_2si(mp,mp,u,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator<<=(const int u) -{ - mpfr_mul_2si(mp,mp,static_cast(u),default_rnd); - return *this; -} - -inline mpreal& mpreal::operator>>=(const unsigned long int u) -{ - mpfr_div_2ui(mp,mp,u,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator>>=(const unsigned int u) -{ - mpfr_div_2ui(mp,mp,static_cast(u),default_rnd); - return *this; -} - -inline mpreal& mpreal::operator>>=(const long int u) -{ - mpfr_div_2si(mp,mp,u,default_rnd); - return *this; -} - -inline mpreal& mpreal::operator>>=(const int u) -{ - mpfr_div_2si(mp,mp,static_cast(u),default_rnd); - 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(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(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(k)); -} - -inline const mpreal operator>>(const mpreal& v, const int k) -{ - return div_2si(v,static_cast(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.mp,v.mp,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.mp,v.mp,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.mp,v.mp,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.mp,v.mp,k,rnd_mode); - return x; -} - -////////////////////////////////////////////////////////////////////////// -//Boolean operators -inline bool operator > (const mpreal& a, const mpreal& b) -{ - return (mpfr_greater_p(a.mp,b.mp)!=0); -} - -inline bool operator > (const mpreal& a, const unsigned long int b) -{ - return a>mpreal(b); -} - -inline bool operator > (const mpreal& a, const unsigned int b) -{ - return a>mpreal(b); -} - -inline bool operator > (const mpreal& a, const long int b) -{ - return a>mpreal(b); -} - -inline bool operator > (const mpreal& a, const int b) -{ - return a>mpreal(b); -} - -inline bool operator > (const mpreal& a, const long double b) -{ - return a>mpreal(b); -} - -inline bool operator > (const mpreal& a, const double b) -{ - return a>mpreal(b); -} - -inline bool operator > (const unsigned long int a, const mpreal& b) -{ - return mpreal(a)>b; -} - -inline bool operator > (const unsigned int a, const mpreal& b) -{ - return mpreal(a)>b; -} - -inline bool operator > (const long int a, const mpreal& b) -{ - return mpreal(a)>b; -} - -inline bool operator > (const int a, const mpreal& b) -{ - return mpreal(a)>b; -} - -inline bool operator > (const long double a, const mpreal& b) -{ - return mpreal(a)>b; -} - -inline bool operator > (const double a, const mpreal& b) -{ - return mpreal(a)>b; -} - -inline bool operator >= (const mpreal& a, const mpreal& b) -{ - return (mpfr_greaterequal_p(a.mp,b.mp)!=0); -} - -inline bool operator >= (const mpreal& a, const unsigned long int b) -{ - return a>=mpreal(b); -} - -inline bool operator >= (const mpreal& a, const unsigned int b) -{ - return a>=mpreal(b); -} - -inline bool operator >= (const mpreal& a, const long int b) -{ - return a>=mpreal(b); -} - -inline bool operator >= (const mpreal& a, const int b) -{ - return a>=mpreal(b); -} - -inline bool operator >= (const mpreal& a, const long double b) -{ - return a>=mpreal(b); -} - -inline bool operator >= (const mpreal& a, const double b) -{ - return a>=mpreal(b); -} - -inline bool operator >= (const unsigned long int a,const mpreal& b) -{ - return mpreal(a)>=b; -} - -inline bool operator >= (const unsigned int a, const mpreal& b) -{ - return mpreal(a)>=b; -} - -inline bool operator >= (const long int a, const mpreal& b) -{ - return mpreal(a)>=b; -} - -inline bool operator >= (const int a, const mpreal& b) -{ - return mpreal(a)>=b; -} - -inline bool operator >= (const long double a, const mpreal& b) -{ - return mpreal(a)>=b; -} - -inline bool operator >= (const double a, const mpreal& b) -{ - return mpreal(a)>=b; -} - -inline bool operator < (const mpreal& a, const mpreal& b) -{ - return (mpfr_less_p(a.mp,b.mp)!=0); -} - -inline bool operator < (const mpreal& a, const unsigned long int b) -{ - return a= MPFR_VERSION_NUM(3,0,0)) -inline bool _isregular(const mpreal& v) -{ - return (mpfr_regular_p(v.mp)); -} -#endif // MPFR 3.0.0 Specifics - -////////////////////////////////////////////////////////////////////////// -// Type Converters -inline mpreal::operator double() const -{ - return mpfr_get_d(mp,default_rnd); -} - -inline mpreal::operator float() const -{ - return (float)mpfr_get_d(mp,default_rnd); -} - -inline mpreal::operator long double() const -{ - return mpfr_get_ld(mp,default_rnd); -} - -inline mpreal::operator unsigned long() const -{ - return mpfr_get_ui(mp,default_rnd); -} - -inline mpreal::operator unsigned int() const -{ - return static_cast(mpfr_get_ui(mp,default_rnd)); -} - -inline mpreal::operator long() const -{ - return mpfr_get_si(mp,default_rnd); -} - -inline mpreal::operator mpfr_ptr() -{ - return mp; -} - -////////////////////////////////////////////////////////////////////////// -// Set/Get number properties -inline int sgn(const mpreal& v) -{ - int r = mpfr_signbit(v.mp); - return (r>0?-1:1); -} - -inline void mpreal::set_sign(int sign, mp_rnd_t rnd_mode) -{ - mpfr_setsign(mp,mp,(sign<0?1:0),rnd_mode); -} - -inline mp_prec_t mpreal::get_prec() const -{ - return mpfr_get_prec(mp); -} - -inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode) -{ - mpfr_prec_round(mp,prec,rnd_mode); -} - -inline void mpreal::set_inf(int sign) -{ - mpfr_set_inf(mp,sign); -} - -inline void mpreal::set_nan() -{ - mpfr_set_nan(mp); -} - -inline mp_exp_t mpreal::get_exp () -{ - return mpfr_get_exp(mp); -} - -inline int mpreal::set_exp (mp_exp_t e) -{ - return mpfr_set_exp(mp,e); -} - -inline const mpreal frexp(const mpreal& v, mp_exp_t* exp) -{ - mpreal x(v); - *exp = x.get_exp(); - x.set_exp(0); - return x; -} - -inline const mpreal ldexp(const mpreal& v, mp_exp_t exp) -{ - mpreal x(v); - - // rounding is not important since we just increasing the exponent - mpfr_mul_2si(x.mp,x.mp,exp,mpreal::default_rnd); - return x; -} - -inline const mpreal machine_epsilon(mp_prec_t prec) -{ - // smallest eps such that 1.0+eps != 1.0 - // depends (of cause) on the precision - mpreal x(1,prec); - return nextabove(x)-x; -} - -inline const mpreal mpreal_min(mp_prec_t prec) -{ - // min = 1/2*2^emin = 2^(emin-1) - - mpreal x(1,prec); - return x <<= mpreal::get_emin()-1; -} - -inline const mpreal mpreal_max(mp_prec_t prec) -{ - // max = (1-eps)*2^emax, assume eps = 0?, - // and use emax-1 to prevent value to be +inf - // max = 2^(emax-1) - - mpreal x(1,prec); - return x <<= mpreal::get_emax()-1; -} - -inline const mpreal modf(const mpreal& v, mpreal& n) -{ - mpreal frac(v); - - // rounding is not important since we are using the same number - mpfr_frac(frac.mp,frac.mp,mpreal::default_rnd); - mpfr_trunc(n.mp,v.mp); - return frac; -} - -inline int mpreal::check_range (int t, mp_rnd_t rnd_mode) -{ - return mpfr_check_range(mp,t,rnd_mode); -} - -inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode) -{ - return mpfr_subnormalize(mp,t,rnd_mode); -} - -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 -////////////////////////////////////////////////////////////////////////// -inline const mpreal sqr(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_sqr(x.mp,x.mp,rnd_mode); - return x; -} - -inline const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_sqrt(x.mp,x.mp,rnd_mode); - return x; -} - -inline const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode) -{ - mpreal x; - mpfr_sqrt_ui(x.mp,v,rnd_mode); - return x; -} - -inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode) -{ - return sqrt(static_cast(v),rnd_mode); -} - -inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode) -{ - if (v>=0) return sqrt(static_cast(v),rnd_mode); - else return mpreal(); // NaN -} - -inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) -{ - if (v>=0) return sqrt(static_cast(v),rnd_mode); - else return mpreal(); // NaN -} - -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& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_cbrt(x.mp,x.mp,rnd_mode); - return x; -} - -inline const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_root(x.mp,x.mp,k,rnd_mode); - return x; -} - -inline const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_abs(x.mp,x.mp,rnd_mode); - return x; -} - -inline const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_abs(x.mp,x.mp,rnd_mode); - return x; -} - -inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode) -{ - mpreal x(a); - mpfr_dim(x.mp,a.mp,b.mp,rnd_mode); - return x; -} - -inline int cmpabs(const mpreal& a,const mpreal& b) -{ - return mpfr_cmpabs(a.mp,b.mp); -} - -inline const mpreal log (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_log(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal log2(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_log2(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_log10(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal exp(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_exp(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal exp2(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_exp2(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_exp10(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_cos(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_sin(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_tan(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_sec(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_csc(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_cot(x.mp,v.mp,rnd_mode); - return x; -} - -inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode) -{ - return mpfr_sin_cos(s.mp,c.mp,v.mp,rnd_mode); -} - -inline const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_acos(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_asin(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_atan(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode) -{ - mpreal a; - mp_prec_t yp, xp; - - yp = y.get_prec(); - xp = x.get_prec(); - - a.set_prec(yp>xp?yp:xp); - - mpfr_atan2(a.mp, y.mp, x.mp, rnd_mode); - - return a; -} - -inline const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_cosh(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_sinh(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_tanh(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_sech(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_csch(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_coth(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_acosh(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_asinh(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_atanh(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal fac_ui (unsigned long int v, mp_rnd_t rnd_mode) -{ - mpreal x; - mpfr_fac_ui(x.mp,v,rnd_mode); - return x; -} - -inline const mpreal log1p (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_log1p(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal expm1 (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_expm1(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_eint(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_gamma(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_lngamma(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_lgamma(x.mp,signp,v.mp,rnd_mode); - return x; -} - -inline const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_zeta(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_erf(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_erfc(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal _j0 (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_j0(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal _j1 (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_j1(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal _jn (long n, const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_jn(x.mp,n,v.mp,rnd_mode); - return x; -} - -inline const mpreal _y0 (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_y0(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal _y1 (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_y1(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal _yn (long n, const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_yn(x.mp,n,v.mp,rnd_mode); - 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) -{ - return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode); -} - -inline const mpreal li2(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_li2(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) -{ - mpreal a; - mp_prec_t yp, xp; - - yp = y.get_prec(); - xp = x.get_prec(); - - 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 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& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_digamma(x.mp,v.mp,rnd_mode); - return x; -} -#endif // MPFR 3.0.0 Specifics - -////////////////////////////////////////////////////////////////////////// -// Constants -inline const mpreal const_log2 (mp_prec_t prec, mp_rnd_t rnd_mode) -{ - mpreal x; - x.set_prec(prec); - mpfr_const_log2(x.mp,rnd_mode); - return x; -} - -inline const mpreal const_pi (mp_prec_t prec, mp_rnd_t rnd_mode) -{ - mpreal x; - x.set_prec(prec); - mpfr_const_pi(x.mp,rnd_mode); - return x; -} - -inline const mpreal const_euler (mp_prec_t prec, mp_rnd_t rnd_mode) -{ - mpreal x; - x.set_prec(prec); - mpfr_const_euler(x.mp,rnd_mode); - return x; -} - -inline const mpreal const_catalan (mp_prec_t prec, mp_rnd_t rnd_mode) -{ - mpreal x; - x.set_prec(prec); - mpfr_const_catalan(x.mp,rnd_mode); - return x; -} - -inline const mpreal const_infinity (int sign, mp_prec_t prec, mp_rnd_t rnd_mode) -{ - mpreal x; - x.set_prec(prec,rnd_mode); - mpfr_set_inf(x.mp, sign); - return x; -} - -////////////////////////////////////////////////////////////////////////// -// Integer Related Functions -inline const mpreal rint(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_rint(x.mp,v.mp,rnd_mode); - return x; -} - -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_ceil (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_rint_ceil(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal rint_floor(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_rint_floor(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal rint_round(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_rint_round(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal rint_trunc(const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_rint_trunc(x.mp,v.mp,rnd_mode); - return x; -} - -inline const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode) -{ - mpreal x(v); - mpfr_frac(x.mp,v.mp,rnd_mode); - return x; -} - -////////////////////////////////////////////////////////////////////////// -// Miscellaneous Functions -inline void swap(mpreal& a, mpreal& b) -{ - mpfr_swap(a.mp,b.mp); -} - -#ifndef max -inline const mpreal max(const mpreal& x, const mpreal& y) -{ - return (x>y?x:y); -} -#endif - -#ifndef min -inline const mpreal min(const mpreal& x, const mpreal& y) -{ - return (x= MPFR_VERSION_NUM(3,0,0)) -// use gmp_randinit_default() to init state, gmp_randclear() to clear -inline const mpreal urandom (gmp_randstate_t& state,mp_rnd_t rnd_mode) -{ - mpreal x; - mpfr_urandom(x.mp,state,rnd_mode); - return x; -} -#endif - -#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) -inline const mpreal random2 (mp_size_t size, mp_exp_t exp) -{ - mpreal x; - mpfr_random2(x.mp,size,exp); - return x; -} -#endif - -////////////////////////////////////////////////////////////////////////// -// Set/Get global properties -inline void mpreal::set_default_prec(mp_prec_t prec) -{ - default_prec = prec; - mpfr_set_default_prec(prec); -} - -inline mp_prec_t mpreal::get_default_prec() -{ - return mpfr_get_default_prec(); -} - -inline void mpreal::set_default_base(int base) -{ - default_base = base; -} - -inline int mpreal::get_default_base() -{ - return default_base; -} - -inline void mpreal::set_default_rnd(mp_rnd_t rnd_mode) -{ - default_rnd = rnd_mode; - mpfr_set_default_rounding_mode(rnd_mode); -} - -inline mp_rnd_t mpreal::get_default_rnd() -{ - return mpfr_get_default_rounding_mode(); -} - -inline void mpreal::set_double_bits(int dbits) -{ - double_bits = dbits; -} - -inline int mpreal::get_double_bits() -{ - return double_bits; -} - -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 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 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 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(b),rnd_mode); -} - -inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode) -{ - 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(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 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(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(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(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(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(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(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(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(a),static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui - else return pow(static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui - else return pow(static_cast(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(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(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(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui - else return pow(mpreal(a),static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui - else return pow(static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui - else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow - }else{ - return pow(mpreal(a),static_cast(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(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(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(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui - else return pow(mpreal(a),static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui - else return pow(static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui - else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow - }else{ - return pow(mpreal(a),static_cast(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(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(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(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(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(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(b),rnd_mode); // mpfr_pow_si -} -} - -// 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 -{ - template <> - inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) - { - return mpfr::swap(x, y); - } -} - -#endif /* __MP_REAL_H__ */ +/* + Multi-precision real number class. C++ wrapper fo MPFR library. + Project homepage: http://www.holoborodko.com/pavel/ + Contact e-mail: pavel@holoborodko.com + + Copyright (c) 2008-2010 Pavel Holoborodko + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + + Contributors: + Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, + Heinz van Saanen, Pere Constans, Dmitriy Gubanov +*/ + +#ifndef __MP_REAL_H__ +#define __MP_REAL_H__ + +#include +#include +#include +#include +#include +#include + +#include + +// Detect compiler using signatures from http://predef.sourceforge.net/ +#if defined(__GNUC__) && defined(__INTEL_COMPILER) + #define IsInf(x) isinf(x) // GNU C/C++ + Intel ICC compiler + +#elif defined(__GNUC__) + #define IsInf(x) std::isinf(x) // GNU C/C++ + +#elif defined(_MSC_VER) + #define IsInf(x) (!_finite(x)) // Microsoft Visual C++ + +#else + #define IsInf(x) std::isinf(x) // C99 conformance +#endif + +namespace mpfr { + +class mpreal { +private: + mpfr_t mp; + +public: + static mp_rnd_t default_rnd; + static mp_prec_t default_prec; + static int default_base; + static int double_bits; + +public: + // Constructors && type conversion + mpreal(); + mpreal(const mpreal& u); + + mpreal(const mpfr_t u); + mpreal(const mpf_t u); + + mpreal(const mpz_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const mpq_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const double u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const long double u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const unsigned long int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const unsigned int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const long int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd); + mpreal(const char* s, mp_prec_t prec = default_prec, int base = default_base, mp_rnd_t mode = default_rnd); + + ~mpreal(); + + // 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 int v); + mpreal& operator=(const long int v); + mpreal& operator=(const int v); + mpreal& operator=(const char* s); + + // + + 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); + 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); + + // Boolean Operators + friend bool operator > (const mpreal& a, const mpreal& b); + friend bool operator >= (const mpreal& a, const mpreal& b); + friend bool operator < (const mpreal& a, const mpreal& b); + friend bool operator <= (const mpreal& a, const mpreal& b); + friend bool operator == (const mpreal& a, const mpreal& b); + friend bool operator != (const mpreal& a, const mpreal& b); + + // Type Conversion operators + inline operator long double() const; + inline operator double() const; + inline operator float() const; + inline operator unsigned long() const; + inline operator unsigned int() const; + inline operator long() const; + operator std::string() const; + inline operator mpfr_ptr(); + + // Math Functions + friend const mpreal sqr(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend int cmpabs(const mpreal& a,const mpreal& b); + + friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + + friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + + friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + + friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fac_ui (unsigned long int v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal log1p (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal expm1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + + friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _j0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _j1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _jn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _y0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _y1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal _yn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend int sgn(const mpreal& v); // -1 or +1 + +// 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 = mpreal::default_rnd); + friend const mpreal li2(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); +#endif + +// MPFR 3.0.0 Specifics +#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) + friend const mpreal digamma(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal urandom (gmp_randstate_t& state,mp_rnd_t rnd_mode = mpreal::default_rnd); // use gmp_randinit_default() to init state, gmp_randclear() to clear + friend bool _isregular(const mpreal& v); +#endif + + // Exponent and mantissa manipulation + friend const mpreal frexp(const mpreal& v, mp_exp_t* exp); + friend const mpreal ldexp(const mpreal& v, mp_exp_t exp); + + // Splits mpreal value into fractional and integer parts. + // Returns fractional part and stores integer part in n. + friend const mpreal modf(const mpreal& v, mpreal& n); + + // 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 = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal const_pi (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal const_euler (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal const_catalan (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + // returns +inf iff sign>=0 otherwise -inf + friend const mpreal const_infinity(int sign = 1, mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd); + + // 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 = mpreal::default_rnd); + 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 = mpreal::default_rnd); + friend const mpreal rint_floor(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal rint_round(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal rint_trunc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); + friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd); + + // 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(const mpreal& v); + friend bool _isinf(const mpreal& v); + friend bool _isnum(const mpreal& v); + friend bool _iszero(const mpreal& v); + friend bool _isint(const mpreal& v); + + // Set/Get instance properties + inline mp_prec_t get_prec() const; + inline void set_prec(mp_prec_t prec, mp_rnd_t rnd_mode = default_rnd); // Change precision with rounding mode + + // Set mpreal to +-inf, NaN + void set_inf(int sign = +1); + void set_nan(); + + // sign = -1 or +1 + void set_sign(int sign, mp_rnd_t rnd_mode = default_rnd); + + //Exponent + mp_exp_t get_exp(); + int set_exp(mp_exp_t e); + int check_range (int t, mp_rnd_t rnd_mode = default_rnd); + int subnormalize (int t,mp_rnd_t rnd_mode = 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 mp_prec_t get_default_prec(); + static void set_default_base(int base); + static int get_default_base(); + static void set_double_bits(int dbits); + static int get_double_bits(); + static void set_default_rnd(mp_rnd_t rnd_mode); + static mp_rnd_t get_default_rnd(); + 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); + + // Get/Set conversions + // Convert mpreal to string with n significant digits in base b + // n = 0 -> convert with the maximum available digits + std::string to_string(size_t n = 0, int b = default_base, mp_rnd_t mode = default_rnd) const; + + // Efficient swapping of two mpreal values + friend void swap(mpreal& x, mpreal& y); + + //Min Max - macros is evil. Needed for systems which defines max and min globally as macros (e.g. Windows) + //Hope that globally defined macros use > < operations only + #ifndef max + friend const mpreal max(const mpreal& x, const mpreal& y); + #endif + + #ifndef min + friend const mpreal min(const mpreal& x, const mpreal& y); + #endif +}; + +////////////////////////////////////////////////////////////////////////// +// Exceptions +class conversion_overflow : public std::exception { +public: + std::string why() { return "inexact conversion from floating point"; } +}; + +////////////////////////////////////////////////////////////////////////// +// + Addition +const mpreal operator+(const mpreal& a, const mpreal& b); + +// + Fast specialized addition - implemented through fast += operations +const mpreal operator+(const mpreal& a, const mpz_t b); +const mpreal operator+(const mpreal& a, const mpq_t b); +const mpreal operator+(const mpreal& a, const long double b); +const mpreal operator+(const mpreal& a, const double b); +const mpreal operator+(const mpreal& a, const unsigned long int b); +const mpreal operator+(const mpreal& a, const unsigned int b); +const mpreal operator+(const mpreal& a, const long int b); +const mpreal operator+(const mpreal& a, const int b); +const mpreal operator+(const mpreal& a, const char* b); +const mpreal operator+(const char* a, const mpreal& b); +const std::string operator+(const mpreal& a, const std::string b); +const std::string operator+(const std::string a, const mpreal& b); + +const mpreal operator+(const mpz_t b, const mpreal& a); +const mpreal operator+(const mpq_t b, const mpreal& a); +const mpreal operator+(const long double b, const mpreal& a); +const mpreal operator+(const double b, const mpreal& a); +const mpreal operator+(const unsigned long int b, const mpreal& a); +const mpreal operator+(const unsigned int b, const mpreal& a); +const mpreal operator+(const long int b, const mpreal& a); +const mpreal operator+(const int b, const mpreal& a); + +////////////////////////////////////////////////////////////////////////// +// - Subtraction +const mpreal operator-(const mpreal& a, const mpreal& b); + +// - Fast specialized subtraction - implemented through fast -= operations +const mpreal operator-(const mpreal& a, const mpz_t b); +const mpreal operator-(const mpreal& a, const mpq_t b); +const mpreal operator-(const mpreal& a, const long double b); +const mpreal operator-(const mpreal& a, const double b); +const mpreal operator-(const mpreal& a, const unsigned long int b); +const mpreal operator-(const mpreal& a, const unsigned int b); +const mpreal operator-(const mpreal& a, const long int b); +const mpreal operator-(const mpreal& a, const int b); +const mpreal operator-(const mpreal& a, const char* b); +const mpreal operator-(const char* a, const mpreal& b); + +const mpreal operator-(const mpz_t b, const mpreal& a); +const mpreal operator-(const mpq_t b, const mpreal& a); +const mpreal operator-(const long double b, const mpreal& a); +//const mpreal operator-(const double b, const mpreal& a); + +////////////////////////////////////////////////////////////////////////// +// * Multiplication +const mpreal operator*(const mpreal& a, const mpreal& b); + +// * Fast specialized multiplication - implemented through fast *= operations +const mpreal operator*(const mpreal& a, const mpz_t b); +const mpreal operator*(const mpreal& a, const mpq_t b); +const mpreal operator*(const mpreal& a, const long double b); +const mpreal operator*(const mpreal& a, const double b); +const mpreal operator*(const mpreal& a, const unsigned long int b); +const mpreal operator*(const mpreal& a, const unsigned int b); +const mpreal operator*(const mpreal& a, const long int b); +const mpreal operator*(const mpreal& a, const int b); + +const mpreal operator*(const mpz_t b, const mpreal& a); +const mpreal operator*(const mpq_t b, const mpreal& a); +const mpreal operator*(const long double b, const mpreal& a); +const mpreal operator*(const double b, const mpreal& a); +const mpreal operator*(const unsigned long int b, const mpreal& a); +const mpreal operator*(const unsigned int b, const mpreal& a); +const mpreal operator*(const long int b, const mpreal& a); +const mpreal operator*(const int b, const mpreal& a); + +////////////////////////////////////////////////////////////////////////// +// / Division +const mpreal operator/(const mpreal& a, const mpreal& b); + +// / Fast specialized division - implemented through fast /= operations +const mpreal operator/(const mpreal& a, const mpz_t b); +const mpreal operator/(const mpreal& a, const mpq_t b); +const mpreal operator/(const mpreal& a, const long double b); +const mpreal operator/(const mpreal& a, const double b); +const mpreal operator/(const mpreal& a, const unsigned long int b); +const mpreal operator/(const mpreal& a, const unsigned int b); +const mpreal operator/(const mpreal& a, const long int b); +const mpreal operator/(const mpreal& a, const int b); + +const mpreal operator/(const long double b, const mpreal& a); + +////////////////////////////////////////////////////////////////////////// +// Shifts operators - Multiplication/Division by a power of 2 +const mpreal operator<<(const mpreal& v, const unsigned long int k); +const mpreal operator<<(const mpreal& v, const unsigned int k); +const mpreal operator<<(const mpreal& v, const long int k); +const mpreal operator<<(const mpreal& v, const int k); + +const mpreal operator>>(const mpreal& v, const unsigned long int k); +const mpreal operator>>(const mpreal& v, const unsigned int k); +const mpreal operator>>(const mpreal& v, const long int k); +const mpreal operator>>(const mpreal& v, const int k); + +////////////////////////////////////////////////////////////////////////// +// Boolean operators +bool operator < (const mpreal& a, const unsigned long int b); +bool operator < (const mpreal& a, const unsigned int b); +bool operator < (const mpreal& a, const long int b); +bool operator < (const mpreal& a, const int b); +bool operator < (const mpreal& a, const long double b); +bool operator < (const mpreal& a, const double b); + +bool operator < (const unsigned long int a,const mpreal& b); +bool operator < (const unsigned int a, const mpreal& b); +bool operator < (const long int a, const mpreal& b); +bool operator < (const int a, const mpreal& b); +bool operator < (const long double a, const mpreal& b); +bool operator < (const double a, const mpreal& b); + +bool operator > (const mpreal& a, const unsigned long int b); +bool operator > (const mpreal& a, const unsigned int b); +bool operator > (const mpreal& a, const long int b); +bool operator > (const mpreal& a, const int b); +bool operator > (const mpreal& a, const long double b); +bool operator > (const mpreal& a, const double b); + +bool operator > (const unsigned long int a,const mpreal& b); +bool operator > (const unsigned int a, const mpreal& b); +bool operator > (const long int a, const mpreal& b); +bool operator > (const int a, const mpreal& b); +bool operator > (const long double a, const mpreal& b); +bool operator > (const double a, const mpreal& b); + +bool operator >= (const mpreal& a, const unsigned long int b); +bool operator >= (const mpreal& a, const unsigned int b); +bool operator >= (const mpreal& a, const long int b); +bool operator >= (const mpreal& a, const int b); +bool operator >= (const mpreal& a, const long double b); +bool operator >= (const mpreal& a, const double b); + +bool operator >= (const unsigned long int a,const mpreal& b); +bool operator >= (const unsigned int a, const mpreal& b); +bool operator >= (const long int a, const mpreal& b); +bool operator >= (const int a, const mpreal& b); +bool operator >= (const long double a, const mpreal& b); +bool operator >= (const double a, const mpreal& b); + +bool operator <= (const mpreal& a, const unsigned long int b); +bool operator <= (const mpreal& a, const unsigned int b); +bool operator <= (const mpreal& a, const long int b); +bool operator <= (const mpreal& a, const int b); +bool operator <= (const mpreal& a, const long double b); +bool operator <= (const mpreal& a, const double b); + +bool operator <= (const unsigned long int a,const mpreal& b); +bool operator <= (const unsigned int a, const mpreal& b); +bool operator <= (const long int a, const mpreal& b); +bool operator <= (const int a, const mpreal& b); +bool operator <= (const long double a, const mpreal& b); +bool operator <= (const double a, const mpreal& b); + +bool operator == (const mpreal& a, const unsigned long int b); +bool operator == (const mpreal& a, const unsigned int b); +bool operator == (const mpreal& a, const long int b); +bool operator == (const mpreal& a, const int b); +bool operator == (const mpreal& a, const long double b); +bool operator == (const mpreal& a, const double b); + +bool operator == (const unsigned long int a,const mpreal& b); +bool operator == (const unsigned int a, const mpreal& b); +bool operator == (const long int a, const mpreal& b); +bool operator == (const int a, const mpreal& b); +bool operator == (const long double a, const mpreal& b); +bool operator == (const double a, const mpreal& b); + +bool operator != (const mpreal& a, const unsigned long int b); +bool operator != (const mpreal& a, const unsigned int b); +bool operator != (const mpreal& a, const long int b); +bool operator != (const mpreal& a, const int b); +bool operator != (const mpreal& a, const long double b); +bool operator != (const mpreal& a, const double b); + +bool operator != (const unsigned long int a,const mpreal& b); +bool operator != (const unsigned int a, const mpreal& b); +bool operator != (const long int a, const mpreal& b); +bool operator != (const int a, const mpreal& b); +bool operator != (const long double a, const mpreal& b); +bool operator != (const double a, const mpreal& b); + +////////////////////////////////////////////////////////////////////////// +// sqrt +const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal sqrt(const long int v, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::default_rnd); + +////////////////////////////////////////////////////////////////////////// +// pow +const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd); +const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd); + +////////////////////////////////////////////////////////////////////////// +// Estimate machine epsilon for the given precision +inline const mpreal machine_epsilon(mp_prec_t prec); +inline const mpreal mpreal_min(mp_prec_t prec); +inline const mpreal mpreal_max(mp_prec_t prec); + +////////////////////////////////////////////////////////////////////////// +// Implementation of inline functions +////////////////////////////////////////////////////////////////////////// + +////////////////////////////////////////////////////////////////////////// +// Operators - Assignment +inline mpreal& mpreal::operator=(const mpreal& v) +{ + if (this!= &v) mpfr_set(mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const mpf_t v) +{ + mpfr_set_f(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const mpz_t v) +{ + mpfr_set_z(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const mpq_t v) +{ + mpfr_set_q(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const long double v) +{ + mpfr_set_ld(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const double v) +{ + if(double_bits == -1 || fits_in_bits(v, double_bits)) + { + mpfr_set_d(mp,v,default_rnd); + } + else + throw conversion_overflow(); + + return *this; +} + +inline mpreal& mpreal::operator=(const unsigned long int v) +{ + mpfr_set_ui(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const unsigned int v) +{ + mpfr_set_ui(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const long int v) +{ + mpfr_set_si(mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator=(const int v) +{ + mpfr_set_si(mp,v,default_rnd); + return *this; +} + +////////////////////////////////////////////////////////////////////////// +// + Addition +inline mpreal& mpreal::operator+=(const mpreal& v) +{ + mpfr_add(mp,mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const mpf_t u) +{ + *this += mpreal(u); + return *this; +} + +inline mpreal& mpreal::operator+=(const mpz_t u) +{ + mpfr_add_z(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const mpq_t u) +{ + mpfr_add_q(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+= (const long double u) +{ + return *this += mpreal(u); +} + +inline mpreal& mpreal::operator+= (const double u) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpfr_add_d(mp,mp,u,default_rnd); + return *this; +#else + return *this += mpreal(u); +#endif +} + +inline mpreal& mpreal::operator+=(const unsigned long int u) +{ + mpfr_add_ui(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const unsigned int u) +{ + mpfr_add_ui(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const long int u) +{ + mpfr_add_si(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator+=(const int u) +{ + mpfr_add_si(mp,mp,u,default_rnd); + return *this; +} + +inline const mpreal mpreal::operator+()const +{ + return mpreal(*this); +} + +inline const mpreal operator+(const mpreal& a, const mpreal& b) +{ + // prec(a+b) = max(prec(a),prec(b)) + if(a.get_prec()>b.get_prec()) return mpreal(a) += b; + else return mpreal(b) += a; +} + +inline const std::string operator+(const mpreal& a, const std::string b) +{ + return (std::string)a+b; +} + +inline const std::string operator+(const std::string a, const mpreal& b) +{ + return a+(std::string)b; +} + +inline const mpreal operator+(const mpreal& a, const mpz_t b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const char* b) +{ + return a+mpreal(b); +} + +inline const mpreal operator+(const char* a, const mpreal& b) +{ + return mpreal(a)+b; + +} + +inline const mpreal operator+(const mpreal& a, const mpq_t b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const long double b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const double b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const unsigned long int b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const unsigned int b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const long int b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpreal& a, const int b) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpz_t b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const mpq_t b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const long double b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const double b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const unsigned long int b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const unsigned int b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const long int b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline const mpreal operator+(const int b, const mpreal& a) +{ + return mpreal(a) += b; +} + +inline mpreal& mpreal::operator++() +{ + *this += 1; + return *this; +} + +inline const mpreal mpreal::operator++ (int) +{ + mpreal x(*this); + *this += 1; + return x; +} + +inline mpreal& mpreal::operator--() +{ + *this -= 1; + return *this; +} + +inline const mpreal mpreal::operator-- (int) +{ + mpreal x(*this); + *this -= 1; + return x; +} + +////////////////////////////////////////////////////////////////////////// +// - Subtraction +inline mpreal& mpreal::operator-= (const mpreal& v) +{ + mpfr_sub(mp,mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const mpz_t v) +{ + mpfr_sub_z(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const mpq_t v) +{ + mpfr_sub_q(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const long double v) +{ + return *this -= mpreal(v); +} + +inline mpreal& mpreal::operator-=(const double v) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpfr_sub_d(mp,mp,v,default_rnd); + return *this; +#else + return *this -= mpreal(v); +#endif +} + +inline mpreal& mpreal::operator-=(const unsigned long int v) +{ + mpfr_sub_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const unsigned int v) +{ + mpfr_sub_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const long int v) +{ + mpfr_sub_si(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator-=(const int v) +{ + mpfr_sub_si(mp,mp,v,default_rnd); + return *this; +} + +inline const mpreal mpreal::operator-()const +{ + mpreal u(*this); + mpfr_neg(u.mp,u.mp,default_rnd); + return u; +} + +inline const mpreal operator-(const mpreal& a, const mpreal& b) +{ + // prec(a-b) = max(prec(a),prec(b)) + if(a.get_prec()>b.get_prec()) return mpreal(a) -= b; + else return -(mpreal(b) -= a); +} + +inline const mpreal operator-(const mpreal& a, const mpz_t b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const mpq_t b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const long double b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const double b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const unsigned long int b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const unsigned int b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const long int b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpreal& a, const int b) +{ + return mpreal(a) -= b; +} + +inline const mpreal operator-(const mpz_t b, const mpreal& a) +{ + return -(mpreal(a) -= b); +} + +inline const mpreal operator-(const mpq_t b, const mpreal& a) +{ + return -(mpreal(a) -= b); +} + +inline const mpreal operator-(const long double b, const mpreal& a) +{ + return -(mpreal(a) -= b); +} + +inline const mpreal operator-(const double b, const mpreal& a) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpreal x(a); + mpfr_d_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +#else + return -(mpreal(a) -= b); +#endif +} + +inline const mpreal operator-(const unsigned long int b, const mpreal& a) +{ + mpreal x(a); + mpfr_ui_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator-(const unsigned int b, const mpreal& a) +{ + mpreal x(a); + mpfr_ui_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator-(const long int b, const mpreal& a) +{ + mpreal x(a); + mpfr_si_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator-(const int b, const mpreal& a) +{ + mpreal x(a); + mpfr_si_sub(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator-(const mpreal& a, const char* b) +{ + return a-mpreal(b); +} + +inline const mpreal operator-(const char* a, const mpreal& b) +{ + return mpreal(a)-b; +} + +////////////////////////////////////////////////////////////////////////// +// * Multiplication +inline mpreal& mpreal::operator*= (const mpreal& v) +{ + mpfr_mul(mp,mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const mpz_t v) +{ + mpfr_mul_z(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const mpq_t v) +{ + mpfr_mul_q(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const long double v) +{ + return *this *= mpreal(v); +} + +inline mpreal& mpreal::operator*=(const double v) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpfr_mul_d(mp,mp,v,default_rnd); + return *this; +#else + return *this *= mpreal(v); +#endif +} + +inline mpreal& mpreal::operator*=(const unsigned long int v) +{ + mpfr_mul_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const unsigned int v) +{ + mpfr_mul_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const long int v) +{ + mpfr_mul_si(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator*=(const int v) +{ + mpfr_mul_si(mp,mp,v,default_rnd); + return *this; +} + +inline const mpreal operator*(const mpreal& a, const mpreal& b) +{ + // prec(a*b) = max(prec(a),prec(b)) + if(a.get_prec()>b.get_prec()) return mpreal(a) *= b; + else return mpreal(b) *= a; +} + +inline const mpreal operator*(const mpreal& a, const mpz_t b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const mpq_t b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const long double b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const double b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const unsigned long int b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const unsigned int b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const long int b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpreal& a, const int b) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpz_t b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const mpq_t b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const long double b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const double b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const unsigned long int b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const unsigned int b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const long int b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +inline const mpreal operator*(const int b, const mpreal& a) +{ + return mpreal(a) *= b; +} + +////////////////////////////////////////////////////////////////////////// +// / Division +inline mpreal& mpreal::operator/=(const mpreal& v) +{ + mpfr_div(mp,mp,v.mp,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const mpz_t v) +{ + mpfr_div_z(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const mpq_t v) +{ + mpfr_div_q(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const long double v) +{ + return *this /= mpreal(v); +} + +inline mpreal& mpreal::operator/=(const double v) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpfr_div_d(mp,mp,v,default_rnd); + return *this; +#else + return *this /= mpreal(v); +#endif +} + +inline mpreal& mpreal::operator/=(const unsigned long int v) +{ + mpfr_div_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const unsigned int v) +{ + mpfr_div_ui(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const long int v) +{ + mpfr_div_si(mp,mp,v,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator/=(const int v) +{ + mpfr_div_si(mp,mp,v,default_rnd); + return *this; +} + +inline const mpreal operator/(const mpreal& a, const mpreal& b) +{ + mpreal x(a); + mp_prec_t pb; + mp_prec_t pa; + + // prec(a/b) = max(prec(a),prec(b)) + pa = a.get_prec(); + pb = b.get_prec(); + if(pb>pa) x.set_prec(pb); + + return x /= b; +} + +inline const mpreal operator/(const mpreal& a, const mpz_t b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const mpq_t b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const long double b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const double b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const unsigned long int b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const unsigned int b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const long int b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const mpreal& a, const int b) +{ + return mpreal(a) /= b; +} + +inline const mpreal operator/(const unsigned long int b, const mpreal& a) +{ + mpreal x(a); + mpfr_ui_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator/(const unsigned int b, const mpreal& a) +{ + mpreal x(a); + mpfr_ui_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator/(const long int b, const mpreal& a) +{ + mpreal x(a); + mpfr_si_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator/(const int b, const mpreal& a) +{ + mpreal x(a); + mpfr_si_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +} + +inline const mpreal operator/(const long double b, const mpreal& a) +{ + mpreal x(b); + return x/a; +} + +inline const mpreal operator/(const double b, const mpreal& a) +{ +#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) + mpreal x(a); + mpfr_d_div(x.mp,b,a.mp,mpreal::default_rnd); + return x; +#else + mpreal x(b); + return x/a; +#endif +} + +////////////////////////////////////////////////////////////////////////// +// Shifts operators - Multiplication/Division by power of 2 +inline mpreal& mpreal::operator<<=(const unsigned long int u) +{ + mpfr_mul_2ui(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator<<=(const unsigned int u) +{ + mpfr_mul_2ui(mp,mp,static_cast(u),default_rnd); + return *this; +} + +inline mpreal& mpreal::operator<<=(const long int u) +{ + mpfr_mul_2si(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator<<=(const int u) +{ + mpfr_mul_2si(mp,mp,static_cast(u),default_rnd); + return *this; +} + +inline mpreal& mpreal::operator>>=(const unsigned long int u) +{ + mpfr_div_2ui(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator>>=(const unsigned int u) +{ + mpfr_div_2ui(mp,mp,static_cast(u),default_rnd); + return *this; +} + +inline mpreal& mpreal::operator>>=(const long int u) +{ + mpfr_div_2si(mp,mp,u,default_rnd); + return *this; +} + +inline mpreal& mpreal::operator>>=(const int u) +{ + mpfr_div_2si(mp,mp,static_cast(u),default_rnd); + 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(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(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(k)); +} + +inline const mpreal operator>>(const mpreal& v, const int k) +{ + return div_2si(v,static_cast(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.mp,v.mp,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.mp,v.mp,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.mp,v.mp,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.mp,v.mp,k,rnd_mode); + return x; +} + +////////////////////////////////////////////////////////////////////////// +//Boolean operators +inline bool operator > (const mpreal& a, const mpreal& b) +{ + return (mpfr_greater_p(a.mp,b.mp)!=0); +} + +inline bool operator > (const mpreal& a, const unsigned long int b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const unsigned int b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const long int b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const int b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const long double b) +{ + return a>mpreal(b); +} + +inline bool operator > (const mpreal& a, const double b) +{ + return a>mpreal(b); +} + +inline bool operator > (const unsigned long int a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const unsigned int a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const long int a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const int a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const long double a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator > (const double a, const mpreal& b) +{ + return mpreal(a)>b; +} + +inline bool operator >= (const mpreal& a, const mpreal& b) +{ + return (mpfr_greaterequal_p(a.mp,b.mp)!=0); +} + +inline bool operator >= (const mpreal& a, const unsigned long int b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const unsigned int b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const long int b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const int b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const long double b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const mpreal& a, const double b) +{ + return a>=mpreal(b); +} + +inline bool operator >= (const unsigned long int a,const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const unsigned int a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const long int a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const int a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const long double a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator >= (const double a, const mpreal& b) +{ + return mpreal(a)>=b; +} + +inline bool operator < (const mpreal& a, const mpreal& b) +{ + return (mpfr_less_p(a.mp,b.mp)!=0); +} + +inline bool operator < (const mpreal& a, const unsigned long int b) +{ + return a= MPFR_VERSION_NUM(3,0,0)) +inline bool _isregular(const mpreal& v) +{ + return (mpfr_regular_p(v.mp)); +} +#endif // MPFR 3.0.0 Specifics + +////////////////////////////////////////////////////////////////////////// +// Type Converters +inline mpreal::operator double() const +{ + return mpfr_get_d(mp,default_rnd); +} + +inline mpreal::operator float() const +{ + return (float)mpfr_get_d(mp,default_rnd); +} + +inline mpreal::operator long double() const +{ + return mpfr_get_ld(mp,default_rnd); +} + +inline mpreal::operator unsigned long() const +{ + return mpfr_get_ui(mp,default_rnd); +} + +inline mpreal::operator unsigned int() const +{ + return static_cast(mpfr_get_ui(mp,default_rnd)); +} + +inline mpreal::operator long() const +{ + return mpfr_get_si(mp,default_rnd); +} + +inline mpreal::operator mpfr_ptr() +{ + return mp; +} + +////////////////////////////////////////////////////////////////////////// +// Set/Get number properties +inline int sgn(const mpreal& v) +{ + int r = mpfr_signbit(v.mp); + return (r>0?-1:1); +} + +inline void mpreal::set_sign(int sign, mp_rnd_t rnd_mode) +{ + mpfr_setsign(mp,mp,(sign<0?1:0),rnd_mode); +} + +inline mp_prec_t mpreal::get_prec() const +{ + return mpfr_get_prec(mp); +} + +inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpfr_prec_round(mp,prec,rnd_mode); +} + +inline void mpreal::set_inf(int sign) +{ + mpfr_set_inf(mp,sign); +} + +inline void mpreal::set_nan() +{ + mpfr_set_nan(mp); +} + +inline mp_exp_t mpreal::get_exp () +{ + return mpfr_get_exp(mp); +} + +inline int mpreal::set_exp (mp_exp_t e) +{ + return mpfr_set_exp(mp,e); +} + +inline const mpreal frexp(const mpreal& v, mp_exp_t* exp) +{ + mpreal x(v); + *exp = x.get_exp(); + x.set_exp(0); + return x; +} + +inline const mpreal ldexp(const mpreal& v, mp_exp_t exp) +{ + mpreal x(v); + + // rounding is not important since we just increasing the exponent + mpfr_mul_2si(x.mp,x.mp,exp,mpreal::default_rnd); + return x; +} + +inline const mpreal machine_epsilon(mp_prec_t prec) +{ + // smallest eps such that 1.0+eps != 1.0 + // depends (of cause) on the precision + mpreal x(1,prec); + return nextabove(x)-x; +} + +inline const mpreal mpreal_min(mp_prec_t prec) +{ + // min = 1/2*2^emin = 2^(emin-1) + + mpreal x(1,prec); + return x <<= mpreal::get_emin()-1; +} + +inline const mpreal mpreal_max(mp_prec_t prec) +{ + // max = (1-eps)*2^emax, assume eps = 0?, + // and use emax-1 to prevent value to be +inf + // max = 2^(emax-1) + + mpreal x(1,prec); + return x <<= mpreal::get_emax()-1; +} + +inline const mpreal modf(const mpreal& v, mpreal& n) +{ + mpreal frac(v); + + // rounding is not important since we are using the same number + mpfr_frac(frac.mp,frac.mp,mpreal::default_rnd); + mpfr_trunc(n.mp,v.mp); + return frac; +} + +inline int mpreal::check_range (int t, mp_rnd_t rnd_mode) +{ + return mpfr_check_range(mp,t,rnd_mode); +} + +inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode) +{ + return mpfr_subnormalize(mp,t,rnd_mode); +} + +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 +////////////////////////////////////////////////////////////////////////// +inline const mpreal sqr(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sqr(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sqrt(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode) +{ + mpreal x; + mpfr_sqrt_ui(x.mp,v,rnd_mode); + return x; +} + +inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode) +{ + return sqrt(static_cast(v),rnd_mode); +} + +inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode) +{ + if (v>=0) return sqrt(static_cast(v),rnd_mode); + else return mpreal(); // NaN +} + +inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) +{ + if (v>=0) return sqrt(static_cast(v),rnd_mode); + else return mpreal(); // NaN +} + +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& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_cbrt(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_root(x.mp,x.mp,k,rnd_mode); + return x; +} + +inline const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_abs(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_abs(x.mp,x.mp,rnd_mode); + return x; +} + +inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode) +{ + mpreal x(a); + mpfr_dim(x.mp,a.mp,b.mp,rnd_mode); + return x; +} + +inline int cmpabs(const mpreal& a,const mpreal& b) +{ + return mpfr_cmpabs(a.mp,b.mp); +} + +inline const mpreal log (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_log(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal log2(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_log2(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_log10(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal exp(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_exp(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal exp2(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_exp2(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_exp10(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_cos(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sin(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_tan(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sec(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_csc(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_cot(x.mp,v.mp,rnd_mode); + return x; +} + +inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode) +{ + return mpfr_sin_cos(s.mp,c.mp,v.mp,rnd_mode); +} + +inline const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_acos(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_asin(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_atan(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode) +{ + mpreal a; + mp_prec_t yp, xp; + + yp = y.get_prec(); + xp = x.get_prec(); + + a.set_prec(yp>xp?yp:xp); + + mpfr_atan2(a.mp, y.mp, x.mp, rnd_mode); + + return a; +} + +inline const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_cosh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sinh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_tanh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_sech(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_csch(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_coth(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_acosh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_asinh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_atanh(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal fac_ui (unsigned long int v, mp_rnd_t rnd_mode) +{ + mpreal x; + mpfr_fac_ui(x.mp,v,rnd_mode); + return x; +} + +inline const mpreal log1p (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_log1p(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal expm1 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_expm1(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_eint(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_gamma(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_lngamma(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_lgamma(x.mp,signp,v.mp,rnd_mode); + return x; +} + +inline const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_zeta(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_erf(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_erfc(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _j0 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_j0(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _j1 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_j1(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _jn (long n, const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_jn(x.mp,n,v.mp,rnd_mode); + return x; +} + +inline const mpreal _y0 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_y0(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _y1 (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_y1(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal _yn (long n, const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_yn(x.mp,n,v.mp,rnd_mode); + 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) +{ + return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode); +} + +inline const mpreal li2(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_li2(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) +{ + mpreal a; + mp_prec_t yp, xp; + + yp = y.get_prec(); + xp = x.get_prec(); + + 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 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& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_digamma(x.mp,v.mp,rnd_mode); + return x; +} +#endif // MPFR 3.0.0 Specifics + +////////////////////////////////////////////////////////////////////////// +// Constants +inline const mpreal const_log2 (mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec); + mpfr_const_log2(x.mp,rnd_mode); + return x; +} + +inline const mpreal const_pi (mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec); + mpfr_const_pi(x.mp,rnd_mode); + return x; +} + +inline const mpreal const_euler (mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec); + mpfr_const_euler(x.mp,rnd_mode); + return x; +} + +inline const mpreal const_catalan (mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec); + mpfr_const_catalan(x.mp,rnd_mode); + return x; +} + +inline const mpreal const_infinity (int sign, mp_prec_t prec, mp_rnd_t rnd_mode) +{ + mpreal x; + x.set_prec(prec,rnd_mode); + mpfr_set_inf(x.mp, sign); + return x; +} + +////////////////////////////////////////////////////////////////////////// +// Integer Related Functions +inline const mpreal rint(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint(x.mp,v.mp,rnd_mode); + return x; +} + +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_ceil (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint_ceil(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal rint_floor(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint_floor(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal rint_round(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint_round(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal rint_trunc(const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_rint_trunc(x.mp,v.mp,rnd_mode); + return x; +} + +inline const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode) +{ + mpreal x(v); + mpfr_frac(x.mp,v.mp,rnd_mode); + return x; +} + +////////////////////////////////////////////////////////////////////////// +// Miscellaneous Functions +inline void swap(mpreal& a, mpreal& b) +{ + mpfr_swap(a.mp,b.mp); +} + +#ifndef max +inline const mpreal max(const mpreal& x, const mpreal& y) +{ + return (x>y?x:y); +} +#endif + +#ifndef min +inline const mpreal min(const mpreal& x, const mpreal& y) +{ + return (x= MPFR_VERSION_NUM(3,0,0)) +// use gmp_randinit_default() to init state, gmp_randclear() to clear +inline const mpreal urandom (gmp_randstate_t& state,mp_rnd_t rnd_mode) +{ + mpreal x; + mpfr_urandom(x.mp,state,rnd_mode); + return x; +} +#endif + +#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) +inline const mpreal random2 (mp_size_t size, mp_exp_t exp) +{ + mpreal x; + mpfr_random2(x.mp,size,exp); + return x; +} +#endif + +////////////////////////////////////////////////////////////////////////// +// Set/Get global properties +inline void mpreal::set_default_prec(mp_prec_t prec) +{ + default_prec = prec; + mpfr_set_default_prec(prec); +} + +inline mp_prec_t mpreal::get_default_prec() +{ + return mpfr_get_default_prec(); +} + +inline void mpreal::set_default_base(int base) +{ + default_base = base; +} + +inline int mpreal::get_default_base() +{ + return default_base; +} + +inline void mpreal::set_default_rnd(mp_rnd_t rnd_mode) +{ + default_rnd = rnd_mode; + mpfr_set_default_rounding_mode(rnd_mode); +} + +inline mp_rnd_t mpreal::get_default_rnd() +{ + return mpfr_get_default_rounding_mode(); +} + +inline void mpreal::set_double_bits(int dbits) +{ + double_bits = dbits; +} + +inline int mpreal::get_double_bits() +{ + return double_bits; +} + +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 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 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 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(b),rnd_mode); +} + +inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode) +{ + 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(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 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(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(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(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(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(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(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(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(a),static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(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(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(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(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(mpreal(a),static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + }else{ + return pow(mpreal(a),static_cast(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(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(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(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(mpreal(a),static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(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(a),static_cast(b),rnd_mode); //mpfr_ui_pow_ui + else return pow(static_cast(a),mpreal(b),rnd_mode); //mpfr_ui_pow + }else{ + return pow(mpreal(a),static_cast(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(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(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(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(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(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(b),rnd_mode); // mpfr_pow_si +} +} + +// 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 +{ + template <> + inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) + { + return mpfr::swap(x, y); + } +} + +#endif /* __MP_REAL_H__ */