Ticket #6296: mpreal.h

File mpreal.h, 98.9 KB (added by anonymous, 11 years ago)
Line 
1/*
2 Multi-precision real number class. C++ interface for MPFR library.
3 Project homepage: http://www.holoborodko.com/pavel/
4 Contact e-mail: pavel@holoborodko.com
5
6 Copyright (c) 2008-2011 Pavel Holoborodko
7
8 Core Developers:
9 Pavel Holoborodko, Dmitriy Gubanov, Konstantin Holoborodko.
10
11 Contributors:
12 Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze,
13 Heinz van Saanen, Pere Constans, Peter van Hoof, Gael Guennebaud,
14 Tsai Chia Cheng, Alexei Zubanov.
15
16 ****************************************************************************
17 This library is free software; you can redistribute it and/or
18 modify it under the terms of the GNU Lesser General Public
19 License as published by the Free Software Foundation; either
20 version 2.1 of the License, or (at your option) any later version.
21
22 This library is distributed in the hope that it will be useful,
23 but WITHOUT ANY WARRANTY; without even the implied warranty of
24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
25 Lesser General Public License for more details.
26
27 You should have received a copy of the GNU Lesser General Public
28 License along with this library; if not, write to the Free Software
29 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
30
31 ****************************************************************************
32 Redistribution and use in source and binary forms, with or without
33 modification, are permitted provided that the following conditions
34 are met:
35
36 1. Redistributions of source code must retain the above copyright
37 notice, this list of conditions and the following disclaimer.
38
39 2. Redistributions in binary form must reproduce the above copyright
40 notice, this list of conditions and the following disclaimer in the
41 documentation and/or other materials provided with the distribution.
42
43 3. The name of the author may be used to endorse or promote products
44 derived from this software without specific prior written permission.
45
46 THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
47 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
48 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
49 ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
50 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
51 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
52 OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
53 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
54 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
55 OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
56 SUCH DAMAGE.
57*/
58
59#ifndef __MPREAL_H__
60#define __MPREAL_H__
61
62#include <string>
63#include <iostream>
64#include <sstream>
65#include <stdexcept>
66#include <cfloat>
67#include <cmath>
68
69// Options
70#define MPREAL_HAVE_INT64_SUPPORT // int64_t support: available only for MSVC 2010 & GCC
71
72// Detect compiler using signatures from http://predef.sourceforge.net/
73#if defined(__GNUC__) && defined(__INTEL_COMPILER)
74 #define IsInf(x) isinf(x) // Intel ICC compiler on Linux
75
76#elif defined(_MSC_VER) // Microsoft Visual C++
77 #define IsInf(x) (!_finite(x))
78
79#elif defined(__GNUC__)
80 #define IsInf(x) std::isinf(x) // GNU C/C++
81
82#else
83 #define IsInf(x) std::isinf(x) // Unknown compiler, just hope for C99 conformance
84#endif
85
86#if defined(MPREAL_HAVE_INT64_SUPPORT)
87
88 #define MPFR_USE_INTMAX_T // should be defined before mpfr.h
89
90 #if defined(_MSC_VER) // <stdint.h> is available only in msvc2010!
91 #if (_MSC_VER >= 1600)
92 #include <stdint.h>
93 #else // MPFR relies on intmax_t which is available only in msvc2010
94 #undef MPREAL_HAVE_INT64_SUPPORT // Besides, MPFR - MPIR have to be compiled with msvc2010
95 #undef MPFR_USE_INTMAX_T // Since we cannot detect this, disable x64 by default
96 // Someone should change this manually if needed.
97 #endif
98 #endif
99
100 #if defined (__GNUC__)
101 #if defined(__amd64__) || defined(__amd64) || defined(__x86_64__) || defined(__x86_64)
102 #undef MPREAL_HAVE_INT64_SUPPORT // remove all shaman dances for x64 builds since
103 #undef MPFR_USE_INTMAX_T // GCC already support x64 as of "long int" is 64-bit integer, nothing left to do
104 #else
105 #include <stdint.h> // use int64_t, uint64_t otherwise.
106 #endif
107 #endif
108
109#endif
110
111#include <mpfr.h>
112
113#if (MPFR_VERSION < MPFR_VERSION_NUM(3,0,0))
114 #include <cstdlib> // needed for random()
115#endif
116
117namespace mpfr {
118
119class mpreal {
120private:
121 mpfr_t mp;
122
123public:
124 static mp_rnd_t default_rnd;
125 static mp_prec_t default_prec;
126 static int default_base;
127 static int double_bits;
128
129public:
130 // Constructors && type conversion
131 mpreal();
132 mpreal(const mpreal& u);
133 mpreal(const mpfr_t u);
134 mpreal(const mpf_t u);
135 mpreal(const mpz_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
136 mpreal(const mpq_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
137 mpreal(const double u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
138 mpreal(const long double u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
139 mpreal(const unsigned long int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
140 mpreal(const unsigned int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
141 mpreal(const long int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
142 mpreal(const int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
143
144#if defined (MPREAL_HAVE_INT64_SUPPORT)
145 mpreal(const uint64_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
146 mpreal(const int64_t u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
147#endif
148
149 mpreal(const char* s, mp_prec_t prec = default_prec, int base = default_base, mp_rnd_t mode = default_rnd);
150 mpreal(const std::string& s, mp_prec_t prec = default_prec, int base = default_base, mp_rnd_t mode = default_rnd);
151
152 ~mpreal();
153
154 // Operations
155 // =
156 // +, -, *, /, ++, --, <<, >>
157 // *=, +=, -=, /=,
158 // <, >, ==, <=, >=
159
160 // =
161 mpreal& operator=(const mpreal& v);
162 mpreal& operator=(const mpf_t v);
163 mpreal& operator=(const mpz_t v);
164 mpreal& operator=(const mpq_t v);
165 mpreal& operator=(const long double v);
166 mpreal& operator=(const double v);
167 mpreal& operator=(const unsigned long int v);
168 mpreal& operator=(const unsigned int v);
169 mpreal& operator=(const long int v);
170 mpreal& operator=(const int v);
171 mpreal& operator=(const char* s);
172
173 // +
174 mpreal& operator+=(const mpreal& v);
175 mpreal& operator+=(const mpf_t v);
176 mpreal& operator+=(const mpz_t v);
177 mpreal& operator+=(const mpq_t v);
178 mpreal& operator+=(const long double u);
179 mpreal& operator+=(const double u);
180 mpreal& operator+=(const unsigned long int u);
181 mpreal& operator+=(const unsigned int u);
182 mpreal& operator+=(const long int u);
183 mpreal& operator+=(const int u);
184
185#if defined (MPREAL_HAVE_INT64_SUPPORT)
186 mpreal& operator+=(const int64_t u);
187 mpreal& operator+=(const uint64_t u);
188 mpreal& operator-=(const int64_t u);
189 mpreal& operator-=(const uint64_t u);
190 mpreal& operator*=(const int64_t u);
191 mpreal& operator*=(const uint64_t u);
192 mpreal& operator/=(const int64_t u);
193 mpreal& operator/=(const uint64_t u);
194#endif
195
196 const mpreal operator+() const;
197 mpreal& operator++ ();
198 const mpreal operator++ (int);
199
200 // -
201 mpreal& operator-=(const mpreal& v);
202 mpreal& operator-=(const mpz_t v);
203 mpreal& operator-=(const mpq_t v);
204 mpreal& operator-=(const long double u);
205 mpreal& operator-=(const double u);
206 mpreal& operator-=(const unsigned long int u);
207 mpreal& operator-=(const unsigned int u);
208 mpreal& operator-=(const long int u);
209 mpreal& operator-=(const int u);
210 const mpreal operator-() const;
211 friend const mpreal operator-(const unsigned long int b, const mpreal& a);
212 friend const mpreal operator-(const unsigned int b, const mpreal& a);
213 friend const mpreal operator-(const long int b, const mpreal& a);
214 friend const mpreal operator-(const int b, const mpreal& a);
215 friend const mpreal operator-(const double b, const mpreal& a);
216 mpreal& operator-- ();
217 const mpreal operator-- (int);
218
219 // *
220 mpreal& operator*=(const mpreal& v);
221 mpreal& operator*=(const mpz_t v);
222 mpreal& operator*=(const mpq_t v);
223 mpreal& operator*=(const long double v);
224 mpreal& operator*=(const double v);
225 mpreal& operator*=(const unsigned long int v);
226 mpreal& operator*=(const unsigned int v);
227 mpreal& operator*=(const long int v);
228 mpreal& operator*=(const int v);
229
230 // /
231 mpreal& operator/=(const mpreal& v);
232 mpreal& operator/=(const mpz_t v);
233 mpreal& operator/=(const mpq_t v);
234 mpreal& operator/=(const long double v);
235 mpreal& operator/=(const double v);
236 mpreal& operator/=(const unsigned long int v);
237 mpreal& operator/=(const unsigned int v);
238 mpreal& operator/=(const long int v);
239 mpreal& operator/=(const int v);
240 friend const mpreal operator/(const unsigned long int b, const mpreal& a);
241 friend const mpreal operator/(const unsigned int b, const mpreal& a);
242 friend const mpreal operator/(const long int b, const mpreal& a);
243 friend const mpreal operator/(const int b, const mpreal& a);
244 friend const mpreal operator/(const double b, const mpreal& a);
245
246 //<<= Fast Multiplication by 2^u
247 mpreal& operator<<=(const unsigned long int u);
248 mpreal& operator<<=(const unsigned int u);
249 mpreal& operator<<=(const long int u);
250 mpreal& operator<<=(const int u);
251
252 //>>= Fast Division by 2^u
253 mpreal& operator>>=(const unsigned long int u);
254 mpreal& operator>>=(const unsigned int u);
255 mpreal& operator>>=(const long int u);
256 mpreal& operator>>=(const int u);
257
258 // Boolean Operators
259 friend bool operator > (const mpreal& a, const mpreal& b);
260 friend bool operator >= (const mpreal& a, const mpreal& b);
261 friend bool operator < (const mpreal& a, const mpreal& b);
262 friend bool operator <= (const mpreal& a, const mpreal& b);
263 friend bool operator == (const mpreal& a, const mpreal& b);
264 friend bool operator != (const mpreal& a, const mpreal& b);
265
266 // Optimized specializations for boolean operators
267 friend bool operator == (const mpreal& a, const unsigned long int b);
268 friend bool operator == (const mpreal& a, const unsigned int b);
269 friend bool operator == (const mpreal& a, const long int b);
270 friend bool operator == (const mpreal& a, const int b);
271 friend bool operator == (const mpreal& a, const long double b);
272 friend bool operator == (const mpreal& a, const double b);
273
274 // Type Conversion operators
275 long toLong() const;
276 unsigned long toULong() const;
277 double toDouble() const;
278 long double toLDouble() const;
279
280#if defined (MPREAL_HAVE_INT64_SUPPORT)
281 int64_t toInt64() const;
282 uint64_t toUInt64() const;
283#endif
284
285 // Convert mpreal to string with n significant digits in base b
286 // n = 0 -> convert with the maximum available digits
287 std::string toString(size_t n = 0, int b = default_base, mp_rnd_t mode = default_rnd) const;
288
289#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
290 std::string toString(const std::string& format) const;
291#endif
292
293 operator std::string() const;
294 inline operator mpfr_ptr();
295 inline operator mpfr_srcptr() const;
296
297 // Math Functions
298 friend const mpreal sqr(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
299 friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
300 friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode = mpreal::default_rnd);
301 friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
302 friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd);
303 friend const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
304 friend const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::default_rnd);
305 friend const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
306 friend const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
307 friend const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
308 friend const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
309 friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
310
311 friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
312 friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
313 friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd);
314 friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::default_rnd);
315 friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd);
316 friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::default_rnd);
317 friend int cmpabs(const mpreal& a,const mpreal& b);
318
319 friend const mpreal log (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
320 friend const mpreal log2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
321 friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
322 friend const mpreal exp (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
323 friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
324 friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
325 friend const mpreal log1p (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
326 friend const mpreal expm1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
327
328 friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
329 friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
330 friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
331 friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
332 friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
333 friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
334 friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
335
336 friend const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
337 friend const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
338 friend const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
339 friend const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::default_rnd);
340 friend const mpreal acot (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
341 friend const mpreal asec (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
342 friend const mpreal acsc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
343
344 friend const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
345 friend const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
346 friend const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
347 friend const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
348 friend const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
349 friend const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
350 friend const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
351 friend const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
352 friend const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
353 friend const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
354 friend const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
355 friend const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
356
357 friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd);
358
359 friend const mpreal fac_ui (unsigned long int v, mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd);
360 friend const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
361
362 friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
363 friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
364 friend const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::default_rnd);
365 friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
366 friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
367 friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
368 friend const mpreal _j0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
369 friend const mpreal _j1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
370 friend const mpreal _jn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
371 friend const mpreal _y0 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
372 friend const mpreal _y1 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
373 friend const mpreal _yn (long n, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
374 friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd);
375 friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd);
376 friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::default_rnd);
377 friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::default_rnd);
378 friend int sgn(const mpreal& v); // -1 or +1
379
380// MPFR 2.4.0 Specifics
381#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
382 friend int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
383 friend const mpreal li2(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
384 friend const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd);
385 friend const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
386#endif
387
388// MPFR 3.0.0 Specifics
389#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
390 friend const mpreal digamma(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
391 friend const mpreal ai(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
392 friend const mpreal urandom (gmp_randstate_t& state,mp_rnd_t rnd_mode = mpreal::default_rnd); // use gmp_randinit_default() to init state, gmp_randclear() to clear
393 friend bool _isregular(const mpreal& v);
394#endif
395
396 // Uniformly distributed random number generation in [0,1] using
397 // Mersenne-Twister algorithm by default.
398 // Use parameter to setup seed, e.g.: random((unsigned)time(NULL))
399 // Check urandom() for more precise control.
400 friend const mpreal random(unsigned int seed = 0);
401
402 // Exponent and mantissa manipulation
403 friend const mpreal frexp(const mpreal& v, mp_exp_t* exp);
404 friend const mpreal ldexp(const mpreal& v, mp_exp_t exp);
405
406 // Splits mpreal value into fractional and integer parts.
407 // Returns fractional part and stores integer part in n.
408 friend const mpreal modf(const mpreal& v, mpreal& n);
409
410 // Constants
411 // don't forget to call mpfr_free_cache() for every thread where you are using const-functions
412 friend const mpreal const_log2 (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd);
413 friend const mpreal const_pi (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd);
414 friend const mpreal const_euler (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd);
415 friend const mpreal const_catalan (mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd);
416 // returns +inf iff sign>=0 otherwise -inf
417 friend const mpreal const_infinity(int sign = 1, mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd);
418
419 // Output/ Input
420 friend std::ostream& operator<<(std::ostream& os, const mpreal& v);
421 friend std::istream& operator>>(std::istream& is, mpreal& v);
422
423 // Integer Related Functions
424 friend const mpreal rint (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
425 friend const mpreal ceil (const mpreal& v);
426 friend const mpreal floor(const mpreal& v);
427 friend const mpreal round(const mpreal& v);
428 friend const mpreal trunc(const mpreal& v);
429 friend const mpreal rint_ceil (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
430 friend const mpreal rint_floor(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
431 friend const mpreal rint_round(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
432 friend const mpreal rint_trunc(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
433 friend const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
434 friend const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd);
435 friend const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd);
436
437 // Miscellaneous Functions
438 friend const mpreal nexttoward (const mpreal& x, const mpreal& y);
439 friend const mpreal nextabove (const mpreal& x);
440 friend const mpreal nextbelow (const mpreal& x);
441
442 // use gmp_randinit_default() to init state, gmp_randclear() to clear
443 friend const mpreal urandomb (gmp_randstate_t& state);
444
445// MPFR < 2.4.2 Specifics
446#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2))
447 friend const mpreal random2 (mp_size_t size, mp_exp_t exp);
448#endif
449
450 // Instance Checkers
451 friend bool isnan (const mpreal& v);
452 friend bool isinf (const mpreal& v);
453 friend bool isfinite(const mpreal& v);
454
455 friend bool _isnum(const mpreal& v);
456 friend bool _iszero(const mpreal& v);
457 friend bool _isint(const mpreal& v);
458
459 // Set/Get instance properties
460 inline mp_prec_t get_prec() const;
461 inline void set_prec(mp_prec_t prec, mp_rnd_t rnd_mode = default_rnd); // Change precision with rounding mode
462
463 // Aliases for get_prec(), set_prec() - needed for compatibility with std::complex<mpreal> interface
464 inline mpreal& setPrecision(int Precision, mp_rnd_t RoundingMode = (mpfr_get_default_rounding_mode)());
465 inline int getPrecision() const;
466
467 // Set mpreal to +/- inf, NaN, +/-0
468 mpreal& setInf (int Sign = +1);
469 mpreal& setNan ();
470 mpreal& setZero (int Sign = +1);
471 mpreal& setSign (int Sign, mp_rnd_t RoundingMode = (mpfr_get_default_rounding_mode)());
472
473 //Exponent
474 mp_exp_t get_exp();
475 int set_exp(mp_exp_t e);
476 int check_range (int t, mp_rnd_t rnd_mode = default_rnd);
477 int subnormalize (int t,mp_rnd_t rnd_mode = default_rnd);
478
479 // Inexact conversion from float
480 inline bool fits_in_bits(double x, int n);
481
482 // Set/Get global properties
483 static void set_default_prec(mp_prec_t prec);
484 static mp_prec_t get_default_prec();
485 static void set_default_base(int base);
486 static int get_default_base();
487 static void set_double_bits(int dbits);
488 static int get_double_bits();
489 static void set_default_rnd(mp_rnd_t rnd_mode);
490 static mp_rnd_t get_default_rnd();
491 static mp_exp_t get_emin (void);
492 static mp_exp_t get_emax (void);
493 static mp_exp_t get_emin_min (void);
494 static mp_exp_t get_emin_max (void);
495 static mp_exp_t get_emax_min (void);
496 static mp_exp_t get_emax_max (void);
497 static int set_emin (mp_exp_t exp);
498 static int set_emax (mp_exp_t exp);
499
500 // Efficient swapping of two mpreal values
501 friend void swap(mpreal& x, mpreal& y);
502
503 //Min Max - macros is evil. Needed for systems which defines max and min globally as macros (e.g. Windows)
504 //Hope that globally defined macros use > < operations only
505 friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = default_rnd);
506 friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = default_rnd);
507
508private:
509 // Optimized dynamic memory allocation/(re-)deallocation.
510 static bool is_custom_malloc;
511 static void *mpreal_allocate (size_t alloc_size);
512 static void *mpreal_reallocate (void *ptr, size_t old_size, size_t new_size);
513 static void mpreal_free (void *ptr, size_t size);
514 inline static void set_custom_malloc (void);
515
516#if defined(_MSC_VER) && defined(_DEBUG)
517 // Human friendly Debug Preview in Visual Studio.
518 // Put one of these lines:
519 //
520 // mpfr::mpreal=<DebugView> ; Show value only
521 // mpfr::mpreal=<DebugView>, <mp[0]._mpfr_prec,u>bits ; Show value & precision
522 //
523 // at the beginning of
524 // [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat
525 double DebugView;
526 inline void setDebugView() { DebugView = this->toDouble(); };
527#endif
528
529};
530
531//////////////////////////////////////////////////////////////////////////
532// Exceptions
533class conversion_overflow : public std::exception {
534public:
535 std::string why() { return "inexact conversion from floating point"; }
536};
537
538//////////////////////////////////////////////////////////////////////////
539// + Addition
540const mpreal operator+(const mpreal& a, const mpreal& b);
541
542// + Fast specialized addition - implemented through fast += operations
543const mpreal operator+(const mpreal& a, const mpz_t b);
544const mpreal operator+(const mpreal& a, const mpq_t b);
545const mpreal operator+(const mpreal& a, const long double b);
546const mpreal operator+(const mpreal& a, const double b);
547const mpreal operator+(const mpreal& a, const unsigned long int b);
548const mpreal operator+(const mpreal& a, const unsigned int b);
549const mpreal operator+(const mpreal& a, const long int b);
550const mpreal operator+(const mpreal& a, const int b);
551const mpreal operator+(const mpreal& a, const char* b);
552const mpreal operator+(const char* a, const mpreal& b);
553const std::string operator+(const mpreal& a, const std::string b);
554const std::string operator+(const std::string a, const mpreal& b);
555
556const mpreal operator+(const mpz_t b, const mpreal& a);
557const mpreal operator+(const mpq_t b, const mpreal& a);
558const mpreal operator+(const long double b, const mpreal& a);
559const mpreal operator+(const double b, const mpreal& a);
560const mpreal operator+(const unsigned long int b, const mpreal& a);
561const mpreal operator+(const unsigned int b, const mpreal& a);
562const mpreal operator+(const long int b, const mpreal& a);
563const mpreal operator+(const int b, const mpreal& a);
564
565#if defined (MPREAL_HAVE_INT64_SUPPORT)
566const mpreal operator+(const mpreal& a, const int64_t b);
567const mpreal operator+(const int64_t b, const mpreal& a);
568const mpreal operator+(const mpreal& a, const uint64_t b);
569const mpreal operator+(const uint64_t b, const mpreal& a);
570const mpreal operator-(const mpreal& a, const int64_t b);
571const mpreal operator-(const int64_t b, const mpreal& a);
572const mpreal operator-(const mpreal& a, const uint64_t b);
573const mpreal operator-(const uint64_t b, const mpreal& a);
574const mpreal operator*(const mpreal& a, const int64_t b);
575const mpreal operator*(const int64_t b, const mpreal& a);
576const mpreal operator*(const mpreal& a, const uint64_t b);
577const mpreal operator*(const uint64_t b, const mpreal& a);
578const mpreal operator/(const mpreal& a, const int64_t b);
579const mpreal operator/(const int64_t b, const mpreal& a);
580const mpreal operator/(const mpreal& a, const uint64_t b);
581const mpreal operator/(const uint64_t b, const mpreal& a);
582#endif
583
584//////////////////////////////////////////////////////////////////////////
585// - Subtraction
586const mpreal operator-(const mpreal& a, const mpreal& b);
587
588// - Fast specialized subtraction - implemented through fast -= operations
589const mpreal operator-(const mpreal& a, const mpz_t b);
590const mpreal operator-(const mpreal& a, const mpq_t b);
591const mpreal operator-(const mpreal& a, const long double b);
592const mpreal operator-(const mpreal& a, const double b);
593const mpreal operator-(const mpreal& a, const unsigned long int b);
594const mpreal operator-(const mpreal& a, const unsigned int b);
595const mpreal operator-(const mpreal& a, const long int b);
596const mpreal operator-(const mpreal& a, const int b);
597const mpreal operator-(const mpreal& a, const char* b);
598const mpreal operator-(const char* a, const mpreal& b);
599
600const mpreal operator-(const mpz_t b, const mpreal& a);
601const mpreal operator-(const mpq_t b, const mpreal& a);
602const mpreal operator-(const long double b, const mpreal& a);
603//const mpreal operator-(const double b, const mpreal& a);
604
605//////////////////////////////////////////////////////////////////////////
606// * Multiplication
607const mpreal operator*(const mpreal& a, const mpreal& b);
608
609// * Fast specialized multiplication - implemented through fast *= operations
610const mpreal operator*(const mpreal& a, const mpz_t b);
611const mpreal operator*(const mpreal& a, const mpq_t b);
612const mpreal operator*(const mpreal& a, const long double b);
613const mpreal operator*(const mpreal& a, const double b);
614const mpreal operator*(const mpreal& a, const unsigned long int b);
615const mpreal operator*(const mpreal& a, const unsigned int b);
616const mpreal operator*(const mpreal& a, const long int b);
617const mpreal operator*(const mpreal& a, const int b);
618
619const mpreal operator*(const mpz_t b, const mpreal& a);
620const mpreal operator*(const mpq_t b, const mpreal& a);
621const mpreal operator*(const long double b, const mpreal& a);
622const mpreal operator*(const double b, const mpreal& a);
623const mpreal operator*(const unsigned long int b, const mpreal& a);
624const mpreal operator*(const unsigned int b, const mpreal& a);
625const mpreal operator*(const long int b, const mpreal& a);
626const mpreal operator*(const int b, const mpreal& a);
627
628//////////////////////////////////////////////////////////////////////////
629// / Division
630const mpreal operator/(const mpreal& a, const mpreal& b);
631
632// / Fast specialized division - implemented through fast /= operations
633const mpreal operator/(const mpreal& a, const mpz_t b);
634const mpreal operator/(const mpreal& a, const mpq_t b);
635const mpreal operator/(const mpreal& a, const long double b);
636const mpreal operator/(const mpreal& a, const double b);
637const mpreal operator/(const mpreal& a, const unsigned long int b);
638const mpreal operator/(const mpreal& a, const unsigned int b);
639const mpreal operator/(const mpreal& a, const long int b);
640const mpreal operator/(const mpreal& a, const int b);
641
642const mpreal operator/(const long double b, const mpreal& a);
643
644//////////////////////////////////////////////////////////////////////////
645// Shifts operators - Multiplication/Division by a power of 2
646const mpreal operator<<(const mpreal& v, const unsigned long int k);
647const mpreal operator<<(const mpreal& v, const unsigned int k);
648const mpreal operator<<(const mpreal& v, const long int k);
649const mpreal operator<<(const mpreal& v, const int k);
650
651const mpreal operator>>(const mpreal& v, const unsigned long int k);
652const mpreal operator>>(const mpreal& v, const unsigned int k);
653const mpreal operator>>(const mpreal& v, const long int k);
654const mpreal operator>>(const mpreal& v, const int k);
655
656//////////////////////////////////////////////////////////////////////////
657// Boolean operators
658bool operator < (const mpreal& a, const unsigned long int b);
659bool operator < (const mpreal& a, const unsigned int b);
660bool operator < (const mpreal& a, const long int b);
661bool operator < (const mpreal& a, const int b);
662bool operator < (const mpreal& a, const long double b);
663bool operator < (const mpreal& a, const double b);
664
665bool operator < (const unsigned long int a,const mpreal& b);
666bool operator < (const unsigned int a, const mpreal& b);
667bool operator < (const long int a, const mpreal& b);
668bool operator < (const int a, const mpreal& b);
669bool operator < (const long double a, const mpreal& b);
670bool operator < (const double a, const mpreal& b);
671
672bool operator > (const mpreal& a, const unsigned long int b);
673bool operator > (const mpreal& a, const unsigned int b);
674bool operator > (const mpreal& a, const long int b);
675bool operator > (const mpreal& a, const int b);
676bool operator > (const mpreal& a, const long double b);
677bool operator > (const mpreal& a, const double b);
678
679bool operator > (const unsigned long int a,const mpreal& b);
680bool operator > (const unsigned int a, const mpreal& b);
681bool operator > (const long int a, const mpreal& b);
682bool operator > (const int a, const mpreal& b);
683bool operator > (const long double a, const mpreal& b);
684bool operator > (const double a, const mpreal& b);
685
686bool operator >= (const mpreal& a, const unsigned long int b);
687bool operator >= (const mpreal& a, const unsigned int b);
688bool operator >= (const mpreal& a, const long int b);
689bool operator >= (const mpreal& a, const int b);
690bool operator >= (const mpreal& a, const long double b);
691bool operator >= (const mpreal& a, const double b);
692
693bool operator >= (const unsigned long int a,const mpreal& b);
694bool operator >= (const unsigned int a, const mpreal& b);
695bool operator >= (const long int a, const mpreal& b);
696bool operator >= (const int a, const mpreal& b);
697bool operator >= (const long double a, const mpreal& b);
698bool operator >= (const double a, const mpreal& b);
699
700bool operator <= (const mpreal& a, const unsigned long int b);
701bool operator <= (const mpreal& a, const unsigned int b);
702bool operator <= (const mpreal& a, const long int b);
703bool operator <= (const mpreal& a, const int b);
704bool operator <= (const mpreal& a, const long double b);
705bool operator <= (const mpreal& a, const double b);
706
707bool operator <= (const unsigned long int a,const mpreal& b);
708bool operator <= (const unsigned int a, const mpreal& b);
709bool operator <= (const long int a, const mpreal& b);
710bool operator <= (const int a, const mpreal& b);
711bool operator <= (const long double a, const mpreal& b);
712bool operator <= (const double a, const mpreal& b);
713
714bool operator == (const unsigned long int a,const mpreal& b);
715bool operator == (const unsigned int a, const mpreal& b);
716bool operator == (const long int a, const mpreal& b);
717bool operator == (const int a, const mpreal& b);
718bool operator == (const long double a, const mpreal& b);
719bool operator == (const double a, const mpreal& b);
720
721bool operator != (const mpreal& a, const unsigned long int b);
722bool operator != (const mpreal& a, const unsigned int b);
723bool operator != (const mpreal& a, const long int b);
724bool operator != (const mpreal& a, const int b);
725bool operator != (const mpreal& a, const long double b);
726bool operator != (const mpreal& a, const double b);
727
728bool operator != (const unsigned long int a,const mpreal& b);
729bool operator != (const unsigned int a, const mpreal& b);
730bool operator != (const long int a, const mpreal& b);
731bool operator != (const int a, const mpreal& b);
732bool operator != (const long double a, const mpreal& b);
733bool operator != (const double a, const mpreal& b);
734
735//////////////////////////////////////////////////////////////////////////
736// sqrt
737const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode = mpreal::default_rnd);
738const mpreal sqrt(const long int v, mp_rnd_t rnd_mode = mpreal::default_rnd);
739const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::default_rnd);
740const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::default_rnd);
741const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::default_rnd);
742
743//////////////////////////////////////////////////////////////////////////
744// pow
745const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
746const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
747const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
748const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
749
750const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
751const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
752const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
753const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
754const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
755
756const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
757const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
758const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
759const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
760const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
761
762const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
763const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
764const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
765const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
766const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
767const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
768
769const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
770const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
771const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
772const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
773const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
774const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
775
776const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
777const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
778const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
779const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
780const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
781const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
782
783const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
784const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
785const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
786const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
787const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
788
789const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::default_rnd);
790const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
791const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
792const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
793const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
794
795//////////////////////////////////////////////////////////////////////////
796// Estimate machine epsilon for the given precision
797// Returns smallest eps such that 1.0 + eps != 1.0
798inline const mpreal machine_epsilon(mp_prec_t prec = mpreal::get_default_prec());
799
800// Returns the positive distance from abs(x) to the next larger in magnitude floating point number of the same precision as x
801inline const mpreal machine_epsilon(const mpreal& x);
802
803inline const mpreal mpreal_min(mp_prec_t prec = mpreal::get_default_prec());
804inline const mpreal mpreal_max(mp_prec_t prec = mpreal::get_default_prec());
805
806//////////////////////////////////////////////////////////////////////////
807// Bits - decimal digits relation
808// bits = ceil(digits*log[2](10))
809// digits = floor(bits*log[10](2))
810
811inline mp_prec_t digits2bits(int d);
812inline int bits2digits(mp_prec_t b);
813
814//////////////////////////////////////////////////////////////////////////
815// min, max
816//const mpreal max(const mpreal& x, const mpreal& y);
817//const mpreal min(const mpreal& x, const mpreal& y);
818
819//////////////////////////////////////////////////////////////////////////
820// Implementation
821//////////////////////////////////////////////////////////////////////////
822
823//////////////////////////////////////////////////////////////////////////
824// Operators - Assignment
825inline mpreal& mpreal::operator=(const mpreal& v)
826{
827 if (this != &v)
828 {
829 mpfr_clear(mp);
830 mpfr_init2(mp,mpfr_get_prec(v.mp));
831 mpfr_set(mp,v.mp,default_rnd);
832
833#if defined(_MSC_VER) && defined(_DEBUG)
834 setDebugView();
835#endif
836
837 }
838 return *this;
839}
840
841inline mpreal& mpreal::operator=(const mpf_t v)
842{
843 mpfr_set_f(mp,v,default_rnd);
844
845#if defined(_MSC_VER) && defined(_DEBUG)
846 setDebugView();
847#endif
848
849 return *this;
850}
851
852inline mpreal& mpreal::operator=(const mpz_t v)
853{
854 mpfr_set_z(mp,v,default_rnd);
855
856#if defined(_MSC_VER) && defined(_DEBUG)
857 setDebugView();
858#endif
859
860 return *this;
861}
862
863inline mpreal& mpreal::operator=(const mpq_t v)
864{
865 mpfr_set_q(mp,v,default_rnd);
866
867#if defined(_MSC_VER) && defined(_DEBUG)
868 setDebugView();
869#endif
870
871 return *this;
872}
873
874inline mpreal& mpreal::operator=(const long double v)
875{
876 mpfr_set_ld(mp,v,default_rnd);
877
878#if defined(_MSC_VER) && defined(_DEBUG)
879 setDebugView();
880#endif
881
882 return *this;
883}
884
885inline mpreal& mpreal::operator=(const double v)
886{
887 if(double_bits == -1 || fits_in_bits(v, double_bits))
888 {
889 mpfr_set_d(mp,v,default_rnd);
890
891#if defined(_MSC_VER) && defined(_DEBUG)
892 setDebugView();
893#endif
894
895 }
896 else
897 throw conversion_overflow();
898
899 return *this;
900}
901
902inline mpreal& mpreal::operator=(const unsigned long int v)
903{
904 mpfr_set_ui(mp,v,default_rnd);
905
906#if defined(_MSC_VER) && defined(_DEBUG)
907 setDebugView();
908#endif
909
910 return *this;
911}
912
913inline mpreal& mpreal::operator=(const unsigned int v)
914{
915 mpfr_set_ui(mp,v,default_rnd);
916
917#if defined(_MSC_VER) && defined(_DEBUG)
918 setDebugView();
919#endif
920
921 return *this;
922}
923
924inline mpreal& mpreal::operator=(const long int v)
925{
926 mpfr_set_si(mp,v,default_rnd);
927
928#if defined(_MSC_VER) && defined(_DEBUG)
929 setDebugView();
930#endif
931
932 return *this;
933}
934
935inline mpreal& mpreal::operator=(const int v)
936{
937 mpfr_set_si(mp,v,default_rnd);
938
939#if defined(_MSC_VER) && defined(_DEBUG)
940 setDebugView();
941#endif
942
943 return *this;
944}
945
946//////////////////////////////////////////////////////////////////////////
947// + Addition
948inline mpreal& mpreal::operator+=(const mpreal& v)
949{
950 mpfr_add(mp,mp,v.mp,default_rnd);
951 return *this;
952}
953
954inline mpreal& mpreal::operator+=(const mpf_t u)
955{
956 *this += mpreal(u);
957 return *this;
958}
959
960inline mpreal& mpreal::operator+=(const mpz_t u)
961{
962 mpfr_add_z(mp,mp,u,default_rnd);
963 return *this;
964}
965
966inline mpreal& mpreal::operator+=(const mpq_t u)
967{
968 mpfr_add_q(mp,mp,u,default_rnd);
969 return *this;
970}
971
972inline mpreal& mpreal::operator+= (const long double u)
973{
974 return *this += mpreal(u);
975}
976
977inline mpreal& mpreal::operator+= (const double u)
978{
979#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
980 mpfr_add_d(mp,mp,u,default_rnd);
981 return *this;
982#else
983 return *this += mpreal(u);
984#endif
985}
986
987inline mpreal& mpreal::operator+=(const unsigned long int u)
988{
989 mpfr_add_ui(mp,mp,u,default_rnd);
990 return *this;
991}
992
993inline mpreal& mpreal::operator+=(const unsigned int u)
994{
995 mpfr_add_ui(mp,mp,u,default_rnd);
996 return *this;
997}
998
999inline mpreal& mpreal::operator+=(const long int u)
1000{
1001 mpfr_add_si(mp,mp,u,default_rnd);
1002 return *this;
1003}
1004
1005inline mpreal& mpreal::operator+=(const int u)
1006{
1007 mpfr_add_si(mp,mp,u,default_rnd);
1008 return *this;
1009}
1010
1011#if defined (MPREAL_HAVE_INT64_SUPPORT)
1012inline mpreal& mpreal::operator+=(const int64_t u){ return *this += mpreal(u);}
1013inline mpreal& mpreal::operator+=(const uint64_t u){ return *this += mpreal(u);}
1014inline mpreal& mpreal::operator-=(const int64_t u){ return *this -= mpreal(u);}
1015inline mpreal& mpreal::operator-=(const uint64_t u){ return *this -= mpreal(u);}
1016inline mpreal& mpreal::operator*=(const int64_t u){ return *this *= mpreal(u);}
1017inline mpreal& mpreal::operator*=(const uint64_t u){ return *this *= mpreal(u);}
1018inline mpreal& mpreal::operator/=(const int64_t u){ return *this /= mpreal(u);}
1019inline mpreal& mpreal::operator/=(const uint64_t u){ return *this /= mpreal(u);}
1020#endif
1021
1022inline const mpreal mpreal::operator+()const
1023{
1024 return mpreal(*this);
1025}
1026
1027inline const mpreal operator+(const mpreal& a, const mpreal& b)
1028{
1029 // prec(a+b) = max(prec(a),prec(b))
1030 if(a.get_prec()>b.get_prec()) return mpreal(a) += b;
1031 else return mpreal(b) += a;
1032}
1033
1034inline const std::string operator+(const mpreal& a, const std::string b)
1035{
1036 return a.toString() + b;
1037}
1038
1039inline const std::string operator+(const std::string a, const mpreal& b)
1040{
1041 return a + b.toString();
1042}
1043
1044inline const mpreal operator+(const mpreal& a, const mpz_t b)
1045{
1046 return mpreal(a) += b;
1047}
1048
1049inline const mpreal operator+(const mpreal& a, const char* b)
1050{
1051 return mpreal(a) += mpreal(b);
1052}
1053
1054inline const mpreal operator+(const char* a, const mpreal& b)
1055{
1056 return mpreal(a) += b;
1057
1058}
1059
1060inline const mpreal operator+(const mpreal& a, const mpq_t b)
1061{
1062 return mpreal(a) += b;
1063}
1064
1065inline const mpreal operator+(const mpreal& a, const long double b)
1066{
1067 return mpreal(a) += b;
1068}
1069
1070inline const mpreal operator+(const mpreal& a, const double b)
1071{
1072 return mpreal(a) += b;
1073}
1074
1075inline const mpreal operator+(const mpreal& a, const unsigned long int b)
1076{
1077 return mpreal(a) += b;
1078}
1079
1080inline const mpreal operator+(const mpreal& a, const unsigned int b)
1081{
1082 return mpreal(a) += b;
1083}
1084
1085inline const mpreal operator+(const mpreal& a, const long int b)
1086{
1087 return mpreal(a) += b;
1088}
1089
1090inline const mpreal operator+(const mpreal& a, const int b)
1091{
1092 return mpreal(a) += b;
1093}
1094
1095inline const mpreal operator+(const mpz_t b, const mpreal& a)
1096{
1097 return mpreal(a) += b;
1098}
1099
1100inline const mpreal operator+(const mpq_t b, const mpreal& a)
1101{
1102 return mpreal(a) += b;
1103}
1104
1105inline const mpreal operator+(const long double b, const mpreal& a)
1106{
1107 return mpreal(a) += b;
1108}
1109
1110inline const mpreal operator+(const double b, const mpreal& a)
1111{
1112 return mpreal(a) += b;
1113}
1114
1115inline const mpreal operator+(const unsigned long int b, const mpreal& a)
1116{
1117 return mpreal(a) += b;
1118}
1119
1120inline const mpreal operator+(const unsigned int b, const mpreal& a)
1121{
1122 return mpreal(a) += b;
1123}
1124
1125inline const mpreal operator+(const long int b, const mpreal& a)
1126{
1127 return mpreal(a) += b;
1128}
1129
1130inline const mpreal operator+(const int b, const mpreal& a)
1131{
1132 return mpreal(a) += b;
1133}
1134
1135#if defined (MPREAL_HAVE_INT64_SUPPORT)
1136
1137inline const mpreal operator+(const mpreal& a, const int64_t b){ return mpreal(a) += b;}
1138inline const mpreal operator+(const int64_t b, const mpreal& a){ return mpreal(a) += b;}
1139inline const mpreal operator+(const mpreal& a, const uint64_t b){ return mpreal(a) += b;}
1140inline const mpreal operator+(const uint64_t b, const mpreal& a){ return mpreal(a) += b;}
1141
1142inline const mpreal operator-(const mpreal& a, const int64_t b){ return mpreal(a) -= b;}
1143inline const mpreal operator-(const int64_t b, const mpreal& a){ return mpreal(b) -= a;}
1144inline const mpreal operator-(const mpreal& a, const uint64_t b){ return mpreal(a) -= b;}
1145inline const mpreal operator-(const uint64_t b, const mpreal& a){ return mpreal(b) -= a;}
1146
1147inline const mpreal operator*(const mpreal& a, const int64_t b){ return mpreal(a) *= b;}
1148inline const mpreal operator*(const int64_t b, const mpreal& a){ return mpreal(a) *= b;}
1149inline const mpreal operator*(const mpreal& a, const uint64_t b){ return mpreal(a) *= b;}
1150inline const mpreal operator*(const uint64_t b, const mpreal& a){ return mpreal(a) *= b;}
1151
1152inline const mpreal operator/(const mpreal& a, const int64_t b){ return mpreal(a) /= b;}
1153inline const mpreal operator/(const int64_t b, const mpreal& a){ return mpreal(b) /= a;}
1154inline const mpreal operator/(const mpreal& a, const uint64_t b){ return mpreal(a) /= b;}
1155inline const mpreal operator/(const uint64_t b, const mpreal& a){ return mpreal(b) /= a;}
1156
1157#endif
1158
1159inline mpreal& mpreal::operator++()
1160{
1161 *this += 1;
1162 return *this;
1163}
1164
1165inline const mpreal mpreal::operator++ (int)
1166{
1167 mpreal x(*this);
1168 *this += 1;
1169 return x;
1170}
1171
1172inline mpreal& mpreal::operator--()
1173{
1174 *this -= 1;
1175 return *this;
1176}
1177
1178inline const mpreal mpreal::operator-- (int)
1179{
1180 mpreal x(*this);
1181 *this -= 1;
1182 return x;
1183}
1184
1185//////////////////////////////////////////////////////////////////////////
1186// - Subtraction
1187inline mpreal& mpreal::operator-= (const mpreal& v)
1188{
1189 mpfr_sub(mp,mp,v.mp,default_rnd);
1190 return *this;
1191}
1192
1193inline mpreal& mpreal::operator-=(const mpz_t v)
1194{
1195 mpfr_sub_z(mp,mp,v,default_rnd);
1196 return *this;
1197}
1198
1199inline mpreal& mpreal::operator-=(const mpq_t v)
1200{
1201 mpfr_sub_q(mp,mp,v,default_rnd);
1202 return *this;
1203}
1204
1205inline mpreal& mpreal::operator-=(const long double v)
1206{
1207 return *this -= mpreal(v);
1208}
1209
1210inline mpreal& mpreal::operator-=(const double v)
1211{
1212#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
1213 mpfr_sub_d(mp,mp,v,default_rnd);
1214 return *this;
1215#else
1216 return *this -= mpreal(v);
1217#endif
1218}
1219
1220inline mpreal& mpreal::operator-=(const unsigned long int v)
1221{
1222 mpfr_sub_ui(mp,mp,v,default_rnd);
1223 return *this;
1224}
1225
1226inline mpreal& mpreal::operator-=(const unsigned int v)
1227{
1228 mpfr_sub_ui(mp,mp,v,default_rnd);
1229 return *this;
1230}
1231
1232inline mpreal& mpreal::operator-=(const long int v)
1233{
1234 mpfr_sub_si(mp,mp,v,default_rnd);
1235 return *this;
1236}
1237
1238inline mpreal& mpreal::operator-=(const int v)
1239{
1240 mpfr_sub_si(mp,mp,v,default_rnd);
1241 return *this;
1242}
1243
1244inline const mpreal mpreal::operator-()const
1245{
1246 mpreal u(*this);
1247 mpfr_neg(u.mp,u.mp,default_rnd);
1248 return u;
1249}
1250
1251inline const mpreal operator-(const mpreal& a, const mpreal& b)
1252{
1253 // prec(a-b) = max(prec(a),prec(b))
1254 if(a.getPrecision() >= b.getPrecision())
1255 {
1256
1257 return mpreal(a) -= b;
1258
1259 }else{
1260
1261 mpreal x(a);
1262
1263 x.setPrecision(b.getPrecision());
1264
1265 return x -= b;
1266 }
1267}
1268
1269inline const mpreal operator-(const mpreal& a, const mpz_t b)
1270{
1271 return mpreal(a) -= b;
1272}
1273
1274inline const mpreal operator-(const mpreal& a, const mpq_t b)
1275{
1276 return mpreal(a) -= b;
1277}
1278
1279inline const mpreal operator-(const mpreal& a, const long double b)
1280{
1281 return mpreal(a) -= b;
1282}
1283
1284inline const mpreal operator-(const mpreal& a, const double b)
1285{
1286 return mpreal(a) -= b;
1287}
1288
1289inline const mpreal operator-(const mpreal& a, const unsigned long int b)
1290{
1291 return mpreal(a) -= b;
1292}
1293
1294inline const mpreal operator-(const mpreal& a, const unsigned int b)
1295{
1296 return mpreal(a) -= b;
1297}
1298
1299inline const mpreal operator-(const mpreal& a, const long int b)
1300{
1301 return mpreal(a) -= b;
1302}
1303
1304inline const mpreal operator-(const mpreal& a, const int b)
1305{
1306 return mpreal(a) -= b;
1307}
1308
1309inline const mpreal operator-(const mpz_t b, const mpreal& a)
1310{
1311 return mpreal(b) -= a;
1312}
1313
1314inline const mpreal operator-(const mpq_t b, const mpreal& a)
1315{
1316 return mpreal(b) -= a;
1317}
1318
1319inline const mpreal operator-(const long double b, const mpreal& a)
1320{
1321 return mpreal(b) -= a;
1322}
1323
1324inline const mpreal operator-(const double b, const mpreal& a)
1325{
1326#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
1327 mpreal x(a);
1328 mpfr_d_sub(x.mp,b,a.mp,mpreal::default_rnd);
1329 return x;
1330#else
1331 return mpreal(b) -= a;
1332#endif
1333}
1334
1335inline const mpreal operator-(const unsigned long int b, const mpreal& a)
1336{
1337 mpreal x(a);
1338 mpfr_ui_sub(x.mp,b,a.mp,mpreal::default_rnd);
1339 return x;
1340}
1341
1342inline const mpreal operator-(const unsigned int b, const mpreal& a)
1343{
1344 mpreal x(a);
1345 mpfr_ui_sub(x.mp,b,a.mp,mpreal::default_rnd);
1346 return x;
1347}
1348
1349inline const mpreal operator-(const long int b, const mpreal& a)
1350{
1351 mpreal x(a);
1352 mpfr_si_sub(x.mp,b,a.mp,mpreal::default_rnd);
1353 return x;
1354}
1355
1356inline const mpreal operator-(const int b, const mpreal& a)
1357{
1358 mpreal x(a);
1359 mpfr_si_sub(x.mp,b,a.mp,mpreal::default_rnd);
1360 return x;
1361}
1362
1363inline const mpreal operator-(const mpreal& a, const char* b)
1364{
1365 return a - mpreal(b);
1366}
1367
1368inline const mpreal operator-(const char* a, const mpreal& b)
1369{
1370 return mpreal(a) - b;
1371}
1372
1373//////////////////////////////////////////////////////////////////////////
1374// * Multiplication
1375inline mpreal& mpreal::operator*= (const mpreal& v)
1376{
1377 mpfr_mul(mp,mp,v.mp,default_rnd);
1378 return *this;
1379}
1380
1381inline mpreal& mpreal::operator*=(const mpz_t v)
1382{
1383 mpfr_mul_z(mp,mp,v,default_rnd);
1384 return *this;
1385}
1386
1387inline mpreal& mpreal::operator*=(const mpq_t v)
1388{
1389 mpfr_mul_q(mp,mp,v,default_rnd);
1390 return *this;
1391}
1392
1393inline mpreal& mpreal::operator*=(const long double v)
1394{
1395 return *this *= mpreal(v);
1396}
1397
1398inline mpreal& mpreal::operator*=(const double v)
1399{
1400#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
1401 mpfr_mul_d(mp,mp,v,default_rnd);
1402 return *this;
1403#else
1404 return *this *= mpreal(v);
1405#endif
1406}
1407
1408inline mpreal& mpreal::operator*=(const unsigned long int v)
1409{
1410 mpfr_mul_ui(mp,mp,v,default_rnd);
1411 return *this;
1412}
1413
1414inline mpreal& mpreal::operator*=(const unsigned int v)
1415{
1416 mpfr_mul_ui(mp,mp,v,default_rnd);
1417 return *this;
1418}
1419
1420inline mpreal& mpreal::operator*=(const long int v)
1421{
1422 mpfr_mul_si(mp,mp,v,default_rnd);
1423 return *this;
1424}
1425
1426inline mpreal& mpreal::operator*=(const int v)
1427{
1428 mpfr_mul_si(mp,mp,v,default_rnd);
1429 return *this;
1430}
1431
1432inline const mpreal operator*(const mpreal& a, const mpreal& b)
1433{
1434 // prec(a*b) = max(prec(a),prec(b))
1435 if(a.getPrecision() >= b.getPrecision()) return mpreal(a) *= b;
1436 else return mpreal(b) *= a;
1437}
1438
1439inline const mpreal operator*(const mpreal& a, const mpz_t b)
1440{
1441 return mpreal(a) *= b;
1442}
1443
1444inline const mpreal operator*(const mpreal& a, const mpq_t b)
1445{
1446 return mpreal(a) *= b;
1447}
1448
1449inline const mpreal operator*(const mpreal& a, const long double b)
1450{
1451 return mpreal(a) *= b;
1452}
1453
1454inline const mpreal operator*(const mpreal& a, const double b)
1455{
1456 return mpreal(a) *= b;
1457}
1458
1459inline const mpreal operator*(const mpreal& a, const unsigned long int b)
1460{
1461 return mpreal(a) *= b;
1462}
1463
1464inline const mpreal operator*(const mpreal& a, const unsigned int b)
1465{
1466 return mpreal(a) *= b;
1467}
1468
1469inline const mpreal operator*(const mpreal& a, const long int b)
1470{
1471 return mpreal(a) *= b;
1472}
1473
1474inline const mpreal operator*(const mpreal& a, const int b)
1475{
1476 return mpreal(a) *= b;
1477}
1478
1479inline const mpreal operator*(const mpz_t b, const mpreal& a)
1480{
1481 return mpreal(a) *= b;
1482}
1483
1484inline const mpreal operator*(const mpq_t b, const mpreal& a)
1485{
1486 return mpreal(a) *= b;
1487}
1488
1489inline const mpreal operator*(const long double b, const mpreal& a)
1490{
1491 return mpreal(a) *= b;
1492}
1493
1494inline const mpreal operator*(const double b, const mpreal& a)
1495{
1496 return mpreal(a) *= b;
1497}
1498
1499inline const mpreal operator*(const unsigned long int b, const mpreal& a)
1500{
1501 return mpreal(a) *= b;
1502}
1503
1504inline const mpreal operator*(const unsigned int b, const mpreal& a)
1505{
1506 return mpreal(a) *= b;
1507}
1508
1509inline const mpreal operator*(const long int b, const mpreal& a)
1510{
1511 return mpreal(a) *= b;
1512}
1513
1514inline const mpreal operator*(const int b, const mpreal& a)
1515{
1516 return mpreal(a) *= b;
1517}
1518
1519//////////////////////////////////////////////////////////////////////////
1520// / Division
1521inline mpreal& mpreal::operator/=(const mpreal& v)
1522{
1523 mpfr_div(mp,mp,v.mp,default_rnd);
1524 return *this;
1525}
1526
1527inline mpreal& mpreal::operator/=(const mpz_t v)
1528{
1529 mpfr_div_z(mp,mp,v,default_rnd);
1530 return *this;
1531}
1532
1533inline mpreal& mpreal::operator/=(const mpq_t v)
1534{
1535 mpfr_div_q(mp,mp,v,default_rnd);
1536 return *this;
1537}
1538
1539inline mpreal& mpreal::operator/=(const long double v)
1540{
1541 return *this /= mpreal(v);
1542}
1543
1544inline mpreal& mpreal::operator/=(const double v)
1545{
1546#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
1547 mpfr_div_d(mp,mp,v,default_rnd);
1548 return *this;
1549#else
1550 return *this /= mpreal(v);
1551#endif
1552}
1553
1554inline mpreal& mpreal::operator/=(const unsigned long int v)
1555{
1556 mpfr_div_ui(mp,mp,v,default_rnd);
1557 return *this;
1558}
1559
1560inline mpreal& mpreal::operator/=(const unsigned int v)
1561{
1562 mpfr_div_ui(mp,mp,v,default_rnd);
1563 return *this;
1564}
1565
1566inline mpreal& mpreal::operator/=(const long int v)
1567{
1568 mpfr_div_si(mp,mp,v,default_rnd);
1569 return *this;
1570}
1571
1572inline mpreal& mpreal::operator/=(const int v)
1573{
1574 mpfr_div_si(mp,mp,v,default_rnd);
1575 return *this;
1576}
1577
1578inline const mpreal operator/(const mpreal& a, const mpreal& b)
1579{
1580 // prec(a/b) = max(prec(a),prec(b))
1581 if(a.getPrecision() >= b.getPrecision())
1582 {
1583 return mpreal(a) /= b;
1584 }else{
1585
1586 mpreal x(a);
1587 x.setPrecision(b.getPrecision());
1588 return x /= b;
1589 }
1590}
1591
1592inline const mpreal operator/(const mpreal& a, const mpz_t b)
1593{
1594 return mpreal(a) /= b;
1595}
1596
1597inline const mpreal operator/(const mpreal& a, const mpq_t b)
1598{
1599 return mpreal(a) /= b;
1600}
1601
1602inline const mpreal operator/(const mpreal& a, const long double b)
1603{
1604 return mpreal(a) /= b;
1605}
1606
1607inline const mpreal operator/(const mpreal& a, const double b)
1608{
1609 return mpreal(a) /= b;
1610}
1611
1612inline const mpreal operator/(const mpreal& a, const unsigned long int b)
1613{
1614 return mpreal(a) /= b;
1615}
1616
1617inline const mpreal operator/(const mpreal& a, const unsigned int b)
1618{
1619 return mpreal(a) /= b;
1620}
1621
1622inline const mpreal operator/(const mpreal& a, const long int b)
1623{
1624 return mpreal(a) /= b;
1625}
1626
1627inline const mpreal operator/(const mpreal& a, const int b)
1628{
1629 return mpreal(a) /= b;
1630}
1631
1632inline const mpreal operator/(const unsigned long int b, const mpreal& a)
1633{
1634 mpreal x(a);
1635 mpfr_ui_div(x.mp,b,a.mp,mpreal::default_rnd);
1636 return x;
1637}
1638
1639inline const mpreal operator/(const unsigned int b, const mpreal& a)
1640{
1641 mpreal x(a);
1642 mpfr_ui_div(x.mp,b,a.mp,mpreal::default_rnd);
1643 return x;
1644}
1645
1646inline const mpreal operator/(const long int b, const mpreal& a)
1647{
1648 mpreal x(a);
1649 mpfr_si_div(x.mp,b,a.mp,mpreal::default_rnd);
1650 return x;
1651}
1652
1653inline const mpreal operator/(const int b, const mpreal& a)
1654{
1655 mpreal x(a);
1656 mpfr_si_div(x.mp,b,a.mp,mpreal::default_rnd);
1657 return x;
1658}
1659
1660inline const mpreal operator/(const long double b, const mpreal& a)
1661{
1662 return mpreal(b) /= a;
1663}
1664
1665inline const mpreal operator/(const double b, const mpreal& a)
1666{
1667#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
1668 mpreal x(a);
1669 mpfr_d_div(x.mp,b,a.mp,mpreal::default_rnd);
1670 return x;
1671#else
1672 return mpreal(b) /= a;
1673#endif
1674}
1675
1676//////////////////////////////////////////////////////////////////////////
1677// Shifts operators - Multiplication/Division by power of 2
1678inline mpreal& mpreal::operator<<=(const unsigned long int u)
1679{
1680 mpfr_mul_2ui(mp,mp,u,default_rnd);
1681 return *this;
1682}
1683
1684inline mpreal& mpreal::operator<<=(const unsigned int u)
1685{
1686 mpfr_mul_2ui(mp,mp,static_cast<unsigned long int>(u),default_rnd);
1687 return *this;
1688}
1689
1690inline mpreal& mpreal::operator<<=(const long int u)
1691{
1692 mpfr_mul_2si(mp,mp,u,default_rnd);
1693 return *this;
1694}
1695
1696inline mpreal& mpreal::operator<<=(const int u)
1697{
1698 mpfr_mul_2si(mp,mp,static_cast<long int>(u),default_rnd);
1699 return *this;
1700}
1701
1702inline mpreal& mpreal::operator>>=(const unsigned long int u)
1703{
1704 mpfr_div_2ui(mp,mp,u,default_rnd);
1705 return *this;
1706}
1707
1708inline mpreal& mpreal::operator>>=(const unsigned int u)
1709{
1710 mpfr_div_2ui(mp,mp,static_cast<unsigned long int>(u),default_rnd);
1711 return *this;
1712}
1713
1714inline mpreal& mpreal::operator>>=(const long int u)
1715{
1716 mpfr_div_2si(mp,mp,u,default_rnd);
1717 return *this;
1718}
1719
1720inline mpreal& mpreal::operator>>=(const int u)
1721{
1722 mpfr_div_2si(mp,mp,static_cast<long int>(u),default_rnd);
1723 return *this;
1724}
1725
1726inline const mpreal operator<<(const mpreal& v, const unsigned long int k)
1727{
1728 return mul_2ui(v,k);
1729}
1730
1731inline const mpreal operator<<(const mpreal& v, const unsigned int k)
1732{
1733 return mul_2ui(v,static_cast<unsigned long int>(k));
1734}
1735
1736inline const mpreal operator<<(const mpreal& v, const long int k)
1737{
1738 return mul_2si(v,k);
1739}
1740
1741inline const mpreal operator<<(const mpreal& v, const int k)
1742{
1743 return mul_2si(v,static_cast<long int>(k));
1744}
1745
1746inline const mpreal operator>>(const mpreal& v, const unsigned long int k)
1747{
1748 return div_2ui(v,k);
1749}
1750
1751inline const mpreal operator>>(const mpreal& v, const long int k)
1752{
1753 return div_2si(v,k);
1754}
1755
1756inline const mpreal operator>>(const mpreal& v, const unsigned int k)
1757{
1758 return div_2ui(v,static_cast<unsigned long int>(k));
1759}
1760
1761inline const mpreal operator>>(const mpreal& v, const int k)
1762{
1763 return div_2si(v,static_cast<long int>(k));
1764}
1765
1766// mul_2ui
1767inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)
1768{
1769 mpreal x(v);
1770 mpfr_mul_2ui(x.mp,v.mp,k,rnd_mode);
1771 return x;
1772}
1773
1774// mul_2si
1775inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)
1776{
1777 mpreal x(v);
1778 mpfr_mul_2si(x.mp,v.mp,k,rnd_mode);
1779 return x;
1780}
1781
1782inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)
1783{
1784 mpreal x(v);
1785 mpfr_div_2ui(x.mp,v.mp,k,rnd_mode);
1786 return x;
1787}
1788
1789inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode)
1790{
1791 mpreal x(v);
1792 mpfr_div_2si(x.mp,v.mp,k,rnd_mode);
1793 return x;
1794}
1795
1796//////////////////////////////////////////////////////////////////////////
1797//Boolean operators
1798inline bool operator > (const mpreal& a, const mpreal& b)
1799{
1800 return (mpfr_greater_p(a.mp,b.mp)!=0);
1801}
1802
1803inline bool operator > (const mpreal& a, const unsigned long int b)
1804{
1805 return a > mpreal(b);
1806}
1807
1808inline bool operator > (const mpreal& a, const unsigned int b)
1809{
1810 return a > mpreal(b);
1811}
1812
1813inline bool operator > (const mpreal& a, const long int b)
1814{
1815 return a > mpreal(b);
1816}
1817
1818inline bool operator > (const mpreal& a, const int b)
1819{
1820 return a > mpreal(b);
1821}
1822
1823inline bool operator > (const mpreal& a, const long double b)
1824{
1825 return a > mpreal(b);
1826}
1827
1828inline bool operator > (const mpreal& a, const double b)
1829{
1830 return a > mpreal(b);
1831}
1832
1833inline bool operator > (const unsigned long int a, const mpreal& b)
1834{
1835 return mpreal(a) > b;
1836}
1837
1838inline bool operator > (const unsigned int a, const mpreal& b)
1839{
1840 return mpreal(a) > b;
1841}
1842
1843inline bool operator > (const long int a, const mpreal& b)
1844{
1845 return mpreal(a) > b;
1846}
1847
1848inline bool operator > (const int a, const mpreal& b)
1849{
1850 return mpreal(a) > b;
1851}
1852
1853inline bool operator > (const long double a, const mpreal& b)
1854{
1855 return mpreal(a) > b;
1856}
1857
1858inline bool operator > (const double a, const mpreal& b)
1859{
1860 return mpreal(a) > b;
1861}
1862
1863inline bool operator >= (const mpreal& a, const mpreal& b)
1864{
1865 return (mpfr_greaterequal_p(a.mp,b.mp)!=0);
1866}
1867
1868inline bool operator >= (const mpreal& a, const unsigned long int b)
1869{
1870 return a >= mpreal(b);
1871}
1872
1873inline bool operator >= (const mpreal& a, const unsigned int b)
1874{
1875 return a >= mpreal(b);
1876}
1877
1878inline bool operator >= (const mpreal& a, const long int b)
1879{
1880 return a >= mpreal(b);
1881}
1882
1883inline bool operator >= (const mpreal& a, const int b)
1884{
1885 return a >= mpreal(b);
1886}
1887
1888inline bool operator >= (const mpreal& a, const long double b)
1889{
1890 return a >= mpreal(b);
1891}
1892
1893inline bool operator >= (const mpreal& a, const double b)
1894{
1895 return a >= mpreal(b);
1896}
1897
1898inline bool operator >= (const unsigned long int a,const mpreal& b)
1899{
1900 return mpreal(a) >= b;
1901}
1902
1903inline bool operator >= (const unsigned int a, const mpreal& b)
1904{
1905 return mpreal(a) >= b;
1906}
1907
1908inline bool operator >= (const long int a, const mpreal& b)
1909{
1910 return mpreal(a) >= b;
1911}
1912
1913inline bool operator >= (const int a, const mpreal& b)
1914{
1915 return mpreal(a) >= b;
1916}
1917
1918inline bool operator >= (const long double a, const mpreal& b)
1919{
1920 return mpreal(a) >= b;
1921}
1922
1923inline bool operator >= (const double a, const mpreal& b)
1924{
1925 return mpreal(a) >= b;
1926}
1927
1928inline bool operator < (const mpreal& a, const mpreal& b)
1929{
1930 return (mpfr_less_p(a.mp,b.mp)!=0);
1931}
1932
1933inline bool operator < (const mpreal& a, const unsigned long int b)
1934{
1935 return a<mpreal(b);
1936}
1937
1938inline bool operator < (const mpreal& a, const unsigned int b)
1939{
1940 return a<mpreal(b);
1941}
1942
1943inline bool operator < (const mpreal& a, const long int b)
1944{
1945 return a<mpreal(b);
1946}
1947
1948inline bool operator < (const mpreal& a, const int b)
1949{
1950 return a<mpreal(b);
1951}
1952
1953inline bool operator < (const mpreal& a, const long double b)
1954{
1955 return a<mpreal(b);
1956}
1957
1958inline bool operator < (const mpreal& a, const double b)
1959{
1960 return a<mpreal(b);
1961}
1962
1963inline bool operator < (const unsigned long int a, const mpreal& b)
1964{
1965 return mpreal(a)<b;
1966}
1967
1968inline bool operator < (const unsigned int a,const mpreal& b)
1969{
1970 return mpreal(a)<b;
1971}
1972
1973inline bool operator < (const long int a,const mpreal& b)
1974{
1975 return mpreal(a)<b;
1976}
1977
1978inline bool operator < (const int a,const mpreal& b)
1979{
1980 return mpreal(a)<b;
1981}
1982
1983inline bool operator < (const long double a,const mpreal& b)
1984{
1985 return mpreal(a)<b;
1986}
1987
1988inline bool operator < (const double a,const mpreal& b)
1989{
1990 return mpreal(a)<b;
1991}
1992
1993inline bool operator <= (const mpreal& a, const mpreal& b)
1994{
1995 return (mpfr_lessequal_p(a.mp,b.mp)!=0);
1996}
1997
1998inline bool operator <= (const mpreal& a, const unsigned long int b)
1999{
2000 return a<=mpreal(b);
2001}
2002
2003inline bool operator <= (const mpreal& a, const unsigned int b)
2004{
2005 return a<=mpreal(b);
2006}
2007
2008inline bool operator <= (const mpreal& a, const long int b)
2009{
2010 return a<=mpreal(b);
2011}
2012
2013inline bool operator <= (const mpreal& a, const int b)
2014{
2015 return a<=mpreal(b);
2016}
2017
2018inline bool operator <= (const mpreal& a, const long double b)
2019{
2020 return a<=mpreal(b);
2021}
2022
2023inline bool operator <= (const mpreal& a, const double b)
2024{
2025 return a<=mpreal(b);
2026}
2027
2028inline bool operator <= (const unsigned long int a,const mpreal& b)
2029{
2030 return mpreal(a)<=b;
2031}
2032
2033inline bool operator <= (const unsigned int a, const mpreal& b)
2034{
2035 return mpreal(a)<=b;
2036}
2037
2038inline bool operator <= (const long int a, const mpreal& b)
2039{
2040 return mpreal(a)<=b;
2041}
2042
2043inline bool operator <= (const int a, const mpreal& b)
2044{
2045 return mpreal(a)<=b;
2046}
2047
2048inline bool operator <= (const long double a, const mpreal& b)
2049{
2050 return mpreal(a)<=b;
2051}
2052
2053inline bool operator <= (const double a, const mpreal& b)
2054{
2055 return mpreal(a)<=b;
2056}
2057
2058inline bool operator == (const mpreal& a, const mpreal& b)
2059{
2060 return (mpfr_equal_p(a.mp,b.mp)!=0);
2061}
2062
2063inline bool operator == (const mpreal& a, const unsigned long int b)
2064{
2065 return (mpfr_cmp_ui(a.mp,b) == 0);
2066}
2067
2068inline bool operator == (const mpreal& a, const unsigned int b)
2069{
2070 return (mpfr_cmp_ui(a.mp,b) == 0);
2071}
2072
2073inline bool operator == (const mpreal& a, const long int b)
2074{
2075 return (mpfr_cmp_si(a.mp,b) == 0);
2076}
2077
2078inline bool operator == (const mpreal& a, const int b)
2079{
2080 return (mpfr_cmp_si(a.mp,b) == 0);
2081}
2082
2083inline bool operator == (const mpreal& a, const long double b)
2084{
2085 return (mpfr_cmp_ld(a.mp,b) == 0);
2086}
2087
2088inline bool operator == (const mpreal& a, const double b)
2089{
2090 return (mpfr_cmp_d(a.mp,b) == 0);
2091}
2092
2093inline bool operator == (const unsigned long int a, const mpreal& b)
2094{
2095 return b == a;
2096}
2097
2098inline bool operator == (const unsigned int a, const mpreal& b)
2099{
2100 return b == a;
2101}
2102
2103inline bool operator == (const long int a, const mpreal& b)
2104{
2105 return b == a;
2106}
2107
2108inline bool operator == (const int a, const mpreal& b)
2109{
2110 return b == a;
2111}
2112
2113inline bool operator == (const long double a, const mpreal& b)
2114{
2115 return b == a;
2116}
2117
2118inline bool operator == (const double a, const mpreal& b)
2119{
2120 return b == a;
2121}
2122
2123inline bool operator != (const mpreal& a, const mpreal& b)
2124{
2125 return (mpfr_lessgreater_p(a.mp,b.mp)!=0);
2126}
2127
2128inline bool operator != (const mpreal& a, const unsigned long int b)
2129{
2130 return !(a == b);
2131}
2132
2133inline bool operator != (const mpreal& a, const unsigned int b)
2134{
2135 return !(a == b);
2136}
2137
2138inline bool operator != (const mpreal& a, const long int b)
2139{
2140 return !(a == b);
2141}
2142
2143inline bool operator != (const mpreal& a, const int b)
2144{
2145 return !(a == b);
2146}
2147
2148inline bool operator != (const mpreal& a, const long double b)
2149{
2150 return !(a == b);
2151}
2152
2153inline bool operator != (const mpreal& a, const double b)
2154{
2155 return !(a == b);
2156}
2157
2158inline bool operator != (const unsigned long int a,const mpreal& b)
2159{
2160 return !(a == b);
2161}
2162
2163inline bool operator != (const unsigned int a, const mpreal& b)
2164{
2165 return !(a == b);
2166}
2167
2168inline bool operator != (const long int a, const mpreal& b)
2169{
2170 return !(a == b);
2171}
2172
2173inline bool operator != (const int a, const mpreal& b)
2174{
2175 return !(a == b);
2176}
2177
2178inline bool operator != (const long double a, const mpreal& b)
2179{
2180 return !(a == b);
2181}
2182
2183inline bool operator != (const double a, const mpreal& b)
2184{
2185 return !(a == b);
2186}
2187
2188inline bool isnan(const mpreal& v)
2189{
2190 return (mpfr_nan_p(v.mp)!=0);
2191}
2192
2193inline bool isinf(const mpreal& v)
2194{
2195 return (mpfr_inf_p(v.mp)!=0);
2196}
2197
2198inline bool isfinite(const mpreal& v)
2199{
2200 return (mpfr_number_p(v.mp)!=0);
2201}
2202
2203inline bool _iszero(const mpreal& v)
2204{
2205 return (mpfr_zero_p(v.mp)!=0);
2206}
2207
2208inline bool _isint(const mpreal& v)
2209{
2210 return (mpfr_integer_p(v.mp)!=0);
2211}
2212
2213#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
2214inline bool _isregular(const mpreal& v)
2215{
2216 return (mpfr_regular_p(v.mp));
2217}
2218#endif // MPFR 3.0.0 Specifics
2219
2220//////////////////////////////////////////////////////////////////////////
2221// Type Converters
2222long inline mpreal::toLong() const
2223{
2224 return mpfr_get_si(mp,GMP_RNDZ);
2225}
2226
2227unsigned long inline mpreal::toULong() const
2228{
2229 return mpfr_get_ui(mp,GMP_RNDZ);
2230}
2231
2232double inline mpreal::toDouble() const
2233{
2234 return mpfr_get_d(mp,default_rnd);
2235}
2236
2237long double inline mpreal::toLDouble() const
2238{
2239 return mpfr_get_ld(mp,default_rnd);
2240}
2241
2242#if defined (MPREAL_HAVE_INT64_SUPPORT)
2243int64_t inline mpreal::toInt64() const
2244{
2245 return mpfr_get_sj(mp,GMP_RNDZ);
2246}
2247
2248uint64_t inline mpreal::toUInt64() const
2249{
2250 return mpfr_get_uj(mp,GMP_RNDZ);
2251}
2252#endif
2253
2254inline mpreal::operator mpfr_ptr()
2255{
2256 return mp;
2257}
2258
2259inline mpreal::operator mpfr_srcptr() const
2260{
2261 return const_cast<mpfr_srcptr>(mp);
2262}
2263
2264//////////////////////////////////////////////////////////////////////////
2265// Bits - decimal digits relation
2266// bits = ceil(digits*log[2](10))
2267// digits = floor(bits*log[10](2))
2268
2269inline mp_prec_t digits2bits(int d)
2270{
2271 const double LOG2_10 = 3.3219280948873624;
2272
2273 d = 10>d?10:d;
2274
2275 return (mp_prec_t)std::ceil((d)*LOG2_10);
2276}
2277
2278inline int bits2digits(mp_prec_t b)
2279{
2280 const double LOG10_2 = 0.30102999566398119;
2281
2282 b = 34>b?34:b;
2283
2284 return (int)std::floor((b)*LOG10_2);
2285}
2286
2287//////////////////////////////////////////////////////////////////////////
2288// Set/Get number properties
2289inline int sgn(const mpreal& v)
2290{
2291 int r = mpfr_signbit(v.mp);
2292 return (r>0?-1:1);
2293}
2294
2295inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode)
2296{
2297 mpfr_setsign(mp,mp,(sign<0?1:0),RoundingMode);
2298 return *this;
2299}
2300
2301inline int mpreal::getPrecision() const
2302{
2303 return mpfr_get_prec(mp);
2304}
2305
2306inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode)
2307{
2308 mpfr_prec_round(mp,Precision, RoundingMode);
2309 return *this;
2310}
2311
2312inline mpreal& mpreal::setInf(int sign)
2313{
2314 mpfr_set_inf(mp,sign);
2315 return *this;
2316}
2317
2318inline mpreal& mpreal::setNan()
2319{
2320 mpfr_set_nan(mp);
2321 return *this;
2322}
2323
2324inline mpreal& mpreal::setZero(int sign)
2325{
2326 mpfr_set_zero(mp,sign);
2327 return *this;
2328}
2329
2330inline mp_prec_t mpreal::get_prec() const
2331{
2332 return mpfr_get_prec(mp);
2333}
2334
2335inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode)
2336{
2337 mpfr_prec_round(mp,prec,rnd_mode);
2338}
2339
2340inline mp_exp_t mpreal::get_exp ()
2341{
2342 return mpfr_get_exp(mp);
2343}
2344
2345inline int mpreal::set_exp (mp_exp_t e)
2346{
2347 return mpfr_set_exp(mp,e);
2348}
2349
2350inline const mpreal frexp(const mpreal& v, mp_exp_t* exp)
2351{
2352 mpreal x(v);
2353 *exp = x.get_exp();
2354 x.set_exp(0);
2355 return x;
2356}
2357
2358inline const mpreal ldexp(const mpreal& v, mp_exp_t exp)
2359{
2360 mpreal x(v);
2361
2362 // rounding is not important since we just increasing the exponent
2363 mpfr_mul_2si(x.mp,x.mp,exp,mpreal::default_rnd);
2364 return x;
2365}
2366
2367inline const mpreal machine_epsilon(mp_prec_t prec)
2368{
2369 // the smallest eps such that 1.0+eps != 1.0
2370 // depends (of cause) on the precision
2371 return machine_epsilon(mpreal(1,prec));
2372}
2373
2374inline const mpreal machine_epsilon(const mpreal& x)
2375{
2376 if( x < 0)
2377 {
2378 return nextabove(-x)+x;
2379 }else{
2380 return nextabove(x)-x;
2381 }
2382}
2383
2384inline const mpreal mpreal_min(mp_prec_t prec)
2385{
2386 // min = 1/2*2^emin = 2^(emin-1)
2387
2388 return mpreal(1,prec) << mpreal::get_emin()-1;
2389}
2390
2391inline const mpreal mpreal_max(mp_prec_t prec)
2392{
2393 // max = (1-eps)*2^emax, assume eps = 0?,
2394 // and use emax-1 to prevent value to be +inf
2395 // max = 2^(emax-1)
2396
2397 return mpreal(1,prec) << mpreal::get_emax()-1;
2398}
2399
2400inline const mpreal modf(const mpreal& v, mpreal& n)
2401{
2402 mpreal frac(v);
2403
2404 // rounding is not important since we are using the same number
2405 mpfr_frac(frac.mp,frac.mp,mpreal::default_rnd);
2406 mpfr_trunc(n.mp,v.mp);
2407 return frac;
2408}
2409
2410inline int mpreal::check_range (int t, mp_rnd_t rnd_mode)
2411{
2412 return mpfr_check_range(mp,t,rnd_mode);
2413}
2414
2415inline int mpreal::subnormalize (int t,mp_rnd_t rnd_mode)
2416{
2417 return mpfr_subnormalize(mp,t,rnd_mode);
2418}
2419
2420inline mp_exp_t mpreal::get_emin (void)
2421{
2422 return mpfr_get_emin();
2423}
2424
2425inline int mpreal::set_emin (mp_exp_t exp)
2426{
2427 return mpfr_set_emin(exp);
2428}
2429
2430inline mp_exp_t mpreal::get_emax (void)
2431{
2432 return mpfr_get_emax();
2433}
2434
2435inline int mpreal::set_emax (mp_exp_t exp)
2436{
2437 return mpfr_set_emax(exp);
2438}
2439
2440inline mp_exp_t mpreal::get_emin_min (void)
2441{
2442 return mpfr_get_emin_min();
2443}
2444
2445inline mp_exp_t mpreal::get_emin_max (void)
2446{
2447 return mpfr_get_emin_max();
2448}
2449
2450inline mp_exp_t mpreal::get_emax_min (void)
2451{
2452 return mpfr_get_emax_min();
2453}
2454
2455inline mp_exp_t mpreal::get_emax_max (void)
2456{
2457 return mpfr_get_emax_max();
2458}
2459
2460//////////////////////////////////////////////////////////////////////////
2461// Mathematical Functions
2462//////////////////////////////////////////////////////////////////////////
2463inline const mpreal sqr(const mpreal& v, mp_rnd_t rnd_mode)
2464{
2465 mpreal x(v);
2466 mpfr_sqr(x.mp,x.mp,rnd_mode);
2467 return x;
2468}
2469
2470inline const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode)
2471{
2472 mpreal x(v);
2473 mpfr_sqrt(x.mp,x.mp,rnd_mode);
2474 return x;
2475}
2476
2477inline const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode)
2478{
2479 mpreal x;
2480 mpfr_sqrt_ui(x.mp,v,rnd_mode);
2481 return x;
2482}
2483
2484inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode)
2485{
2486 return sqrt(static_cast<unsigned long int>(v),rnd_mode);
2487}
2488
2489inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode)
2490{
2491 if (v>=0) return sqrt(static_cast<unsigned long int>(v),rnd_mode);
2492 else return mpreal().setNan(); // NaN
2493}
2494
2495inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode)
2496{
2497 if (v>=0) return sqrt(static_cast<unsigned long int>(v),rnd_mode);
2498 else return mpreal().setNan(); // NaN
2499}
2500
2501inline const mpreal sqrt(const long double v, mp_rnd_t rnd_mode)
2502{
2503 return sqrt(mpreal(v),rnd_mode);
2504}
2505
2506inline const mpreal sqrt(const double v, mp_rnd_t rnd_mode)
2507{
2508 return sqrt(mpreal(v),rnd_mode);
2509}
2510
2511inline const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode)
2512{
2513 mpreal x(v);
2514 mpfr_cbrt(x.mp,x.mp,rnd_mode);
2515 return x;
2516}
2517
2518inline const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode)
2519{
2520 mpreal x(v);
2521 mpfr_root(x.mp,x.mp,k,rnd_mode);
2522 return x;
2523}
2524
2525inline const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode)
2526{
2527 mpreal x(v);
2528 mpfr_abs(x.mp,x.mp,rnd_mode);
2529 return x;
2530}
2531
2532inline const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode)
2533{
2534 mpreal x(v);
2535 mpfr_abs(x.mp,x.mp,rnd_mode);
2536 return x;
2537}
2538
2539inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode)
2540{
2541 mpreal x(a);
2542 mpfr_dim(x.mp,a.mp,b.mp,rnd_mode);
2543 return x;
2544}
2545
2546inline int cmpabs(const mpreal& a,const mpreal& b)
2547{
2548 return mpfr_cmpabs(a.mp,b.mp);
2549}
2550
2551inline const mpreal log (const mpreal& v, mp_rnd_t rnd_mode)
2552{
2553 mpreal x(v);
2554 mpfr_log(x.mp,v.mp,rnd_mode);
2555 return x;
2556}
2557
2558inline const mpreal log2(const mpreal& v, mp_rnd_t rnd_mode)
2559{
2560 mpreal x(v);
2561 mpfr_log2(x.mp,v.mp,rnd_mode);
2562 return x;
2563}
2564
2565inline const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode)
2566{
2567 mpreal x(v);
2568 mpfr_log10(x.mp,v.mp,rnd_mode);
2569 return x;
2570}
2571
2572inline const mpreal exp(const mpreal& v, mp_rnd_t rnd_mode)
2573{
2574 mpreal x(v);
2575 mpfr_exp(x.mp,v.mp,rnd_mode);
2576 return x;
2577}
2578
2579inline const mpreal exp2(const mpreal& v, mp_rnd_t rnd_mode)
2580{
2581 mpreal x(v);
2582 mpfr_exp2(x.mp,v.mp,rnd_mode);
2583 return x;
2584}
2585
2586inline const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode)
2587{
2588 mpreal x(v);
2589 mpfr_exp10(x.mp,v.mp,rnd_mode);
2590 return x;
2591}
2592
2593inline const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode)
2594{
2595 mpreal x(v);
2596 mpfr_cos(x.mp,v.mp,rnd_mode);
2597 return x;
2598}
2599
2600inline const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode)
2601{
2602 mpreal x(v);
2603 mpfr_sin(x.mp,v.mp,rnd_mode);
2604 return x;
2605}
2606
2607inline const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode)
2608{
2609 mpreal x(v);
2610 mpfr_tan(x.mp,v.mp,rnd_mode);
2611 return x;
2612}
2613
2614inline const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode)
2615{
2616 mpreal x(v);
2617 mpfr_sec(x.mp,v.mp,rnd_mode);
2618 return x;
2619}
2620
2621inline const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode)
2622{
2623 mpreal x(v);
2624 mpfr_csc(x.mp,v.mp,rnd_mode);
2625 return x;
2626}
2627
2628inline const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode)
2629{
2630 mpreal x(v);
2631 mpfr_cot(x.mp,v.mp,rnd_mode);
2632 return x;
2633}
2634
2635inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)
2636{
2637 return mpfr_sin_cos(s.mp,c.mp,v.mp,rnd_mode);
2638}
2639
2640inline const mpreal acos (const mpreal& v, mp_rnd_t rnd_mode)
2641{
2642 mpreal x(v);
2643 mpfr_acos(x.mp,v.mp,rnd_mode);
2644 return x;
2645}
2646
2647inline const mpreal asin (const mpreal& v, mp_rnd_t rnd_mode)
2648{
2649 mpreal x(v);
2650 mpfr_asin(x.mp,v.mp,rnd_mode);
2651 return x;
2652}
2653
2654inline const mpreal atan (const mpreal& v, mp_rnd_t rnd_mode)
2655{
2656 mpreal x(v);
2657 mpfr_atan(x.mp,v.mp,rnd_mode);
2658 return x;
2659}
2660
2661inline const mpreal acot (const mpreal& v, mp_rnd_t rnd_mode)
2662{
2663 return atan(1/v, rnd_mode);
2664}
2665
2666inline const mpreal asec (const mpreal& v, mp_rnd_t rnd_mode)
2667{
2668 return acos(1/v, rnd_mode);
2669}
2670
2671inline const mpreal acsc (const mpreal& v, mp_rnd_t rnd_mode)
2672{
2673 return asin(1/v, rnd_mode);
2674}
2675
2676inline const mpreal acoth (const mpreal& v, mp_rnd_t rnd_mode)
2677{
2678 return atanh(1/v, rnd_mode);
2679}
2680
2681inline const mpreal asech (const mpreal& v, mp_rnd_t rnd_mode)
2682{
2683 return acosh(1/v, rnd_mode);
2684}
2685
2686inline const mpreal acsch (const mpreal& v, mp_rnd_t rnd_mode)
2687{
2688 return asinh(1/v, rnd_mode);
2689}
2690
2691inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode)
2692{
2693 mpreal a;
2694 mp_prec_t yp, xp;
2695
2696 yp = y.get_prec();
2697 xp = x.get_prec();
2698
2699 a.set_prec(yp>xp?yp:xp);
2700
2701 mpfr_atan2(a.mp, y.mp, x.mp, rnd_mode);
2702
2703 return a;
2704}
2705
2706inline const mpreal cosh (const mpreal& v, mp_rnd_t rnd_mode)
2707{
2708 mpreal x(v);
2709 mpfr_cosh(x.mp,v.mp,rnd_mode);
2710 return x;
2711}
2712
2713inline const mpreal sinh (const mpreal& v, mp_rnd_t rnd_mode)
2714{
2715 mpreal x(v);
2716 mpfr_sinh(x.mp,v.mp,rnd_mode);
2717 return x;
2718}
2719
2720inline const mpreal tanh (const mpreal& v, mp_rnd_t rnd_mode)
2721{
2722 mpreal x(v);
2723 mpfr_tanh(x.mp,v.mp,rnd_mode);
2724 return x;
2725}
2726
2727inline const mpreal sech (const mpreal& v, mp_rnd_t rnd_mode)
2728{
2729 mpreal x(v);
2730 mpfr_sech(x.mp,v.mp,rnd_mode);
2731 return x;
2732}
2733
2734inline const mpreal csch (const mpreal& v, mp_rnd_t rnd_mode)
2735{
2736 mpreal x(v);
2737 mpfr_csch(x.mp,v.mp,rnd_mode);
2738 return x;
2739}
2740
2741inline const mpreal coth (const mpreal& v, mp_rnd_t rnd_mode)
2742{
2743 mpreal x(v);
2744 mpfr_coth(x.mp,v.mp,rnd_mode);
2745 return x;
2746}
2747
2748inline const mpreal acosh (const mpreal& v, mp_rnd_t rnd_mode)
2749{
2750 mpreal x(v);
2751 mpfr_acosh(x.mp,v.mp,rnd_mode);
2752 return x;
2753}
2754
2755inline const mpreal asinh (const mpreal& v, mp_rnd_t rnd_mode)
2756{
2757 mpreal x(v);
2758 mpfr_asinh(x.mp,v.mp,rnd_mode);
2759 return x;
2760}
2761
2762inline const mpreal atanh (const mpreal& v, mp_rnd_t rnd_mode)
2763{
2764 mpreal x(v);
2765 mpfr_atanh(x.mp,v.mp,rnd_mode);
2766 return x;
2767}
2768
2769inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
2770{
2771 mpreal a;
2772 mp_prec_t yp, xp;
2773
2774 yp = y.get_prec();
2775 xp = x.get_prec();
2776
2777 a.set_prec(yp>xp?yp:xp);
2778
2779 mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode);
2780
2781 return a;
2782}
2783
2784inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
2785{
2786 mpreal a;
2787 mp_prec_t yp, xp;
2788
2789 yp = y.get_prec();
2790 xp = x.get_prec();
2791
2792 a.set_prec(yp>xp?yp:xp);
2793
2794 mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode);
2795
2796 return a;
2797}
2798
2799inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode)
2800{
2801 mpreal x(0,prec);
2802 mpfr_fac_ui(x.mp,v,rnd_mode);
2803 return x;
2804}
2805
2806inline const mpreal log1p (const mpreal& v, mp_rnd_t rnd_mode)
2807{
2808 mpreal x(v);
2809 mpfr_log1p(x.mp,v.mp,rnd_mode);
2810 return x;
2811}
2812
2813inline const mpreal expm1 (const mpreal& v, mp_rnd_t rnd_mode)
2814{
2815 mpreal x(v);
2816 mpfr_expm1(x.mp,v.mp,rnd_mode);
2817 return x;
2818}
2819
2820inline const mpreal eint (const mpreal& v, mp_rnd_t rnd_mode)
2821{
2822 mpreal x(v);
2823 mpfr_eint(x.mp,v.mp,rnd_mode);
2824 return x;
2825}
2826
2827inline const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode)
2828{
2829 mpreal x(v);
2830 mpfr_gamma(x.mp,v.mp,rnd_mode);
2831 return x;
2832}
2833
2834inline const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode)
2835{
2836 mpreal x(v);
2837 mpfr_lngamma(x.mp,v.mp,rnd_mode);
2838 return x;
2839}
2840
2841inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode)
2842{
2843 mpreal x(v);
2844 int tsignp;
2845
2846 if(signp)
2847 mpfr_lgamma(x.mp,signp,v.mp,rnd_mode);
2848 else
2849 mpfr_lgamma(x.mp,&tsignp,v.mp,rnd_mode);
2850
2851 return x;
2852}
2853
2854inline const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode)
2855{
2856 mpreal x(v);
2857 mpfr_zeta(x.mp,v.mp,rnd_mode);
2858 return x;
2859}
2860
2861inline const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode)
2862{
2863 mpreal x(v);
2864 mpfr_erf(x.mp,v.mp,rnd_mode);
2865 return x;
2866}
2867
2868inline const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode)
2869{
2870 mpreal x(v);
2871 mpfr_erfc(x.mp,v.mp,rnd_mode);
2872 return x;
2873}
2874
2875inline const mpreal _j0 (const mpreal& v, mp_rnd_t rnd_mode)
2876{
2877 mpreal x(v);
2878 mpfr_j0(x.mp,v.mp,rnd_mode);
2879 return x;
2880}
2881
2882inline const mpreal _j1 (const mpreal& v, mp_rnd_t rnd_mode)
2883{
2884 mpreal x(v);
2885 mpfr_j1(x.mp,v.mp,rnd_mode);
2886 return x;
2887}
2888
2889inline const mpreal _jn (long n, const mpreal& v, mp_rnd_t rnd_mode)
2890{
2891 mpreal x(v);
2892 mpfr_jn(x.mp,n,v.mp,rnd_mode);
2893 return x;
2894}
2895
2896inline const mpreal _y0 (const mpreal& v, mp_rnd_t rnd_mode)
2897{
2898 mpreal x(v);
2899 mpfr_y0(x.mp,v.mp,rnd_mode);
2900 return x;
2901}
2902
2903inline const mpreal _y1 (const mpreal& v, mp_rnd_t rnd_mode)
2904{
2905 mpreal x(v);
2906 mpfr_y1(x.mp,v.mp,rnd_mode);
2907 return x;
2908}
2909
2910inline const mpreal _yn (long n, const mpreal& v, mp_rnd_t rnd_mode)
2911{
2912 mpreal x(v);
2913 mpfr_yn(x.mp,n,v.mp,rnd_mode);
2914 return x;
2915}
2916
2917//////////////////////////////////////////////////////////////////////////
2918// MPFR 2.4.0 Specifics
2919#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
2920
2921inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode)
2922{
2923 return mpfr_sinh_cosh(s.mp,c.mp,v.mp,rnd_mode);
2924}
2925
2926inline const mpreal li2(const mpreal& v, mp_rnd_t rnd_mode)
2927{
2928 mpreal x(v);
2929 mpfr_li2(x.mp,v.mp,rnd_mode);
2930 return x;
2931}
2932
2933inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
2934{
2935 mpreal a;
2936 mp_prec_t yp, xp;
2937
2938 yp = y.get_prec();
2939 xp = x.get_prec();
2940
2941 a.set_prec(yp>xp?yp:xp);
2942
2943 mpfr_fmod(a.mp, x.mp, y.mp, rnd_mode);
2944
2945 return a;
2946}
2947
2948inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode)
2949{
2950 mpreal x(v);
2951 mpfr_rec_sqrt(x.mp,v.mp,rnd_mode);
2952 return x;
2953}
2954#endif // MPFR 2.4.0 Specifics
2955
2956//////////////////////////////////////////////////////////////////////////
2957// MPFR 3.0.0 Specifics
2958#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
2959
2960inline const mpreal digamma(const mpreal& v, mp_rnd_t rnd_mode)
2961{
2962 mpreal x(v);
2963 mpfr_digamma(x.mp,v.mp,rnd_mode);
2964 return x;
2965}
2966
2967inline const mpreal ai(const mpreal& v, mp_rnd_t rnd_mode)
2968{
2969 mpreal x(v);
2970 mpfr_ai(x.mp,v.mp,rnd_mode);
2971 return x;
2972}
2973
2974#endif // MPFR 3.0.0 Specifics
2975
2976//////////////////////////////////////////////////////////////////////////
2977// Constants
2978inline const mpreal const_log2 (mp_prec_t prec, mp_rnd_t rnd_mode)
2979{
2980 mpreal x;
2981 x.set_prec(prec);
2982 mpfr_const_log2(x.mp,rnd_mode);
2983 return x;
2984}
2985
2986inline const mpreal const_pi (mp_prec_t prec, mp_rnd_t rnd_mode)
2987{
2988 mpreal x;
2989 x.set_prec(prec);
2990 mpfr_const_pi(x.mp,rnd_mode);
2991 return x;
2992}
2993
2994inline const mpreal const_euler (mp_prec_t prec, mp_rnd_t rnd_mode)
2995{
2996 mpreal x;
2997 x.set_prec(prec);
2998 mpfr_const_euler(x.mp,rnd_mode);
2999 return x;
3000}
3001
3002inline const mpreal const_catalan (mp_prec_t prec, mp_rnd_t rnd_mode)
3003{
3004 mpreal x;
3005 x.set_prec(prec);
3006 mpfr_const_catalan(x.mp,rnd_mode);
3007 return x;
3008}
3009
3010inline const mpreal const_infinity (int sign, mp_prec_t prec, mp_rnd_t rnd_mode)
3011{
3012 mpreal x;
3013 x.set_prec(prec,rnd_mode);
3014 mpfr_set_inf(x.mp, sign);
3015 return x;
3016}
3017
3018//////////////////////////////////////////////////////////////////////////
3019// Integer Related Functions
3020inline const mpreal rint(const mpreal& v, mp_rnd_t rnd_mode)
3021{
3022 mpreal x(v);
3023 mpfr_rint(x.mp,v.mp,rnd_mode);
3024 return x;
3025}
3026
3027inline const mpreal ceil(const mpreal& v)
3028{
3029 mpreal x(v);
3030 mpfr_ceil(x.mp,v.mp);
3031 return x;
3032
3033}
3034
3035inline const mpreal floor(const mpreal& v)
3036{
3037 mpreal x(v);
3038 mpfr_floor(x.mp,v.mp);
3039 return x;
3040}
3041
3042inline const mpreal round(const mpreal& v)
3043{
3044 mpreal x(v);
3045 mpfr_round(x.mp,v.mp);
3046 return x;
3047}
3048
3049inline const mpreal trunc(const mpreal& v)
3050{
3051 mpreal x(v);
3052 mpfr_trunc(x.mp,v.mp);
3053 return x;
3054}
3055
3056inline const mpreal rint_ceil (const mpreal& v, mp_rnd_t rnd_mode)
3057{
3058 mpreal x(v);
3059 mpfr_rint_ceil(x.mp,v.mp,rnd_mode);
3060 return x;
3061}
3062
3063inline const mpreal rint_floor(const mpreal& v, mp_rnd_t rnd_mode)
3064{
3065 mpreal x(v);
3066 mpfr_rint_floor(x.mp,v.mp,rnd_mode);
3067 return x;
3068}
3069
3070inline const mpreal rint_round(const mpreal& v, mp_rnd_t rnd_mode)
3071{
3072 mpreal x(v);
3073 mpfr_rint_round(x.mp,v.mp,rnd_mode);
3074 return x;
3075}
3076
3077inline const mpreal rint_trunc(const mpreal& v, mp_rnd_t rnd_mode)
3078{
3079 mpreal x(v);
3080 mpfr_rint_trunc(x.mp,v.mp,rnd_mode);
3081 return x;
3082}
3083
3084inline const mpreal frac (const mpreal& v, mp_rnd_t rnd_mode)
3085{
3086 mpreal x(v);
3087 mpfr_frac(x.mp,v.mp,rnd_mode);
3088 return x;
3089}
3090
3091//////////////////////////////////////////////////////////////////////////
3092// Miscellaneous Functions
3093inline void swap(mpreal& a, mpreal& b)
3094{
3095 mpfr_swap(a.mp,b.mp);
3096}
3097
3098inline const mpreal (max)(const mpreal& x, const mpreal& y)
3099{
3100 return (x>y?x:y);
3101}
3102
3103inline const mpreal (min)(const mpreal& x, const mpreal& y)
3104{
3105 return (x<y?x:y);
3106}
3107
3108inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
3109{
3110 mpreal a;
3111 mpfr_max(a.mp,x.mp,y.mp,rnd_mode);
3112 return a;
3113}
3114
3115inline const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
3116{
3117 mpreal a;
3118 mpfr_min(a.mp,x.mp,y.mp,rnd_mode);
3119 return a;
3120}
3121
3122inline const mpreal nexttoward (const mpreal& x, const mpreal& y)
3123{
3124 mpreal a(x);
3125 mpfr_nexttoward(a.mp,y.mp);
3126 return a;
3127}
3128
3129inline const mpreal nextabove (const mpreal& x)
3130{
3131 mpreal a(x);
3132 mpfr_nextabove(a.mp);
3133 return a;
3134}
3135
3136inline const mpreal nextbelow (const mpreal& x)
3137{
3138 mpreal a(x);
3139 mpfr_nextbelow(a.mp);
3140 return a;
3141}
3142
3143inline const mpreal urandomb (gmp_randstate_t& state)
3144{
3145 mpreal x;
3146 mpfr_urandomb(x.mp,state);
3147 return x;
3148}
3149
3150#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
3151// use gmp_randinit_default() to init state, gmp_randclear() to clear
3152inline const mpreal urandom (gmp_randstate_t& state, mp_rnd_t rnd_mode)
3153{
3154 mpreal x;
3155 mpfr_urandom(x.mp,state,rnd_mode);
3156 return x;
3157}
3158#endif
3159
3160#if (MPFR_VERSION <= MPFR_VERSION_NUM(2,4,2))
3161inline const mpreal random2 (mp_size_t size, mp_exp_t exp)
3162{
3163 mpreal x;
3164 mpfr_random2(x.mp,size,exp);
3165 return x;
3166}
3167#endif
3168
3169// Uniformly distributed random number generation
3170// a = random(seed); <- initialization & first random number generation
3171// a = random(); <- next random numbers generation
3172// seed != 0
3173inline const mpreal random(unsigned int seed)
3174{
3175
3176#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
3177 static gmp_randstate_t state;
3178 static bool isFirstTime = true;
3179
3180 if(isFirstTime)
3181 {
3182 gmp_randinit_default(state);
3183 gmp_randseed_ui(state,0);
3184 isFirstTime = false;
3185 }
3186
3187 if(seed != 0) gmp_randseed_ui(state,seed);
3188
3189 return mpfr::urandom(state);
3190#else
3191 if(seed != 0) std::srand(seed);
3192 return mpfr::mpreal(std::rand()/(double)RAND_MAX);
3193#endif
3194
3195}
3196
3197//////////////////////////////////////////////////////////////////////////
3198// Set/Get global properties
3199inline void mpreal::set_default_prec(mp_prec_t prec)
3200{
3201 default_prec = prec;
3202 mpfr_set_default_prec(prec);
3203}
3204
3205inline mp_prec_t mpreal::get_default_prec()
3206{
3207 return (mpfr_get_default_prec)();
3208}
3209
3210inline void mpreal::set_default_base(int base)
3211{
3212 default_base = base;
3213}
3214
3215inline int mpreal::get_default_base()
3216{
3217 return default_base;
3218}
3219
3220inline void mpreal::set_default_rnd(mp_rnd_t rnd_mode)
3221{
3222 default_rnd = rnd_mode;
3223 mpfr_set_default_rounding_mode(rnd_mode);
3224}
3225
3226inline mp_rnd_t mpreal::get_default_rnd()
3227{
3228 return static_cast<mp_rnd_t>((mpfr_get_default_rounding_mode)());
3229}
3230
3231inline void mpreal::set_double_bits(int dbits)
3232{
3233 double_bits = dbits;
3234}
3235
3236inline int mpreal::get_double_bits()
3237{
3238 return double_bits;
3239}
3240
3241inline bool mpreal::fits_in_bits(double x, int n)
3242{
3243 int i;
3244 double t;
3245 return IsInf(x) || (std::modf ( std::ldexp ( std::frexp ( x, &i ), n ), &t ) == 0.0);
3246}
3247
3248inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode)
3249{
3250 mpreal x(a);
3251 mpfr_pow(x.mp,x.mp,b.mp,rnd_mode);
3252 return x;
3253}
3254
3255inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode)
3256{
3257 mpreal x(a);
3258 mpfr_pow_z(x.mp,x.mp,b,rnd_mode);
3259 return x;
3260}
3261
3262inline const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode)
3263{
3264 mpreal x(a);
3265 mpfr_pow_ui(x.mp,x.mp,b,rnd_mode);
3266 return x;
3267}
3268
3269inline const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode)
3270{
3271 return pow(a,static_cast<unsigned long int>(b),rnd_mode);
3272}
3273
3274inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode)
3275{
3276 mpreal x(a);
3277 mpfr_pow_si(x.mp,x.mp,b,rnd_mode);
3278 return x;
3279}
3280
3281inline const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode)
3282{
3283 return pow(a,static_cast<long int>(b),rnd_mode);
3284}
3285
3286inline const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode)
3287{
3288 return pow(a,mpreal(b),rnd_mode);
3289}
3290
3291inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode)
3292{
3293 return pow(a,mpreal(b),rnd_mode);
3294}
3295
3296inline const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode)
3297{
3298 mpreal x(a);
3299 mpfr_ui_pow(x.mp,a,b.mp,rnd_mode);
3300 return x;
3301}
3302
3303inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode)
3304{
3305 return pow(static_cast<unsigned long int>(a),b,rnd_mode);
3306}
3307
3308inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode)
3309{
3310 if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode);
3311 else return pow(mpreal(a),b,rnd_mode);
3312}
3313
3314inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode)
3315{
3316 if (a>=0) return pow(static_cast<unsigned long int>(a),b,rnd_mode);
3317 else return pow(mpreal(a),b,rnd_mode);
3318}
3319
3320inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode)
3321{
3322 return pow(mpreal(a),b,rnd_mode);
3323}
3324
3325inline const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode)
3326{
3327 return pow(mpreal(a),b,rnd_mode);
3328}
3329
3330// pow unsigned long int
3331inline const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode)
3332{
3333 mpreal x(a);
3334 mpfr_ui_pow_ui(x.mp,a,b,rnd_mode);
3335 return x;
3336}
3337
3338inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode)
3339{
3340 return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3341}
3342
3343inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode)
3344{
3345 if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3346 else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
3347}
3348
3349inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode)
3350{
3351 if(b>0) return pow(a,static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3352 else return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
3353}
3354
3355inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode)
3356{
3357 return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
3358}
3359
3360inline const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode)
3361{
3362 return pow(a,mpreal(b),rnd_mode); //mpfr_ui_pow
3363}
3364
3365// pow unsigned int
3366inline const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode)
3367{
3368 return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui
3369}
3370
3371inline const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode)
3372{
3373 return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3374}
3375
3376inline const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode)
3377{
3378 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3379 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3380}
3381
3382inline const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode)
3383{
3384 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3385 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3386}
3387
3388inline const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode)
3389{
3390 return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3391}
3392
3393inline const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode)
3394{
3395 return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3396}
3397
3398// pow long int
3399inline const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode)
3400{
3401 if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui
3402 else return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui
3403}
3404
3405inline const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode)
3406{
3407 if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3408 else return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui
3409}
3410
3411inline const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode)
3412{
3413 if (a>0)
3414 {
3415 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3416 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3417 }else{
3418 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si
3419 }
3420}
3421
3422inline const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode)
3423{
3424 if (a>0)
3425 {
3426 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3427 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3428 }else{
3429 return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si
3430 }
3431}
3432
3433inline const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode)
3434{
3435 if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3436 else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow
3437}
3438
3439inline const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode)
3440{
3441 if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3442 else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow
3443}
3444
3445// pow int
3446inline const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode)
3447{
3448 if (a>0) return pow(static_cast<unsigned long int>(a),b,rnd_mode); //mpfr_ui_pow_ui
3449 else return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui
3450}
3451
3452inline const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode)
3453{
3454 if (a>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3455 else return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui
3456}
3457
3458inline const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode)
3459{
3460 if (a>0)
3461 {
3462 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3463 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3464 }else{
3465 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si
3466 }
3467}
3468
3469inline const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode)
3470{
3471 if (a>0)
3472 {
3473 if(b>0) return pow(static_cast<unsigned long int>(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_ui_pow_ui
3474 else return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3475 }else{
3476 return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si
3477 }
3478}
3479
3480inline const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode)
3481{
3482 if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3483 else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow
3484}
3485
3486inline const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode)
3487{
3488 if (a>=0) return pow(static_cast<unsigned long int>(a),mpreal(b),rnd_mode); //mpfr_ui_pow
3489 else return pow(mpreal(a),mpreal(b),rnd_mode); //mpfr_pow
3490}
3491
3492// pow long double
3493inline const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode)
3494{
3495 return pow(mpreal(a),mpreal(b),rnd_mode);
3496}
3497
3498inline const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode)
3499{
3500 return pow(mpreal(a),b,rnd_mode); //mpfr_pow_ui
3501}
3502
3503inline const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode)
3504{
3505 return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); //mpfr_pow_ui
3506}
3507
3508inline const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode)
3509{
3510 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si
3511}
3512
3513inline const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode)
3514{
3515 return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si
3516}
3517
3518inline const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode)
3519{
3520 return pow(mpreal(a),mpreal(b),rnd_mode);
3521}
3522
3523inline const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode)
3524{
3525 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_ui
3526}
3527
3528inline const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode)
3529{
3530 return pow(mpreal(a),static_cast<unsigned long int>(b),rnd_mode); // mpfr_pow_ui
3531}
3532
3533inline const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode)
3534{
3535 return pow(mpreal(a),b,rnd_mode); // mpfr_pow_si
3536}
3537
3538inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode)
3539{
3540 return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si
3541}
3542
3543} // End of mpfr namespace
3544
3545// Explicit specialization of std::swap for mpreal numbers
3546// Thus standard algorithms will use efficient version of swap (due to Koenig lookup)
3547// Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap
3548namespace std
3549{
3550 template <>
3551 inline void swap(mpfr::mpreal& x, mpfr::mpreal& y)
3552 {
3553 return mpfr::swap(x, y);
3554 }
3555}
3556
3557#endif /* __MPREAL_H__ */