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