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