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