mpreal.h
Go to the documentation of this file.
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
149  inline static mp_rnd_t get_default_rnd() { return (mp_rnd_t)(mpfr_get_default_rounding_mode()); }
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
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
549 
550  // "Smart" resources deallocation. Checks if instance initialized before deletion.
551  void clear(::mpfr_ptr);
552 };
553 
555 // Exceptions
556 class conversion_overflow : public std::exception {
557 public:
558  std::string why() { return "inexact conversion from floating point"; }
559 };
560 
562 // Constructors & converters
563 // Default constructor: creates mp number and initializes it to 0.
565 {
566  mpfr_init2(mpfr_ptr(), mpreal::get_default_prec());
568 
570 }
571 
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 
578 }
579 
580 #ifdef MPREAL_HAVE_MOVE_SUPPORT
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 
587 }
588 
589 inline mpreal& mpreal::operator=(mpreal&& other)
590 {
591  mpfr_swap(mpfr_ptr(), other.mpfr_ptr());
592 
594  return *this;
595 }
596 #endif
597 
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 
611 }
612 
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 
619 }
620 
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 
627 }
628 
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 
635 }
636 
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)
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 
652 }
653 
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 
660 }
661 
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 
668 }
669 
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 
676 }
677 
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 
684 }
685 
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 
692 }
693 
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 
700 }
701 
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 
708 }
709 
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 
716 }
717 
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 
724 }
725 
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 
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 
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 
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 
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 
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 
897 // min, max
898 const mpreal (max)(const mpreal& x, const mpreal& y);
899 const mpreal (min)(const mpreal& x, const mpreal& y);
900 
902 // Implementation
904 
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 
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 
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 
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 
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 
955  return *this;
956 }
957 
958 inline mpreal& mpreal::operator=(const double v)
959 {
960 #if (MPREAL_DOUBLE_BITS_OVERFLOW > -1)
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 
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 
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 
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 
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 
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 
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 
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());
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());
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 
1075 // + Addition
1077 {
1080  return *this;
1081 }
1082 
1083 inline mpreal& mpreal::operator+=(const mpf_t u)
1084 {
1085  *this += mpreal(u);
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());
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());
1101  return *this;
1102 }
1103 
1104 inline mpreal& mpreal::operator+= (const long double u)
1105 {
1106  *this += mpreal(u);
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 
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());
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());
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());
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());
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 
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 
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 
1194 // - Subtraction
1196 {
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());
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());
1213  return *this;
1214 }
1215 
1216 inline mpreal& mpreal::operator-=(const long double v)
1217 {
1218  *this -= mpreal(v);
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 
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());
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());
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());
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());
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 
1319 // * Multiplication
1321 {
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());
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());
1338  return *this;
1339 }
1340 
1341 inline mpreal& mpreal::operator*=(const long double v)
1342 {
1343  *this *= mpreal(v);
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
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());
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());
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());
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());
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 
1395 // / Division
1397 {
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());
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());
1414  return *this;
1415 }
1416 
1417 inline mpreal& mpreal::operator/=(const long double v)
1418 {
1419  *this /= mpreal(v);
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
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());
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());
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());
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());
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 
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());
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());
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());
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());
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());
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());
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());
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());
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 
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 
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 
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(),
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 
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 
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);
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);
1955  return *this;
1956 }
1957 
1959 {
1960  mpfr_set_inf(mpfr_ptr(), sign);
1962  return *this;
1963 }
1964 
1966 {
1967  mpfr_set_nan(mpfr_ptr());
1969  return *this;
1970 }
1971 
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 
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);
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);
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 
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 
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);
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 
2157 // Mathematical Functions
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())
2166 
2167 inline const mpreal sqrt (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
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 ); }
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 ); }
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 ); }
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 ); }
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 
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 {
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 
2464 // MPFR 3.0.0 Specifics
2465 #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
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 
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 
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 ); }
2542 inline const mpreal frac (const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(frac ); }
2543 
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 
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__ */
friend const mpreal nexttoward(const mpreal &x, const mpreal &y)
Definition: mpreal.h:2564
mp_exp_t get_exp()
Definition: mpreal.h:1996
Matrix3f m
static const mpfr::mpreal infinity()
Definition: mpreal.h:3031
bool operator>=(const mpreal &a, const mpreal &b)
Definition: mpreal.h:1660
friend const mpreal tan(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2232
mpreal & setInf(int Sign=+1)
Definition: mpreal.h:1958
friend const mpreal exp2(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2228
friend const mpreal sum(const mpreal tab[], const unsigned long int n, int &status, mp_rnd_t rnd_mode)
Definition: mpreal.h:2381
friend int sin_cos(mpreal &s, mpreal &c, const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2213
friend const mpreal acosh(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2255
friend bool() isinf(const mpreal &v)
Definition: mpreal.h:1701
friend const mpreal div_2si(const mpreal &v, long int k, mp_rnd_t rnd_mode)
Definition: mpreal.h:1632
friend const mpreal rint_trunc(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2541
friend const mpreal acoth(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2245
friend std::istream & operator>>(std::istream &is, mpreal &v)
Definition: mpreal.h:1904
friend const mpreal root(const mpreal &v, unsigned long int k, mp_rnd_t rnd_mode)
Definition: mpreal.h:2194
friend const mpreal remquo(long *q, const mpreal &x, const mpreal &y, mp_rnd_t rnd_mode)
Definition: mpreal.h:2294
friend const mpreal hypot(const mpreal &x, const mpreal &y, mp_rnd_t rnd_mode)
Definition: mpreal.h:2280
friend const mpreal modf(const mpreal &v, mpreal &n)
Definition: mpreal.h:2094
friend const mpreal urandomb(gmp_randstate_t &state)
Definition: mpreal.h:2585
static int set_emax(mp_exp_t exp)
Definition: mpreal.h:2131
#define EIGEN_NOT_A_MACRO
Definition: Macros.h:327
friend const mpreal operator/(const unsigned long int b, const mpreal &a)
Definition: mpreal.h:1470
mp_prec_t digits2bits(int d)
Definition: mpreal.h:1918
Scalar * y
const mpreal scalbn(const mpreal &v, mp_exp_t exp)
Definition: mpreal.h:2029
friend bool isint(const mpreal &v)
Definition: mpreal.h:1704
Vector v2
friend int sgn(const mpreal &v)
Definition: mpreal.h:1934
Scalar * b
Definition: benchVecAdd.cpp:17
static int max_digits10(mp_prec_t precision=mpfr::mpreal::get_default_prec())
Definition: mpreal.h:3077
static Point3 p3
static mp_prec_t get_default_prec()
Definition: mpreal.h:150
return int(ret)+1
static mpfr::mpreal epsilon(mp_prec_t precision=mpfr::mpreal::get_default_prec())
Definition: mpreal.h:3018
Vector v1
Vector3f p1
void set_prec(mp_prec_t prec, mp_rnd_t rnd_mode=get_default_rnd())
Definition: mpreal.h:1990
mp_prec_t get_prec() const
Definition: mpreal.h:1985
friend const mpreal log(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2224
friend int cmpabs(const mpreal &a, const mpreal &b)
Definition: mpreal.h:2208
friend const mpreal rint_ceil(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2538
mpreal & operator>>=(const unsigned long int u)
Definition: mpreal.h:1541
friend bool isnum(const mpreal &v)
const mpreal ldexp(const mpreal &v, mp_exp_t exp)
Definition: mpreal.h:2020
ArrayXcf v
Definition: Cwise_arg.cpp:1
static const mpfr::mpreal denorm_min()
Definition: mpreal.h:3034
friend const mpreal mod(const mpreal &x, const mpreal &y, mp_rnd_t rnd_mode)
Definition: mpreal.h:2415
static mpfr::mpreal epsilon(const mpfr::mpreal &x)
Definition: mpreal.h:3021
friend const mpreal rint_round(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2540
friend const mpreal fmin(const mpreal &x, const mpreal &y, mp_rnd_t rnd_mode)
Definition: mpreal.h:2557
int n
static float_round_style round_style()
Definition: mpreal.h:3050
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
double toDouble(mp_rnd_t mode=GMP_RNDN) const
Definition: mpreal.h:1716
#define MPREAL_DOUBLE_BITS_OVERFLOW
Definition: mpreal.h:125
#define MPREAL_MSVC_DEBUGVIEW_CODE
Definition: mpreal.h:114
bool signbit(const mpreal &x)
Definition: mpreal.h:2089
friend const mpreal eint(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2261
friend const mpreal sqrt(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2167
static mp_exp_t get_emin(void)
Definition: mpreal.h:2116
Definition: Half.h:150
std::ostream & output(std::ostream &os) const
Definition: mpreal.h:1873
friend const mpreal rint_floor(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2539
Definition: mpreal.h:140
friend const mpreal random2(mp_size_t size, mp_exp_t exp)
Definition: mpreal.h:2602
mpreal & setPrecision(int Precision, mp_rnd_t RoundingMode=get_default_rnd())
Definition: mpreal.h:1951
bool operator==(const mpreal &a, const mpreal &b)
Definition: mpreal.h:1684
friend const mpreal atan(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2238
const mpreal() max(const mpreal &x, const mpreal &y)
Definition: mpreal.h:2547
mpreal & setSign(int Sign, mp_rnd_t RoundingMode=get_default_rnd())
Definition: mpreal.h:1939
std::string why()
Definition: mpreal.h:558
bool operator!=(const mpreal &a, const mpreal &b)
Definition: mpreal.h:1692
friend const mpreal urandom(gmp_randstate_t &state, mp_rnd_t rnd_mode)
Definition: mpreal.h:2593
friend const mpreal random(unsigned int seed)
Definition: mpreal.h:2614
friend const mpreal tanh(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2251
friend const mpreal digamma(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2466
friend const mpreal div_2ui(const mpreal &v, unsigned long int k, mp_rnd_t rnd_mode)
Definition: mpreal.h:1625
#define MPREAL_PERMISSIVE_EXPR
Definition: mpreal.h:137
friend const mpreal agm(const mpreal &v1, const mpreal &v2, mp_rnd_t rnd_mode)
Definition: mpreal.h:2366
mpreal & operator--()
Definition: mpreal.h:1181
friend const mpreal log10(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2226
friend const mpreal lngamma(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2264
friend const mpreal nextabove(const mpreal &x)
Definition: mpreal.h:2571
friend const mpreal const_pi(mp_prec_t prec, mp_rnd_t rnd_mode)
Definition: mpreal.h:2479
static int digits10(const mpfr::mpreal &x)
Definition: mpreal.h:3072
#define IsInf(x)
Definition: mpreal.h:80
friend const mpreal mul_2si(const mpreal &v, long int k, mp_rnd_t rnd_mode)
Definition: mpreal.h:1618
friend const mpreal grandom(gmp_randstate_t &state, mp_rnd_t rnd_mode)
Definition: mpreal.h:2639
static const mpfr::mpreal signaling_NaN()
Definition: mpreal.h:3033
friend const mpreal besseljn(long n, const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2322
friend const mpreal asec(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2243
Array33i a
friend bool() isnan(const mpreal &v)
Definition: mpreal.h:1700
mpfr_t mp
Definition: mpreal.h:144
friend const mpreal remainder(const mpreal &x, const mpreal &y, mp_rnd_t rnd_mode)
Definition: mpreal.h:2287
long long toLLong(mp_rnd_t mode=GMP_RNDZ) const
Definition: mpreal.h:1718
long double toLDouble(mp_rnd_t mode=GMP_RNDN) const
Definition: mpreal.h:1717
friend const mpreal fma(const mpreal &v1, const mpreal &v2, const mpreal &v3, mp_rnd_t rnd_mode)
Definition: mpreal.h:2336
friend const mpreal sinh(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2250
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
void swap(mpreal &a, mpreal &b)
Definition: mpreal.h:2546
friend bool iszero(const mpreal &v)
Definition: mpreal.h:1703
std::string toString(int n=-1, int b=10, mp_rnd_t mode=mpreal::get_default_rnd()) const
Definition: mpreal.h:1755
EIGEN_DEVICE_FUNC const CeilReturnType ceil() const
friend const mpreal sech(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2252
friend void swap(mpreal &x, mpreal &y)
Definition: mpreal.h:2546
friend const mpreal acsch(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2247
Vector v3
static mp_exp_t get_emin_min(void)
Definition: mpreal.h:2136
mpreal copysign(const mpreal &x, const mpreal &y, mp_rnd_t rnd_mode=mpreal::get_default_rnd())
Definition: mpreal.h:2082
const mpreal modf(const mpreal &v, mpreal &n)
Definition: mpreal.h:2094
static mp_exp_t get_emax(void)
Definition: mpreal.h:2126
const mpreal operator-() const
Definition: mpreal.h:1263
friend const mpreal abs(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2223
const mpreal urandom(gmp_randstate_t &state, mp_rnd_t rnd_mode=mpreal::get_default_rnd())
Definition: mpreal.h:2593
EIGEN_DEVICE_FUNC const SignReturnType sign() const
friend const mpreal lgamma(const mpreal &v, int *signp, mp_rnd_t rnd_mode)
Definition: mpreal.h:2310
friend const mpreal coth(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2254
friend const mpreal fmax(const mpreal &x, const mpreal &y, mp_rnd_t rnd_mode)
Definition: mpreal.h:2550
friend const mpreal gamma(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2262
static mp_exp_t get_emax_max(void)
Definition: mpreal.h:2151
const mpreal() min(const mpreal &x, const mpreal &y)
Definition: mpreal.h:2548
friend const mpreal csch(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2253
friend const mpreal exp(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2227
float * ptr
bool operator>(const mpreal &a, const mpreal &b)
Definition: mpreal.h:1652
friend bool isregular(const mpreal &v)
Definition: mpreal.h:1707
static const mpfr::mpreal quiet_NaN()
Definition: mpreal.h:3032
const mpreal grandom(gmp_randstate_t &state, mp_rnd_t rnd_mode=mpreal::get_default_rnd())
Definition: mpreal.h:2639
friend const mpreal floor(const mpreal &v)
Definition: mpreal.h:2516
mpreal maxval(mp_prec_t prec=mpreal::get_default_prec())
Definition: mpreal.h:2059
friend const mpreal tgamma(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2263
friend const mpreal fmod(const mpreal &x, const mpreal &y, mp_rnd_t rnd_mode)
Definition: mpreal.h:2440
friend const mpreal besselj1(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2269
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
unsigned long toULong(mp_rnd_t mode=GMP_RNDZ) const
Definition: mpreal.h:1714
friend const mpreal acos(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2236
static void set_default_prec(mp_prec_t prec)
Definition: mpreal.h:2666
Array< double, 1, 3 > e(1./3., 0.5, 2.)
RealScalar s
#define MPREAL_UNARY_MATH_FUNCTION_BODY(f)
Definition: mpreal.h:2159
friend const mpreal besselyn(long n, const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2329
EIGEN_DEVICE_FUNC const Scalar & q
friend const mpreal asech(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2246
mpreal & operator=(const mpreal &v)
Definition: mpreal.h:907
mpreal minval(mp_prec_t prec=mpreal::get_default_prec())
Definition: mpreal.h:2052
friend const mpreal const_euler(mp_prec_t prec, mp_rnd_t rnd_mode)
Definition: mpreal.h:2486
friend const mpreal expm1(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2260
mpreal & operator<<=(const unsigned long int u)
Definition: mpreal.h:1513
friend const mpreal zeta(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2265
bool operator<=(const mpreal &a, const mpreal &b)
Definition: mpreal.h:1676
friend std::ostream & operator<<(std::ostream &os, const mpreal &v)
Definition: mpreal.h:1899
#define MPREAL_MSVC_DEBUGVIEW_DATA
Definition: mpreal.h:115
#define NULL
Definition: ccolamd.c:609
MPREAL_MSVC_DEBUGVIEW_DATA void clear(::mpfr_ptr)
Definition: mpreal.h:726
static int digits(const mpfr::mpreal &x)
Definition: mpreal.h:3065
::mpfr_srcptr mpfr_srcptr() const
Definition: mpreal.h:1723
friend const mpreal ai(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2467
::mpfr_ptr mpfr_ptr()
Definition: mpreal.h:1721
friend const mpreal asin(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2237
const mpreal frexp(const mpreal &x, mp_exp_t *exp, mp_rnd_t mode=mpreal::get_default_rnd())
Definition: mpreal.h:2008
friend const mpreal fabs(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2222
mpreal & operator++()
Definition: mpreal.h:1169
friend const mpreal trunc(const mpreal &v)
Definition: mpreal.h:2530
friend const mpreal atanh(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2257
friend const mpreal asinh(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2256
mpreal machine_epsilon(mp_prec_t prec=mpreal::get_default_prec())
Definition: mpreal.h:2034
mpreal & operator*=(const long long int u)
Definition: mpreal.h:1155
unsigned long long toULLong(mp_rnd_t mode=GMP_RNDZ) const
Definition: mpreal.h:1719
bool fits_in_bits(double x, int n)
Definition: mpreal.h:2676
#define mpfr_set_zero_fast(x)
Definition: mpreal.h:132
friend const mpreal mul_2ui(const mpreal &v, unsigned long int k, mp_rnd_t rnd_mode)
Definition: mpreal.h:1610
friend bool() isfinite(const mpreal &v)
Definition: mpreal.h:1702
int set_exp(mp_exp_t e)
Definition: mpreal.h:2001
friend const mpreal cot(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2235
bool toBool() const
Definition: mpreal.h:1712
int subnormalize(int t, mp_rnd_t rnd_mode=get_default_rnd())
Definition: mpreal.h:2109
EIGEN_DEVICE_FUNC const FloorReturnType floor() const
ofstream os("timeSchurFactors.csv")
const mpreal operator+() const
Definition: mpreal.h:1160
mpreal & setNan()
Definition: mpreal.h:1965
friend const mpreal round(const mpreal &v)
Definition: mpreal.h:2523
const internal::result_type< Rhs >::type operator*(const mpreal &lhs, const Rhs &rhs)
Definition: mpreal.h:781
friend const mpreal erf(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2266
int bits2digits(mp_prec_t b)
Definition: mpreal.h:1925
bool operator<(const mpreal &a, const mpreal &b)
Definition: mpreal.h:1668
friend const mpreal acsc(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2244
mpreal & setZero(int Sign=+1)
Definition: mpreal.h:1972
static int set_emin(mp_exp_t exp)
Definition: mpreal.h:2121
static mpfr::mpreal round_error(mp_prec_t precision=mpfr::mpreal::get_default_prec())
Definition: mpreal.h:3023
cout precision(2)
friend const mpreal bessely0(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2270
static mp_exp_t get_emax_min(void)
Definition: mpreal.h:2146
friend const mpreal erfc(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2267
float * p
friend const mpreal pow(const mpreal &a, const mpreal &b, mp_rnd_t rnd_mode)
Definition: mpreal.h:2683
friend const mpreal const_infinity(int sign, mp_prec_t prec)
Definition: mpreal.h:2500
mpreal & operator+=(const mpreal &v)
Definition: mpreal.h:1076
friend const mpreal log2(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2225
mpreal & operator-=(const long long int u)
Definition: mpreal.h:1153
static Point3 p2
friend const mpreal cos(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2230
friend const mpreal sqr(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2164
int getPrecision() const
Definition: mpreal.h:1946
friend const mpreal exp10(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2229
friend const mpreal logb(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2240
friend const mpreal const_log2(mp_prec_t prec, mp_rnd_t rnd_mode)
Definition: mpreal.h:2472
friend const mpreal sec(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2233
const mpreal const_infinity(int sign=1, mp_prec_t p=mpreal::get_default_prec())
Definition: mpreal.h:2500
friend const mpreal rint(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2537
Annotation indicating that a class derives from another given type.
Definition: attr.h:42
friend const mpreal acot(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2242
static mpfr::mpreal lowest(mp_prec_t precision=mpfr::mpreal::get_default_prec())
Definition: mpreal.h:3015
friend const mpreal log1p(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2259
friend const mpreal nextbelow(const mpreal &x)
Definition: mpreal.h:2578
friend const mpreal rec_sqrt(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2455
mpreal & operator/=(const long long int u)
Definition: mpreal.h:1157
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
friend const mpreal fms(const mpreal &v1, const mpreal &v2, const mpreal &v3, mp_rnd_t rnd_mode)
Definition: mpreal.h:2351
friend const mpreal frac(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2542
friend const mpreal sin(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2231
friend const mpreal dim(const mpreal &a, const mpreal &b, mp_rnd_t rnd_mode)
Definition: mpreal.h:2201
Values initialize(const NonlinearFactorGraph &graph, bool useOdometricPath)
Definition: lago.cpp:338
friend const mpreal rem(const mpreal &x, const mpreal &y, mp_rnd_t rnd_mode)
Definition: mpreal.h:2409
bool isEqualFuzzy(const mpreal &a, const mpreal &b, const mpreal &eps)
Definition: mpreal.h:2070
friend const mpreal const_catalan(mp_prec_t prec, mp_rnd_t rnd_mode)
Definition: mpreal.h:2493
friend const mpreal ceil(const mpreal &v)
Definition: mpreal.h:2509
friend const mpreal li2(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2404
friend const mpreal besselj0(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2268
static mp_rnd_t get_default_rnd()
Definition: mpreal.h:149
float toFloat(mp_rnd_t mode=GMP_RNDN) const
Definition: mpreal.h:1715
static mp_exp_t get_emin_max(void)
Definition: mpreal.h:2141
long toLong(mp_rnd_t mode=GMP_RNDZ) const
Definition: mpreal.h:1713
friend const mpreal cbrt(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2221
static void set_default_rnd(mp_rnd_t rnd_mode)
Definition: mpreal.h:2671
friend const mpreal bessely1(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2271
bool isEqualUlps(const mpreal &a, const mpreal &b, int maxUlps)
Definition: mpreal.h:2065
Point2 t(10, 10)
friend const mpreal fac_ui(unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode)
Definition: mpreal.h:2301
friend const mpreal atan2(const mpreal &y, const mpreal &x, mp_rnd_t rnd_mode)
Definition: mpreal.h:2273
static int digits10(mp_prec_t precision=mpfr::mpreal::get_default_prec())
Definition: mpreal.h:3067
friend int sinh_cosh(mpreal &s, mpreal &c, const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2399
friend const mpreal csc(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2234
friend const mpreal cosh(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2249
int check_range(int t, mp_rnd_t rnd_mode=get_default_rnd())
Definition: mpreal.h:2104


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:43:01