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