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