source: pacpussensors/trunk/Vislab/lib3dv/eigen/unsupported/test/mpreal/mpreal.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

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