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