/* | |
Multi-precision floating point number class for C++. | |
Based on MPFR library: http://mpfr.org | |
Project homepage: http://www.holoborodko.com/pavel/mpfr | |
Contact e-mail: pavel@holoborodko.com | |
Copyright (c) 2008-2012 Pavel Holoborodko | |
Contributors: | |
Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman, | |
Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen, | |
Pere Constans, Peter van Hoof, Gael Guennebaud, Tsai Chia Cheng, | |
Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood. | |
**************************************************************************** | |
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 | |
**************************************************************************** | |
Redistribution and use in source and binary forms, with or without | |
modification, are permitted provided that the following conditions | |
are met: | |
1. Redistributions of source code must retain the above copyright | |
notice, this list of conditions and the following disclaimer. | |
2. Redistributions in binary form must reproduce the above copyright | |
notice, this list of conditions and the following disclaimer in the | |
documentation and/or other materials provided with the distribution. | |
3. The name of the author may not be used to endorse or promote products | |
derived from this software without specific prior written permission. | |
THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND | |
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | |
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | |
ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE | |
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | |
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | |
OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | |
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | |
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | |
OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | |
SUCH DAMAGE. | |
*/ | |
#ifndef __MPREAL_H__ | |
#define __MPREAL_H__ | |
#include <string> | |
#include <iostream> | |
#include <sstream> | |
#include <stdexcept> | |
#include <cfloat> | |
#include <cmath> | |
#include <limits> | |
// Options | |
#define MPREAL_HAVE_INT64_SUPPORT // Enable int64_t support if possible. Available only for MSVC 2010 & GCC. | |
#define MPREAL_HAVE_MSVC_DEBUGVIEW // Enable Debugger Visualizer for "Debug" builds in MSVC. | |
// Detect compiler using signatures from http://predef.sourceforge.net/ | |
#if defined(__GNUC__) && defined(__INTEL_COMPILER) | |
#define IsInf(x) isinf(x) // Intel ICC compiler on Linux | |
#elif defined(_MSC_VER) // Microsoft Visual C++ | |
#define IsInf(x) (!_finite(x)) | |
#else | |
#define IsInf(x) std::isinf(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance | |
#endif | |
#if defined(MPREAL_HAVE_INT64_SUPPORT) | |
#define MPFR_USE_INTMAX_T // Should be defined before mpfr.h | |
#if defined(_MSC_VER) // MSVC + Windows | |
#if (_MSC_VER >= 1600) | |
#include <stdint.h> // <stdint.h> is available only in msvc2010! | |
#else // MPFR relies on intmax_t which is available only in msvc2010 | |
#undef MPREAL_HAVE_INT64_SUPPORT // Besides, MPFR & MPIR have to be compiled with msvc2010 | |
#undef MPFR_USE_INTMAX_T // Since we cannot detect this, disable x64 by default | |
// Someone should change this manually if needed. | |
#endif | |
#elif defined (__GNUC__) && defined(__linux__) | |
#if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64) || defined(__ia64) || defined(__itanium__) || defined(_M_IA64) | |
#undef MPREAL_HAVE_INT64_SUPPORT // Remove all shaman dances for x64 builds since | |
#undef MPFR_USE_INTMAX_T // GCC already supports x64 as of "long int" is 64-bit integer, nothing left to do | |
#else | |
#include <stdint.h> // use int64_t, uint64_t otherwise | |
#endif | |
#else | |
#include <stdint.h> // rely on int64_t, uint64_t in all other cases, Mac OSX, etc. | |
#endif | |
#endif | |
#if defined(MPREAL_HAVE_MSVC_DEBUGVIEW) && defined(_MSC_VER) && defined(_DEBUG) | |
#define MPREAL_MSVC_DEBUGVIEW_CODE DebugView = toString(); | |
#define MPREAL_MSVC_DEBUGVIEW_DATA std::string DebugView; | |
#else | |
#define MPREAL_MSVC_DEBUGVIEW_CODE | |
#define MPREAL_MSVC_DEBUGVIEW_DATA | |
#endif | |
#include <mpfr.h> | |
#if (MPFR_VERSION < MPFR_VERSION_NUM(3,0,0)) | |
#include <cstdlib> // Needed for random() | |
#endif | |
// Less important options | |
#define MPREAL_DOUBLE_BITS_OVERFLOW -1 // Triggers overflow exception during conversion to double if mpreal | |
// cannot fit in MPREAL_DOUBLE_BITS_OVERFLOW bits | |
// = -1 disables overflow checks (default) | |
#if defined(__GNUC__) | |
#define MPREAL_PERMISSIVE_EXPR __extension__ | |
#else | |
#define MPREAL_PERMISSIVE_EXPR | |
#endif | |
namespace mpfr { | |
class mpreal { | |
private: | |
mpfr_t mp; | |
public: | |
// Get default rounding mode & precision | |
inline static mp_rnd_t get_default_rnd() { return (mp_rnd_t)(mpfr_get_default_rounding_mode()); } | |
inline static mp_prec_t get_default_prec() { return mpfr_get_default_prec(); } | |
// Constructors && type conversions | |
mpreal(); | |
mpreal(const mpreal& u); | |
mpreal(const mpfr_t u); | |
mpreal(const mpf_t u); | |
mpreal(const mpz_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
mpreal(const mpq_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
mpreal(const double u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
mpreal(const long double u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
mpreal(const unsigned long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
mpreal(const unsigned int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
mpreal(const long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
mpreal(const int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
#if defined (MPREAL_HAVE_INT64_SUPPORT) | |
mpreal(const uint64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
mpreal(const int64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); | |
#endif | |
mpreal(const char* s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, mp_rnd_t mode = mpreal::get_default_rnd()); | |
mpreal(const std::string& s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, mp_rnd_t mode = mpreal::get_default_rnd()); | |
~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 std::string& 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); | |
#if defined (MPREAL_HAVE_INT64_SUPPORT) | |
mpreal& operator+=(const int64_t u); | |
mpreal& operator+=(const uint64_t u); | |
mpreal& operator-=(const int64_t u); | |
mpreal& operator-=(const uint64_t u); | |
mpreal& operator*=(const int64_t u); | |
mpreal& operator*=(const uint64_t u); | |
mpreal& operator/=(const int64_t u); | |
mpreal& operator/=(const uint64_t u); | |
#endif | |
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); | |
// Optimized specializations for boolean operators | |
friend bool operator == (const mpreal& a, const unsigned long int b); | |
friend bool operator == (const mpreal& a, const unsigned int b); | |
friend bool operator == (const mpreal& a, const long int b); | |
friend bool operator == (const mpreal& a, const int b); | |
friend bool operator == (const mpreal& a, const long double b); | |
friend bool operator == (const mpreal& a, const double b); | |
// Type Conversion operators | |
long toLong (mp_rnd_t mode = GMP_RNDZ) const; | |
unsigned long toULong (mp_rnd_t mode = GMP_RNDZ) const; | |
double toDouble (mp_rnd_t mode = GMP_RNDN) const; | |
long double toLDouble (mp_rnd_t mode = GMP_RNDN) const; | |
#if defined (MPREAL_HAVE_INT64_SUPPORT) | |
int64_t toInt64 (mp_rnd_t mode = GMP_RNDZ) const; | |
uint64_t toUInt64 (mp_rnd_t mode = GMP_RNDZ) const; | |
#endif | |
// Get raw pointers so that mpreal can be directly used in raw mpfr_* functions | |
::mpfr_ptr mpfr_ptr(); | |
::mpfr_srcptr mpfr_ptr() const; | |
::mpfr_srcptr mpfr_srcptr() const; | |
// Convert mpreal to string with n significant digits in base b | |
// n = 0 -> convert with the maximum available digits | |
std::string toString(int n = 0, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const; | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) | |
std::string toString(const std::string& format) const; | |
#endif | |
// Math Functions | |
friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend int cmpabs(const mpreal& a,const mpreal& b); | |
friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal acot (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal asec (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal acsc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal fac_ui (unsigned long int v, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal besseljn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal bessely0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal bessely1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend int sgn(const mpreal& v); // returns -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::get_default_rnd()); | |
friend const mpreal li2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal rec_sqrt (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
// MATLAB's semantic equivalents | |
friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // Remainder after division | |
friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // Modulus after division | |
#endif | |
// MPFR 3.0.0 Specifics | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) | |
friend const mpreal digamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal ai (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); // use gmp_randinit_default() to init state, gmp_randclear() to clear | |
#endif | |
// Uniformly distributed random number generation in [0,1] using | |
// Mersenne-Twister algorithm by default. | |
// Use parameter to setup seed, e.g.: random((unsigned)time(NULL)) | |
// Check urandom() for more precise control. | |
friend const mpreal random(unsigned int seed = 0); | |
// 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::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal const_pi (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal const_euler (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal const_catalan (mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
// returns +inf iff sign>=0 otherwise -inf | |
friend const mpreal const_infinity(int sign = 1, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t rnd_mode = mpreal::get_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::get_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::get_default_rnd()); | |
friend const mpreal rint_floor (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal rint_round (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal rint_trunc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal remainder ( const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_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 isfinite (const mpreal& v); | |
friend bool isnum (const mpreal& v); | |
friend bool iszero (const mpreal& v); | |
friend bool isint (const mpreal& v); | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) | |
friend bool isregular(const mpreal& v); | |
#endif | |
// Set/Get instance properties | |
inline mp_prec_t get_prec() const; | |
inline void set_prec(mp_prec_t prec, mp_rnd_t rnd_mode = get_default_rnd()); // Change precision with rounding mode | |
// Aliases for get_prec(), set_prec() - needed for compatibility with std::complex<mpreal> interface | |
inline mpreal& setPrecision(int Precision, mp_rnd_t RoundingMode = get_default_rnd()); | |
inline int getPrecision() const; | |
// Set mpreal to +/- inf, NaN, +/-0 | |
mpreal& setInf (int Sign = +1); | |
mpreal& setNan (); | |
mpreal& setZero (int Sign = +1); | |
mpreal& setSign (int Sign, mp_rnd_t RoundingMode = get_default_rnd()); | |
//Exponent | |
mp_exp_t get_exp(); | |
int set_exp(mp_exp_t e); | |
int check_range (int t, mp_rnd_t rnd_mode = get_default_rnd()); | |
int subnormalize (int t,mp_rnd_t rnd_mode = get_default_rnd()); | |
// Inexact conversion from float | |
inline bool fits_in_bits(double x, int n); | |
// Set/Get global properties | |
static void set_default_prec(mp_prec_t prec); | |
static void set_default_rnd(mp_rnd_t rnd_mode); | |
static mp_exp_t get_emin (void); | |
static mp_exp_t get_emax (void); | |
static mp_exp_t get_emin_min (void); | |
static mp_exp_t get_emin_max (void); | |
static mp_exp_t get_emax_min (void); | |
static mp_exp_t get_emax_max (void); | |
static int set_emin (mp_exp_t exp); | |
static int set_emax (mp_exp_t exp); | |
// Efficient swapping of two mpreal values - needed for std algorithms | |
friend void swap(mpreal& x, mpreal& y); | |
friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
private: | |
// Human friendly Debug Preview in Visual Studio. | |
// Put one of these lines: | |
// | |
// mpfr::mpreal=<DebugView> ; Show value only | |
// mpfr::mpreal=<DebugView>, <mp[0]._mpfr_prec,u>bits ; Show value & precision | |
// | |
// at the beginning of | |
// [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat | |
MPREAL_MSVC_DEBUGVIEW_DATA | |
}; | |
////////////////////////////////////////////////////////////////////////// | |
// Exceptions | |
class conversion_overflow : public std::exception { | |
public: | |
std::string why() { return "inexact conversion from floating point"; } | |
}; | |
////////////////////////////////////////////////////////////////////////// | |
// Constructors & converters | |
// Default constructor: creates mp number and initializes it to 0. | |
inline mpreal::mpreal() | |
{ | |
mpfr_init2(mp,mpreal::get_default_prec()); | |
mpfr_set_ui(mp,0,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const mpreal& u) | |
{ | |
mpfr_init2(mp,mpfr_get_prec(u.mp)); | |
mpfr_set(mp,u.mp,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const mpfr_t u) | |
{ | |
mpfr_init2(mp,mpfr_get_prec(u)); | |
mpfr_set(mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const mpf_t u) | |
{ | |
mpfr_init2(mp,(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t) | |
mpfr_set_f(mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp,prec); | |
mpfr_set_z(mp,u,mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp,prec); | |
mpfr_set_q(mp,u,mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp, prec); | |
#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1) | |
if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW)) | |
{ | |
mpfr_set_d(mp, u, mode); | |
}else | |
throw conversion_overflow(); | |
#else | |
mpfr_set_d(mp, u, mode); | |
#endif | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp,prec); | |
mpfr_set_ld(mp,u,mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp,prec); | |
mpfr_set_ui(mp,u,mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp,prec); | |
mpfr_set_ui(mp,u,mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp,prec); | |
mpfr_set_si(mp,u,mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp,prec); | |
mpfr_set_si(mp,u,mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
#if defined (MPREAL_HAVE_INT64_SUPPORT) | |
inline mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp,prec); | |
mpfr_set_uj(mp, u, mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp,prec); | |
mpfr_set_sj(mp, u, mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
#endif | |
inline mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp, prec); | |
mpfr_set_str(mp, s, base, mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode) | |
{ | |
mpfr_init2(mp, prec); | |
mpfr_set_str(mp, s.c_str(), base, mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mpreal::~mpreal() | |
{ | |
mpfr_clear(mp); | |
} | |
// internal namespace needed for template magic | |
namespace internal{ | |
// Use SFINAE to restrict arithmetic operations instantiation only for numeric types | |
// This is needed for smooth integration with libraries based on expression templates, like Eigen. | |
// TODO: Do the same for boolean operators. | |
template <typename ArgumentType> struct result_type {}; | |
template <> struct result_type<mpreal> {typedef mpreal type;}; | |
template <> struct result_type<mpz_t> {typedef mpreal type;}; | |
template <> struct result_type<mpq_t> {typedef mpreal type;}; | |
template <> struct result_type<long double> {typedef mpreal type;}; | |
template <> struct result_type<double> {typedef mpreal type;}; | |
template <> struct result_type<unsigned long int> {typedef mpreal type;}; | |
template <> struct result_type<unsigned int> {typedef mpreal type;}; | |
template <> struct result_type<long int> {typedef mpreal type;}; | |
template <> struct result_type<int> {typedef mpreal type;}; | |
#if defined (MPREAL_HAVE_INT64_SUPPORT) | |
template <> struct result_type<int64_t > {typedef mpreal type;}; | |
template <> struct result_type<uint64_t > {typedef mpreal type;}; | |
#endif | |
} | |
// + Addition | |
template <typename Rhs> | |
inline const typename internal::result_type<Rhs>::type | |
operator+(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) += rhs; } | |
template <typename Lhs> | |
inline const typename internal::result_type<Lhs>::type | |
operator+(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) += lhs; } | |
// - Subtraction | |
template <typename Rhs> | |
inline const typename internal::result_type<Rhs>::type | |
operator-(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) -= rhs; } | |
template <typename Lhs> | |
inline const typename internal::result_type<Lhs>::type | |
operator-(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) -= rhs; } | |
// * Multiplication | |
template <typename Rhs> | |
inline const typename internal::result_type<Rhs>::type | |
operator*(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) *= rhs; } | |
template <typename Lhs> | |
inline const typename internal::result_type<Lhs>::type | |
operator*(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) *= lhs; } | |
// / Division | |
template <typename Rhs> | |
inline const typename internal::result_type<Rhs>::type | |
operator/(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) /= rhs; } | |
template <typename Lhs> | |
inline const typename internal::result_type<Lhs>::type | |
operator/(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) /= rhs; } | |
////////////////////////////////////////////////////////////////////////// | |
// sqrt | |
const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal sqrt(const long int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
////////////////////////////////////////////////////////////////////////// | |
// pow | |
const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); | |
////////////////////////////////////////////////////////////////////////// | |
// Estimate machine epsilon for the given precision | |
// Returns smallest eps such that 1.0 + eps != 1.0 | |
inline mpreal machine_epsilon(mp_prec_t prec = mpreal::get_default_prec()); | |
// Returns smallest eps such that x + eps != x (relative machine epsilon) | |
inline mpreal machine_epsilon(const mpreal& x); | |
// Gives max & min values for the required precision, | |
// minval is 'safe' meaning 1 / minval does not overflow | |
// maxval is 'safe' meaning 1 / maxval does not underflow | |
inline mpreal minval(mp_prec_t prec = mpreal::get_default_prec()); | |
inline mpreal maxval(mp_prec_t prec = mpreal::get_default_prec()); | |
// 'Dirty' equality check 1: |a-b| < min{|a|,|b|} * eps | |
inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps); | |
// 'Dirty' equality check 2: |a-b| < min{|a|,|b|} * eps( min{|a|,|b|} ) | |
inline bool isEqualFuzzy(const mpreal& a, const mpreal& b); | |
// 'Bitwise' equality check | |
// maxUlps - a and b can be apart by maxUlps binary numbers. | |
inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps); | |
////////////////////////////////////////////////////////////////////////// | |
// Convert precision in 'bits' to decimal digits and vice versa. | |
// bits = ceil(digits*log[2](10)) | |
// digits = floor(bits*log[10](2)) | |
inline mp_prec_t digits2bits(int d); | |
inline int bits2digits(mp_prec_t b); | |
////////////////////////////////////////////////////////////////////////// | |
// min, max | |
const mpreal (max)(const mpreal& x, const mpreal& y); | |
const mpreal (min)(const mpreal& x, const mpreal& y); | |
////////////////////////////////////////////////////////////////////////// | |
// Implementation | |
////////////////////////////////////////////////////////////////////////// | |
////////////////////////////////////////////////////////////////////////// | |
// Operators - Assignment | |
inline mpreal& mpreal::operator=(const mpreal& v) | |
{ | |
if (this != &v) | |
{ | |
mp_prec_t tp = mpfr_get_prec(mp); | |
mp_prec_t vp = mpfr_get_prec(v.mp); | |
if(tp != vp){ | |
mpfr_clear(mp); | |
mpfr_init2(mp, vp); | |
} | |
mpfr_set(mp, v.mp, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const mpf_t v) | |
{ | |
mpfr_set_f(mp, v, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const mpz_t v) | |
{ | |
mpfr_set_z(mp, v, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const mpq_t v) | |
{ | |
mpfr_set_q(mp, v, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const long double v) | |
{ | |
mpfr_set_ld(mp, v, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const double v) | |
{ | |
#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1) | |
if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW)) | |
{ | |
mpfr_set_d(mp,v,mpreal::get_default_rnd()); | |
}else | |
throw conversion_overflow(); | |
#else | |
mpfr_set_d(mp,v,mpreal::get_default_rnd()); | |
#endif | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const unsigned long int v) | |
{ | |
mpfr_set_ui(mp, v, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const unsigned int v) | |
{ | |
mpfr_set_ui(mp, v, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const long int v) | |
{ | |
mpfr_set_si(mp, v, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const int v) | |
{ | |
mpfr_set_si(mp, v, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const char* s) | |
{ | |
// Use other converters for more precise control on base & precision & rounding: | |
// | |
// mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) | |
// mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode) | |
// | |
// Here we assume base = 10 and we use precision of target variable. | |
mpfr_t t; | |
mpfr_init2(t, mpfr_get_prec(mp)); | |
if(0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd())) | |
{ | |
mpfr_set(mp, t, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
mpfr_clear(t); | |
return *this; | |
} | |
inline mpreal& mpreal::operator=(const std::string& s) | |
{ | |
// Use other converters for more precise control on base & precision & rounding: | |
// | |
// mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) | |
// mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode) | |
// | |
// Here we assume base = 10 and we use precision of target variable. | |
mpfr_t t; | |
mpfr_init2(t, mpfr_get_prec(mp)); | |
if(0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd())) | |
{ | |
mpfr_set(mp, t, mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
mpfr_clear(t); | |
return *this; | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// + Addition | |
inline mpreal& mpreal::operator+=(const mpreal& v) | |
{ | |
mpfr_add(mp,mp,v.mp,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator+=(const mpf_t u) | |
{ | |
*this += mpreal(u); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator+=(const mpz_t u) | |
{ | |
mpfr_add_z(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator+=(const mpq_t u) | |
{ | |
mpfr_add_q(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator+= (const long double u) | |
{ | |
*this += mpreal(u); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator+= (const double u) | |
{ | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) | |
mpfr_add_d(mp,mp,u,mpreal::get_default_rnd()); | |
#else | |
*this += mpreal(u); | |
#endif | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator+=(const unsigned long int u) | |
{ | |
mpfr_add_ui(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator+=(const unsigned int u) | |
{ | |
mpfr_add_ui(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator+=(const long int u) | |
{ | |
mpfr_add_si(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator+=(const int u) | |
{ | |
mpfr_add_si(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
#if defined (MPREAL_HAVE_INT64_SUPPORT) | |
inline mpreal& mpreal::operator+=(const int64_t u){ *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } | |
inline mpreal& mpreal::operator+=(const uint64_t u){ *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } | |
inline mpreal& mpreal::operator-=(const int64_t u){ *this -= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } | |
inline mpreal& mpreal::operator-=(const uint64_t u){ *this -= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } | |
inline mpreal& mpreal::operator*=(const int64_t u){ *this *= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } | |
inline mpreal& mpreal::operator*=(const uint64_t u){ *this *= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } | |
inline mpreal& mpreal::operator/=(const int64_t u){ *this /= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } | |
inline mpreal& mpreal::operator/=(const uint64_t u){ *this /= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } | |
#endif | |
inline const mpreal mpreal::operator+()const { return mpreal(*this); } | |
inline const mpreal operator+(const mpreal& a, const mpreal& b) | |
{ | |
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); | |
mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return c; | |
} | |
inline mpreal& mpreal::operator++() | |
{ | |
return *this += 1; | |
} | |
inline const mpreal mpreal::operator++ (int) | |
{ | |
mpreal x(*this); | |
*this += 1; | |
return x; | |
} | |
inline mpreal& mpreal::operator--() | |
{ | |
return *this -= 1; | |
} | |
inline const mpreal mpreal::operator-- (int) | |
{ | |
mpreal x(*this); | |
*this -= 1; | |
return x; | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// - Subtraction | |
inline mpreal& mpreal::operator-=(const mpreal& v) | |
{ | |
mpfr_sub(mp,mp,v.mp,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator-=(const mpz_t v) | |
{ | |
mpfr_sub_z(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator-=(const mpq_t v) | |
{ | |
mpfr_sub_q(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator-=(const long double v) | |
{ | |
*this -= mpreal(v); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator-=(const double v) | |
{ | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) | |
mpfr_sub_d(mp,mp,v,mpreal::get_default_rnd()); | |
#else | |
*this -= mpreal(v); | |
#endif | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator-=(const unsigned long int v) | |
{ | |
mpfr_sub_ui(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator-=(const unsigned int v) | |
{ | |
mpfr_sub_ui(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator-=(const long int v) | |
{ | |
mpfr_sub_si(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator-=(const int v) | |
{ | |
mpfr_sub_si(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline const mpreal mpreal::operator-()const | |
{ | |
mpreal u(*this); | |
mpfr_neg(u.mp,u.mp,mpreal::get_default_rnd()); | |
return u; | |
} | |
inline const mpreal operator-(const mpreal& a, const mpreal& b) | |
{ | |
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); | |
mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return c; | |
} | |
inline const mpreal operator-(const double b, const mpreal& a) | |
{ | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_d_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return x; | |
#else | |
mpreal x(b, mpfr_get_prec(a.mpfr_ptr())); | |
x -= a; | |
return x; | |
#endif | |
} | |
inline const mpreal operator-(const unsigned long int b, const mpreal& a) | |
{ | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return x; | |
} | |
inline const mpreal operator-(const unsigned int b, const mpreal& a) | |
{ | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return x; | |
} | |
inline const mpreal operator-(const long int b, const mpreal& a) | |
{ | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return x; | |
} | |
inline const mpreal operator-(const int b, const mpreal& a) | |
{ | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return x; | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// * Multiplication | |
inline mpreal& mpreal::operator*= (const mpreal& v) | |
{ | |
mpfr_mul(mp,mp,v.mp,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator*=(const mpz_t v) | |
{ | |
mpfr_mul_z(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator*=(const mpq_t v) | |
{ | |
mpfr_mul_q(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator*=(const long double v) | |
{ | |
*this *= mpreal(v); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator*=(const double v) | |
{ | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) | |
mpfr_mul_d(mp,mp,v,mpreal::get_default_rnd()); | |
#else | |
*this *= mpreal(v); | |
#endif | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator*=(const unsigned long int v) | |
{ | |
mpfr_mul_ui(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator*=(const unsigned int v) | |
{ | |
mpfr_mul_ui(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator*=(const long int v) | |
{ | |
mpfr_mul_si(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator*=(const int v) | |
{ | |
mpfr_mul_si(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline const mpreal operator*(const mpreal& a, const mpreal& b) | |
{ | |
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); | |
mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return c; | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// / Division | |
inline mpreal& mpreal::operator/=(const mpreal& v) | |
{ | |
mpfr_div(mp,mp,v.mp,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator/=(const mpz_t v) | |
{ | |
mpfr_div_z(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator/=(const mpq_t v) | |
{ | |
mpfr_div_q(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator/=(const long double v) | |
{ | |
*this /= mpreal(v); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator/=(const double v) | |
{ | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) | |
mpfr_div_d(mp,mp,v,mpreal::get_default_rnd()); | |
#else | |
*this /= mpreal(v); | |
#endif | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator/=(const unsigned long int v) | |
{ | |
mpfr_div_ui(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator/=(const unsigned int v) | |
{ | |
mpfr_div_ui(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator/=(const long int v) | |
{ | |
mpfr_div_si(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator/=(const int v) | |
{ | |
mpfr_div_si(mp,mp,v,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline const mpreal operator/(const mpreal& a, const mpreal& b) | |
{ | |
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); | |
mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return c; | |
} | |
inline const mpreal operator/(const unsigned long int b, const mpreal& a) | |
{ | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return x; | |
} | |
inline const mpreal operator/(const unsigned int b, const mpreal& a) | |
{ | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); | |
return x; | |
} | |
inline const mpreal operator/(const long int b, const mpreal& a) | |
{ | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd()); | |
return x; | |
} | |
inline const mpreal operator/(const int b, const mpreal& a) | |
{ | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd()); | |
return x; | |
} | |
inline const mpreal operator/(const double b, const mpreal& a) | |
{ | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(),mpreal::get_default_rnd()); | |
return x; | |
#else | |
mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); | |
x /= a; | |
return x; | |
#endif | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// Shifts operators - Multiplication/Division by power of 2 | |
inline mpreal& mpreal::operator<<=(const unsigned long int u) | |
{ | |
mpfr_mul_2ui(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator<<=(const unsigned int u) | |
{ | |
mpfr_mul_2ui(mp,mp,static_cast<unsigned long int>(u),mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator<<=(const long int u) | |
{ | |
mpfr_mul_2si(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator<<=(const int u) | |
{ | |
mpfr_mul_2si(mp,mp,static_cast<long int>(u),mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator>>=(const unsigned long int u) | |
{ | |
mpfr_div_2ui(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator>>=(const unsigned int u) | |
{ | |
mpfr_div_2ui(mp,mp,static_cast<unsigned long int>(u),mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator>>=(const long int u) | |
{ | |
mpfr_div_2si(mp,mp,u,mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::operator>>=(const int u) | |
{ | |
mpfr_div_2si(mp,mp,static_cast<long int>(u),mpreal::get_default_rnd()); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline const mpreal operator<<(const mpreal& v, const unsigned long int k) | |
{ | |
return mul_2ui(v,k); | |
} | |
inline const mpreal operator<<(const mpreal& v, const unsigned int k) | |
{ | |
return mul_2ui(v,static_cast<unsigned long int>(k)); | |
} | |
inline const mpreal operator<<(const mpreal& v, const long int k) | |
{ | |
return mul_2si(v,k); | |
} | |
inline const mpreal operator<<(const mpreal& v, const int k) | |
{ | |
return mul_2si(v,static_cast<long int>(k)); | |
} | |
inline const mpreal operator>>(const mpreal& v, const unsigned long int k) | |
{ | |
return div_2ui(v,k); | |
} | |
inline const mpreal operator>>(const mpreal& v, const long int k) | |
{ | |
return div_2si(v,k); | |
} | |
inline const mpreal operator>>(const mpreal& v, const unsigned int k) | |
{ | |
return div_2ui(v,static_cast<unsigned long int>(k)); | |
} | |
inline const mpreal operator>>(const mpreal& v, const int k) | |
{ | |
return div_2si(v,static_cast<long int>(k)); | |
} | |
// mul_2ui | |
inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) | |
{ | |
mpreal x(v); | |
mpfr_mul_2ui(x.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 mpreal& b){ return (mpfr_greaterequal_p(a.mp,b.mp) !=0); } | |
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 mpreal& b){ return (mpfr_lessequal_p(a.mp,b.mp) !=0); } | |
inline bool operator == (const mpreal& a, const mpreal& b){ return (mpfr_equal_p(a.mp,b.mp) !=0); } | |
inline bool operator != (const mpreal& a, const mpreal& b){ return (mpfr_lessgreater_p(a.mp,b.mp) !=0); } | |
inline bool operator == (const mpreal& a, const unsigned long int b ){ return (mpfr_cmp_ui(a.mp,b) == 0); } | |
inline bool operator == (const mpreal& a, const unsigned int b ){ return (mpfr_cmp_ui(a.mp,b) == 0); } | |
inline bool operator == (const mpreal& a, const long int b ){ return (mpfr_cmp_si(a.mp,b) == 0); } | |
inline bool operator == (const mpreal& a, const int b ){ return (mpfr_cmp_si(a.mp,b) == 0); } | |
inline bool operator == (const mpreal& a, const long double b ){ return (mpfr_cmp_ld(a.mp,b) == 0); } | |
inline bool operator == (const mpreal& a, const double b ){ return (mpfr_cmp_d(a.mp,b) == 0); } | |
inline bool isnan (const mpreal& v){ return (mpfr_nan_p(v.mp) != 0); } | |
inline bool isinf (const mpreal& v){ return (mpfr_inf_p(v.mp) != 0); } | |
inline bool isfinite (const mpreal& v){ return (mpfr_number_p(v.mp) != 0); } | |
inline bool iszero (const mpreal& v){ return (mpfr_zero_p(v.mp) != 0); } | |
inline bool isint (const mpreal& v){ return (mpfr_integer_p(v.mp) != 0); } | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) | |
inline bool isregular(const mpreal& v){ return (mpfr_regular_p(v.mp));} | |
#endif | |
////////////////////////////////////////////////////////////////////////// | |
// Type Converters | |
inline long mpreal::toLong (mp_rnd_t mode) const { return mpfr_get_si(mp, mode); } | |
inline unsigned long mpreal::toULong (mp_rnd_t mode) const { return mpfr_get_ui(mp, mode); } | |
inline double mpreal::toDouble (mp_rnd_t mode) const { return mpfr_get_d (mp, mode); } | |
inline long double mpreal::toLDouble(mp_rnd_t mode) const { return mpfr_get_ld(mp, mode); } | |
#if defined (MPREAL_HAVE_INT64_SUPPORT) | |
inline int64_t mpreal::toInt64 (mp_rnd_t mode) const{ return mpfr_get_sj(mp, mode); } | |
inline uint64_t mpreal::toUInt64(mp_rnd_t mode) const{ return mpfr_get_uj(mp, mode); } | |
#endif | |
inline ::mpfr_ptr mpreal::mpfr_ptr() { return mp; } | |
inline ::mpfr_srcptr mpreal::mpfr_ptr() const { return mp; } | |
inline ::mpfr_srcptr mpreal::mpfr_srcptr() const { return mp; } | |
template <class T> | |
inline std::string toString(T t, std::ios_base & (*f)(std::ios_base&)) | |
{ | |
std::ostringstream oss; | |
oss << f << t; | |
return oss.str(); | |
} | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) | |
inline std::string mpreal::toString(const std::string& format) const | |
{ | |
char *s = NULL; | |
std::string out; | |
if( !format.empty() ) | |
{ | |
if(!(mpfr_asprintf(&s,format.c_str(),mp) < 0)) | |
{ | |
out = std::string(s); | |
mpfr_free_str(s); | |
} | |
} | |
return out; | |
} | |
#endif | |
inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const | |
{ | |
(void)b; | |
(void)mode; | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) | |
// Use MPFR native function for output | |
char format[128]; | |
int digits; | |
digits = n > 0 ? n : bits2digits(mpfr_get_prec(mp)); | |
sprintf(format,"%%.%dRNg",digits); // Default format | |
return toString(std::string(format)); | |
#else | |
char *s, *ns = NULL; | |
size_t slen, nslen; | |
mp_exp_t exp; | |
std::string out; | |
if(mpfr_inf_p(mp)) | |
{ | |
if(mpfr_sgn(mp)>0) return "+Inf"; | |
else return "-Inf"; | |
} | |
if(mpfr_zero_p(mp)) return "0"; | |
if(mpfr_nan_p(mp)) return "NaN"; | |
s = mpfr_get_str(NULL,&exp,b,0,mp,mode); | |
ns = mpfr_get_str(NULL,&exp,b,n,mp,mode); | |
if(s!=NULL && ns!=NULL) | |
{ | |
slen = strlen(s); | |
nslen = strlen(ns); | |
if(nslen<=slen) | |
{ | |
mpfr_free_str(s); | |
s = ns; | |
slen = nslen; | |
} | |
else { | |
mpfr_free_str(ns); | |
} | |
// Make human eye-friendly formatting if possible | |
if (exp>0 && static_cast<size_t>(exp)<slen) | |
{ | |
if(s[0]=='-') | |
{ | |
// Remove zeros starting from right end | |
char* ptr = s+slen-1; | |
while (*ptr=='0' && ptr>s+exp) ptr--; | |
if(ptr==s+exp) out = std::string(s,exp+1); | |
else out = std::string(s,exp+1)+'.'+std::string(s+exp+1,ptr-(s+exp+1)+1); | |
//out = string(s,exp+1)+'.'+string(s+exp+1); | |
} | |
else | |
{ | |
// Remove zeros starting from right end | |
char* ptr = s+slen-1; | |
while (*ptr=='0' && ptr>s+exp-1) ptr--; | |
if(ptr==s+exp-1) out = std::string(s,exp); | |
else out = std::string(s,exp)+'.'+std::string(s+exp,ptr-(s+exp)+1); | |
//out = string(s,exp)+'.'+string(s+exp); | |
} | |
}else{ // exp<0 || exp>slen | |
if(s[0]=='-') | |
{ | |
// Remove zeros starting from right end | |
char* ptr = s+slen-1; | |
while (*ptr=='0' && ptr>s+1) ptr--; | |
if(ptr==s+1) out = std::string(s,2); | |
else out = std::string(s,2)+'.'+std::string(s+2,ptr-(s+2)+1); | |
//out = string(s,2)+'.'+string(s+2); | |
} | |
else | |
{ | |
// Remove zeros starting from right end | |
char* ptr = s+slen-1; | |
while (*ptr=='0' && ptr>s) ptr--; | |
if(ptr==s) out = std::string(s,1); | |
else out = std::string(s,1)+'.'+std::string(s+1,ptr-(s+1)+1); | |
//out = string(s,1)+'.'+string(s+1); | |
} | |
// Make final string | |
if(--exp) | |
{ | |
if(exp>0) out += "e+"+mpfr::toString<mp_exp_t>(exp,std::dec); | |
else out += "e"+mpfr::toString<mp_exp_t>(exp,std::dec); | |
} | |
} | |
mpfr_free_str(s); | |
return out; | |
}else{ | |
return "conversion error!"; | |
} | |
#endif | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// I/O | |
inline std::ostream& operator<<(std::ostream& os, const mpreal& v) | |
{ | |
return os<<v.toString(static_cast<int>(os.precision())); | |
} | |
inline std::istream& operator>>(std::istream &is, mpreal& v) | |
{ | |
// ToDo, use cout::hexfloat and other flags to setup base | |
std::string tmp; | |
is >> tmp; | |
mpfr_set_str(v.mp, tmp.c_str(), 10, mpreal::get_default_rnd()); | |
return is; | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// Bits - decimal digits relation | |
// bits = ceil(digits*log[2](10)) | |
// digits = floor(bits*log[10](2)) | |
inline mp_prec_t digits2bits(int d) | |
{ | |
const double LOG2_10 = 3.3219280948873624; | |
return (mp_prec_t) std::ceil( d * LOG2_10 ); | |
} | |
inline int bits2digits(mp_prec_t b) | |
{ | |
const double LOG10_2 = 0.30102999566398119; | |
return (int) std::floor( b * LOG10_2 ); | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// Set/Get number properties | |
inline int sgn(const mpreal& v) | |
{ | |
int r = mpfr_signbit(v.mp); | |
return (r>0?-1:1); | |
} | |
inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode) | |
{ | |
mpfr_setsign(mp,mp,(sign < 0 ? 1 : 0),RoundingMode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline int mpreal::getPrecision() const | |
{ | |
return mpfr_get_prec(mp); | |
} | |
inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode) | |
{ | |
mpfr_prec_round(mp, Precision, RoundingMode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::setInf(int sign) | |
{ | |
mpfr_set_inf(mp,sign); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::setNan() | |
{ | |
mpfr_set_nan(mp); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mpreal& mpreal::setZero(int sign) | |
{ | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) | |
mpfr_set_zero(mp, sign); | |
#else | |
mpfr_set_si(mp, 0, (mpfr_get_default_rounding_mode)()); | |
setSign(sign); | |
#endif | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return *this; | |
} | |
inline mp_prec_t mpreal::get_prec() const | |
{ | |
return mpfr_get_prec(mp); | |
} | |
inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode) | |
{ | |
mpfr_prec_round(mp,prec,rnd_mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
} | |
inline mp_exp_t mpreal::get_exp () | |
{ | |
return mpfr_get_exp(mp); | |
} | |
inline int mpreal::set_exp (mp_exp_t e) | |
{ | |
int x = mpfr_set_exp(mp, e); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return x; | |
} | |
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::get_default_rnd()); | |
return x; | |
} | |
inline mpreal machine_epsilon(mp_prec_t prec) | |
{ | |
/* the smallest eps such that 1 + eps != 1 */ | |
return machine_epsilon(mpreal(1, prec)); | |
} | |
inline mpreal machine_epsilon(const mpreal& x) | |
{ | |
/* the smallest eps such that x + eps != x */ | |
if( x < 0) | |
{ | |
return nextabove(-x)+x; | |
}else{ | |
return nextabove(x)-x; | |
} | |
} | |
// minval is 'safe' meaning 1 / minval does not overflow | |
inline mpreal minval(mp_prec_t prec) | |
{ | |
/* min = 1/2 * 2^emin = 2^(emin - 1) */ | |
return mpreal(1, prec) << mpreal::get_emin()-1; | |
} | |
// maxval is 'safe' meaning 1 / maxval does not underflow | |
inline mpreal maxval(mp_prec_t prec) | |
{ | |
/* max = (1 - eps) * 2^emax, eps is machine epsilon */ | |
return (mpreal(1, prec) - machine_epsilon(prec)) << mpreal::get_emax(); | |
} | |
inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps) | |
{ | |
return abs(a - b) <= machine_epsilon((max)(abs(a), abs(b))) * maxUlps; | |
} | |
inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps) | |
{ | |
return abs(a - b) <= (min)(abs(a), abs(b)) * eps; | |
} | |
inline bool isEqualFuzzy(const mpreal& a, const mpreal& b) | |
{ | |
return isEqualFuzzy(a, b, machine_epsilon((min)(abs(a), abs(b)))); | |
} | |
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::get_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) | |
{ | |
int r = mpfr_subnormalize(mp,t,rnd_mode); | |
MPREAL_MSVC_DEBUGVIEW_CODE; | |
return r; | |
} | |
inline mp_exp_t mpreal::get_emin (void) | |
{ | |
return mpfr_get_emin(); | |
} | |
inline int mpreal::set_emin (mp_exp_t exp) | |
{ | |
return mpfr_set_emin(exp); | |
} | |
inline mp_exp_t mpreal::get_emax (void) | |
{ | |
return mpfr_get_emax(); | |
} | |
inline int mpreal::set_emax (mp_exp_t exp) | |
{ | |
return mpfr_set_emax(exp); | |
} | |
inline mp_exp_t mpreal::get_emin_min (void) | |
{ | |
return mpfr_get_emin_min(); | |
} | |
inline mp_exp_t mpreal::get_emin_max (void) | |
{ | |
return mpfr_get_emin_max(); | |
} | |
inline mp_exp_t mpreal::get_emax_min (void) | |
{ | |
return mpfr_get_emax_min(); | |
} | |
inline mp_exp_t mpreal::get_emax_max (void) | |
{ | |
return mpfr_get_emax_max(); | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// Mathematical Functions | |
////////////////////////////////////////////////////////////////////////// | |
#define MPREAL_UNARY_MATH_FUNCTION_BODY(f) \ | |
mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); \ | |
mpfr_##f(y.mpfr_ptr(), x.mpfr_srcptr(), r); \ | |
return y; | |
inline const mpreal sqr (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sqr ); } | |
inline const mpreal sqrt (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt); } | |
inline const mpreal sqrt(const unsigned long int x, mp_rnd_t r) | |
{ | |
mpreal y; | |
mpfr_sqrt_ui(y.mpfr_ptr(), x, r); | |
return y; | |
} | |
inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode) | |
{ | |
return sqrt(static_cast<unsigned long int>(v),rnd_mode); | |
} | |
inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode) | |
{ | |
if (v>=0) return sqrt(static_cast<unsigned long int>(v),rnd_mode); | |
else return mpreal().setNan(); // NaN | |
} | |
inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) | |
{ | |
if (v>=0) return sqrt(static_cast<unsigned long int>(v),rnd_mode); | |
else return mpreal().setNan(); // NaN | |
} | |
inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r) | |
{ | |
mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); | |
mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r); | |
return y; | |
} | |
inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r) | |
{ | |
mpreal y(0, mpfr_get_prec(a.mpfr_srcptr())); | |
mpfr_dim(y.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), r); | |
return y; | |
} | |
inline int cmpabs(const mpreal& a,const mpreal& b) | |
{ | |
return mpfr_cmpabs(a.mp,b.mp); | |
} | |
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 sqrt (const long double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v),rnd_mode); } | |
inline const mpreal sqrt (const double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v),rnd_mode); } | |
inline const mpreal cbrt (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt ); } | |
inline const mpreal fabs (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); } | |
inline const mpreal abs (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); } | |
inline const mpreal log (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log ); } | |
inline const mpreal log2 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log2 ); } | |
inline const mpreal log10 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log10); } | |
inline const mpreal exp (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp ); } | |
inline const mpreal exp2 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 ); } | |
inline const mpreal exp10 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp10); } | |
inline const mpreal cos (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cos ); } | |
inline const mpreal sin (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sin ); } | |
inline const mpreal tan (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(tan ); } | |
inline const mpreal sec (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sec ); } | |
inline const mpreal csc (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(csc ); } | |
inline const mpreal cot (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cot ); } | |
inline const mpreal acos (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(acos); } | |
inline const mpreal asin (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(asin); } | |
inline const mpreal atan (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(atan); } | |
inline const mpreal acot (const mpreal& v, mp_rnd_t r) { return atan (1/v, r); } | |
inline const mpreal asec (const mpreal& v, mp_rnd_t r) { return acos (1/v, r); } | |
inline const mpreal acsc (const mpreal& v, mp_rnd_t r) { return asin (1/v, r); } | |
inline const mpreal acoth (const mpreal& v, mp_rnd_t r) { return atanh(1/v, r); } | |
inline const mpreal asech (const mpreal& v, mp_rnd_t r) { return acosh(1/v, r); } | |
inline const mpreal acsch (const mpreal& v, mp_rnd_t r) { return asinh(1/v, r); } | |
inline const mpreal cosh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(cosh ); } | |
inline const mpreal sinh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sinh ); } | |
inline const mpreal tanh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(tanh ); } | |
inline const mpreal sech (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(sech ); } | |
inline const mpreal csch (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(csch ); } | |
inline const mpreal coth (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(coth ); } | |
inline const mpreal acosh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(acosh); } | |
inline const mpreal asinh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(asinh); } | |
inline const mpreal atanh (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(atanh); } | |
inline const mpreal log1p (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(log1p ); } | |
inline const mpreal expm1 (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(expm1 ); } | |
inline const mpreal eint (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(eint ); } | |
inline const mpreal gamma (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(gamma ); } | |
inline const mpreal lngamma (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma); } | |
inline const mpreal zeta (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(zeta ); } | |
inline const mpreal erf (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(erf ); } | |
inline const mpreal erfc (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(erfc ); } | |
inline const mpreal besselj0(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(j0 ); } | |
inline const mpreal besselj1(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(j1 ); } | |
inline const mpreal bessely0(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(y0 ); } | |
inline const mpreal bessely1(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(y1 ); } | |
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 hypot (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_hypot(a.mp, x.mp, y.mp, rnd_mode); | |
return a; | |
} | |
inline const mpreal remainder (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_remainder(a.mp, x.mp, y.mp, rnd_mode); | |
return a; | |
} | |
inline const mpreal remquo (long* q, 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_remquo(a.mp,q, x.mp, y.mp, rnd_mode); | |
return a; | |
} | |
inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode) | |
{ | |
mpreal x(0, prec); | |
mpfr_fac_ui(x.mp,v,rnd_mode); | |
return x; | |
} | |
inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode) | |
{ | |
mpreal x(v); | |
int tsignp; | |
if(signp) mpfr_lgamma(x.mp,signp,v.mp,rnd_mode); | |
else mpfr_lgamma(x.mp,&tsignp,v.mp,rnd_mode); | |
return x; | |
} | |
inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r) | |
{ | |
mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); | |
mpfr_jn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r); | |
return y; | |
} | |
inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r) | |
{ | |
mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); | |
mpfr_yn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r); | |
return y; | |
} | |
inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode) | |
{ | |
mpreal a; | |
mp_prec_t p1, p2, p3; | |
p1 = v1.get_prec(); | |
p2 = v2.get_prec(); | |
p3 = v3.get_prec(); | |
a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); | |
mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); | |
return a; | |
} | |
inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode) | |
{ | |
mpreal a; | |
mp_prec_t p1, p2, p3; | |
p1 = v1.get_prec(); | |
p2 = v2.get_prec(); | |
p3 = v3.get_prec(); | |
a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); | |
mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); | |
return a; | |
} | |
inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode) | |
{ | |
mpreal a; | |
mp_prec_t p1, p2; | |
p1 = v1.get_prec(); | |
p2 = v2.get_prec(); | |
a.set_prec(p1>p2?p1:p2); | |
mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode); | |
return a; | |
} | |
inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode) | |
{ | |
mpreal x; | |
mpfr_ptr* t; | |
unsigned long int i; | |
t = new mpfr_ptr[n]; | |
for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp; | |
mpfr_sum(x.mp,t,n,rnd_mode); | |
delete[] t; | |
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& x, mp_rnd_t r) | |
{ | |
MPREAL_UNARY_MATH_FUNCTION_BODY(li2); | |
} | |
inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) | |
{ | |
/* R = rem(X,Y) if Y != 0, returns X - n * Y where n = trunc(X/Y). */ | |
return fmod(x, y, rnd_mode); | |
} | |
inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) | |
{ | |
(void)rnd_mode; | |
/* | |
m = mod(x,y) if y != 0, returns x - n*y where n = floor(x/y) | |
The following are true by convention: | |
- mod(x,0) is x | |
- mod(x,x) is 0 | |
- mod(x,y) for x != y and y != 0 has the same sign as y. | |
*/ | |
if(iszero(y)) return x; | |
if(x == y) return 0; | |
mpreal m = x - floor(x / y) * y; | |
m.setSign(sgn(y)); // make sure result has the same sign as Y | |
return m; | |
} | |
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& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(digamma); } | |
inline const mpreal ai (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(ai); } | |
#endif // MPFR 3.0.0 Specifics | |
////////////////////////////////////////////////////////////////////////// | |
// Constants | |
inline const mpreal const_log2 (mp_prec_t p, mp_rnd_t r) | |
{ | |
mpreal x(0, p); | |
mpfr_const_log2(x.mpfr_ptr(), r); | |
return x; | |
} | |
inline const mpreal const_pi (mp_prec_t p, mp_rnd_t r) | |
{ | |
mpreal x(0, p); | |
mpfr_const_pi(x.mpfr_ptr(), r); | |
return x; | |
} | |
inline const mpreal const_euler (mp_prec_t p, mp_rnd_t r) | |
{ | |
mpreal x(0, p); | |
mpfr_const_euler(x.mpfr_ptr(), r); | |
return x; | |
} | |
inline const mpreal const_catalan (mp_prec_t p, mp_rnd_t r) | |
{ | |
mpreal x(0, p); | |
mpfr_const_catalan(x.mpfr_ptr(), r); | |
return x; | |
} | |
inline const mpreal const_infinity (int sign, mp_prec_t p, mp_rnd_t /*r*/) | |
{ | |
mpreal x(0, p); | |
mpfr_set_inf(x.mpfr_ptr(), sign); | |
return x; | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// Integer Related Functions | |
inline const mpreal ceil(const mpreal& v) | |
{ | |
mpreal x(v); | |
mpfr_ceil(x.mp,v.mp); | |
return x; | |
} | |
inline const mpreal floor(const mpreal& v) | |
{ | |
mpreal x(v); | |
mpfr_floor(x.mp,v.mp); | |
return x; | |
} | |
inline const mpreal round(const mpreal& v) | |
{ | |
mpreal x(v); | |
mpfr_round(x.mp,v.mp); | |
return x; | |
} | |
inline const mpreal trunc(const mpreal& v) | |
{ | |
mpreal x(v); | |
mpfr_trunc(x.mp,v.mp); | |
return x; | |
} | |
inline const mpreal rint (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint ); } | |
inline const mpreal rint_ceil (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil ); } | |
inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor); } | |
inline const mpreal rint_round (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round); } | |
inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc); } | |
inline const mpreal frac (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(frac ); } | |
////////////////////////////////////////////////////////////////////////// | |
// Miscellaneous Functions | |
inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mp,b.mp); } | |
inline const mpreal (max)(const mpreal& x, const mpreal& y){ return (x>y?x:y); } | |
inline const mpreal (min)(const mpreal& x, const mpreal& y){ return (x<y?x:y); } | |
inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) | |
{ | |
mpreal a; | |
mpfr_max(a.mp,x.mp,y.mp,rnd_mode); | |
return a; | |
} | |
inline const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) | |
{ | |
mpreal a; | |
mpfr_min(a.mp,x.mp,y.mp,rnd_mode); | |
return a; | |
} | |
inline const mpreal nexttoward (const mpreal& x, const mpreal& y) | |
{ | |
mpreal a(x); | |
mpfr_nexttoward(a.mp,y.mp); | |
return a; | |
} | |
inline const mpreal nextabove (const mpreal& x) | |
{ | |
mpreal a(x); | |
mpfr_nextabove(a.mp); | |
return a; | |
} | |
inline const mpreal nextbelow (const mpreal& x) | |
{ | |
mpreal a(x); | |
mpfr_nextbelow(a.mp); | |
return a; | |
} | |
inline const mpreal urandomb (gmp_randstate_t& state) | |
{ | |
mpreal x; | |
mpfr_urandomb(x.mp,state); | |
return x; | |
} | |
#if (MPFR_VERSION >= 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 | |
// Uniformly distributed random number generation | |
// a = random(seed); <- initialization & first random number generation | |
// a = random(); <- next random numbers generation | |
// seed != 0 | |
inline const mpreal random(unsigned int seed) | |
{ | |
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) | |
static gmp_randstate_t state; | |
static bool isFirstTime = true; | |
if(isFirstTime) | |
{ | |
gmp_randinit_default(state); | |
gmp_randseed_ui(state,0); | |
isFirstTime = false; | |
} | |
if(seed != 0) gmp_randseed_ui(state,seed); | |
return mpfr::urandom(state); | |
#else | |
if(seed != 0) std::srand(seed); | |
return mpfr::mpreal(std::rand()/(double)RAND_MAX); | |
#endif | |
} | |
////////////////////////////////////////////////////////////////////////// | |
// Set/Get global properties | |
inline void mpreal::set_default_prec(mp_prec_t prec) | |
{ | |
mpfr_set_default_prec(prec); | |
} | |
inline void mpreal::set_default_rnd(mp_rnd_t rnd_mode) | |
{ | |
mpfr_set_default_rounding_mode(rnd_mode); | |
} | |
inline bool mpreal::fits_in_bits(double x, int n) | |
{ | |
int i; | |
double t; | |
return IsInf(x) || (std::modf ( std::ldexp ( std::frexp ( x, &i ), n ), &t ) == 0.0); | |
} | |
inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode) | |
{ | |
mpreal 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<unsigned long int>(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<long int>(b),rnd_mode); | |
} | |
inline const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode) | |
{ | |
return pow(a,mpreal(b),rnd_mode); | |
} | |
inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode) | |
{ | |
return pow(a,mpreal(b),rnd_mode); | |
} | |
inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode) | |
{ | |
mpreal x(a); | |
mpfr_ui_pow(x.mp,a,b.mp,rnd_mode); | |
return x; | |
} | |
inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode) | |
{ | |
return pow(static_cast<unsigned long int>(a),b,rnd_mode); | |
} | |
inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode) | |
{ | |
if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); | |
else return pow(mpreal(a),b,rnd_mode); | |
} | |
inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode) | |
{ | |
if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); | |
else return pow(mpreal(a),b,rnd_mode); | |
} | |
inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),b,rnd_mode); | |
} | |
inline const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),b,rnd_mode); | |
} | |
// pow unsigned long int | |
inline const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode) | |
{ | |
mpreal x(a); | |
mpfr_ui_pow_ui(x.mp,a,b,rnd_mode); | |
return x; | |
} | |
inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
} | |
inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode) | |
{ | |
if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow | |
} | |
inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode) | |
{ | |
if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow | |
} | |
inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode) | |
{ | |
return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow | |
} | |
inline const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode) | |
{ | |
return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow | |
} | |
// pow unsigned int | |
inline const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui | |
} | |
inline const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
} | |
inline const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode) | |
{ | |
if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
} | |
inline const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode) | |
{ | |
if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
} | |
inline const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode) | |
{ | |
return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
} | |
inline const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode) | |
{ | |
return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
} | |
// pow long int | |
inline const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode) | |
{ | |
if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui | |
else return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui | |
} | |
inline const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode) | |
{ | |
if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui | |
} | |
inline const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode) | |
{ | |
if (a>0) | |
{ | |
if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
}else{ | |
return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si | |
} | |
} | |
inline const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode) | |
{ | |
if (a>0) | |
{ | |
if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
}else{ | |
return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si | |
} | |
} | |
inline const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode) | |
{ | |
if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow | |
} | |
inline const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode) | |
{ | |
if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow | |
} | |
// pow int | |
inline const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode) | |
{ | |
if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui | |
else return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui | |
} | |
inline const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode) | |
{ | |
if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui | |
} | |
inline const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode) | |
{ | |
if (a>0) | |
{ | |
if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
}else{ | |
return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si | |
} | |
} | |
inline const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode) | |
{ | |
if (a>0) | |
{ | |
if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui | |
else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
}else{ | |
return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si | |
} | |
} | |
inline const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode) | |
{ | |
if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow | |
} | |
inline const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode) | |
{ | |
if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow | |
else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow | |
} | |
// pow long double | |
inline const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),mpreal(b),rnd_mode); | |
} | |
inline const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui | |
} | |
inline const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui | |
} | |
inline const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si | |
} | |
inline const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si | |
} | |
inline const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),mpreal(b),rnd_mode); | |
} | |
inline const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),b,rnd_mode); // mpfr_pow_ui | |
} | |
inline const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); // mpfr_pow_ui | |
} | |
inline const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si | |
} | |
inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode) | |
{ | |
return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si | |
} | |
} // End of mpfr namespace | |
// Explicit specialization of std::swap for mpreal numbers | |
// Thus standard algorithms will use efficient version of swap (due to Koenig lookup) | |
// Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap | |
namespace std | |
{ | |
// only allowed to extend namespace std with specializations | |
template <> | |
inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) | |
{ | |
return mpfr::swap(x, y); | |
} | |
template<> | |
class numeric_limits<mpfr::mpreal> | |
{ | |
public: | |
static const bool is_specialized = true; | |
static const bool is_signed = true; | |
static const bool is_integer = false; | |
static const bool is_exact = false; | |
static const int radix = 2; | |
static const bool has_infinity = true; | |
static const bool has_quiet_NaN = true; | |
static const bool has_signaling_NaN = true; | |
static const bool is_iec559 = true; // = IEEE 754 | |
static const bool is_bounded = true; | |
static const bool is_modulo = false; | |
static const bool traps = true; | |
static const bool tinyness_before = true; | |
static const float_denorm_style has_denorm = denorm_absent; | |
inline static float_round_style round_style() | |
{ | |
mp_rnd_t r = mpfr::mpreal::get_default_rnd(); | |
switch (r) | |
{ | |
case MPFR_RNDN: return round_to_nearest; | |
case MPFR_RNDZ: return round_toward_zero; | |
case MPFR_RNDU: return round_toward_infinity; | |
case MPFR_RNDD: return round_toward_neg_infinity; | |
default: return round_indeterminate; | |
} | |
} | |
inline static mpfr::mpreal (min) (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::minval(precision); } | |
inline static mpfr::mpreal (max) (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(precision); } | |
inline static mpfr::mpreal lowest (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(precision); } | |
// Returns smallest eps such that 1 + eps != 1 (classic machine epsilon) | |
inline static mpfr::mpreal epsilon(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(precision); } | |
// Returns smallest eps such that x + eps != x (relative machine epsilon) | |
inline static mpfr::mpreal epsilon(const mpfr::mpreal& x) { return mpfr::machine_epsilon(x); } | |
inline static mpfr::mpreal round_error(mp_prec_t precision = mpfr::mpreal::get_default_prec()) | |
{ | |
mp_rnd_t r = mpfr::mpreal::get_default_rnd(); | |
if(r == MPFR_RNDN) return mpfr::mpreal(0.5, precision); | |
else return mpfr::mpreal(1.0, precision); | |
} | |
inline static const mpfr::mpreal infinity() { return mpfr::const_infinity(); } | |
inline static const mpfr::mpreal quiet_NaN() { return mpfr::mpreal().setNan(); } | |
inline static const mpfr::mpreal signaling_NaN() { return mpfr::mpreal().setNan(); } | |
inline static const mpfr::mpreal denorm_min() { return (min)(); } | |
// Please note, exponent range is not fixed in MPFR | |
static const int min_exponent = MPFR_EMIN_DEFAULT; | |
static const int max_exponent = MPFR_EMAX_DEFAULT; | |
MPREAL_PERMISSIVE_EXPR static const int min_exponent10 = (int) (MPFR_EMIN_DEFAULT * 0.3010299956639811); | |
MPREAL_PERMISSIVE_EXPR static const int max_exponent10 = (int) (MPFR_EMAX_DEFAULT * 0.3010299956639811); | |
// Should be constant according to standard, but 'digits' depends on precision in MPFR | |
inline static int digits() { return mpfr::mpreal::get_default_prec(); } | |
inline static int digits(const mpfr::mpreal& x) { return x.getPrecision(); } | |
inline static int digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec()) | |
{ | |
return mpfr::bits2digits(precision); | |
} | |
inline static int digits10(const mpfr::mpreal& x) | |
{ | |
return mpfr::bits2digits(x.getPrecision()); | |
} | |
inline static int max_digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec()) | |
{ | |
return digits10(precision); | |
} | |
}; | |
} | |
#endif /* __MPREAL_H__ */ |