1/* 2 MPFR C++: Multi-precision floating point number class for C++. 3 Based on MPFR library: http://mpfr.org 4 5 Project homepage: http://www.holoborodko.com/pavel/mpfr 6 Contact e-mail: pavel@holoborodko.com 7 8 Copyright (c) 2008-2014 Pavel Holoborodko 9 10 Contributors: 11 Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman, 12 Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen, 13 Pere Constans, Peter van Hoof, Gael Guennebaud, Tsai Chia Cheng, 14 Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood, 15 Petr Aleksandrov, Orion Poplawski, Charles Karney. 16 17 Licensing: 18 (A) MPFR C++ is under GNU General Public License ("GPL"). 19 20 (B) Non-free licenses may also be purchased from the author, for users who 21 do not want their programs protected by the GPL. 22 23 The non-free licenses are for users that wish to use MPFR C++ in 24 their products but are unwilling to release their software 25 under the GPL (which would require them to release source code 26 and allow free redistribution). 27 28 Such users can purchase an unlimited-use license from the author. 29 Contact us for more details. 30 31 GNU General Public License ("GPL") copyright permissions statement: 32 ************************************************************************** 33 This program is free software: you can redistribute it and/or modify 34 it under the terms of the GNU General Public License as published by 35 the Free Software Foundation, either version 3 of the License, or 36 (at your option) any later version. 37 38 This program is distributed in the hope that it will be useful, 39 but WITHOUT ANY WARRANTY; without even the implied warranty of 40 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 41 GNU General Public License for more details. 42 43 You should have received a copy of the GNU General Public License 44 along with this program. If not, see <http://www.gnu.org/licenses/>. 45*/ 46 47#ifndef __MPREAL_H__ 48#define __MPREAL_H__ 49 50#include <string> 51#include <iostream> 52#include <sstream> 53#include <stdexcept> 54#include <cfloat> 55#include <cmath> 56#include <cstring> 57#include <limits> 58 59// Options 60// FIXME HAVE_INT64_SUPPORT leads to clashes with long int and int64_t on some systems. 61//#define MPREAL_HAVE_INT64_SUPPORT // Enable int64_t support if possible. Available only for MSVC 2010 & GCC. 62#define MPREAL_HAVE_MSVC_DEBUGVIEW // Enable Debugger Visualizer for "Debug" builds in MSVC. 63#define MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS // Enable extended std::numeric_limits<mpfr::mpreal> specialization. 64 // Meaning that "digits", "round_style" and similar members are defined as functions, not constants. 65 // See std::numeric_limits<mpfr::mpreal> at the end of the file for more information. 66 67// Library version 68#define MPREAL_VERSION_MAJOR 3 69#define MPREAL_VERSION_MINOR 5 70#define MPREAL_VERSION_PATCHLEVEL 9 71#define MPREAL_VERSION_STRING "3.5.9" 72 73// Detect compiler using signatures from http://predef.sourceforge.net/ 74#if defined(__GNUC__) && defined(__INTEL_COMPILER) 75 #define IsInf(x) isinf(x) // Intel ICC compiler on Linux 76 77#elif defined(_MSC_VER) // Microsoft Visual C++ 78 #define IsInf(x) (!_finite(x)) 79 80#else 81 #define IsInf(x) std::isinf(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance 82#endif 83 84// A Clang feature extension to determine compiler features. 85#ifndef __has_feature 86 #define __has_feature(x) 0 87#endif 88 89// Detect support for r-value references (move semantic). Borrowed from Eigen. 90#if (__has_feature(cxx_rvalue_references) || \ 91 defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \ 92 (defined(_MSC_VER) && _MSC_VER >= 1600)) 93 94 #define MPREAL_HAVE_MOVE_SUPPORT 95 96 // Use fields in mpfr_t structure to check if it was initialized / set dummy initialization 97 #define mpfr_is_initialized(x) (0 != (x)->_mpfr_d) 98 #define mpfr_set_uninitialized(x) ((x)->_mpfr_d = 0 ) 99#endif 100 101// Detect support for explicit converters. 102#if (__has_feature(cxx_explicit_conversions) || \ 103 defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \ 104 (defined(_MSC_VER) && _MSC_VER >= 1800)) 105 106 #define MPREAL_HAVE_EXPLICIT_CONVERTERS 107#endif 108 109// Detect available 64-bit capabilities 110#if defined(MPREAL_HAVE_INT64_SUPPORT) 111 112 #define MPFR_USE_INTMAX_T // Should be defined before mpfr.h 113 114 #if defined(_MSC_VER) // MSVC + Windows 115 #if (_MSC_VER >= 1600) 116 #include <stdint.h> // <stdint.h> is available only in msvc2010! 117 118 #else // MPFR relies on intmax_t which is available only in msvc2010 119 #undef MPREAL_HAVE_INT64_SUPPORT // Besides, MPFR & MPIR have to be compiled with msvc2010 120 #undef MPFR_USE_INTMAX_T // Since we cannot detect this, disable x64 by default 121 // Someone should change this manually if needed. 122 #endif 123 124 #elif defined (__GNUC__) && defined(__linux__) 125 #if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64) || defined(__ia64) || defined(__itanium__) || defined(_M_IA64) || defined (__PPC64__) 126 #undef MPREAL_HAVE_INT64_SUPPORT // Remove all shaman dances for x64 builds since 127 #undef MPFR_USE_INTMAX_T // GCC already supports x64 as of "long int" is 64-bit integer, nothing left to do 128 #else 129 #include <stdint.h> // use int64_t, uint64_t otherwise 130 #endif 131 132 #else 133 #include <stdint.h> // rely on int64_t, uint64_t in all other cases, Mac OSX, etc. 134 #endif 135 136#endif 137 138#if defined(MPREAL_HAVE_MSVC_DEBUGVIEW) && defined(_MSC_VER) && defined(_DEBUG) 139 #define MPREAL_MSVC_DEBUGVIEW_CODE DebugView = toString(); 140 #define MPREAL_MSVC_DEBUGVIEW_DATA std::string DebugView; 141#else 142 #define MPREAL_MSVC_DEBUGVIEW_CODE 143 #define MPREAL_MSVC_DEBUGVIEW_DATA 144#endif 145 146#include <mpfr.h> 147 148#if (MPFR_VERSION < MPFR_VERSION_NUM(3,0,0)) 149 #include <cstdlib> // Needed for random() 150#endif 151 152// Less important options 153#define MPREAL_DOUBLE_BITS_OVERFLOW -1 // Triggers overflow exception during conversion to double if mpreal 154 // cannot fit in MPREAL_DOUBLE_BITS_OVERFLOW bits 155 // = -1 disables overflow checks (default) 156#if defined(__GNUC__) 157 #define MPREAL_PERMISSIVE_EXPR __extension__ 158#else 159 #define MPREAL_PERMISSIVE_EXPR 160#endif 161 162namespace mpfr { 163 164class mpreal { 165private: 166 mpfr_t mp; 167 168public: 169 170 // Get default rounding mode & precision 171 inline static mp_rnd_t get_default_rnd() { return (mp_rnd_t)(mpfr_get_default_rounding_mode()); } 172 inline static mp_prec_t get_default_prec() { return mpfr_get_default_prec(); } 173 174 // Constructors && type conversions 175 mpreal(); 176 mpreal(const mpreal& u); 177 mpreal(const mpf_t u); 178 mpreal(const mpz_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 179 mpreal(const mpq_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 180 mpreal(const double u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 181 mpreal(const long double u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 182 mpreal(const unsigned long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 183 mpreal(const unsigned int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 184 mpreal(const long int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 185 mpreal(const int u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 186 187 // Construct mpreal from mpfr_t structure. 188 // shared = true allows to avoid deep copy, so that mpreal and 'u' share the same data & pointers. 189 mpreal(const mpfr_t u, bool shared = false); 190 191#if defined (MPREAL_HAVE_INT64_SUPPORT) 192 mpreal(const uint64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 193 mpreal(const int64_t u, mp_prec_t prec = mpreal::get_default_prec(), mp_rnd_t mode = mpreal::get_default_rnd()); 194#endif 195 196 mpreal(const char* s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, mp_rnd_t mode = mpreal::get_default_rnd()); 197 mpreal(const std::string& s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, mp_rnd_t mode = mpreal::get_default_rnd()); 198 199 ~mpreal(); 200 201#ifdef MPREAL_HAVE_MOVE_SUPPORT 202 mpreal& operator=(mpreal&& v); 203 mpreal(mpreal&& u); 204#endif 205 206 // Operations 207 // = 208 // +, -, *, /, ++, --, <<, >> 209 // *=, +=, -=, /=, 210 // <, >, ==, <=, >= 211 212 // = 213 mpreal& operator=(const mpreal& v); 214 mpreal& operator=(const mpf_t v); 215 mpreal& operator=(const mpz_t v); 216 mpreal& operator=(const mpq_t v); 217 mpreal& operator=(const long double v); 218 mpreal& operator=(const double v); 219 mpreal& operator=(const unsigned long int v); 220 mpreal& operator=(const unsigned int v); 221 mpreal& operator=(const long int v); 222 mpreal& operator=(const int v); 223 mpreal& operator=(const char* s); 224 mpreal& operator=(const std::string& s); 225 226 // + 227 mpreal& operator+=(const mpreal& v); 228 mpreal& operator+=(const mpf_t v); 229 mpreal& operator+=(const mpz_t v); 230 mpreal& operator+=(const mpq_t v); 231 mpreal& operator+=(const long double u); 232 mpreal& operator+=(const double u); 233 mpreal& operator+=(const unsigned long int u); 234 mpreal& operator+=(const unsigned int u); 235 mpreal& operator+=(const long int u); 236 mpreal& operator+=(const int u); 237 238#if defined (MPREAL_HAVE_INT64_SUPPORT) 239 mpreal& operator+=(const int64_t u); 240 mpreal& operator+=(const uint64_t u); 241 mpreal& operator-=(const int64_t u); 242 mpreal& operator-=(const uint64_t u); 243 mpreal& operator*=(const int64_t u); 244 mpreal& operator*=(const uint64_t u); 245 mpreal& operator/=(const int64_t u); 246 mpreal& operator/=(const uint64_t u); 247#endif 248 249 const mpreal operator+() const; 250 mpreal& operator++ (); 251 const mpreal operator++ (int); 252 253 // - 254 mpreal& operator-=(const mpreal& v); 255 mpreal& operator-=(const mpz_t v); 256 mpreal& operator-=(const mpq_t v); 257 mpreal& operator-=(const long double u); 258 mpreal& operator-=(const double u); 259 mpreal& operator-=(const unsigned long int u); 260 mpreal& operator-=(const unsigned int u); 261 mpreal& operator-=(const long int u); 262 mpreal& operator-=(const int u); 263 const mpreal operator-() const; 264 friend const mpreal operator-(const unsigned long int b, const mpreal& a); 265 friend const mpreal operator-(const unsigned int b, const mpreal& a); 266 friend const mpreal operator-(const long int b, const mpreal& a); 267 friend const mpreal operator-(const int b, const mpreal& a); 268 friend const mpreal operator-(const double b, const mpreal& a); 269 mpreal& operator-- (); 270 const mpreal operator-- (int); 271 272 // * 273 mpreal& operator*=(const mpreal& v); 274 mpreal& operator*=(const mpz_t v); 275 mpreal& operator*=(const mpq_t v); 276 mpreal& operator*=(const long double v); 277 mpreal& operator*=(const double v); 278 mpreal& operator*=(const unsigned long int v); 279 mpreal& operator*=(const unsigned int v); 280 mpreal& operator*=(const long int v); 281 mpreal& operator*=(const int v); 282 283 // / 284 mpreal& operator/=(const mpreal& v); 285 mpreal& operator/=(const mpz_t v); 286 mpreal& operator/=(const mpq_t v); 287 mpreal& operator/=(const long double v); 288 mpreal& operator/=(const double v); 289 mpreal& operator/=(const unsigned long int v); 290 mpreal& operator/=(const unsigned int v); 291 mpreal& operator/=(const long int v); 292 mpreal& operator/=(const int v); 293 friend const mpreal operator/(const unsigned long int b, const mpreal& a); 294 friend const mpreal operator/(const unsigned int b, const mpreal& a); 295 friend const mpreal operator/(const long int b, const mpreal& a); 296 friend const mpreal operator/(const int b, const mpreal& a); 297 friend const mpreal operator/(const double b, const mpreal& a); 298 299 //<<= Fast Multiplication by 2^u 300 mpreal& operator<<=(const unsigned long int u); 301 mpreal& operator<<=(const unsigned int u); 302 mpreal& operator<<=(const long int u); 303 mpreal& operator<<=(const int u); 304 305 //>>= Fast Division by 2^u 306 mpreal& operator>>=(const unsigned long int u); 307 mpreal& operator>>=(const unsigned int u); 308 mpreal& operator>>=(const long int u); 309 mpreal& operator>>=(const int u); 310 311 // Boolean Operators 312 friend bool operator > (const mpreal& a, const mpreal& b); 313 friend bool operator >= (const mpreal& a, const mpreal& b); 314 friend bool operator < (const mpreal& a, const mpreal& b); 315 friend bool operator <= (const mpreal& a, const mpreal& b); 316 friend bool operator == (const mpreal& a, const mpreal& b); 317 friend bool operator != (const mpreal& a, const mpreal& b); 318 319 // Optimized specializations for boolean operators 320 friend bool operator == (const mpreal& a, const unsigned long int b); 321 friend bool operator == (const mpreal& a, const unsigned int b); 322 friend bool operator == (const mpreal& a, const long int b); 323 friend bool operator == (const mpreal& a, const int b); 324 friend bool operator == (const mpreal& a, const long double b); 325 friend bool operator == (const mpreal& a, const double b); 326 327 // Type Conversion operators 328 bool toBool (mp_rnd_t mode = GMP_RNDZ) const; 329 long toLong (mp_rnd_t mode = GMP_RNDZ) const; 330 unsigned long toULong (mp_rnd_t mode = GMP_RNDZ) const; 331 float toFloat (mp_rnd_t mode = GMP_RNDN) const; 332 double toDouble (mp_rnd_t mode = GMP_RNDN) const; 333 long double toLDouble (mp_rnd_t mode = GMP_RNDN) const; 334 335#if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS) 336 explicit operator bool () const { return toBool(); } 337 explicit operator int () const { return toLong(); } 338 explicit operator long () const { return toLong(); } 339 explicit operator long long () const { return toLong(); } 340 explicit operator unsigned () const { return toULong(); } 341 explicit operator unsigned long () const { return toULong(); } 342 explicit operator unsigned long long () const { return toULong(); } 343 explicit operator float () const { return toFloat(); } 344 explicit operator double () const { return toDouble(); } 345 explicit operator long double () const { return toLDouble(); } 346#endif 347 348#if defined (MPREAL_HAVE_INT64_SUPPORT) 349 int64_t toInt64 (mp_rnd_t mode = GMP_RNDZ) const; 350 uint64_t toUInt64 (mp_rnd_t mode = GMP_RNDZ) const; 351 352 #if defined (MPREAL_HAVE_EXPLICIT_CONVERTERS) 353 explicit operator int64_t () const { return toInt64(); } 354 explicit operator uint64_t () const { return toUInt64(); } 355 #endif 356#endif 357 358 // Get raw pointers so that mpreal can be directly used in raw mpfr_* functions 359 ::mpfr_ptr mpfr_ptr(); 360 ::mpfr_srcptr mpfr_ptr() const; 361 ::mpfr_srcptr mpfr_srcptr() const; 362 363 // Convert mpreal to string with n significant digits in base b 364 // n = -1 -> convert with the maximum available digits 365 std::string toString(int n = -1, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const; 366 367#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 368 std::string toString(const std::string& format) const; 369#endif 370 371 std::ostream& output(std::ostream& os) const; 372 373 // Math Functions 374 friend const mpreal sqr (const mpreal& v, mp_rnd_t rnd_mode); 375 friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode); 376 friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode); 377 friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode); 378 friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode); 379 friend const mpreal pow (const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode); 380 friend const mpreal pow (const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode); 381 friend const mpreal pow (const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode); 382 friend const mpreal pow (const mpreal& a, const long int b, mp_rnd_t rnd_mode); 383 friend const mpreal pow (const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode); 384 friend const mpreal pow (const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode); 385 friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode); 386 387 friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode); 388 friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode); 389 friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode); 390 friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode); 391 friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode); 392 friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode); 393 friend int cmpabs(const mpreal& a,const mpreal& b); 394 395 friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode); 396 friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode); 397 friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode); 398 friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode); 399 friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode); 400 friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode); 401 friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode); 402 friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode); 403 404 friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode); 405 friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode); 406 friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode); 407 friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode); 408 friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode); 409 friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode); 410 friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode); 411 412 friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode); 413 friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode); 414 friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode); 415 friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode); 416 friend const mpreal acot (const mpreal& v, mp_rnd_t rnd_mode); 417 friend const mpreal asec (const mpreal& v, mp_rnd_t rnd_mode); 418 friend const mpreal acsc (const mpreal& v, mp_rnd_t rnd_mode); 419 420 friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode); 421 friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode); 422 friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode); 423 friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode); 424 friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode); 425 friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode); 426 friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode); 427 friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode); 428 friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode); 429 friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode); 430 friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode); 431 friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode); 432 433 friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); 434 435 friend const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode); 436 friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode); 437 438 friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode); 439 friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode); 440 friend const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode); 441 friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode); 442 friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode); 443 friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode); 444 friend const mpreal besselj0 (const mpreal& v, mp_rnd_t rnd_mode); 445 friend const mpreal besselj1 (const mpreal& v, mp_rnd_t rnd_mode); 446 friend const mpreal besseljn (long n, const mpreal& v, mp_rnd_t rnd_mode); 447 friend const mpreal bessely0 (const mpreal& v, mp_rnd_t rnd_mode); 448 friend const mpreal bessely1 (const mpreal& v, mp_rnd_t rnd_mode); 449 friend const mpreal besselyn (long n, const mpreal& v, mp_rnd_t rnd_mode); 450 friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode); 451 friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode); 452 friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode); 453 friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode); 454 friend int sgn(const mpreal& v); // returns -1 or +1 455 456// MPFR 2.4.0 Specifics 457#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 458 friend int sinh_cosh (mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode); 459 friend const mpreal li2 (const mpreal& v, mp_rnd_t rnd_mode); 460 friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); 461 friend const mpreal rec_sqrt (const mpreal& v, mp_rnd_t rnd_mode); 462 463 // MATLAB's semantic equivalents 464 friend const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Remainder after division 465 friend const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Modulus after division 466#endif 467 468// MPFR 3.0.0 Specifics 469#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) 470 friend const mpreal digamma (const mpreal& v, mp_rnd_t rnd_mode); 471 friend const mpreal ai (const mpreal& v, mp_rnd_t rnd_mode); 472 friend const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear 473 friend const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear 474 friend const mpreal grandom (unsigned int seed); 475#endif 476 477 // Uniformly distributed random number generation in [0,1] using 478 // Mersenne-Twister algorithm by default. 479 // Use parameter to setup seed, e.g.: random((unsigned)time(NULL)) 480 // Check urandom() for more precise control. 481 friend const mpreal random(unsigned int seed); 482 483 // Exponent and mantissa manipulation 484 friend const mpreal frexp(const mpreal& v, mp_exp_t* exp); 485 friend const mpreal ldexp(const mpreal& v, mp_exp_t exp); 486 487 // Splits mpreal value into fractional and integer parts. 488 // Returns fractional part and stores integer part in n. 489 friend const mpreal modf(const mpreal& v, mpreal& n); 490 491 // Constants 492 // don't forget to call mpfr_free_cache() for every thread where you are using const-functions 493 friend const mpreal const_log2 (mp_prec_t prec, mp_rnd_t rnd_mode); 494 friend const mpreal const_pi (mp_prec_t prec, mp_rnd_t rnd_mode); 495 friend const mpreal const_euler (mp_prec_t prec, mp_rnd_t rnd_mode); 496 friend const mpreal const_catalan (mp_prec_t prec, mp_rnd_t rnd_mode); 497 498 // returns +inf iff sign>=0 otherwise -inf 499 friend const mpreal const_infinity(int sign, mp_prec_t prec); 500 501 // Output/ Input 502 friend std::ostream& operator<<(std::ostream& os, const mpreal& v); 503 friend std::istream& operator>>(std::istream& is, mpreal& v); 504 505 // Integer Related Functions 506 friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode); 507 friend const mpreal ceil (const mpreal& v); 508 friend const mpreal floor(const mpreal& v); 509 friend const mpreal round(const mpreal& v); 510 friend const mpreal trunc(const mpreal& v); 511 friend const mpreal rint_ceil (const mpreal& v, mp_rnd_t rnd_mode); 512 friend const mpreal rint_floor (const mpreal& v, mp_rnd_t rnd_mode); 513 friend const mpreal rint_round (const mpreal& v, mp_rnd_t rnd_mode); 514 friend const mpreal rint_trunc (const mpreal& v, mp_rnd_t rnd_mode); 515 friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode); 516 friend const mpreal remainder ( const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); 517 friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); 518 519 // Miscellaneous Functions 520 friend const mpreal nexttoward (const mpreal& x, const mpreal& y); 521 friend const mpreal nextabove (const mpreal& x); 522 friend const mpreal nextbelow (const mpreal& x); 523 524 // use gmp_randinit_default() to init state, gmp_randclear() to clear 525 friend const mpreal urandomb (gmp_randstate_t& state); 526 527// MPFR < 2.4.2 Specifics 528#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) 529 friend const mpreal random2 (mp_size_t size, mp_exp_t exp); 530#endif 531 532 // Instance Checkers 533 friend bool isnan (const mpreal& v); 534 friend bool isinf (const mpreal& v); 535 friend bool isfinite (const mpreal& v); 536 537 friend bool isnum (const mpreal& v); 538 friend bool iszero (const mpreal& v); 539 friend bool isint (const mpreal& v); 540 541#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) 542 friend bool isregular(const mpreal& v); 543#endif 544 545 // Set/Get instance properties 546 inline mp_prec_t get_prec() const; 547 inline void set_prec(mp_prec_t prec, mp_rnd_t rnd_mode = get_default_rnd()); // Change precision with rounding mode 548 549 // Aliases for get_prec(), set_prec() - needed for compatibility with std::complex<mpreal> interface 550 inline mpreal& setPrecision(int Precision, mp_rnd_t RoundingMode = get_default_rnd()); 551 inline int getPrecision() const; 552 553 // Set mpreal to +/- inf, NaN, +/-0 554 mpreal& setInf (int Sign = +1); 555 mpreal& setNan (); 556 mpreal& setZero (int Sign = +1); 557 mpreal& setSign (int Sign, mp_rnd_t RoundingMode = get_default_rnd()); 558 559 //Exponent 560 mp_exp_t get_exp(); 561 int set_exp(mp_exp_t e); 562 int check_range (int t, mp_rnd_t rnd_mode = get_default_rnd()); 563 int subnormalize (int t,mp_rnd_t rnd_mode = get_default_rnd()); 564 565 // Inexact conversion from float 566 inline bool fits_in_bits(double x, int n); 567 568 // Set/Get global properties 569 static void set_default_prec(mp_prec_t prec); 570 static void set_default_rnd(mp_rnd_t rnd_mode); 571 572 static mp_exp_t get_emin (void); 573 static mp_exp_t get_emax (void); 574 static mp_exp_t get_emin_min (void); 575 static mp_exp_t get_emin_max (void); 576 static mp_exp_t get_emax_min (void); 577 static mp_exp_t get_emax_max (void); 578 static int set_emin (mp_exp_t exp); 579 static int set_emax (mp_exp_t exp); 580 581 // Efficient swapping of two mpreal values - needed for std algorithms 582 friend void swap(mpreal& x, mpreal& y); 583 584 friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); 585 friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); 586 587private: 588 // Human friendly Debug Preview in Visual Studio. 589 // Put one of these lines: 590 // 591 // mpfr::mpreal=<DebugView> ; Show value only 592 // mpfr::mpreal=<DebugView>, <mp[0]._mpfr_prec,u>bits ; Show value & precision 593 // 594 // at the beginning of 595 // [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat 596 MPREAL_MSVC_DEBUGVIEW_DATA 597 598 // "Smart" resources deallocation. Checks if instance initialized before deletion. 599 void clear(::mpfr_ptr); 600}; 601 602////////////////////////////////////////////////////////////////////////// 603// Exceptions 604class conversion_overflow : public std::exception { 605public: 606 std::string why() { return "inexact conversion from floating point"; } 607}; 608 609////////////////////////////////////////////////////////////////////////// 610// Constructors & converters 611// Default constructor: creates mp number and initializes it to 0. 612inline mpreal::mpreal() 613{ 614 mpfr_init2 (mpfr_ptr(), mpreal::get_default_prec()); 615 mpfr_set_ui(mpfr_ptr(), 0, mpreal::get_default_rnd()); 616 617 MPREAL_MSVC_DEBUGVIEW_CODE; 618} 619 620inline mpreal::mpreal(const mpreal& u) 621{ 622 mpfr_init2(mpfr_ptr(),mpfr_get_prec(u.mpfr_srcptr())); 623 mpfr_set (mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd()); 624 625 MPREAL_MSVC_DEBUGVIEW_CODE; 626} 627 628#ifdef MPREAL_HAVE_MOVE_SUPPORT 629inline mpreal::mpreal(mpreal&& other) 630{ 631 mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds no pinter to actual data 632 mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); 633 634 MPREAL_MSVC_DEBUGVIEW_CODE; 635} 636 637inline mpreal& mpreal::operator=(mpreal&& other) 638{ 639 mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); 640 641 MPREAL_MSVC_DEBUGVIEW_CODE; 642 return *this; 643} 644#endif 645 646inline mpreal::mpreal(const mpfr_t u, bool shared) 647{ 648 if(shared) 649 { 650 std::memcpy(mpfr_ptr(), u, sizeof(mpfr_t)); 651 } 652 else 653 { 654 mpfr_init2(mpfr_ptr(), mpfr_get_prec(u)); 655 mpfr_set (mpfr_ptr(), u, mpreal::get_default_rnd()); 656 } 657 658 MPREAL_MSVC_DEBUGVIEW_CODE; 659} 660 661inline mpreal::mpreal(const mpf_t u) 662{ 663 mpfr_init2(mpfr_ptr(),(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t) 664 mpfr_set_f(mpfr_ptr(),u,mpreal::get_default_rnd()); 665 666 MPREAL_MSVC_DEBUGVIEW_CODE; 667} 668 669inline mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode) 670{ 671 mpfr_init2(mpfr_ptr(), prec); 672 mpfr_set_z(mpfr_ptr(), u, mode); 673 674 MPREAL_MSVC_DEBUGVIEW_CODE; 675} 676 677inline mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode) 678{ 679 mpfr_init2(mpfr_ptr(), prec); 680 mpfr_set_q(mpfr_ptr(), u, mode); 681 682 MPREAL_MSVC_DEBUGVIEW_CODE; 683} 684 685inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) 686{ 687 mpfr_init2(mpfr_ptr(), prec); 688 689#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1) 690 if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW)) 691 { 692 mpfr_set_d(mpfr_ptr(), u, mode); 693 }else 694 throw conversion_overflow(); 695#else 696 mpfr_set_d(mpfr_ptr(), u, mode); 697#endif 698 699 MPREAL_MSVC_DEBUGVIEW_CODE; 700} 701 702inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) 703{ 704 mpfr_init2 (mpfr_ptr(), prec); 705 mpfr_set_ld(mpfr_ptr(), u, mode); 706 707 MPREAL_MSVC_DEBUGVIEW_CODE; 708} 709 710inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) 711{ 712 mpfr_init2 (mpfr_ptr(), prec); 713 mpfr_set_ui(mpfr_ptr(), u, mode); 714 715 MPREAL_MSVC_DEBUGVIEW_CODE; 716} 717 718inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) 719{ 720 mpfr_init2 (mpfr_ptr(), prec); 721 mpfr_set_ui(mpfr_ptr(), u, mode); 722 723 MPREAL_MSVC_DEBUGVIEW_CODE; 724} 725 726inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) 727{ 728 mpfr_init2 (mpfr_ptr(), prec); 729 mpfr_set_si(mpfr_ptr(), u, mode); 730 731 MPREAL_MSVC_DEBUGVIEW_CODE; 732} 733 734inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) 735{ 736 mpfr_init2 (mpfr_ptr(), prec); 737 mpfr_set_si(mpfr_ptr(), u, mode); 738 739 MPREAL_MSVC_DEBUGVIEW_CODE; 740} 741 742#if defined (MPREAL_HAVE_INT64_SUPPORT) 743inline mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode) 744{ 745 mpfr_init2 (mpfr_ptr(), prec); 746 mpfr_set_uj(mpfr_ptr(), u, mode); 747 748 MPREAL_MSVC_DEBUGVIEW_CODE; 749} 750 751inline mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode) 752{ 753 mpfr_init2 (mpfr_ptr(), prec); 754 mpfr_set_sj(mpfr_ptr(), u, mode); 755 756 MPREAL_MSVC_DEBUGVIEW_CODE; 757} 758#endif 759 760inline mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) 761{ 762 mpfr_init2 (mpfr_ptr(), prec); 763 mpfr_set_str(mpfr_ptr(), s, base, mode); 764 765 MPREAL_MSVC_DEBUGVIEW_CODE; 766} 767 768inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode) 769{ 770 mpfr_init2 (mpfr_ptr(), prec); 771 mpfr_set_str(mpfr_ptr(), s.c_str(), base, mode); 772 773 MPREAL_MSVC_DEBUGVIEW_CODE; 774} 775 776inline void mpreal::clear(::mpfr_ptr x) 777{ 778#ifdef MPREAL_HAVE_MOVE_SUPPORT 779 if(mpfr_is_initialized(x)) 780#endif 781 mpfr_clear(x); 782} 783 784inline mpreal::~mpreal() 785{ 786 clear(mpfr_ptr()); 787} 788 789// internal namespace needed for template magic 790namespace internal{ 791 792 // Use SFINAE to restrict arithmetic operations instantiation only for numeric types 793 // This is needed for smooth integration with libraries based on expression templates, like Eigen. 794 // TODO: Do the same for boolean operators. 795 template <typename ArgumentType> struct result_type {}; 796 797 template <> struct result_type<mpreal> {typedef mpreal type;}; 798 template <> struct result_type<mpz_t> {typedef mpreal type;}; 799 template <> struct result_type<mpq_t> {typedef mpreal type;}; 800 template <> struct result_type<long double> {typedef mpreal type;}; 801 template <> struct result_type<double> {typedef mpreal type;}; 802 template <> struct result_type<unsigned long int> {typedef mpreal type;}; 803 template <> struct result_type<unsigned int> {typedef mpreal type;}; 804 template <> struct result_type<long int> {typedef mpreal type;}; 805 template <> struct result_type<int> {typedef mpreal type;}; 806 807#if defined (MPREAL_HAVE_INT64_SUPPORT) 808 template <> struct result_type<int64_t > {typedef mpreal type;}; 809 template <> struct result_type<uint64_t > {typedef mpreal type;}; 810#endif 811} 812 813// + Addition 814template <typename Rhs> 815inline const typename internal::result_type<Rhs>::type 816 operator+(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) += rhs; } 817 818template <typename Lhs> 819inline const typename internal::result_type<Lhs>::type 820 operator+(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) += lhs; } 821 822// - Subtraction 823template <typename Rhs> 824inline const typename internal::result_type<Rhs>::type 825 operator-(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) -= rhs; } 826 827template <typename Lhs> 828inline const typename internal::result_type<Lhs>::type 829 operator-(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) -= rhs; } 830 831// * Multiplication 832template <typename Rhs> 833inline const typename internal::result_type<Rhs>::type 834 operator*(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) *= rhs; } 835 836template <typename Lhs> 837inline const typename internal::result_type<Lhs>::type 838 operator*(const Lhs& lhs, const mpreal& rhs){ return mpreal(rhs) *= lhs; } 839 840// / Division 841template <typename Rhs> 842inline const typename internal::result_type<Rhs>::type 843 operator/(const mpreal& lhs, const Rhs& rhs){ return mpreal(lhs) /= rhs; } 844 845template <typename Lhs> 846inline const typename internal::result_type<Lhs>::type 847 operator/(const Lhs& lhs, const mpreal& rhs){ return mpreal(lhs) /= rhs; } 848 849////////////////////////////////////////////////////////////////////////// 850// sqrt 851const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 852const mpreal sqrt(const long int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 853const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 854const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 855const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 856 857// abs 858inline const mpreal abs(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()); 859 860////////////////////////////////////////////////////////////////////////// 861// pow 862const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 863const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 864const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 865const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 866 867const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 868const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 869const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 870const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 871const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 872 873const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 874const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 875const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 876const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 877const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 878 879const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 880const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 881const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 882const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 883const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 884const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 885 886const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 887const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 888const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 889const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 890const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 891const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 892 893const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 894const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 895const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 896const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 897const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 898const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 899 900const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 901const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 902const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 903const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 904const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 905 906const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 907const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 908const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 909const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 910const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 911 912inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 913inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 914inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 915inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); 916 917////////////////////////////////////////////////////////////////////////// 918// Estimate machine epsilon for the given precision 919// Returns smallest eps such that 1.0 + eps != 1.0 920inline mpreal machine_epsilon(mp_prec_t prec = mpreal::get_default_prec()); 921 922// Returns smallest eps such that x + eps != x (relative machine epsilon) 923inline mpreal machine_epsilon(const mpreal& x); 924 925// Gives max & min values for the required precision, 926// minval is 'safe' meaning 1 / minval does not overflow 927// maxval is 'safe' meaning 1 / maxval does not underflow 928inline mpreal minval(mp_prec_t prec = mpreal::get_default_prec()); 929inline mpreal maxval(mp_prec_t prec = mpreal::get_default_prec()); 930 931// 'Dirty' equality check 1: |a-b| < min{|a|,|b|} * eps 932inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps); 933 934// 'Dirty' equality check 2: |a-b| < min{|a|,|b|} * eps( min{|a|,|b|} ) 935inline bool isEqualFuzzy(const mpreal& a, const mpreal& b); 936 937// 'Bitwise' equality check 938// maxUlps - a and b can be apart by maxUlps binary numbers. 939inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps); 940 941////////////////////////////////////////////////////////////////////////// 942// Convert precision in 'bits' to decimal digits and vice versa. 943// bits = ceil(digits*log[2](10)) 944// digits = floor(bits*log[10](2)) 945 946inline mp_prec_t digits2bits(int d); 947inline int bits2digits(mp_prec_t b); 948 949////////////////////////////////////////////////////////////////////////// 950// min, max 951const mpreal (max)(const mpreal& x, const mpreal& y); 952const mpreal (min)(const mpreal& x, const mpreal& y); 953 954////////////////////////////////////////////////////////////////////////// 955// Implementation 956////////////////////////////////////////////////////////////////////////// 957 958////////////////////////////////////////////////////////////////////////// 959// Operators - Assignment 960inline mpreal& mpreal::operator=(const mpreal& v) 961{ 962 if (this != &v) 963 { 964 mp_prec_t tp = mpfr_get_prec( mpfr_srcptr()); 965 mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr()); 966 967 if(tp != vp){ 968 clear(mpfr_ptr()); 969 mpfr_init2(mpfr_ptr(), vp); 970 } 971 972 mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd()); 973 974 MPREAL_MSVC_DEBUGVIEW_CODE; 975 } 976 return *this; 977} 978 979inline mpreal& mpreal::operator=(const mpf_t v) 980{ 981 mpfr_set_f(mpfr_ptr(), v, mpreal::get_default_rnd()); 982 983 MPREAL_MSVC_DEBUGVIEW_CODE; 984 return *this; 985} 986 987inline mpreal& mpreal::operator=(const mpz_t v) 988{ 989 mpfr_set_z(mpfr_ptr(), v, mpreal::get_default_rnd()); 990 991 MPREAL_MSVC_DEBUGVIEW_CODE; 992 return *this; 993} 994 995inline mpreal& mpreal::operator=(const mpq_t v) 996{ 997 mpfr_set_q(mpfr_ptr(), v, mpreal::get_default_rnd()); 998 999 MPREAL_MSVC_DEBUGVIEW_CODE; 1000 return *this; 1001} 1002 1003inline mpreal& mpreal::operator=(const long double v) 1004{ 1005 mpfr_set_ld(mpfr_ptr(), v, mpreal::get_default_rnd()); 1006 1007 MPREAL_MSVC_DEBUGVIEW_CODE; 1008 return *this; 1009} 1010 1011inline mpreal& mpreal::operator=(const double v) 1012{ 1013#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1) 1014 if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW)) 1015 { 1016 mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd()); 1017 }else 1018 throw conversion_overflow(); 1019#else 1020 mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd()); 1021#endif 1022 1023 MPREAL_MSVC_DEBUGVIEW_CODE; 1024 return *this; 1025} 1026 1027inline mpreal& mpreal::operator=(const unsigned long int v) 1028{ 1029 mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); 1030 1031 MPREAL_MSVC_DEBUGVIEW_CODE; 1032 return *this; 1033} 1034 1035inline mpreal& mpreal::operator=(const unsigned int v) 1036{ 1037 mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); 1038 1039 MPREAL_MSVC_DEBUGVIEW_CODE; 1040 return *this; 1041} 1042 1043inline mpreal& mpreal::operator=(const long int v) 1044{ 1045 mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); 1046 1047 MPREAL_MSVC_DEBUGVIEW_CODE; 1048 return *this; 1049} 1050 1051inline mpreal& mpreal::operator=(const int v) 1052{ 1053 mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); 1054 1055 MPREAL_MSVC_DEBUGVIEW_CODE; 1056 return *this; 1057} 1058 1059inline mpreal& mpreal::operator=(const char* s) 1060{ 1061 // Use other converters for more precise control on base & precision & rounding: 1062 // 1063 // mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) 1064 // mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode) 1065 // 1066 // Here we assume base = 10 and we use precision of target variable. 1067 1068 mpfr_t t; 1069 1070 mpfr_init2(t, mpfr_get_prec(mpfr_srcptr())); 1071 1072 if(0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd())) 1073 { 1074 mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); 1075 MPREAL_MSVC_DEBUGVIEW_CODE; 1076 } 1077 1078 clear(t); 1079 return *this; 1080} 1081 1082inline mpreal& mpreal::operator=(const std::string& s) 1083{ 1084 // Use other converters for more precise control on base & precision & rounding: 1085 // 1086 // mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) 1087 // mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode) 1088 // 1089 // Here we assume base = 10 and we use precision of target variable. 1090 1091 mpfr_t t; 1092 1093 mpfr_init2(t, mpfr_get_prec(mpfr_srcptr())); 1094 1095 if(0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd())) 1096 { 1097 mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); 1098 MPREAL_MSVC_DEBUGVIEW_CODE; 1099 } 1100 1101 clear(t); 1102 return *this; 1103} 1104 1105 1106////////////////////////////////////////////////////////////////////////// 1107// + Addition 1108inline mpreal& mpreal::operator+=(const mpreal& v) 1109{ 1110 mpfr_add(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd()); 1111 MPREAL_MSVC_DEBUGVIEW_CODE; 1112 return *this; 1113} 1114 1115inline mpreal& mpreal::operator+=(const mpf_t u) 1116{ 1117 *this += mpreal(u); 1118 MPREAL_MSVC_DEBUGVIEW_CODE; 1119 return *this; 1120} 1121 1122inline mpreal& mpreal::operator+=(const mpz_t u) 1123{ 1124 mpfr_add_z(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1125 MPREAL_MSVC_DEBUGVIEW_CODE; 1126 return *this; 1127} 1128 1129inline mpreal& mpreal::operator+=(const mpq_t u) 1130{ 1131 mpfr_add_q(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1132 MPREAL_MSVC_DEBUGVIEW_CODE; 1133 return *this; 1134} 1135 1136inline mpreal& mpreal::operator+= (const long double u) 1137{ 1138 *this += mpreal(u); 1139 MPREAL_MSVC_DEBUGVIEW_CODE; 1140 return *this; 1141} 1142 1143inline mpreal& mpreal::operator+= (const double u) 1144{ 1145#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 1146 mpfr_add_d(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1147#else 1148 *this += mpreal(u); 1149#endif 1150 1151 MPREAL_MSVC_DEBUGVIEW_CODE; 1152 return *this; 1153} 1154 1155inline mpreal& mpreal::operator+=(const unsigned long int u) 1156{ 1157 mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1158 MPREAL_MSVC_DEBUGVIEW_CODE; 1159 return *this; 1160} 1161 1162inline mpreal& mpreal::operator+=(const unsigned int u) 1163{ 1164 mpfr_add_ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1165 MPREAL_MSVC_DEBUGVIEW_CODE; 1166 return *this; 1167} 1168 1169inline mpreal& mpreal::operator+=(const long int u) 1170{ 1171 mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1172 MPREAL_MSVC_DEBUGVIEW_CODE; 1173 return *this; 1174} 1175 1176inline mpreal& mpreal::operator+=(const int u) 1177{ 1178 mpfr_add_si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1179 MPREAL_MSVC_DEBUGVIEW_CODE; 1180 return *this; 1181} 1182 1183#if defined (MPREAL_HAVE_INT64_SUPPORT) 1184inline mpreal& mpreal::operator+=(const int64_t u){ *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } 1185inline mpreal& mpreal::operator+=(const uint64_t u){ *this += mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } 1186inline mpreal& mpreal::operator-=(const int64_t u){ *this -= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } 1187inline mpreal& mpreal::operator-=(const uint64_t u){ *this -= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } 1188inline mpreal& mpreal::operator*=(const int64_t u){ *this *= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } 1189inline mpreal& mpreal::operator*=(const uint64_t u){ *this *= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } 1190inline mpreal& mpreal::operator/=(const int64_t u){ *this /= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } 1191inline mpreal& mpreal::operator/=(const uint64_t u){ *this /= mpreal(u); MPREAL_MSVC_DEBUGVIEW_CODE; return *this; } 1192#endif 1193 1194inline const mpreal mpreal::operator+()const { return mpreal(*this); } 1195 1196inline const mpreal operator+(const mpreal& a, const mpreal& b) 1197{ 1198 mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); 1199 mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); 1200 return c; 1201} 1202 1203inline mpreal& mpreal::operator++() 1204{ 1205 return *this += 1; 1206} 1207 1208inline const mpreal mpreal::operator++ (int) 1209{ 1210 mpreal x(*this); 1211 *this += 1; 1212 return x; 1213} 1214 1215inline mpreal& mpreal::operator--() 1216{ 1217 return *this -= 1; 1218} 1219 1220inline const mpreal mpreal::operator-- (int) 1221{ 1222 mpreal x(*this); 1223 *this -= 1; 1224 return x; 1225} 1226 1227////////////////////////////////////////////////////////////////////////// 1228// - Subtraction 1229inline mpreal& mpreal::operator-=(const mpreal& v) 1230{ 1231 mpfr_sub(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd()); 1232 MPREAL_MSVC_DEBUGVIEW_CODE; 1233 return *this; 1234} 1235 1236inline mpreal& mpreal::operator-=(const mpz_t v) 1237{ 1238 mpfr_sub_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1239 MPREAL_MSVC_DEBUGVIEW_CODE; 1240 return *this; 1241} 1242 1243inline mpreal& mpreal::operator-=(const mpq_t v) 1244{ 1245 mpfr_sub_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1246 MPREAL_MSVC_DEBUGVIEW_CODE; 1247 return *this; 1248} 1249 1250inline mpreal& mpreal::operator-=(const long double v) 1251{ 1252 *this -= mpreal(v); 1253 MPREAL_MSVC_DEBUGVIEW_CODE; 1254 return *this; 1255} 1256 1257inline mpreal& mpreal::operator-=(const double v) 1258{ 1259#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 1260 mpfr_sub_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1261#else 1262 *this -= mpreal(v); 1263#endif 1264 1265 MPREAL_MSVC_DEBUGVIEW_CODE; 1266 return *this; 1267} 1268 1269inline mpreal& mpreal::operator-=(const unsigned long int v) 1270{ 1271 mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1272 MPREAL_MSVC_DEBUGVIEW_CODE; 1273 return *this; 1274} 1275 1276inline mpreal& mpreal::operator-=(const unsigned int v) 1277{ 1278 mpfr_sub_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1279 MPREAL_MSVC_DEBUGVIEW_CODE; 1280 return *this; 1281} 1282 1283inline mpreal& mpreal::operator-=(const long int v) 1284{ 1285 mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1286 MPREAL_MSVC_DEBUGVIEW_CODE; 1287 return *this; 1288} 1289 1290inline mpreal& mpreal::operator-=(const int v) 1291{ 1292 mpfr_sub_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1293 MPREAL_MSVC_DEBUGVIEW_CODE; 1294 return *this; 1295} 1296 1297inline const mpreal mpreal::operator-()const 1298{ 1299 mpreal u(*this); 1300 mpfr_neg(u.mpfr_ptr(),u.mpfr_srcptr(),mpreal::get_default_rnd()); 1301 return u; 1302} 1303 1304inline const mpreal operator-(const mpreal& a, const mpreal& b) 1305{ 1306 mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); 1307 mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); 1308 return c; 1309} 1310 1311inline const mpreal operator-(const double b, const mpreal& a) 1312{ 1313#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 1314 mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); 1315 mpfr_d_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1316 return x; 1317#else 1318 mpreal x(b, mpfr_get_prec(a.mpfr_ptr())); 1319 x -= a; 1320 return x; 1321#endif 1322} 1323 1324inline const mpreal operator-(const unsigned long int b, const mpreal& a) 1325{ 1326 mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); 1327 mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1328 return x; 1329} 1330 1331inline const mpreal operator-(const unsigned int b, const mpreal& a) 1332{ 1333 mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); 1334 mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1335 return x; 1336} 1337 1338inline const mpreal operator-(const long int b, const mpreal& a) 1339{ 1340 mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); 1341 mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1342 return x; 1343} 1344 1345inline const mpreal operator-(const int b, const mpreal& a) 1346{ 1347 mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); 1348 mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1349 return x; 1350} 1351 1352////////////////////////////////////////////////////////////////////////// 1353// * Multiplication 1354inline mpreal& mpreal::operator*= (const mpreal& v) 1355{ 1356 mpfr_mul(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd()); 1357 MPREAL_MSVC_DEBUGVIEW_CODE; 1358 return *this; 1359} 1360 1361inline mpreal& mpreal::operator*=(const mpz_t v) 1362{ 1363 mpfr_mul_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1364 MPREAL_MSVC_DEBUGVIEW_CODE; 1365 return *this; 1366} 1367 1368inline mpreal& mpreal::operator*=(const mpq_t v) 1369{ 1370 mpfr_mul_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1371 MPREAL_MSVC_DEBUGVIEW_CODE; 1372 return *this; 1373} 1374 1375inline mpreal& mpreal::operator*=(const long double v) 1376{ 1377 *this *= mpreal(v); 1378 MPREAL_MSVC_DEBUGVIEW_CODE; 1379 return *this; 1380} 1381 1382inline mpreal& mpreal::operator*=(const double v) 1383{ 1384#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 1385 mpfr_mul_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1386#else 1387 *this *= mpreal(v); 1388#endif 1389 MPREAL_MSVC_DEBUGVIEW_CODE; 1390 return *this; 1391} 1392 1393inline mpreal& mpreal::operator*=(const unsigned long int v) 1394{ 1395 mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1396 MPREAL_MSVC_DEBUGVIEW_CODE; 1397 return *this; 1398} 1399 1400inline mpreal& mpreal::operator*=(const unsigned int v) 1401{ 1402 mpfr_mul_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1403 MPREAL_MSVC_DEBUGVIEW_CODE; 1404 return *this; 1405} 1406 1407inline mpreal& mpreal::operator*=(const long int v) 1408{ 1409 mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1410 MPREAL_MSVC_DEBUGVIEW_CODE; 1411 return *this; 1412} 1413 1414inline mpreal& mpreal::operator*=(const int v) 1415{ 1416 mpfr_mul_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1417 MPREAL_MSVC_DEBUGVIEW_CODE; 1418 return *this; 1419} 1420 1421inline const mpreal operator*(const mpreal& a, const mpreal& b) 1422{ 1423 mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); 1424 mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); 1425 return c; 1426} 1427 1428////////////////////////////////////////////////////////////////////////// 1429// / Division 1430inline mpreal& mpreal::operator/=(const mpreal& v) 1431{ 1432 mpfr_div(mpfr_ptr(),mpfr_srcptr(),v.mpfr_srcptr(),mpreal::get_default_rnd()); 1433 MPREAL_MSVC_DEBUGVIEW_CODE; 1434 return *this; 1435} 1436 1437inline mpreal& mpreal::operator/=(const mpz_t v) 1438{ 1439 mpfr_div_z(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1440 MPREAL_MSVC_DEBUGVIEW_CODE; 1441 return *this; 1442} 1443 1444inline mpreal& mpreal::operator/=(const mpq_t v) 1445{ 1446 mpfr_div_q(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1447 MPREAL_MSVC_DEBUGVIEW_CODE; 1448 return *this; 1449} 1450 1451inline mpreal& mpreal::operator/=(const long double v) 1452{ 1453 *this /= mpreal(v); 1454 MPREAL_MSVC_DEBUGVIEW_CODE; 1455 return *this; 1456} 1457 1458inline mpreal& mpreal::operator/=(const double v) 1459{ 1460#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 1461 mpfr_div_d(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1462#else 1463 *this /= mpreal(v); 1464#endif 1465 MPREAL_MSVC_DEBUGVIEW_CODE; 1466 return *this; 1467} 1468 1469inline mpreal& mpreal::operator/=(const unsigned long int v) 1470{ 1471 mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1472 MPREAL_MSVC_DEBUGVIEW_CODE; 1473 return *this; 1474} 1475 1476inline mpreal& mpreal::operator/=(const unsigned int v) 1477{ 1478 mpfr_div_ui(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1479 MPREAL_MSVC_DEBUGVIEW_CODE; 1480 return *this; 1481} 1482 1483inline mpreal& mpreal::operator/=(const long int v) 1484{ 1485 mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1486 MPREAL_MSVC_DEBUGVIEW_CODE; 1487 return *this; 1488} 1489 1490inline mpreal& mpreal::operator/=(const int v) 1491{ 1492 mpfr_div_si(mpfr_ptr(),mpfr_srcptr(),v,mpreal::get_default_rnd()); 1493 MPREAL_MSVC_DEBUGVIEW_CODE; 1494 return *this; 1495} 1496 1497inline const mpreal operator/(const mpreal& a, const mpreal& b) 1498{ 1499 mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_srcptr()), mpfr_get_prec(b.mpfr_srcptr()))); 1500 mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); 1501 return c; 1502} 1503 1504inline const mpreal operator/(const unsigned long int b, const mpreal& a) 1505{ 1506 mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); 1507 mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1508 return x; 1509} 1510 1511inline const mpreal operator/(const unsigned int b, const mpreal& a) 1512{ 1513 mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); 1514 mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1515 return x; 1516} 1517 1518inline const mpreal operator/(const long int b, const mpreal& a) 1519{ 1520 mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); 1521 mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1522 return x; 1523} 1524 1525inline const mpreal operator/(const int b, const mpreal& a) 1526{ 1527 mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); 1528 mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1529 return x; 1530} 1531 1532inline const mpreal operator/(const double b, const mpreal& a) 1533{ 1534#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 1535 mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); 1536 mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); 1537 return x; 1538#else 1539 mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); 1540 x /= a; 1541 return x; 1542#endif 1543} 1544 1545////////////////////////////////////////////////////////////////////////// 1546// Shifts operators - Multiplication/Division by power of 2 1547inline mpreal& mpreal::operator<<=(const unsigned long int u) 1548{ 1549 mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1550 MPREAL_MSVC_DEBUGVIEW_CODE; 1551 return *this; 1552} 1553 1554inline mpreal& mpreal::operator<<=(const unsigned int u) 1555{ 1556 mpfr_mul_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd()); 1557 MPREAL_MSVC_DEBUGVIEW_CODE; 1558 return *this; 1559} 1560 1561inline mpreal& mpreal::operator<<=(const long int u) 1562{ 1563 mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1564 MPREAL_MSVC_DEBUGVIEW_CODE; 1565 return *this; 1566} 1567 1568inline mpreal& mpreal::operator<<=(const int u) 1569{ 1570 mpfr_mul_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd()); 1571 MPREAL_MSVC_DEBUGVIEW_CODE; 1572 return *this; 1573} 1574 1575inline mpreal& mpreal::operator>>=(const unsigned long int u) 1576{ 1577 mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1578 MPREAL_MSVC_DEBUGVIEW_CODE; 1579 return *this; 1580} 1581 1582inline mpreal& mpreal::operator>>=(const unsigned int u) 1583{ 1584 mpfr_div_2ui(mpfr_ptr(),mpfr_srcptr(),static_cast<unsigned long int>(u),mpreal::get_default_rnd()); 1585 MPREAL_MSVC_DEBUGVIEW_CODE; 1586 return *this; 1587} 1588 1589inline mpreal& mpreal::operator>>=(const long int u) 1590{ 1591 mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),u,mpreal::get_default_rnd()); 1592 MPREAL_MSVC_DEBUGVIEW_CODE; 1593 return *this; 1594} 1595 1596inline mpreal& mpreal::operator>>=(const int u) 1597{ 1598 mpfr_div_2si(mpfr_ptr(),mpfr_srcptr(),static_cast<long int>(u),mpreal::get_default_rnd()); 1599 MPREAL_MSVC_DEBUGVIEW_CODE; 1600 return *this; 1601} 1602 1603inline const mpreal operator<<(const mpreal& v, const unsigned long int k) 1604{ 1605 return mul_2ui(v,k); 1606} 1607 1608inline const mpreal operator<<(const mpreal& v, const unsigned int k) 1609{ 1610 return mul_2ui(v,static_cast<unsigned long int>(k)); 1611} 1612 1613inline const mpreal operator<<(const mpreal& v, const long int k) 1614{ 1615 return mul_2si(v,k); 1616} 1617 1618inline const mpreal operator<<(const mpreal& v, const int k) 1619{ 1620 return mul_2si(v,static_cast<long int>(k)); 1621} 1622 1623inline const mpreal operator>>(const mpreal& v, const unsigned long int k) 1624{ 1625 return div_2ui(v,k); 1626} 1627 1628inline const mpreal operator>>(const mpreal& v, const long int k) 1629{ 1630 return div_2si(v,k); 1631} 1632 1633inline const mpreal operator>>(const mpreal& v, const unsigned int k) 1634{ 1635 return div_2ui(v,static_cast<unsigned long int>(k)); 1636} 1637 1638inline const mpreal operator>>(const mpreal& v, const int k) 1639{ 1640 return div_2si(v,static_cast<long int>(k)); 1641} 1642 1643// mul_2ui 1644inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) 1645{ 1646 mpreal x(v); 1647 mpfr_mul_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode); 1648 return x; 1649} 1650 1651// mul_2si 1652inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode) 1653{ 1654 mpreal x(v); 1655 mpfr_mul_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode); 1656 return x; 1657} 1658 1659inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) 1660{ 1661 mpreal x(v); 1662 mpfr_div_2ui(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode); 1663 return x; 1664} 1665 1666inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode) 1667{ 1668 mpreal x(v); 1669 mpfr_div_2si(x.mpfr_ptr(),v.mpfr_srcptr(),k,rnd_mode); 1670 return x; 1671} 1672 1673////////////////////////////////////////////////////////////////////////// 1674//Boolean operators 1675inline bool operator > (const mpreal& a, const mpreal& b){ return (mpfr_greater_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); } 1676inline bool operator >= (const mpreal& a, const mpreal& b){ return (mpfr_greaterequal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); } 1677inline bool operator < (const mpreal& a, const mpreal& b){ return (mpfr_less_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); } 1678inline bool operator <= (const mpreal& a, const mpreal& b){ return (mpfr_lessequal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); } 1679inline bool operator == (const mpreal& a, const mpreal& b){ return (mpfr_equal_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); } 1680inline bool operator != (const mpreal& a, const mpreal& b){ return (mpfr_lessgreater_p (a.mpfr_srcptr(),b.mpfr_srcptr()) !=0 ); } 1681 1682inline bool operator == (const mpreal& a, const unsigned long int b ){ return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); } 1683inline bool operator == (const mpreal& a, const unsigned int b ){ return (mpfr_cmp_ui(a.mpfr_srcptr(),b) == 0 ); } 1684inline bool operator == (const mpreal& a, const long int b ){ return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); } 1685inline bool operator == (const mpreal& a, const int b ){ return (mpfr_cmp_si(a.mpfr_srcptr(),b) == 0 ); } 1686inline bool operator == (const mpreal& a, const long double b ){ return (mpfr_cmp_ld(a.mpfr_srcptr(),b) == 0 ); } 1687inline bool operator == (const mpreal& a, const double b ){ return (mpfr_cmp_d (a.mpfr_srcptr(),b) == 0 ); } 1688 1689 1690inline bool isnan (const mpreal& op){ return (mpfr_nan_p (op.mpfr_srcptr()) != 0 ); } 1691inline bool isinf (const mpreal& op){ return (mpfr_inf_p (op.mpfr_srcptr()) != 0 ); } 1692inline bool isfinite (const mpreal& op){ return (mpfr_number_p (op.mpfr_srcptr()) != 0 ); } 1693inline bool iszero (const mpreal& op){ return (mpfr_zero_p (op.mpfr_srcptr()) != 0 ); } 1694inline bool isint (const mpreal& op){ return (mpfr_integer_p(op.mpfr_srcptr()) != 0 ); } 1695 1696#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) 1697inline bool isregular(const mpreal& op){ return (mpfr_regular_p(op.mpfr_srcptr()));} 1698#endif 1699 1700////////////////////////////////////////////////////////////////////////// 1701// Type Converters 1702inline bool mpreal::toBool (mp_rnd_t /*mode*/) const { return mpfr_zero_p (mpfr_srcptr()) == 0; } 1703inline long mpreal::toLong (mp_rnd_t mode) const { return mpfr_get_si (mpfr_srcptr(), mode); } 1704inline unsigned long mpreal::toULong (mp_rnd_t mode) const { return mpfr_get_ui (mpfr_srcptr(), mode); } 1705inline float mpreal::toFloat (mp_rnd_t mode) const { return mpfr_get_flt(mpfr_srcptr(), mode); } 1706inline double mpreal::toDouble (mp_rnd_t mode) const { return mpfr_get_d (mpfr_srcptr(), mode); } 1707inline long double mpreal::toLDouble(mp_rnd_t mode) const { return mpfr_get_ld (mpfr_srcptr(), mode); } 1708 1709#if defined (MPREAL_HAVE_INT64_SUPPORT) 1710inline int64_t mpreal::toInt64 (mp_rnd_t mode) const{ return mpfr_get_sj(mpfr_srcptr(), mode); } 1711inline uint64_t mpreal::toUInt64(mp_rnd_t mode) const{ return mpfr_get_uj(mpfr_srcptr(), mode); } 1712#endif 1713 1714inline ::mpfr_ptr mpreal::mpfr_ptr() { return mp; } 1715inline ::mpfr_srcptr mpreal::mpfr_ptr() const { return mp; } 1716inline ::mpfr_srcptr mpreal::mpfr_srcptr() const { return mp; } 1717 1718template <class T> 1719inline std::string toString(T t, std::ios_base & (*f)(std::ios_base&)) 1720{ 1721 std::ostringstream oss; 1722 oss << f << t; 1723 return oss.str(); 1724} 1725 1726#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 1727 1728inline std::string mpreal::toString(const std::string& format) const 1729{ 1730 char *s = NULL; 1731 std::string out; 1732 1733 if( !format.empty() ) 1734 { 1735 if(!(mpfr_asprintf(&s, format.c_str(), mpfr_srcptr()) < 0)) 1736 { 1737 out = std::string(s); 1738 1739 mpfr_free_str(s); 1740 } 1741 } 1742 1743 return out; 1744} 1745 1746#endif 1747 1748inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const 1749{ 1750 // TODO: Add extended format specification (f, e, rounding mode) as it done in output operator 1751 (void)b; 1752 (void)mode; 1753 1754#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 1755 1756 std::ostringstream format; 1757 1758 int digits = (n >= 0) ? n : bits2digits(mpfr_get_prec(mpfr_srcptr())); 1759 1760 format << "%." << digits << "RNg"; 1761 1762 return toString(format.str()); 1763 1764#else 1765 1766 char *s, *ns = NULL; 1767 size_t slen, nslen; 1768 mp_exp_t exp; 1769 std::string out; 1770 1771 if(mpfr_inf_p(mp)) 1772 { 1773 if(mpfr_sgn(mp)>0) return "+Inf"; 1774 else return "-Inf"; 1775 } 1776 1777 if(mpfr_zero_p(mp)) return "0"; 1778 if(mpfr_nan_p(mp)) return "NaN"; 1779 1780 s = mpfr_get_str(NULL, &exp, b, 0, mp, mode); 1781 ns = mpfr_get_str(NULL, &exp, b, (std::max)(0,n), mp, mode); 1782 1783 if(s!=NULL && ns!=NULL) 1784 { 1785 slen = strlen(s); 1786 nslen = strlen(ns); 1787 if(nslen<=slen) 1788 { 1789 mpfr_free_str(s); 1790 s = ns; 1791 slen = nslen; 1792 } 1793 else { 1794 mpfr_free_str(ns); 1795 } 1796 1797 // Make human eye-friendly formatting if possible 1798 if (exp>0 && static_cast<size_t>(exp)<slen) 1799 { 1800 if(s[0]=='-') 1801 { 1802 // Remove zeros starting from right end 1803 char* ptr = s+slen-1; 1804 while (*ptr=='0' && ptr>s+exp) ptr--; 1805 1806 if(ptr==s+exp) out = std::string(s,exp+1); 1807 else out = std::string(s,exp+1)+'.'+std::string(s+exp+1,ptr-(s+exp+1)+1); 1808 1809 //out = string(s,exp+1)+'.'+string(s+exp+1); 1810 } 1811 else 1812 { 1813 // Remove zeros starting from right end 1814 char* ptr = s+slen-1; 1815 while (*ptr=='0' && ptr>s+exp-1) ptr--; 1816 1817 if(ptr==s+exp-1) out = std::string(s,exp); 1818 else out = std::string(s,exp)+'.'+std::string(s+exp,ptr-(s+exp)+1); 1819 1820 //out = string(s,exp)+'.'+string(s+exp); 1821 } 1822 1823 }else{ // exp<0 || exp>slen 1824 if(s[0]=='-') 1825 { 1826 // Remove zeros starting from right end 1827 char* ptr = s+slen-1; 1828 while (*ptr=='0' && ptr>s+1) ptr--; 1829 1830 if(ptr==s+1) out = std::string(s,2); 1831 else out = std::string(s,2)+'.'+std::string(s+2,ptr-(s+2)+1); 1832 1833 //out = string(s,2)+'.'+string(s+2); 1834 } 1835 else 1836 { 1837 // Remove zeros starting from right end 1838 char* ptr = s+slen-1; 1839 while (*ptr=='0' && ptr>s) ptr--; 1840 1841 if(ptr==s) out = std::string(s,1); 1842 else out = std::string(s,1)+'.'+std::string(s+1,ptr-(s+1)+1); 1843 1844 //out = string(s,1)+'.'+string(s+1); 1845 } 1846 1847 // Make final string 1848 if(--exp) 1849 { 1850 if(exp>0) out += "e+"+mpfr::toString<mp_exp_t>(exp,std::dec); 1851 else out += "e"+mpfr::toString<mp_exp_t>(exp,std::dec); 1852 } 1853 } 1854 1855 mpfr_free_str(s); 1856 return out; 1857 }else{ 1858 return "conversion error!"; 1859 } 1860#endif 1861} 1862 1863 1864////////////////////////////////////////////////////////////////////////// 1865// I/O 1866inline std::ostream& mpreal::output(std::ostream& os) const 1867{ 1868 std::ostringstream format; 1869 const std::ios::fmtflags flags = os.flags(); 1870 1871 format << ((flags & std::ios::showpos) ? "%+" : "%"); 1872 if (os.precision() >= 0) 1873 format << '.' << os.precision() << "R*" 1874 << ((flags & std::ios::floatfield) == std::ios::fixed ? 'f' : 1875 (flags & std::ios::floatfield) == std::ios::scientific ? 'e' : 1876 'g'); 1877 else 1878 format << "R*e"; 1879 1880 char *s = NULL; 1881 if(!(mpfr_asprintf(&s, format.str().c_str(), 1882 mpfr::mpreal::get_default_rnd(), 1883 mpfr_srcptr()) 1884 < 0)) 1885 { 1886 os << std::string(s); 1887 mpfr_free_str(s); 1888 } 1889 return os; 1890} 1891 1892inline std::ostream& operator<<(std::ostream& os, const mpreal& v) 1893{ 1894 return v.output(os); 1895} 1896 1897inline std::istream& operator>>(std::istream &is, mpreal& v) 1898{ 1899 // TODO: use cout::hexfloat and other flags to setup base 1900 std::string tmp; 1901 is >> tmp; 1902 mpfr_set_str(v.mpfr_ptr(), tmp.c_str(), 10, mpreal::get_default_rnd()); 1903 return is; 1904} 1905 1906////////////////////////////////////////////////////////////////////////// 1907// Bits - decimal digits relation 1908// bits = ceil(digits*log[2](10)) 1909// digits = floor(bits*log[10](2)) 1910 1911inline mp_prec_t digits2bits(int d) 1912{ 1913 const double LOG2_10 = 3.3219280948873624; 1914 1915 return mp_prec_t(std::ceil( d * LOG2_10 )); 1916} 1917 1918inline int bits2digits(mp_prec_t b) 1919{ 1920 const double LOG10_2 = 0.30102999566398119; 1921 1922 return int(std::floor( b * LOG10_2 )); 1923} 1924 1925////////////////////////////////////////////////////////////////////////// 1926// Set/Get number properties 1927inline int sgn(const mpreal& op) 1928{ 1929 int r = mpfr_signbit(op.mpfr_srcptr()); 1930 return (r > 0? -1 : 1); 1931} 1932 1933inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode) 1934{ 1935 mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), (sign < 0 ? 1 : 0), RoundingMode); 1936 MPREAL_MSVC_DEBUGVIEW_CODE; 1937 return *this; 1938} 1939 1940inline int mpreal::getPrecision() const 1941{ 1942 return int(mpfr_get_prec(mpfr_srcptr())); 1943} 1944 1945inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode) 1946{ 1947 mpfr_prec_round(mpfr_ptr(), Precision, RoundingMode); 1948 MPREAL_MSVC_DEBUGVIEW_CODE; 1949 return *this; 1950} 1951 1952inline mpreal& mpreal::setInf(int sign) 1953{ 1954 mpfr_set_inf(mpfr_ptr(), sign); 1955 MPREAL_MSVC_DEBUGVIEW_CODE; 1956 return *this; 1957} 1958 1959inline mpreal& mpreal::setNan() 1960{ 1961 mpfr_set_nan(mpfr_ptr()); 1962 MPREAL_MSVC_DEBUGVIEW_CODE; 1963 return *this; 1964} 1965 1966inline mpreal& mpreal::setZero(int sign) 1967{ 1968 1969#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) 1970 mpfr_set_zero(mpfr_ptr(), sign); 1971#else 1972 mpfr_set_si(mpfr_ptr(), 0, (mpfr_get_default_rounding_mode)()); 1973 setSign(sign); 1974#endif 1975 1976 MPREAL_MSVC_DEBUGVIEW_CODE; 1977 return *this; 1978} 1979 1980inline mp_prec_t mpreal::get_prec() const 1981{ 1982 return mpfr_get_prec(mpfr_srcptr()); 1983} 1984 1985inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode) 1986{ 1987 mpfr_prec_round(mpfr_ptr(),prec,rnd_mode); 1988 MPREAL_MSVC_DEBUGVIEW_CODE; 1989} 1990 1991inline mp_exp_t mpreal::get_exp () 1992{ 1993 return mpfr_get_exp(mpfr_srcptr()); 1994} 1995 1996inline int mpreal::set_exp (mp_exp_t e) 1997{ 1998 int x = mpfr_set_exp(mpfr_ptr(), e); 1999 MPREAL_MSVC_DEBUGVIEW_CODE; 2000 return x; 2001} 2002 2003inline const mpreal frexp(const mpreal& v, mp_exp_t* exp) 2004{ 2005 mpreal x(v); 2006 *exp = x.get_exp(); 2007 x.set_exp(0); 2008 return x; 2009} 2010 2011inline const mpreal ldexp(const mpreal& v, mp_exp_t exp) 2012{ 2013 mpreal x(v); 2014 2015 // rounding is not important since we just increasing the exponent 2016 mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd()); 2017 return x; 2018} 2019 2020inline mpreal machine_epsilon(mp_prec_t prec) 2021{ 2022 /* the smallest eps such that 1 + eps != 1 */ 2023 return machine_epsilon(mpreal(1, prec)); 2024} 2025 2026inline mpreal machine_epsilon(const mpreal& x) 2027{ 2028 /* the smallest eps such that x + eps != x */ 2029 if( x < 0) 2030 { 2031 return nextabove(-x) + x; 2032 }else{ 2033 return nextabove( x) - x; 2034 } 2035} 2036 2037// minval is 'safe' meaning 1 / minval does not overflow 2038inline mpreal minval(mp_prec_t prec) 2039{ 2040 /* min = 1/2 * 2^emin = 2^(emin - 1) */ 2041 return mpreal(1, prec) << mpreal::get_emin()-1; 2042} 2043 2044// maxval is 'safe' meaning 1 / maxval does not underflow 2045inline mpreal maxval(mp_prec_t prec) 2046{ 2047 /* max = (1 - eps) * 2^emax, eps is machine epsilon */ 2048 return (mpreal(1, prec) - machine_epsilon(prec)) << mpreal::get_emax(); 2049} 2050 2051inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps) 2052{ 2053 return abs(a - b) <= machine_epsilon((max)(abs(a), abs(b))) * maxUlps; 2054} 2055 2056inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps) 2057{ 2058 return abs(a - b) <= eps; 2059} 2060 2061inline bool isEqualFuzzy(const mpreal& a, const mpreal& b) 2062{ 2063 return isEqualFuzzy(a, b, machine_epsilon((max)(1, (min)(abs(a), abs(b))))); 2064} 2065 2066inline const mpreal modf(const mpreal& v, mpreal& n) 2067{ 2068 mpreal f(v); 2069 2070 // rounding is not important since we are using the same number 2071 mpfr_frac (f.mpfr_ptr(),f.mpfr_srcptr(),mpreal::get_default_rnd()); 2072 mpfr_trunc(n.mpfr_ptr(),v.mpfr_srcptr()); 2073 return f; 2074} 2075 2076inline int mpreal::check_range (int t, mp_rnd_t rnd_mode) 2077{ 2078 return mpfr_check_range(mpfr_ptr(),t,rnd_mode); 2079} 2080 2081inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode) 2082{ 2083 int r = mpfr_subnormalize(mpfr_ptr(),t,rnd_mode); 2084 MPREAL_MSVC_DEBUGVIEW_CODE; 2085 return r; 2086} 2087 2088inline mp_exp_t mpreal::get_emin (void) 2089{ 2090 return mpfr_get_emin(); 2091} 2092 2093inline int mpreal::set_emin (mp_exp_t exp) 2094{ 2095 return mpfr_set_emin(exp); 2096} 2097 2098inline mp_exp_t mpreal::get_emax (void) 2099{ 2100 return mpfr_get_emax(); 2101} 2102 2103inline int mpreal::set_emax (mp_exp_t exp) 2104{ 2105 return mpfr_set_emax(exp); 2106} 2107 2108inline mp_exp_t mpreal::get_emin_min (void) 2109{ 2110 return mpfr_get_emin_min(); 2111} 2112 2113inline mp_exp_t mpreal::get_emin_max (void) 2114{ 2115 return mpfr_get_emin_max(); 2116} 2117 2118inline mp_exp_t mpreal::get_emax_min (void) 2119{ 2120 return mpfr_get_emax_min(); 2121} 2122 2123inline mp_exp_t mpreal::get_emax_max (void) 2124{ 2125 return mpfr_get_emax_max(); 2126} 2127 2128////////////////////////////////////////////////////////////////////////// 2129// Mathematical Functions 2130////////////////////////////////////////////////////////////////////////// 2131#define MPREAL_UNARY_MATH_FUNCTION_BODY(f) \ 2132 mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); \ 2133 mpfr_##f(y.mpfr_ptr(), x.mpfr_srcptr(), r); \ 2134 return y; 2135 2136inline const mpreal sqr (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) 2137{ MPREAL_UNARY_MATH_FUNCTION_BODY(sqr ); } 2138 2139inline const mpreal sqrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) 2140{ MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt); } 2141 2142inline const mpreal sqrt(const unsigned long int x, mp_rnd_t r) 2143{ 2144 mpreal y; 2145 mpfr_sqrt_ui(y.mpfr_ptr(), x, r); 2146 return y; 2147} 2148 2149inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode) 2150{ 2151 return sqrt(static_cast<unsigned long int>(v),rnd_mode); 2152} 2153 2154inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode) 2155{ 2156 if (v>=0) return sqrt(static_cast<unsigned long int>(v),rnd_mode); 2157 else return mpreal().setNan(); // NaN 2158} 2159 2160inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) 2161{ 2162 if (v>=0) return sqrt(static_cast<unsigned long int>(v),rnd_mode); 2163 else return mpreal().setNan(); // NaN 2164} 2165 2166inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r = mpreal::get_default_rnd()) 2167{ 2168 mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); 2169 mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r); 2170 return y; 2171} 2172 2173inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r = mpreal::get_default_rnd()) 2174{ 2175 mpreal y(0, mpfr_get_prec(a.mpfr_srcptr())); 2176 mpfr_dim(y.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), r); 2177 return y; 2178} 2179 2180inline int cmpabs(const mpreal& a,const mpreal& b) 2181{ 2182 return mpfr_cmpabs(a.mpfr_ptr(), b.mpfr_srcptr()); 2183} 2184 2185inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2186{ 2187 return mpfr_sin_cos(s.mpfr_ptr(), c.mpfr_ptr(), v.mpfr_srcptr(), rnd_mode); 2188} 2189 2190inline const mpreal sqrt (const long double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v),rnd_mode); } 2191inline const mpreal sqrt (const double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v),rnd_mode); } 2192 2193inline const mpreal cbrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt ); } 2194inline const mpreal fabs (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); } 2195inline const mpreal abs (const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs ); } 2196inline const mpreal log (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log ); } 2197inline const mpreal log2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log2 ); } 2198inline const mpreal log10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log10); } 2199inline const mpreal exp (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp ); } 2200inline const mpreal exp2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp2 ); } 2201inline const mpreal exp10 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(exp10); } 2202inline const mpreal cos (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cos ); } 2203inline const mpreal sin (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sin ); } 2204inline const mpreal tan (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(tan ); } 2205inline const mpreal sec (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sec ); } 2206inline const mpreal csc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(csc ); } 2207inline const mpreal cot (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cot ); } 2208inline const mpreal acos (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(acos ); } 2209inline const mpreal asin (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(asin ); } 2210inline const mpreal atan (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(atan ); } 2211 2212inline const mpreal acot (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return atan (1/v, r); } 2213inline const mpreal asec (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return acos (1/v, r); } 2214inline const mpreal acsc (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return asin (1/v, r); } 2215inline const mpreal acoth (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return atanh(1/v, r); } 2216inline const mpreal asech (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return acosh(1/v, r); } 2217inline const mpreal acsch (const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return asinh(1/v, r); } 2218 2219inline const mpreal cosh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(cosh ); } 2220inline const mpreal sinh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sinh ); } 2221inline const mpreal tanh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(tanh ); } 2222inline const mpreal sech (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(sech ); } 2223inline const mpreal csch (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(csch ); } 2224inline const mpreal coth (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(coth ); } 2225inline const mpreal acosh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(acosh); } 2226inline const mpreal asinh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(asinh); } 2227inline const mpreal atanh (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(atanh); } 2228 2229inline const mpreal log1p (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(log1p ); } 2230inline const mpreal expm1 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(expm1 ); } 2231inline const mpreal eint (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(eint ); } 2232inline const mpreal gamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(gamma ); } 2233inline const mpreal lngamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma); } 2234inline const mpreal zeta (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(zeta ); } 2235inline const mpreal erf (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(erf ); } 2236inline const mpreal erfc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(erfc ); } 2237inline const mpreal besselj0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(j0 ); } 2238inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(j1 ); } 2239inline const mpreal bessely0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y0 ); } 2240inline const mpreal bessely1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y1 ); } 2241 2242inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2243{ 2244 mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision())); 2245 mpfr_atan2(a.mpfr_ptr(), y.mpfr_srcptr(), x.mpfr_srcptr(), rnd_mode); 2246 return a; 2247} 2248 2249inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2250{ 2251 mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision())); 2252 mpfr_hypot(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode); 2253 return a; 2254} 2255 2256inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2257{ 2258 mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision())); 2259 mpfr_remainder(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode); 2260 return a; 2261} 2262 2263inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2264{ 2265 mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision())); 2266 mpfr_remquo(a.mpfr_ptr(),q, x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode); 2267 return a; 2268} 2269 2270inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec = mpreal::get_default_prec(), 2271 mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2272{ 2273 mpreal x(0, prec); 2274 mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode); 2275 return x; 2276} 2277 2278 2279inline const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2280{ 2281 mpreal x(v); 2282 int tsignp; 2283 2284 if(signp) mpfr_lgamma(x.mpfr_ptr(), signp,v.mpfr_srcptr(),rnd_mode); 2285 else mpfr_lgamma(x.mpfr_ptr(),&tsignp,v.mpfr_srcptr(),rnd_mode); 2286 2287 return x; 2288} 2289 2290 2291inline const mpreal besseljn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) 2292{ 2293 mpreal y(0, x.getPrecision()); 2294 mpfr_jn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r); 2295 return y; 2296} 2297 2298inline const mpreal besselyn (long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) 2299{ 2300 mpreal y(0, x.getPrecision()); 2301 mpfr_yn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r); 2302 return y; 2303} 2304 2305inline const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2306{ 2307 mpreal a; 2308 mp_prec_t p1, p2, p3; 2309 2310 p1 = v1.get_prec(); 2311 p2 = v2.get_prec(); 2312 p3 = v3.get_prec(); 2313 2314 a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); 2315 2316 mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); 2317 return a; 2318} 2319 2320inline const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2321{ 2322 mpreal a; 2323 mp_prec_t p1, p2, p3; 2324 2325 p1 = v1.get_prec(); 2326 p2 = v2.get_prec(); 2327 p3 = v3.get_prec(); 2328 2329 a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); 2330 2331 mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); 2332 return a; 2333} 2334 2335inline const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2336{ 2337 mpreal a; 2338 mp_prec_t p1, p2; 2339 2340 p1 = v1.get_prec(); 2341 p2 = v2.get_prec(); 2342 2343 a.set_prec(p1>p2?p1:p2); 2344 2345 mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode); 2346 2347 return a; 2348} 2349 2350inline const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2351{ 2352 mpreal x; 2353 mpfr_ptr* t; 2354 unsigned long int i; 2355 2356 t = new mpfr_ptr[n]; 2357 for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp; 2358 mpfr_sum(x.mp,t,n,rnd_mode); 2359 delete[] t; 2360 return x; 2361} 2362 2363////////////////////////////////////////////////////////////////////////// 2364// MPFR 2.4.0 Specifics 2365#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 2366 2367inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2368{ 2369 return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode); 2370} 2371 2372inline const mpreal li2 (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) 2373{ 2374 MPREAL_UNARY_MATH_FUNCTION_BODY(li2); 2375} 2376 2377inline const mpreal rem (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2378{ 2379 /* R = rem(X,Y) if Y != 0, returns X - n * Y where n = trunc(X/Y). */ 2380 return fmod(x, y, rnd_mode); 2381} 2382 2383inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2384{ 2385 (void)rnd_mode; 2386 2387 /* 2388 2389 m = mod(x,y) if y != 0, returns x - n*y where n = floor(x/y) 2390 2391 The following are true by convention: 2392 - mod(x,0) is x 2393 - mod(x,x) is 0 2394 - mod(x,y) for x != y and y != 0 has the same sign as y. 2395 2396 */ 2397 2398 if(iszero(y)) return x; 2399 if(x == y) return 0; 2400 2401 mpreal m = x - floor(x / y) * y; 2402 2403 m.setSign(sgn(y)); // make sure result has the same sign as Y 2404 2405 return m; 2406} 2407 2408inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2409{ 2410 mpreal a; 2411 mp_prec_t yp, xp; 2412 2413 yp = y.get_prec(); 2414 xp = x.get_prec(); 2415 2416 a.set_prec(yp>xp?yp:xp); 2417 2418 mpfr_fmod(a.mp, x.mp, y.mp, rnd_mode); 2419 2420 return a; 2421} 2422 2423inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2424{ 2425 mpreal x(v); 2426 mpfr_rec_sqrt(x.mp,v.mp,rnd_mode); 2427 return x; 2428} 2429#endif // MPFR 2.4.0 Specifics 2430 2431////////////////////////////////////////////////////////////////////////// 2432// MPFR 3.0.0 Specifics 2433#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) 2434inline const mpreal digamma (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(digamma); } 2435inline const mpreal ai (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(ai); } 2436#endif // MPFR 3.0.0 Specifics 2437 2438////////////////////////////////////////////////////////////////////////// 2439// Constants 2440inline const mpreal const_log2 (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd()) 2441{ 2442 mpreal x(0, p); 2443 mpfr_const_log2(x.mpfr_ptr(), r); 2444 return x; 2445} 2446 2447inline const mpreal const_pi (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd()) 2448{ 2449 mpreal x(0, p); 2450 mpfr_const_pi(x.mpfr_ptr(), r); 2451 return x; 2452} 2453 2454inline const mpreal const_euler (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd()) 2455{ 2456 mpreal x(0, p); 2457 mpfr_const_euler(x.mpfr_ptr(), r); 2458 return x; 2459} 2460 2461inline const mpreal const_catalan (mp_prec_t p = mpreal::get_default_prec(), mp_rnd_t r = mpreal::get_default_rnd()) 2462{ 2463 mpreal x(0, p); 2464 mpfr_const_catalan(x.mpfr_ptr(), r); 2465 return x; 2466} 2467 2468inline const mpreal const_infinity (int sign = 1, mp_prec_t p = mpreal::get_default_prec()) 2469{ 2470 mpreal x(0, p); 2471 mpfr_set_inf(x.mpfr_ptr(), sign); 2472 return x; 2473} 2474 2475////////////////////////////////////////////////////////////////////////// 2476// Integer Related Functions 2477inline const mpreal ceil(const mpreal& v) 2478{ 2479 mpreal x(v); 2480 mpfr_ceil(x.mp,v.mp); 2481 return x; 2482} 2483 2484inline const mpreal floor(const mpreal& v) 2485{ 2486 mpreal x(v); 2487 mpfr_floor(x.mp,v.mp); 2488 return x; 2489} 2490 2491inline const mpreal round(const mpreal& v) 2492{ 2493 mpreal x(v); 2494 mpfr_round(x.mp,v.mp); 2495 return x; 2496} 2497 2498inline const mpreal trunc(const mpreal& v) 2499{ 2500 mpreal x(v); 2501 mpfr_trunc(x.mp,v.mp); 2502 return x; 2503} 2504 2505inline const mpreal rint (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint ); } 2506inline const mpreal rint_ceil (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil ); } 2507inline const mpreal rint_floor (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor); } 2508inline const mpreal rint_round (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round); } 2509inline const mpreal rint_trunc (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc); } 2510inline const mpreal frac (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(frac ); } 2511 2512////////////////////////////////////////////////////////////////////////// 2513// Miscellaneous Functions 2514inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mp,b.mp); } 2515inline const mpreal (max)(const mpreal& x, const mpreal& y){ return (x>y?x:y); } 2516inline const mpreal (min)(const mpreal& x, const mpreal& y){ return (x<y?x:y); } 2517 2518inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2519{ 2520 mpreal a; 2521 mpfr_max(a.mp,x.mp,y.mp,rnd_mode); 2522 return a; 2523} 2524 2525inline const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2526{ 2527 mpreal a; 2528 mpfr_min(a.mp,x.mp,y.mp,rnd_mode); 2529 return a; 2530} 2531 2532inline const mpreal nexttoward (const mpreal& x, const mpreal& y) 2533{ 2534 mpreal a(x); 2535 mpfr_nexttoward(a.mp,y.mp); 2536 return a; 2537} 2538 2539inline const mpreal nextabove (const mpreal& x) 2540{ 2541 mpreal a(x); 2542 mpfr_nextabove(a.mp); 2543 return a; 2544} 2545 2546inline const mpreal nextbelow (const mpreal& x) 2547{ 2548 mpreal a(x); 2549 mpfr_nextbelow(a.mp); 2550 return a; 2551} 2552 2553inline const mpreal urandomb (gmp_randstate_t& state) 2554{ 2555 mpreal x; 2556 mpfr_urandomb(x.mp,state); 2557 return x; 2558} 2559 2560#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0)) 2561// use gmp_randinit_default() to init state, gmp_randclear() to clear 2562inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2563{ 2564 mpreal x; 2565 mpfr_urandom(x.mp,state,rnd_mode); 2566 return x; 2567} 2568 2569inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2570{ 2571 mpreal x; 2572 mpfr_grandom(x.mp, NULL, state, rnd_mode); 2573 return x; 2574} 2575 2576#endif 2577 2578#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2)) 2579inline const mpreal random2 (mp_size_t size, mp_exp_t exp) 2580{ 2581 mpreal x; 2582 mpfr_random2(x.mp,size,exp); 2583 return x; 2584} 2585#endif 2586 2587// Uniformly distributed random number generation 2588// a = random(seed); <- initialization & first random number generation 2589// a = random(); <- next random numbers generation 2590// seed != 0 2591inline const mpreal random(unsigned int seed = 0) 2592{ 2593 2594#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) 2595 static gmp_randstate_t state; 2596 static bool isFirstTime = true; 2597 2598 if(isFirstTime) 2599 { 2600 gmp_randinit_default(state); 2601 gmp_randseed_ui(state,0); 2602 isFirstTime = false; 2603 } 2604 2605 if(seed != 0) gmp_randseed_ui(state,seed); 2606 2607 return mpfr::urandom(state); 2608#else 2609 if(seed != 0) std::srand(seed); 2610 return mpfr::mpreal(std::rand()/(double)RAND_MAX); 2611#endif 2612 2613} 2614 2615#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0)) 2616inline const mpreal grandom(unsigned int seed = 0) 2617{ 2618 static gmp_randstate_t state; 2619 static bool isFirstTime = true; 2620 2621 if(isFirstTime) 2622 { 2623 gmp_randinit_default(state); 2624 gmp_randseed_ui(state,0); 2625 isFirstTime = false; 2626 } 2627 2628 if(seed != 0) gmp_randseed_ui(state,seed); 2629 2630 return mpfr::grandom(state); 2631} 2632#endif 2633 2634////////////////////////////////////////////////////////////////////////// 2635// Set/Get global properties 2636inline void mpreal::set_default_prec(mp_prec_t prec) 2637{ 2638 mpfr_set_default_prec(prec); 2639} 2640 2641inline void mpreal::set_default_rnd(mp_rnd_t rnd_mode) 2642{ 2643 mpfr_set_default_rounding_mode(rnd_mode); 2644} 2645 2646inline bool mpreal::fits_in_bits(double x, int n) 2647{ 2648 int i; 2649 double t; 2650 return IsInf(x) || (std::modf ( std::ldexp ( std::frexp ( x, &i ), n ), &t ) == 0.0); 2651} 2652 2653inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2654{ 2655 mpreal x(a); 2656 mpfr_pow(x.mp,x.mp,b.mp,rnd_mode); 2657 return x; 2658} 2659 2660inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2661{ 2662 mpreal x(a); 2663 mpfr_pow_z(x.mp,x.mp,b,rnd_mode); 2664 return x; 2665} 2666 2667inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2668{ 2669 mpreal x(a); 2670 mpfr_pow_ui(x.mp,x.mp,b,rnd_mode); 2671 return x; 2672} 2673 2674inline const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode) 2675{ 2676 return pow(a,static_cast<unsigned long int>(b),rnd_mode); 2677} 2678 2679inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2680{ 2681 mpreal x(a); 2682 mpfr_pow_si(x.mp,x.mp,b,rnd_mode); 2683 return x; 2684} 2685 2686inline const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode) 2687{ 2688 return pow(a,static_cast<long int>(b),rnd_mode); 2689} 2690 2691inline const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode) 2692{ 2693 return pow(a,mpreal(b),rnd_mode); 2694} 2695 2696inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode) 2697{ 2698 return pow(a,mpreal(b),rnd_mode); 2699} 2700 2701inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) 2702{ 2703 mpreal x(a); 2704 mpfr_ui_pow(x.mp,a,b.mp,rnd_mode); 2705 return x; 2706} 2707 2708inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode) 2709{ 2710 return pow(static_cast<unsigned long int>(a),b,rnd_mode); 2711} 2712 2713inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode) 2714{ 2715 if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); 2716 else return pow(mpreal(a),b,rnd_mode); 2717} 2718 2719inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode) 2720{ 2721 if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); 2722 else return pow(mpreal(a),b,rnd_mode); 2723} 2724 2725inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode) 2726{ 2727 return pow(mpreal(a),b,rnd_mode); 2728} 2729 2730inline const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode) 2731{ 2732 return pow(mpreal(a),b,rnd_mode); 2733} 2734 2735// pow unsigned long int 2736inline const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode) 2737{ 2738 mpreal x(a); 2739 mpfr_ui_pow_ui(x.mp,a,b,rnd_mode); 2740 return x; 2741} 2742 2743inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode) 2744{ 2745 return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2746} 2747 2748inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode) 2749{ 2750 if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2751 else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow 2752} 2753 2754inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode) 2755{ 2756 if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2757 else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow 2758} 2759 2760inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode) 2761{ 2762 return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow 2763} 2764 2765inline const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode) 2766{ 2767 return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow 2768} 2769 2770// pow unsigned int 2771inline const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode) 2772{ 2773 return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui 2774} 2775 2776inline const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode) 2777{ 2778 return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2779} 2780 2781inline const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode) 2782{ 2783 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2784 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2785} 2786 2787inline const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode) 2788{ 2789 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2790 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2791} 2792 2793inline const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode) 2794{ 2795 return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2796} 2797 2798inline const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode) 2799{ 2800 return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2801} 2802 2803// pow long int 2804inline const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode) 2805{ 2806 if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui 2807 else return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui 2808} 2809 2810inline const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode) 2811{ 2812 if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2813 else return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui 2814} 2815 2816inline const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode) 2817{ 2818 if (a>0) 2819 { 2820 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2821 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2822 }else{ 2823 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si 2824 } 2825} 2826 2827inline const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode) 2828{ 2829 if (a>0) 2830 { 2831 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2832 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2833 }else{ 2834 return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si 2835 } 2836} 2837 2838inline const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode) 2839{ 2840 if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2841 else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow 2842} 2843 2844inline const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode) 2845{ 2846 if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2847 else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow 2848} 2849 2850// pow int 2851inline const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode) 2852{ 2853 if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui 2854 else return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui 2855} 2856 2857inline const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode) 2858{ 2859 if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2860 else return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui 2861} 2862 2863inline const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode) 2864{ 2865 if (a>0) 2866 { 2867 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2868 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2869 }else{ 2870 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si 2871 } 2872} 2873 2874inline const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode) 2875{ 2876 if (a>0) 2877 { 2878 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui 2879 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2880 }else{ 2881 return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si 2882 } 2883} 2884 2885inline const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode) 2886{ 2887 if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2888 else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow 2889} 2890 2891inline const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode) 2892{ 2893 if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow 2894 else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow 2895} 2896 2897// pow long double 2898inline const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode) 2899{ 2900 return pow(mpreal(a),mpreal(b),rnd_mode); 2901} 2902 2903inline const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode) 2904{ 2905 return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui 2906} 2907 2908inline const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode) 2909{ 2910 return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui 2911} 2912 2913inline const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode) 2914{ 2915 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si 2916} 2917 2918inline const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode) 2919{ 2920 return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si 2921} 2922 2923inline const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode) 2924{ 2925 return pow(mpreal(a),mpreal(b),rnd_mode); 2926} 2927 2928inline const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode) 2929{ 2930 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_ui 2931} 2932 2933inline const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode) 2934{ 2935 return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); // mpfr_pow_ui 2936} 2937 2938inline const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode) 2939{ 2940 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si 2941} 2942 2943inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode) 2944{ 2945 return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si 2946} 2947} // End of mpfr namespace 2948 2949// Explicit specialization of std::swap for mpreal numbers 2950// Thus standard algorithms will use efficient version of swap (due to Koenig lookup) 2951// Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap 2952namespace std 2953{ 2954 // we are allowed to extend namespace std with specializations only 2955 template <> 2956 inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) 2957 { 2958 return mpfr::swap(x, y); 2959 } 2960 2961 template<> 2962 class numeric_limits<mpfr::mpreal> 2963 { 2964 public: 2965 static const bool is_specialized = true; 2966 static const bool is_signed = true; 2967 static const bool is_integer = false; 2968 static const bool is_exact = false; 2969 static const int radix = 2; 2970 2971 static const bool has_infinity = true; 2972 static const bool has_quiet_NaN = true; 2973 static const bool has_signaling_NaN = true; 2974 2975 static const bool is_iec559 = true; // = IEEE 754 2976 static const bool is_bounded = true; 2977 static const bool is_modulo = false; 2978 static const bool traps = true; 2979 static const bool tinyness_before = true; 2980 2981 static const float_denorm_style has_denorm = denorm_absent; 2982 2983 inline static mpfr::mpreal (min) (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::minval(precision); } 2984 inline static mpfr::mpreal (max) (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(precision); } 2985 inline static mpfr::mpreal lowest (mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(precision); } 2986 2987 // Returns smallest eps such that 1 + eps != 1 (classic machine epsilon) 2988 inline static mpfr::mpreal epsilon(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(precision); } 2989 2990 // Returns smallest eps such that x + eps != x (relative machine epsilon) 2991 inline static mpfr::mpreal epsilon(const mpfr::mpreal& x) { return mpfr::machine_epsilon(x); } 2992 2993 inline static mpfr::mpreal round_error(mp_prec_t precision = mpfr::mpreal::get_default_prec()) 2994 { 2995 mp_rnd_t r = mpfr::mpreal::get_default_rnd(); 2996 2997 if(r == GMP_RNDN) return mpfr::mpreal(0.5, precision); 2998 else return mpfr::mpreal(1.0, precision); 2999 } 3000 3001 inline static const mpfr::mpreal infinity() { return mpfr::const_infinity(); } 3002 inline static const mpfr::mpreal quiet_NaN() { return mpfr::mpreal().setNan(); } 3003 inline static const mpfr::mpreal signaling_NaN() { return mpfr::mpreal().setNan(); } 3004 inline static const mpfr::mpreal denorm_min() { return (min)(); } 3005 3006 // Please note, exponent range is not fixed in MPFR 3007 static const int min_exponent = MPFR_EMIN_DEFAULT; 3008 static const int max_exponent = MPFR_EMAX_DEFAULT; 3009 MPREAL_PERMISSIVE_EXPR static const int min_exponent10 = (int) (MPFR_EMIN_DEFAULT * 0.3010299956639811); 3010 MPREAL_PERMISSIVE_EXPR static const int max_exponent10 = (int) (MPFR_EMAX_DEFAULT * 0.3010299956639811); 3011 3012#ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS 3013 3014 // Following members should be constant according to standard, but they can be variable in MPFR 3015 // So we define them as functions here. 3016 // 3017 // This is preferable way for std::numeric_limits<mpfr::mpreal> specialization. 3018 // But it is incompatible with standard std::numeric_limits and might not work with other libraries, e.g. boost. 3019 // See below for compatible implementation. 3020 inline static float_round_style round_style() 3021 { 3022 mp_rnd_t r = mpfr::mpreal::get_default_rnd(); 3023 3024 switch (r) 3025 { 3026 case GMP_RNDN: return round_to_nearest; 3027 case GMP_RNDZ: return round_toward_zero; 3028 case GMP_RNDU: return round_toward_infinity; 3029 case GMP_RNDD: return round_toward_neg_infinity; 3030 default: return round_indeterminate; 3031 } 3032 } 3033 3034 inline static int digits() { return int(mpfr::mpreal::get_default_prec()); } 3035 inline static int digits(const mpfr::mpreal& x) { return x.getPrecision(); } 3036 3037 inline static int digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec()) 3038 { 3039 return mpfr::bits2digits(precision); 3040 } 3041 3042 inline static int digits10(const mpfr::mpreal& x) 3043 { 3044 return mpfr::bits2digits(x.getPrecision()); 3045 } 3046 3047 inline static int max_digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec()) 3048 { 3049 return digits10(precision); 3050 } 3051#else 3052 // Digits and round_style are NOT constants when it comes to mpreal. 3053 // If possible, please use functions digits() and round_style() defined above. 3054 // 3055 // These (default) values are preserved for compatibility with existing libraries, e.g. boost. 3056 // Change them accordingly to your application. 3057 // 3058 // For example, if you use 256 bits of precision uniformly in your program, then: 3059 // digits = 256 3060 // digits10 = 77 3061 // max_digits10 = 78 3062 // 3063 // Approximate formula for decimal digits is: digits10 = floor(log10(2) * digits). See bits2digits() for more details. 3064 3065 static const std::float_round_style round_style = round_to_nearest; 3066 static const int digits = 53; 3067 static const int digits10 = 15; 3068 static const int max_digits10 = 16; 3069#endif 3070 }; 3071 3072} 3073 3074#endif /* __MPREAL_H__ */ 3075