?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge commit 'crater/master'
[dragonfly.git] / contrib / mpfr / mpfr-impl.h
blobfd5ce58a1cef872116240bc4c09ec6dfc6311fc5
1 /* Utilities for MPFR developers, not exported.
2
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 2.1 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
22
23 #ifndef __MPFR_IMPL_H__
24 #define __MPFR_IMPL_H__
25
26 /* Let's include some standard headers unconditionally as they are
27 already needed by several source files or when some options are
28 enabled/disabled, and it is easy to forget them (some configure
29 options may hide the error).
30 Note: If some source file must not have such a header included
31 (which is very unlikely and probably means something broken in
32 this source file), we should do that with some macro (that would
33 also force to disable incompatible features). */
34 #if defined (__cplusplus)
35 #include <cstdio>
36 #else
37 #include <stdio.h>
38 #endif
39 #include <limits.h>
40
41 /*
42 * DragonFly: needed to build gcc44 on DragonFly-2.2, since the
43 * headers are confused and define INTMAX_C et al, but not intmax_t.
44 * So we include stdint.h explicitly to have these types for sure.
45 */
46 #include <stdint.h>
47
48 /* Check if we are inside a build of MPFR or inside the test suite.
49 This is needed in mpfr.h to export or import the functions.
50 It matters only for Windows DLL */
51 #ifndef __MPFR_TEST_H__
52 # define __MPFR_WITHIN_MPFR 1
53 #endif
54
55 /******************************************************
56 ****************** Include files *********************
57 ******************************************************/
58
59 /* Include 'config.h' before using ANY configure macros if needed
60 NOTE: It isn't MPFR 'config.h', but GMP's one! */
61 #ifdef HAVE_CONFIG_H
62 # include "config.h"
63 #endif
64
65 #ifdef MPFR_HAVE_GMP_IMPL /* Build with gmp internals*/
66
67 # ifndef __GMP_H__
68 # include "gmp.h"
69 # endif
70 # ifndef __GMP_IMPL_H__
71 # include "gmp-impl.h"
72 # endif
73 # ifdef MPFR_NEED_LONGLONG_H
74 # include "longlong.h"
75 # endif
76 # ifndef __MPFR_H
77 # include "mpfr.h"
78 # endif
79
80 #else /* Build without gmp internals */
81
82 # ifndef __GMP_H__
83 # include "gmp.h"
84 # endif
85 # ifndef __MPFR_H
86 # include "mpfr.h"
87 # endif
88 # ifndef __GMPFR_GMP_H__
89 # include "mpfr-gmp.h"
90 # endif
91 # ifdef MPFR_NEED_LONGLONG_H
92 # include "mpfr-longlong.h"
93 # endif
94
95 #endif
96 #undef MPFR_NEED_LONGLONG_H
97
98 /* For the definition of MPFR_THREAD_ATTR. GCC/ICC detection macros are
99 no longer used, as they sometimes gave incorrect information about
100 the support of thread-local variables. A configure check is now done.
101 If the use of detection macros is needed in the future, this could be
102 moved below (after the detection macros are defined). */
103 #include "mpfr-thread.h"
104
105
106 /******************************************************
107 ***************** Detection macros *******************
108 ******************************************************/
109
110 /* Macros to detect STDC, GCC, GLIBC, GMP and ICC version */
111 #if defined(__STDC_VERSION__)
112 # define __MPFR_STDC(version) (__STDC_VERSION__>=(version))
113 #elif defined(__STDC__)
114 # define __MPFR_STDC(version) (0 == (version))
115 #else
116 # define __MPFR_STDC(version) 0
117 #endif
118
119 #if defined(__ICC)
120 # define __MPFR_ICC(a,b,c) (__ICC >= (a)*100+(b)*10+(c))
121 #elif defined(__INTEL_COMPILER)
122 # define __MPFR_ICC(a,b,c) (__INTEL_COMPILER >= (a)*100+(b)*10+(c))
123 #else
124 # define __MPFR_ICC(a,b,c) 0
125 #endif
126
127 #if defined(__GNUC__) && defined(__GNUC_MINOR__) && ! __MPFR_ICC(0,0,0)
128 # define __MPFR_GNUC(a,i) \
129 (MPFR_VERSION_NUM(__GNUC__,__GNUC_MINOR__,0) >= MPFR_VERSION_NUM(a,i,0))
130 #else
131 # define __MPFR_GNUC(a,i) 0
132 #endif
133
134 #if defined(__GLIBC__) && defined(__GLIBC_MINOR__)
135 # define __MPFR_GLIBC(a,i) \
136 (MPFR_VERSION_NUM(__GLIBC__,__GLIBC_MINOR__,0) >= MPFR_VERSION_NUM(a,i,0))
137 #else
138 # define __MPFR_GLIBC(a,i) 0
139 #endif
140
141 #if defined(__GNU_MP_VERSION) && \
142 defined(__GNU_MP_VERSION_MINOR) && \
143 defined(__GNU_MP_VERSION_PATCHLEVEL)
144 # define __MPFR_GMP(a,b,c) \
145 (MPFR_VERSION_NUM(__GNU_MP_VERSION,__GNU_MP_VERSION_MINOR,__GNU_MP_VERSION_PATCHLEVEL) >= MPFR_VERSION_NUM(a,b,c))
146 #else
147 # define __MPFR_GMP(a,b,c) 0
148 #endif
149
150
151
152 /******************************************************
153 ******************** Check GMP ***********************
154 ******************************************************/
155
156 #if !__MPFR_GMP(4,1,0)
157 # error "GMP 4.1.0 or newer needed"
158 #endif
159
160 #if GMP_NAIL_BITS != 0
161 # error "MPFR doesn't support nonzero values of GMP_NAIL_BITS"
162 #endif
163
164 #if (BITS_PER_MP_LIMB<32) || (BITS_PER_MP_LIMB & (BITS_PER_MP_LIMB - 1))
165 # error "BITS_PER_MP_LIMB must be a power of 2, and >= 32"
166 #endif
167
168 #if BITS_PER_MP_LIMB == 16
169 # define MPFR_LOG2_BITS_PER_MP_LIMB 4
170 #elif BITS_PER_MP_LIMB == 32
171 # define MPFR_LOG2_BITS_PER_MP_LIMB 5
172 #elif BITS_PER_MP_LIMB == 64
173 # define MPFR_LOG2_BITS_PER_MP_LIMB 6
174 #elif BITS_PER_MP_LIMB == 128
175 # define MPFR_LOG2_BITS_PER_MP_LIMB 7
176 #elif BITS_PER_MP_LIMB == 256
177 # define MPFR_LOG2_BITS_PER_MP_LIMB 8
178 #else
179 # error "Can't compute log2(BITS_PER_MP_LIMB)"
180 #endif
181
182 #if __MPFR_GNUC(3,0) || __MPFR_ICC(8,1,0)
183 # define MPFR_NORETURN_ATTR __attribute__ ((noreturn))
184 # define MPFR_CONST_ATTR __attribute__ ((const))
185 #else
186 # define MPFR_NORETURN_ATTR
187 # define MPFR_CONST_ATTR
188 #endif
189
190 /******************************************************
191 ************* Global Internal Variables **************
192 ******************************************************/
193
194 /* Cache struct */
195 struct __gmpfr_cache_s {
196 mpfr_t x;
197 int inexact;
198 int (*func)(mpfr_ptr, mpfr_rnd_t);
199 };
200 typedef struct __gmpfr_cache_s mpfr_cache_t[1];
201
202 #if defined (__cplusplus)
203 extern "C" {
204 #endif
205
206 __MPFR_DECLSPEC extern MPFR_THREAD_ATTR unsigned int __gmpfr_flags;
207 __MPFR_DECLSPEC extern MPFR_THREAD_ATTR mp_exp_t __gmpfr_emin;
208 __MPFR_DECLSPEC extern MPFR_THREAD_ATTR mp_exp_t __gmpfr_emax;
209 __MPFR_DECLSPEC extern MPFR_THREAD_ATTR mp_prec_t __gmpfr_default_fp_bit_precision;
210 __MPFR_DECLSPEC extern MPFR_THREAD_ATTR mpfr_rnd_t __gmpfr_default_rounding_mode;
211 __MPFR_DECLSPEC extern MPFR_THREAD_ATTR mpfr_cache_t __gmpfr_cache_const_pi;
212 __MPFR_DECLSPEC extern MPFR_THREAD_ATTR mpfr_cache_t __gmpfr_cache_const_log2;
213 __MPFR_DECLSPEC extern MPFR_THREAD_ATTR mpfr_cache_t __gmpfr_cache_const_euler;
214 __MPFR_DECLSPEC extern MPFR_THREAD_ATTR mpfr_cache_t __gmpfr_cache_const_catalan;
215
216 #define BASE_MAX 36
217 __MPFR_DECLSPEC extern const __mpfr_struct __gmpfr_l2b[BASE_MAX-1][2];
218
219 /* Note: do not use the following values when they can be outside the
220 current exponent range, e.g. when the exponent range has not been
221 extended yet; under such a condition, they can be used only in
222 mpfr_cmpabs. */
223 __MPFR_DECLSPEC extern const mpfr_t __gmpfr_one;
224 __MPFR_DECLSPEC extern const mpfr_t __gmpfr_two;
225 __MPFR_DECLSPEC extern const mpfr_t __gmpfr_four;
226
227
228 #if defined (__cplusplus)
229 }
230 #endif
231
232 /* Flags of __gmpfr_flags */
233 #define MPFR_FLAGS_UNDERFLOW 1
234 #define MPFR_FLAGS_OVERFLOW 2
235 #define MPFR_FLAGS_NAN 4
236 #define MPFR_FLAGS_INEXACT 8
237 #define MPFR_FLAGS_ERANGE 16
238 #define MPFR_FLAGS_ALL 31
239
240 /* Replace some commun functions for direct access to the global vars */
241 #define mpfr_get_emin() (__gmpfr_emin + 0)
242 #define mpfr_get_emax() (__gmpfr_emax + 0)
243 #define mpfr_get_default_rounding_mode() (__gmpfr_default_rounding_mode + 0)
244 #define mpfr_get_default_prec() (__gmpfr_default_fp_bit_precision + 0)
245
246 #define mpfr_clear_flags() \
247 ((void) (__gmpfr_flags = 0))
248 #define mpfr_clear_underflow() \
249 ((void) (__gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_UNDERFLOW))
250 #define mpfr_clear_overflow() \
251 ((void) (__gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_OVERFLOW))
252 #define mpfr_clear_nanflag() \
253 ((void) (__gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_NAN))
254 #define mpfr_clear_inexflag() \
255 ((void) (__gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_INEXACT))
256 #define mpfr_clear_erangeflag() \
257 ((void) (__gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE))
258 #define mpfr_underflow_p() \
259 ((int) (__gmpfr_flags & MPFR_FLAGS_UNDERFLOW))
260 #define mpfr_overflow_p() \
261 ((int) (__gmpfr_flags & MPFR_FLAGS_OVERFLOW))
262 #define mpfr_nanflag_p() \
263 ((int) (__gmpfr_flags & MPFR_FLAGS_NAN))
264 #define mpfr_inexflag_p() \
265 ((int) (__gmpfr_flags & MPFR_FLAGS_INEXACT))
266 #define mpfr_erangeflag_p() \
267 ((int) (__gmpfr_flags & MPFR_FLAGS_ERANGE))
268
269 /* Testing an exception flag correctly is tricky. There are mainly two
270 pitfalls: First, one needs to remember to clear the corresponding
271 flag, in case it was set before the function call or during some
272 intermediate computations (in practice, one can clear all the flags).
273 Secondly, one needs to test the flag early enough, i.e. before it
274 can be modified by another function. Moreover, it is quite difficult
275 (if not impossible) to reliably check problems with "make check". To
276 avoid these pitfalls, it is recommended to use the following macros.
277 Other use of the exception-flag predicate functions/macros will be
278 detected by mpfrlint.
279 Note: _op can be either a statement or an expression.
280 MPFR_BLOCK_EXCEP should be used only inside a block; it is useful to
281 detect some exception in order to exit the block as soon as possible. */
282 #define MPFR_BLOCK_DECL(_flags) unsigned int _flags
283 #define MPFR_BLOCK(_flags,_op) \
284 do \
285 { \
286 mpfr_clear_flags (); \
287 _op; \
288 (_flags) = __gmpfr_flags; \
289 } \
290 while (0)
291 #define MPFR_BLOCK_TEST(_flags,_f) MPFR_UNLIKELY ((_flags) & (_f))
292 #define MPFR_BLOCK_EXCEP (__gmpfr_flags & (MPFR_FLAGS_UNDERFLOW | \
293 MPFR_FLAGS_OVERFLOW | \
294 MPFR_FLAGS_NAN))
295 /* Let's use a MPFR_ prefix, because e.g. OVERFLOW is defined by glibc's
296 math.h, though this is not a reserved identifier! */
297 #define MPFR_UNDERFLOW(_flags) MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_UNDERFLOW)
298 #define MPFR_OVERFLOW(_flags) MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_OVERFLOW)
299 #define MPFR_NANFLAG(_flags) MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_NAN)
300 #define MPFR_INEXFLAG(_flags) MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_INEXACT)
301 #define MPFR_ERANGEFLAG(_flags) MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_ERANGE)
302
303
304 /******************************************************
305 ******************** Assertions **********************
306 ******************************************************/
307
308 /* Compile with -DWANT_ASSERT to check all assert statements */
309
310 /* Note: do not use GMP macros ASSERT_ALWAYS and ASSERT as they are not
311 expressions, and as a consequence, they cannot be used in a for(),
312 with a comma operator and so on. */
313
314 /* MPFR_ASSERTN(expr): assertions that should always be checked */
315 #define MPFR_ASSERTN(expr) \
316 ((void) ((MPFR_UNLIKELY(expr)) || MPFR_UNLIKELY( (ASSERT_FAIL(expr),0) )))
317
318 /* MPFR_ASSERTD(expr): assertions that should be checked when testing */
319 #ifdef WANT_ASSERT
320 # define MPFR_EXP_CHECK 1
321 # define MPFR_ASSERTD(expr) MPFR_ASSERTN (expr)
322 #else
323 # define MPFR_ASSERTD(expr) ((void) 0)
324 #endif
325
326 /* Code to deal with impossible
327 WARNING: It doesn't use do { } while (0) for Insure++*/
328 #define MPFR_RET_NEVER_GO_HERE() {MPFR_ASSERTN(0); return 0;}
329
330
331 /******************************************************
332 ******************** Warnings ************************
333 ******************************************************/
334
335 /* MPFR_WARNING is no longer useful, but let's keep the macro in case
336 it needs to be used again in the future. */
337
338 #ifdef MPFR_USE_WARNINGS
339 # include <stdlib.h>
340 # define MPFR_WARNING(W) \
341 do \
342 { \
343 char *q = getenv ("MPFR_QUIET"); \
344 if (q == NULL || *q == 0) \
345 fprintf (stderr, "MPFR: %s\n", W); \
346 } \
347 while (0)
348 #else
349 # define MPFR_WARNING(W) ((void) 0)
350 #endif
351
352
353 /******************************************************
354 ****************** double macros *********************
355 ******************************************************/
356
357 /* Definition of constants */
358 #define LOG2 0.69314718055994528622 /* log(2) rounded to zero on 53 bits */
359 #define ALPHA 4.3191365662914471407 /* a+2 = a*log(a), rounded to +infinity */
360 #define EXPM1 0.36787944117144227851 /* exp(-1), rounded to zero */
361
362 /* MPFR_DOUBLE_SPEC = 1 if the C type 'double' corresponds to IEEE-754
363 double precision, 0 if it doesn't, and undefined if one doesn't know.
364 On all the tested machines, MPFR_DOUBLE_SPEC = 1. To have this macro
365 defined here, #include <float.h> is needed. If need be, other values
366 could be defined for other specs (once they are known). */
367 #if !defined(MPFR_DOUBLE_SPEC) && defined(FLT_RADIX) && \
368 defined(DBL_MANT_DIG) && defined(DBL_MIN_EXP) && defined(DBL_MAX_EXP)
369 # if FLT_RADIX == 2 && DBL_MANT_DIG == 53 && \
370 DBL_MIN_EXP == -1021 && DBL_MAX_EXP == 1024
371 # define MPFR_DOUBLE_SPEC 1
372 # else
373 # define MPFR_DOUBLE_SPEC 0
374 # endif
375 #endif
376
377 /* Debug non IEEE floats */
378 #ifdef XDEBUG
379 # undef _GMP_IEEE_FLOATS
380 #endif
381 #ifndef _GMP_IEEE_FLOATS
382 # define _GMP_IEEE_FLOATS 0
383 #endif
384
385 #ifndef IEEE_DBL_MANT_DIG
386 #define IEEE_DBL_MANT_DIG 53
387 #endif
388 #define MPFR_LIMBS_PER_DOUBLE ((IEEE_DBL_MANT_DIG-1)/BITS_PER_MP_LIMB+1)
389
390 /* Visual C++ doesn't support +1.0/.00, -1.0/0.0 and 0.0/0.0
391 at compile time. */
392 #if defined(_MSC_VER) && defined(_WIN32) && (_MSC_VER >= 1200)
393 static double double_zero = 0.0;
394 # define DBL_NAN (double_zero/double_zero)
395 # define DBL_POS_INF ((double) 1.0/double_zero)
396 # define DBL_NEG_INF ((double)-1.0/double_zero)
397 # define DBL_NEG_ZERO (-double_zero)
398 #else
399 # define DBL_POS_INF ((double) 1.0/0.0)
400 # define DBL_NEG_INF ((double)-1.0/0.0)
401 # define DBL_NAN ((double) 0.0/0.0)
402 # define DBL_NEG_ZERO (-0.0)
403 #endif
404
405 /* Note: the argument x must be a lvalue of type double. */
406 #if _GMP_IEEE_FLOATS
407 typedef union ieee_double_extract Ieee_double_extract;
408
409 # define DOUBLE_ISNANorINF(x) (((Ieee_double_extract *)&(x))->s.exp == 0x7ff)
410 # define DOUBLE_ISINF(x) (DOUBLE_ISNANorINF(x) && \
411 (((Ieee_double_extract *)&(x))->s.manl == 0) && \
412 (((Ieee_double_extract *)&(x))->s.manh == 0))
413 # define DOUBLE_ISNAN(x) (DOUBLE_ISNANorINF(x) && \
414 ((((Ieee_double_extract *)&(x))->s.manl != 0) || \
415 (((Ieee_double_extract *)&(x))->s.manh != 0)))
416 #else
417 /* Below, the &(x) == &(x) || &(x) != &(x) allows to make sure that x
418 is a lvalue without (probably) any warning from the compiler. The
419 &(x) != &(x) is needed to avoid a failure under Mac OS X 10.4.11
420 (with Xcode 2.4.1, i.e. the latest one). */
421 # define LVALUE(x) (&(x) == &(x) || &(x) != &(x))
422 # define DOUBLE_ISINF(x) (LVALUE(x) && ((x) > DBL_MAX || (x) < -DBL_MAX))
423 # ifdef MPFR_NANISNAN
424 /* Avoid MIPSpro / IRIX64 / gcc -ffast-math (incorrect) optimizations.
425 The + must not be replaced by a ||. With gcc -ffast-math, NaN is
426 regarded as a positive number or something like that; the second
427 test catches this case. */
428 # define DOUBLE_ISNAN(x) \
429 (LVALUE(x) && !((((x) >= 0.0) + ((x) <= 0.0)) && -(x)*(x) <= 0.0))
430 # else
431 # define DOUBLE_ISNAN(x) (LVALUE(x) && (x) != (x))
432 # endif
433 #endif
434
435 /******************************************************
436 *************** Long double macros *******************
437 ******************************************************/
438
439 /* We try to get the exact value of the precision of long double
440 (provided by the implementation) in order to provide correct
441 rounding in this case (not guaranteed if the C implementation
442 does not have an adequate long double arithmetic). Note that
443 it may be lower than the precision of some numbers that can
444 be represented in a long double; e.g. on FreeBSD/x86, it is
445 53 because the processor is configured to round in double
446 precision (even when using the long double type -- this is a
447 limitation of the x87 arithmetic), and on Mac OS X, it is 106
448 because the implementation is a double-double arithmetic.
449 Otherwise (e.g. in base 10), we get an upper bound of the
450 precision, and correct rounding isn't currently provided.
451 */
452 #if defined(LDBL_MANT_DIG) && FLT_RADIX == 2
453 # define MPFR_LDBL_MANT_DIG LDBL_MANT_DIG
454 #else
455 # define MPFR_LDBL_MANT_DIG \
456 (sizeof(long double)*BITS_PER_MP_LIMB/sizeof(mp_limb_t))
457 #endif
458 #define MPFR_LIMBS_PER_LONG_DOUBLE \
459 ((sizeof(long double)-1)/sizeof(mp_limb_t)+1)
460
461 /* LONGDOUBLE_NAN_ACTION executes the code "action" if x is a NaN. */
462
463 /* On hppa2.0n-hp-hpux10 with the unbundled HP cc, the test x!=x on a NaN
464 has been seen false, meaning NaNs are not detected. This seemed to
465 happen only after other comparisons, not sure what's really going on. In
466 any case we can pick apart the bytes to identify a NaN. */
467 #ifdef HAVE_LDOUBLE_IEEE_QUAD_BIG
468 # define LONGDOUBLE_NAN_ACTION(x, action) \
469 do { \
470 union { \
471 long double ld; \
472 struct { \
473 unsigned int sign : 1; \
474 unsigned int exp : 15; \
475 unsigned int man3 : 16; \
476 unsigned int man2 : 32; \
477 unsigned int man1 : 32; \
478 unsigned int man0 : 32; \
479 } s; \
480 } u; \
481 u.ld = (x); \
482 if (u.s.exp == 0x7FFFL \
483 && (u.s.man0 | u.s.man1 | u.s.man2 | u.s.man3) != 0) \
484 { action; } \
485 } while (0)
486 #endif
487
488 /* Under IEEE rules, NaN is not equal to anything, including itself.
489 "volatile" here stops "cc" on mips64-sgi-irix6.5 from optimizing away
490 x!=x. */
491 #ifndef LONGDOUBLE_NAN_ACTION
492 # define LONGDOUBLE_NAN_ACTION(x, action) \
493 do { \
494 volatile long double __x = LONGDOUBLE_VOLATILE (x); \
495 if ((x) != __x) \
496 { action; } \
497 } while (0)
498 # define WANT_LONGDOUBLE_VOLATILE 1
499 #endif
500
501 /* If we don't have a proper "volatile" then volatile is #defined to empty,
502 in this case call through an external function to stop the compiler
503 optimizing anything. */
504 #ifdef WANT_LONGDOUBLE_VOLATILE
505 # ifdef volatile
506 __MPFR_DECLSPEC long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long double)) MPFR_CONST_ATTR;
507 # define LONGDOUBLE_VOLATILE(x) (__gmpfr_longdouble_volatile (x))
508 # define WANT_GMPFR_LONGDOUBLE_VOLATILE 1
509 # else
510 # define LONGDOUBLE_VOLATILE(x) (x)
511 # endif
512 #endif
513
514 /* Some special case for IEEE_EXT Litle Endian */
515 #if HAVE_LDOUBLE_IEEE_EXT_LITTLE
516
517 typedef union {
518 long double ld;
519 struct {
520 unsigned int manl : 32;
521 unsigned int manh : 32;
522 unsigned int expl : 8 ;
523 unsigned int exph : 7;
524 unsigned int sign : 1;
525 } s;
526 } mpfr_long_double_t;
527
528 /* #undef MPFR_LDBL_MANT_DIG */
529 #undef MPFR_LIMBS_PER_LONG_DOUBLE
530 /* #define MPFR_LDBL_MANT_DIG 64 */
531 #define MPFR_LIMBS_PER_LONG_DOUBLE ((64-1)/BITS_PER_MP_LIMB+1)
532
533 #endif
534
535 /******************************************************
536 *************** _Decimal64 support *******************
537 ******************************************************/
538
539 #ifdef MPFR_WANT_DECIMAL_FLOATS
540 /* to cast between binary64 and decimal64 */
541 union ieee_double_decimal64 { double d; _Decimal64 d64; };
542 #endif
543
544 /******************************************************
545 **************** mpfr_t properties *******************
546 ******************************************************/
547
548 #define MPFR_PREC(x) ((x)->_mpfr_prec)
549 #define MPFR_EXP(x) ((x)->_mpfr_exp)
550 #define MPFR_MANT(x) ((x)->_mpfr_d)
551 #define MPFR_LIMB_SIZE(x) ((MPFR_PREC((x))-1)/BITS_PER_MP_LIMB+1)
552
553 #if _MPFR_PREC_FORMAT == 1
554 # define MPFR_INTPREC_MAX (USHRT_MAX & ~(unsigned int) (BITS_PER_MP_LIMB - 1))
555 #elif _MPFR_PREC_FORMAT == 2
556 # define MPFR_INTPREC_MAX (UINT_MAX & ~(unsigned int) (BITS_PER_MP_LIMB - 1))
557 #elif _MPFR_PREC_FORMAT == 3
558 # define MPFR_INTPREC_MAX (ULONG_MAX & ~(unsigned long) (BITS_PER_MP_LIMB - 1))
559 #else
560 # error "Invalid MPFR Prec format"
561 #endif
562
563
564
565 /******************************************************
566 ***************** exponent limits ********************
567 ******************************************************/
568
569 /* Define limits and unsigned type of exponent. The following definitions
570 * depend on mp_exp_t; if this type changes in GMP, these definitions will
571 * need to be modified (alternatively, a mpfr_exp_t type could be defined).
572 */
573 #if __GMP_MP_SIZE_T_INT == 1
574 typedef unsigned int mpfr_uexp_t;
575 # define MPFR_EXP_MAX (INT_MAX)
576 # define MPFR_EXP_MIN (INT_MIN)
577 #else
578 typedef unsigned long int mpfr_uexp_t;
579 # define MPFR_EXP_MAX (LONG_MAX)
580 # define MPFR_EXP_MIN (LONG_MIN)
581 #endif
582 #ifndef mp_exp_unsigned_t
583 # define mp_exp_unsigned_t mpfr_uexp_t
584 #endif
585
586 #if MPFR_EXP_MIN >= LONG_MIN && MPFR_EXP_MAX <= LONG_MAX
587 typedef long int mpfr_eexp_t;
588 # define mpfr_get_exp_t(x,r) mpfr_get_si((x),(r))
589 # define mpfr_set_exp_t(x,e,r) mpfr_set_si((x),(e),(r))
590 #elif defined (_MPFR_H_HAVE_INTMAX_T)
591 typedef intmax_t mpfr_eexp_t;
592 # define mpfr_get_exp_t(x,r) mpfr_get_sj((x),(r))
593 # define mpfr_set_exp_t(x,e,r) mpfr_set_sj((x),(e),(r))
594 #else
595 # error "Cannot define mpfr_get_exp_t and mpfr_set_exp_t"
596 #endif
597
598 /* Invalid exponent value (to track bugs...) */
599 #define MPFR_EXP_INVALID \
600 ((mp_exp_t) 1 << (BITS_PER_MP_LIMB*sizeof(mp_exp_t)/sizeof(mp_limb_t)-2))
601
602 /* Definition of the exponent limits for MPFR numbers.
603 * These limits are chosen so that if e is such an exponent, then 2e-1 and
604 * 2e+1 are representable. This is useful for intermediate computations,
605 * in particular the multiplication.
606 */
607 #undef MPFR_EMIN_MIN
608 #undef MPFR_EMIN_MAX
609 #undef MPFR_EMAX_MIN
610 #undef MPFR_EMAX_MAX
611 #define MPFR_EMIN_MIN (1-MPFR_EXP_INVALID)
612 #define MPFR_EMIN_MAX (MPFR_EXP_INVALID-1)
613 #define MPFR_EMAX_MIN (1-MPFR_EXP_INVALID)
614 #define MPFR_EMAX_MAX (MPFR_EXP_INVALID-1)
615
616 /* Use MPFR_GET_EXP and MPFR_SET_EXP instead of MPFR_EXP directly,
617 unless when the exponent may be out-of-range, for instance when
618 setting the exponent before calling mpfr_check_range.
619 MPFR_EXP_CHECK is defined when WANT_ASSERT is defined, but if you
620 don't use WANT_ASSERT (for speed reasons), you can still define
621 MPFR_EXP_CHECK by setting -DMPFR_EXP_CHECK in $CFLAGS. */
622
623 #ifdef MPFR_EXP_CHECK
624 # define MPFR_GET_EXP(x) (mpfr_get_exp) (x)
625 # define MPFR_SET_EXP(x, exp) MPFR_ASSERTN (!mpfr_set_exp ((x), (exp)))
626 # define MPFR_SET_INVALID_EXP(x) ((void) (MPFR_EXP (x) = MPFR_EXP_INVALID))
627 #else
628 # define MPFR_GET_EXP(x) MPFR_EXP (x)
629 # define MPFR_SET_EXP(x, exp) ((void) (MPFR_EXP (x) = (exp)))
630 # define MPFR_SET_INVALID_EXP(x) ((void) 0)
631 #endif
632
633
634
635 /******************************************************
636 ********** Singular Values (NAN, INF, ZERO) **********
637 ******************************************************/
638
639 /*
640 * Clear flags macros are still defined and should be still used
641 * since the functions must not assume the internal format.
642 * How to deal with special values ?
643 * 1. Check if is a special value (Zero, Nan, Inf) wiht MPFR_IS_SINGULAR
644 * 2. Deal with the special value with MPFR_IS_NAN, MPFR_IS_INF, etc
645 * 3. Else clear the flags of the dest (it must be done after since src
646 * may be also the dest!)
647 * MPFR_SET_INF, MPFR_SET_NAN, MPFR_SET_ZERO must clear by
648 * themselves the other flags.
649 */
650
651 /* Enum special value of exponent.*/
652 # define MPFR_EXP_ZERO (MPFR_EXP_MIN+1)
653 # define MPFR_EXP_NAN (MPFR_EXP_MIN+2)
654 # define MPFR_EXP_INF (MPFR_EXP_MIN+3)
655
656 #define MPFR_CLEAR_FLAGS(x)
657
658 #define MPFR_IS_NAN(x) (MPFR_EXP(x) == MPFR_EXP_NAN)
659 #define MPFR_SET_NAN(x) (MPFR_EXP(x) = MPFR_EXP_NAN)
660 #define MPFR_IS_INF(x) (MPFR_EXP(x) == MPFR_EXP_INF)
661 #define MPFR_SET_INF(x) (MPFR_EXP(x) = MPFR_EXP_INF)
662 #define MPFR_IS_ZERO(x) (MPFR_EXP(x) == MPFR_EXP_ZERO)
663 #define MPFR_SET_ZERO(x) (MPFR_EXP(x) = MPFR_EXP_ZERO)
664 #define MPFR_NOTZERO(x) (MPFR_EXP(x) != MPFR_EXP_ZERO)
665
666 #define MPFR_IS_FP(x) (!MPFR_IS_NAN(x) && !MPFR_IS_INF(x))
667 #define MPFR_IS_SINGULAR(x) (MPFR_EXP(x) <= MPFR_EXP_INF)
668 #define MPFR_IS_PURE_FP(x) (!MPFR_IS_SINGULAR(x))
669
670 #define MPFR_ARE_SINGULAR(x,y) \
671 (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)) || MPFR_UNLIKELY(MPFR_IS_SINGULAR(y)))
672
673
674
675 /******************************************************
676 ********************* Sign Macros ********************
677 ******************************************************/
678
679 #define MPFR_SIGN_POS (1)
680 #define MPFR_SIGN_NEG (-1)
681
682 #define MPFR_IS_STRICTPOS(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x) > 0)
683 #define MPFR_IS_STRICTNEG(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x) < 0)
684
685 #define MPFR_IS_NEG(x) (MPFR_SIGN(x) < 0)
686 #define MPFR_IS_POS(x) (MPFR_SIGN(x) > 0)
687
688 #define MPFR_SET_POS(x) (MPFR_SIGN(x) = MPFR_SIGN_POS)
689 #define MPFR_SET_NEG(x) (MPFR_SIGN(x) = MPFR_SIGN_NEG)
690
691 #define MPFR_CHANGE_SIGN(x) (MPFR_SIGN(x) = -MPFR_SIGN(x))
692 #define MPFR_SET_SAME_SIGN(x, y) (MPFR_SIGN(x) = MPFR_SIGN(y))
693 #define MPFR_SET_OPPOSITE_SIGN(x, y) (MPFR_SIGN(x) = -MPFR_SIGN(y))
694 #define MPFR_ASSERT_SIGN(s) \
695 (MPFR_ASSERTD((s) == MPFR_SIGN_POS || (s) == MPFR_SIGN_NEG))
696 #define MPFR_SET_SIGN(x, s) \
697 (MPFR_ASSERT_SIGN(s), MPFR_SIGN(x) = s)
698 #define MPFR_IS_POS_SIGN(s1) (s1 > 0)
699 #define MPFR_IS_NEG_SIGN(s1) (s1 < 0)
700 #define MPFR_MULT_SIGN(s1, s2) ((s1) * (s2))
701 /* Transform a sign to 1 or -1 */
702 #define MPFR_FROM_SIGN_TO_INT(s) (s)
703 #define MPFR_INT_SIGN(x) MPFR_FROM_SIGN_TO_INT(MPFR_SIGN(x))
704
705
706
707 /******************************************************
708 ***************** Ternary Value Macros ***************
709 ******************************************************/
710
711 /* Special inexact value */
712 #define MPFR_EVEN_INEX 2
713
714 /* When returning the ternary inexact value, ALWAYS use one of the
715 following two macros, unless the flag comes from another function
716 returning the ternary inexact value */
717 #define MPFR_RET(I) return \
718 (I) ? ((__gmpfr_flags |= MPFR_FLAGS_INEXACT), (I)) : 0
719 #define MPFR_RET_NAN return (__gmpfr_flags |= MPFR_FLAGS_NAN), 0
720
721 #define MPFR_SET_ERANGE() (__gmpfr_flags |= MPFR_FLAGS_ERANGE)
722
723 #define SIGN(I) ((I) < 0 ? -1 : (I) > 0)
724 #define SAME_SIGN(I1,I2) (SIGN (I1) == SIGN (I2))
725
726
727
728 /******************************************************
729 ************** Rounding mode macros *****************
730 ******************************************************/
731
732 /* We want to test this :
733 * (rnd == GMP_RNDU && test) || (rnd == RNDD && !test)
734 * ie it transforms RNDU or RNDD to Away or Zero according to the sign */
735 #define MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd, test) \
736 (((rnd) + (test)) == GMP_RNDD)
737
738 /* We want to test if rnd = Zero, or Away.
739 'test' is true iff negative. */
740 #define MPFR_IS_LIKE_RNDZ(rnd, test) \
741 ((rnd==GMP_RNDZ) || MPFR_IS_RNDUTEST_OR_RNDDNOTTEST (rnd, test))
742
743 /* Invert a rounding mode */
744 #define MPFR_INVERT_RND(rnd) ((rnd == GMP_RNDU) ? GMP_RNDD : \
745 ((rnd == GMP_RNDD) ? GMP_RNDU : rnd))
746
747 /* Transform RNDU and RNDD to RNDA or RNDZ */
748 #define MPFR_UPDATE_RND_MODE(rnd, test) \
749 do { \
750 if (MPFR_UNLIKELY(MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd, test))) \
751 rnd = GMP_RNDZ; \
752 } while (0)
753
754
755
756 /******************************************************
757 ******************* Limb Macros **********************
758 ******************************************************/
759
760 /* Definition of MPFR_LIMB_HIGHBIT */
761 #if defined(GMP_LIMB_HIGHBIT)
762 # define MPFR_LIMB_HIGHBIT GMP_LIMB_HIGHBIT
763 #elif defined(MP_LIMB_T_HIGHBIT)
764 # define MPFR_LIMB_HIGHBIT MP_LIMB_T_HIGHBIT
765 #else
766 # error "Neither GMP_LIMB_HIGHBIT nor MP_LIMB_T_HIGHBIT defined in GMP"
767 #endif
768
769 /* Mask to get the Most Significant Bit of a limb */
770 #define MPFR_LIMB_MSB(l) ((l)&MPFR_LIMB_HIGHBIT)
771
772 /* Definition of MPFR_LIMB_ONE & MPFR_LIMB_ZERO*/
773 #ifdef CNST_LIMB
774 # define MPFR_LIMB_ONE CNST_LIMB(1)
775 # define MPFR_LIMB_ZERO CNST_LIMB(0)
776 #else
777 # define MPFR_LIMB_ONE ((mp_limb_t) 1L)
778 # define MPFR_LIMB_ZERO ((mp_limb_t) 0L)
779 #endif
780
781 /* Mask for the low 's' bits of a limb */
782 #define MPFR_LIMB_MASK(s) ((MPFR_LIMB_ONE<<(s))-MPFR_LIMB_ONE)
783
784
785
786 /******************************************************
787 ********************** Memory ************************
788 ******************************************************/
789
790 /* Heap Memory gestion */
791 typedef union { mp_size_t s; mp_limb_t l; } mpfr_size_limb_t;
792 #define MPFR_GET_ALLOC_SIZE(x) \
793 ( ((mp_size_t*) MPFR_MANT(x))[-1] + 0)
794 #define MPFR_SET_ALLOC_SIZE(x, n) \
795 ( ((mp_size_t*) MPFR_MANT(x))[-1] = n)
796 #define MPFR_MALLOC_SIZE(s) \
797 ( sizeof(mpfr_size_limb_t) + BYTES_PER_MP_LIMB * ((size_t) s) )
798 #define MPFR_SET_MANT_PTR(x,p) \
799 (MPFR_MANT(x) = (mp_limb_t*) ((mpfr_size_limb_t*) p + 1))
800 #define MPFR_GET_REAL_PTR(x) \
801 ((mp_limb_t*) ((mpfr_size_limb_t*) MPFR_MANT(x) - 1))
802
803 /* Temporary memory gestion */
804 #ifndef TMP_SALLOC
805 /* GMP 4.1.x or below or internals */
806 #define MPFR_TMP_DECL TMP_DECL
807 #define MPFR_TMP_MARK TMP_MARK
808 #define MPFR_TMP_ALLOC TMP_ALLOC
809 #define MPFR_TMP_FREE TMP_FREE
810 #else
811 #define MPFR_TMP_DECL(x) TMP_DECL
812 #define MPFR_TMP_MARK(x) TMP_MARK
813 #define MPFR_TMP_ALLOC(s) TMP_ALLOC(s)
814 #define MPFR_TMP_FREE(x) TMP_FREE
815 #endif
816
817 /* This code is experimental: don't use it */
818 #ifdef MPFR_USE_OWN_MPFR_TMP_ALLOC
819 extern unsigned char *mpfr_stack;
820 #undef MPFR_TMP_DECL
821 #undef MPFR_TMP_MARK
822 #undef MPFR_TMP_ALLOC
823 #undef MPFR_TMP_FREE
824 #define MPFR_TMP_DECL(_x) unsigned char *(_x)
825 #define MPFR_TMP_MARK(_x) ((_x) = mpfr_stack)
826 #define MPFR_TMP_ALLOC(_s) (mpfr_stack += (_s), mpfr_stack - (_s))
827 #define MPFR_TMP_FREE(_x) (mpfr_stack = (_x))
828 #endif
829
830 /* temporary allocate 1 limb at xp, and initialize mpfr variable x */
831 /* The temporary var doesn't have any size field, but it doesn't matter
832 * since only functions dealing with the Heap care about it */
833 #define MPFR_TMP_INIT1(xp, x, p) \
834 ( MPFR_PREC(x) = (p), \
835 MPFR_MANT(x) = (xp), \
836 MPFR_SET_POS(x), \
837 MPFR_SET_INVALID_EXP(x))
838
839 #define MPFR_TMP_INIT(xp, x, p, s) \
840 (xp = (mp_ptr) MPFR_TMP_ALLOC(BYTES_PER_MP_LIMB * ((size_t) s)), \
841 MPFR_TMP_INIT1(xp, x, p))
842
843 #define MPFR_TMP_INIT_ABS(d, s) \
844 ( MPFR_PREC(d) = MPFR_PREC(s), \
845 MPFR_MANT(d) = MPFR_MANT(s), \
846 MPFR_SET_POS(d), \
847 MPFR_EXP(d) = MPFR_EXP(s))
848
849
850
851 /******************************************************
852 ***************** Cache macros **********************
853 ******************************************************/
854
855 #define mpfr_const_pi(_d,_r) mpfr_cache(_d, __gmpfr_cache_const_pi,_r)
856 #define mpfr_const_log2(_d,_r) mpfr_cache(_d, __gmpfr_cache_const_log2, _r)
857 #define mpfr_const_euler(_d,_r) mpfr_cache(_d, __gmpfr_cache_const_euler, _r)
858 #define mpfr_const_catalan(_d,_r) mpfr_cache(_d,__gmpfr_cache_const_catalan,_r)
859
860 #define MPFR_DECL_INIT_CACHE(_cache,_func) \
861 mpfr_cache_t MPFR_THREAD_ATTR _cache = \
862 {{{{0,MPFR_SIGN_POS,0,(mp_limb_t*)0}},0,_func}}
863
864
865
866 /******************************************************
867 ******************* Threshold ***********************
868 ******************************************************/
869
870 #include "mparam.h"
871
872 /******************************************************
873 ***************** Useful macros *********************
874 ******************************************************/
875
876 /* Theses macros help the compiler to determine if a test is
877 * likely or unlikely. */
878 #if __MPFR_GNUC(3,0) || __MPFR_ICC(8,1,0)
879 # define MPFR_LIKELY(x) (__builtin_expect(!!(x),1))
880 # define MPFR_UNLIKELY(x) (__builtin_expect((x),0))
881 #else
882 # define MPFR_LIKELY(x) (x)
883 # define MPFR_UNLIKELY(x) (x)
884 #endif
885
886 /* Declare that some variable is initialized before being used (without a
887 dummy initialization) in order to avoid some compiler warnings. Use the
888 VAR = VAR trick (see http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36296)
889 only with gcc as this is undefined behavior, and we don't know what
890 other compilers do (they may also be smarter). This trick could be
891 disabled with future gcc versions. */
892 #if defined(__GNUC__)
893 # define INITIALIZED(VAR) VAR = VAR
894 #else
895 # define INITIALIZED(VAR) VAR
896 #endif
897
898 /* Ceil log 2: If GCC, uses a GCC extension, otherwise calls a function */
899 /* Warning:
900 * Needs to define MPFR_NEED_LONGLONG.
901 * Computes ceil(log2(x)) only for x integer (unsigned long)
902 * Undefined if x is 0 */
903 #if __MPFR_GNUC(2,95) || __MPFR_ICC(8,1,0)
904 # define MPFR_INT_CEIL_LOG2(x) \
905 (MPFR_UNLIKELY ((x) == 1) ? 0 : \
906 __extension__ ({ int _b; mp_limb_t _limb; \
907 MPFR_ASSERTN ((x) > 1); \
908 _limb = (x) - 1; \
909 MPFR_ASSERTN (_limb == (x) - 1); \
910 count_leading_zeros (_b, _limb); \
911 (BITS_PER_MP_LIMB - _b); }))
912 #else
913 # define MPFR_INT_CEIL_LOG2(x) (__gmpfr_int_ceil_log2(x))
914 #endif
915
916 /* Add two integers with overflow handling */
917 /* Example: MPFR_SADD_OVERFLOW (c, a, b, long, unsigned long,
918 * LONG_MIN, LONG_MAX,
919 * goto overflow, goto underflow); */
920 #define MPFR_UADD_OVERFLOW(c,a,b,ACTION_IF_OVERFLOW) \
921 do { \
922 (c) = (a) + (b); \
923 if ((c) < (a)) ACTION_IF_OVERFLOW; \
924 } while (0)
925
926 #define MPFR_SADD_OVERFLOW(c,a,b,STYPE,UTYPE,MIN,MAX,ACTION_IF_POS_OVERFLOW,ACTION_IF_NEG_OVERFLOW) \
927 do { \
928 if ((a) >= 0 && (b) >= 0) { \
929 UTYPE uc,ua,ub; \
930 ua = (UTYPE) a; ub = (UTYPE) b; \
931 MPFR_UADD_OVERFLOW (uc, ua, ub, ACTION_IF_POS_OVERFLOW); \
932 if (uc > (UTYPE)(MAX)) ACTION_IF_POS_OVERFLOW; \
933 else (c) = (STYPE) uc; \
934 } else if ((a) < 0 && (b) < 0) { \
935 UTYPE uc,ua,ub; \
936 ua = -(UTYPE) a; ub = -(UTYPE) b; \
937 MPFR_UADD_OVERFLOW (uc, ua, ub, ACTION_IF_NEG_OVERFLOW); \
938 if (uc >= -(UTYPE)(MIN) || uc > (UTYPE)(MAX)) { \
939 if (uc == -(UTYPE)(MIN)) (c) = (MIN); \
940 else ACTION_IF_NEG_OVERFLOW; } \
941 else (c) = -(STYPE) uc; \
942 } else (c) = (a) + (b); \
943 } while (0)
944
945
946 /* Set a number to 1 (Fast) - It doesn't check if 1 is in the exponent range */
947 #define MPFR_SET_ONE(x) \
948 do { \
949 mp_size_t _size = MPFR_LIMB_SIZE(x) - 1; \
950 MPFR_SET_POS(x); \
951 MPFR_EXP(x) = 1; \
952 MPN_ZERO ( MPFR_MANT(x), _size); \
953 MPFR_MANT(x)[_size] = MPFR_LIMB_HIGHBIT; \
954 } while (0)
955
956 /* Compute s = (-a) % BITS_PER_MP_LIMB
957 * a is unsigned! Check if it works,
958 * otherwise tries another way to compute it */
959 #define MPFR_UNSIGNED_MINUS_MODULO(s, a) \
960 do \
961 { \
962 if (IS_POW2 (BITS_PER_MP_LIMB)) \
963 (s) = (-(a)) % BITS_PER_MP_LIMB; \
964 else \
965 { \
966 (s) = (a) % BITS_PER_MP_LIMB; \
967 if ((s) != 0) \
968 (s) = BITS_PER_MP_LIMB - (s); \
969 } \
970 MPFR_ASSERTD ((s) >= 0 && (s) < BITS_PER_MP_LIMB); \
971 } \
972 while (0)
973
974 /* Use it only for debug reasons */
975 /* MPFR_TRACE (operation) : execute operation iff DEBUG flag is set */
976 /* MPFR_DUMP (x) : print x (a mpfr_t) on stdout */
977 #ifdef DEBUG
978 # define MPFR_TRACE(x) x
979 #else
980 # define MPFR_TRACE(x) (void) 0
981 #endif
982 #define MPFR_DUMP(x) ( printf(#x"="), mpfr_dump(x) )
983
984 /* Test if X (positive) is a power of 2 */
985 #define IS_POW2(X) (((X) & ((X) - 1)) == 0)
986 #define NOT_POW2(X) (((X) & ((X) - 1)) != 0)
987
988 /* Safe absolute value (to avoid possible integer overflow) */
989 /* type is the target (unsigned) type */
990 #define SAFE_ABS(type,x) ((x) >= 0 ? (type)(x) : -(type)(x))
991
992 #define mpfr_get_d1(x) mpfr_get_d(x,__gmpfr_default_rounding_mode)
993
994 /* Store in r the size in bits of the mpz_t z */
995 #define MPFR_MPZ_SIZEINBASE2(r, z) \
996 do { \
997 int _cnt; \
998 mp_size_t _size; \
999 MPFR_ASSERTD (mpz_sgn (z) != 0); \
1000 _size = ABSIZ(z); \
1001 count_leading_zeros (_cnt, PTR(z)[_size-1]); \
1002 (r) = _size * BITS_PER_MP_LIMB - _cnt; \
1003 } while (0)
1004
1005 /* Needs <locale.h> */
1006 #ifdef HAVE_LOCALE_H
1007 #include <locale.h>
1008 /* Warning! In case of signed char, the value of MPFR_DECIMAL_POINT may
1009 be negative (the ISO C99 does not seem to forbid negative values). */
1010 #define MPFR_DECIMAL_POINT (localeconv()->decimal_point[0])
1011 #define MPFR_THOUSANDS_SEPARATOR (localeconv()->thousands_sep[0])
1012 #else
1013 #define MPFR_DECIMAL_POINT ((char) '.')
1014 #define MPFR_THOUSANDS_SEPARATOR ('\0')
1015 #endif
1016
1017
1018 /* Set y to s*significand(x)*2^e, for example MPFR_ALIAS(y,x,1,MPFR_EXP(x))
1019 sets y to |x|, and MPFR_ALIAS(y,x,MPFR_SIGN(x),0) sets y to x*2^f such
1020 that 1/2 <= |y| < 1. Does not check y is in the valid exponent range.
1021 WARNING! x and y share the same mantissa. So, some operations are
1022 not valid if x has been provided via an argument, e.g., trying to
1023 modify the mantissa of y, even temporarily, or calling mpfr_clear on y.
1024 */
1025 #define MPFR_ALIAS(y,x,s,e) \
1026 do \
1027 { \
1028 MPFR_PREC(y) = MPFR_PREC(x); \
1029 MPFR_SIGN(y) = (s); \
1030 MPFR_EXP(y) = (e); \
1031 MPFR_MANT(y) = MPFR_MANT(x); \
1032 } while (0)
1033
1034
1035 /******************************************************
1036 ************** Save exponent macros ****************
1037 ******************************************************/
1038
1039 /* See README.dev for details on how to use the macros.
1040 They are used to set the exponent range to the maximum
1041 temporarily */
1042
1043 typedef struct {
1044 unsigned int saved_flags;
1045 mp_exp_t saved_emin;
1046 mp_exp_t saved_emax;
1047 } mpfr_save_expo_t;
1048
1049 #define MPFR_SAVE_EXPO_DECL(x) mpfr_save_expo_t x
1050 #define MPFR_SAVE_EXPO_MARK(x) \
1051 ((x).saved_flags = __gmpfr_flags, \
1052 (x).saved_emin = __gmpfr_emin, \
1053 (x).saved_emax = __gmpfr_emax, \
1054 __gmpfr_emin = MPFR_EMIN_MIN, \
1055 __gmpfr_emax = MPFR_EMAX_MAX)
1056 #define MPFR_SAVE_EXPO_FREE(x) \
1057 (__gmpfr_flags = (x).saved_flags, \
1058 __gmpfr_emin = (x).saved_emin, \
1059 __gmpfr_emax = (x).saved_emax)
1060 #define MPFR_SAVE_EXPO_UPDATE_FLAGS(x, flags) \
1061 (x).saved_flags |= (flags)
1062
1063 /* Speed up final checking */
1064 #define mpfr_check_range(x,t,r) \
1065 (MPFR_LIKELY (MPFR_EXP (x) >= __gmpfr_emin && MPFR_EXP (x) <= __gmpfr_emax) \
1066 ? ((t) ? (__gmpfr_flags |= MPFR_FLAGS_INEXACT, (t)) : 0) \
1067 : mpfr_check_range(x,t,r))
1068
1069
1070 /******************************************************
1071 ***************** Inline Rounding *******************
1072 ******************************************************/
1073
1074 /*
1075 * Note: due to the labels, one cannot use a macro MPFR_RNDRAW* more than
1076 * once in a function (otherwise these labels would not be unique).
1077 */
1078
1079 /*
1080 * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
1081 * assuming dest's sign is sign.
1082 * In rounding to nearest mode, execute MIDDLE_HANDLER when the value
1083 * is the middle of two consecutive numbers in dest precision.
1084 * Execute OVERFLOW_HANDLER in case of overflow when rounding.
1085 */
1086 #define MPFR_RNDRAW_GEN(inexact, dest, srcp, sprec, rnd, sign, \
1087 MIDDLE_HANDLER, OVERFLOW_HANDLER) \
1088 do { \
1089 mp_size_t _dests, _srcs; \
1090 mp_limb_t *_destp; \
1091 mp_prec_t _destprec, _srcprec; \
1092 \
1093 /* Check Trivial Case when Dest Mantissa has more bits than source */ \
1094 _srcprec = sprec; \
1095 _destprec = MPFR_PREC (dest); \
1096 _destp = MPFR_MANT (dest); \
1097 if (MPFR_UNLIKELY (_destprec >= _srcprec)) \
1098 { \
1099 _srcs = (_srcprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; \
1100 _dests = (_destprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB - _srcs; \
1101 MPN_COPY (_destp + _dests, srcp, _srcs); \
1102 MPN_ZERO (_destp, _dests); \
1103 inexact = 0; \
1104 } \
1105 else \
1106 { \
1107 /* Non trivial case: rounding needed */ \
1108 mp_prec_t _sh; \
1109 mp_limb_t *_sp; \
1110 mp_limb_t _rb, _sb, _ulp; \
1111 \
1112 /* Compute Position and shift */ \
1113 _srcs = (_srcprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; \
1114 _dests = (_destprec + BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB; \
1115 MPFR_UNSIGNED_MINUS_MODULO (_sh, _destprec); \
1116 _sp = srcp + _srcs - _dests; \
1117 \
1118 /* General case when prec % BITS_PER_MP_LIMB != 0 */ \
1119 if (MPFR_LIKELY (_sh != 0)) \
1120 { \
1121 mp_limb_t _mask; \
1122 /* Compute Rounding Bit and Sticky Bit */ \
1123 _mask = MPFR_LIMB_ONE << (_sh - 1); \
1124 _rb = _sp[0] & _mask; \
1125 _sb = _sp[0] & (_mask - 1); \
1126 if (MPFR_UNLIKELY (_sb == 0)) \
1127 { /* TODO: Improve it */ \
1128 mp_limb_t *_tmp; \
1129 mp_size_t _n; \
1130 for (_tmp = _sp, _n = _srcs - _dests ; \
1131 _n != 0 && _sb == 0 ; _n--) \
1132 _sb = *--_tmp; \
1133 } \
1134 _ulp = 2 * _mask; \
1135 } \
1136 else /* _sh == 0 */ \
1137 { \
1138 MPFR_ASSERTD (_dests < _srcs); \
1139 /* Compute Rounding Bit and Sticky Bit */ \
1140 _rb = _sp[-1] & MPFR_LIMB_HIGHBIT; \
1141 _sb = _sp[-1] & (MPFR_LIMB_HIGHBIT-1); \
1142 if (MPFR_UNLIKELY (_sb == 0)) \
1143 { \
1144 mp_limb_t *_tmp; \
1145 mp_size_t _n; \
1146 for (_tmp = _sp - 1, _n = _srcs - _dests - 1 ; \
1147 _n != 0 && _sb == 0 ; _n--) \
1148 _sb = *--_tmp; \
1149 } \
1150 _ulp = MPFR_LIMB_ONE; \
1151 } \
1152 /* Rounding */ \
1153 if (MPFR_LIKELY (rnd == GMP_RNDN)) \
1154 { \
1155 if (_rb == 0) \
1156 { \
1157 trunc: \
1158 inexact = MPFR_LIKELY ((_sb | _rb) != 0) ? -sign : 0; \
1159 trunc_doit: \
1160 MPN_COPY (_destp, _sp, _dests); \
1161 _destp[0] &= ~(_ulp - 1); \
1162 } \
1163 else if (MPFR_UNLIKELY (_sb == 0)) \
1164 { /* Middle of two consecutive representable numbers */ \
1165 MIDDLE_HANDLER; \
1166 } \
1167 else \
1168 { \
1169 if (0) \
1170 goto addoneulp_doit; /* dummy code to avoid warning */ \
1171 addoneulp: \
1172 inexact = sign; \
1173 addoneulp_doit: \
1174 if (MPFR_UNLIKELY (mpn_add_1 (_destp, _sp, _dests, _ulp))) \
1175 { \
1176 _destp[_dests - 1] = MPFR_LIMB_HIGHBIT; \
1177 OVERFLOW_HANDLER; \
1178 } \
1179 _destp[0] &= ~(_ulp - 1); \
1180 } \
1181 } \
1182 else \
1183 { /* Directed rounding mode */ \
1184 if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd, \
1185 MPFR_IS_NEG_SIGN (sign)))) \
1186 goto trunc; \
1187 else if (MPFR_UNLIKELY ((_sb | _rb) == 0)) \
1188 { \
1189 inexact = 0; \
1190 goto trunc_doit; \
1191 } \
1192 else \
1193 goto addoneulp; \
1194 } \
1195 } \
1196 } while (0)
1197
1198 /*
1199 * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
1200 * assuming dest's sign is sign.
1201 * Execute OVERFLOW_HANDLER in case of overflow when rounding.
1202 */
1203 #define MPFR_RNDRAW(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER) \
1204 MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign, \
1205 if ((_sp[0] & _ulp) == 0) \
1206 { \
1207 inexact = -sign; \
1208 goto trunc_doit; \
1209 } \
1210 else \
1211 goto addoneulp; \
1212 , OVERFLOW_HANDLER)
1213
1214 /*
1215 * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
1216 * assuming dest's sign is sign.
1217 * Execute OVERFLOW_HANDLER in case of overflow when rounding.
1218 * Set inexact to +/- MPFR_EVEN_INEX in case of even rounding.
1219 */
1220 #define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, \
1221 OVERFLOW_HANDLER) \
1222 MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign, \
1223 if ((_sp[0] & _ulp) == 0) \
1224 { \
1225 inexact = -MPFR_EVEN_INEX * sign; \
1226 goto trunc_doit; \
1227 } \
1228 else \
1229 { \
1230 inexact = MPFR_EVEN_INEX * sign; \
1231 goto addoneulp_doit; \
1232 } \
1233 , OVERFLOW_HANDLER)
1234
1235 /* Return TRUE if b is non singular and we can round it to precision 'prec'
1236 and determine the ternary value, with rounding mode 'rnd', and with
1237 error at most 'error' */
1238 #define MPFR_CAN_ROUND(b,err,prec,rnd) \
1239 (!MPFR_IS_SINGULAR (b) && mpfr_round_p (MPFR_MANT (b), MPFR_LIMB_SIZE (b), \
1240 (err), (prec) + ((rnd)==GMP_RNDN)))
1241
1242 /* TODO: fix this description (see round_near_x.c). */
1243 /* Assuming that the function has a Taylor expansion which looks like:
1244 y=o(f(x)) = o(v + g(x)) with |g(x)| <= 2^(EXP(v)-err)
1245 we can quickly set y to v if x is small (ie err > prec(y)+1) in most
1246 cases. It assumes that f(x) is not representable exactly as a FP number.
1247 v must not be a singular value (NAN, INF or ZERO); usual values are
1248 v=1 or v=x.
1249
1250 y is the destination (a mpfr_t), v the value to set (a mpfr_t),
1251 err1+err2 with err2 <= 3 the error term (mp_exp_t's), dir (an int) is
1252 the direction of the committed error (if dir = 0, it rounds towards 0,
1253 if dir=1, it rounds away from 0), rnd the rounding mode.
1254
1255 It returns from the function a ternary value in case of success.
1256 If you want to free something, you must fill the "extra" field
1257 in consequences, otherwise put nothing in it.
1258
1259 The test is less restrictive than necessary, but the function
1260 will finish the check itself.
1261
1262 Note: err1 + err2 is allowed to overflow as mp_exp_t, but it must give
1263 its real value as mpfr_uexp_t.
1264 */
1265 #define MPFR_FAST_COMPUTE_IF_SMALL_INPUT(y,v,err1,err2,dir,rnd,extra) \
1266 do { \
1267 mpfr_ptr _y = (y); \
1268 mp_exp_t _err1 = (err1); \
1269 mp_exp_t _err2 = (err2); \
1270 if (_err1 > 0) \
1271 { \
1272 mpfr_uexp_t _err = (mpfr_uexp_t) _err1 + _err2; \
1273 if (MPFR_UNLIKELY (_err > MPFR_PREC (_y) + 1)) \
1274 { \
1275 int _inexact = mpfr_round_near_x (_y,(v),_err,(dir),(rnd)); \
1276 if (_inexact != 0) \
1277 { \
1278 extra; \
1279 return _inexact; \
1280 } \
1281 } \
1282 } \
1283 } while (0)
1284
1285 /* Variant, to be called somewhere after MPFR_SAVE_EXPO_MARK. This variant
1286 is needed when there are some computations before or when some non-zero
1287 real constant is used, such as __gmpfr_one for mpfr_cos. */
1288 #define MPFR_SMALL_INPUT_AFTER_SAVE_EXPO(y,v,err1,err2,dir,rnd,expo,extra) \
1289 do { \
1290 mpfr_ptr _y = (y); \
1291 mp_exp_t _err1 = (err1); \
1292 mp_exp_t _err2 = (err2); \
1293 if (_err1 > 0) \
1294 { \
1295 mpfr_uexp_t _err = (mpfr_uexp_t) _err1 + _err2; \
1296 if (MPFR_UNLIKELY (_err > MPFR_PREC (_y) + 1)) \
1297 { \
1298 int _inexact; \
1299 mpfr_clear_flags (); \
1300 _inexact = mpfr_round_near_x (_y,(v),_err,(dir),(rnd)); \
1301 if (_inexact != 0) \
1302 { \
1303 extra; \
1304 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); \
1305 MPFR_SAVE_EXPO_FREE (expo); \
1306 return mpfr_check_range (_y, _inexact, (rnd)); \
1307 } \
1308 } \
1309 } \
1310 } while (0)
1311
1312 /******************************************************
1313 *************** Ziv Loop Macro *********************
1314 ******************************************************/
1315
1316 #ifndef MPFR_USE_LOGGING
1317
1318 #define MPFR_ZIV_DECL(_x) mp_prec_t _x
1319 #define MPFR_ZIV_INIT(_x, _p) (_x) = BITS_PER_MP_LIMB
1320 #define MPFR_ZIV_NEXT(_x, _p) ((_p) += (_x), (_x) = (_p)/2)
1321 #define MPFR_ZIV_FREE(x)
1322
1323 #else
1324
1325 /* The following test on glibc is there mainly for Darwin (Mac OS X), to
1326 obtain a better error message. The real test should have been a test
1327 concerning nested functions in gcc, which are disabled by default on
1328 Darwin; but it is not possible to do that without a configure test. */
1329 # if defined (__cplusplus) || !(__MPFR_GNUC(3,0) && __MPFR_GLIBC(2,0))
1330 # error "Logging not supported (needs gcc >= 3.0 and GNU C Library >= 2.0)."
1331 # endif
1332
1333 /* Use LOGGING */
1334 #define MPFR_ZIV_DECL(_x) \
1335 mp_prec_t _x; \
1336 int _x ## _cpt = 1; \
1337 static unsigned long _x ## _loop = 0, _x ## _bad = 0; \
1338 static const char *_x ## _fname = __func__; \
1339 auto void __attribute__ ((destructor)) x ## _f (void); \
1340 void __attribute__ ((destructor)) x ## _f (void) { \
1341 if (_x ## _loop != 0 && MPFR_LOG_STAT_F&mpfr_log_type) \
1342 fprintf (mpfr_log_file, \
1343 "%s: Ziv failed %2.2f%% (%lu bad cases / %lu calls)\n", _x ## _fname, \
1344 (double) 100.0 * _x ## _bad / _x ## _loop, _x ## _bad, _x ## _loop ); }
1345
1346 #define MPFR_ZIV_INIT(_x, _p) ((_x) = BITS_PER_MP_LIMB, _x ## _loop ++); \
1347 if (MPFR_LOG_BADCASE_F&mpfr_log_type && mpfr_log_current<=mpfr_log_level) \
1348 fprintf (mpfr_log_file, "%s:ZIV 1st prec=%lu\n", __func__, \
1349 (unsigned long) (_p))
1350
1351 #define MPFR_ZIV_NEXT(_x, _p) \
1352 ((_p)+=(_x),(_x)=(_p)/2, _x ## _bad += (_x ## _cpt == 1), _x ## _cpt ++); \
1353 if (MPFR_LOG_BADCASE_F&mpfr_log_type && mpfr_log_current<=mpfr_log_level) \
1354 fprintf (mpfr_log_file, "%s:ZIV new prec=%lu\n", __func__, \
1355 (unsigned long) (_p))
1356
1357 #define MPFR_ZIV_FREE(_x) \
1358 if (MPFR_LOG_BADCASE_F&mpfr_log_type && _x##_cpt>1 \
1359 && mpfr_log_current<=mpfr_log_level) \
1360 fprintf (mpfr_log_file, "%s:ZIV %d loops\n", __func__, _x ## _cpt)
1361
1362 #endif
1363
1364
1365 /******************************************************
1366 *************** Logging Macros *********************
1367 ******************************************************/
1368
1369 /* The different kind of LOG */
1370 #define MPFR_LOG_INPUT_F 1
1371 #define MPFR_LOG_OUTPUT_F 2
1372 #define MPFR_LOG_INTERNAL_F 4
1373 #define MPFR_LOG_TIME_F 8
1374 #define MPFR_LOG_BADCASE_F 16
1375 #define MPFR_LOG_MSG_F 32
1376 #define MPFR_LOG_STAT_F 64
1377
1378 #ifdef MPFR_USE_LOGGING
1379
1380 /* Check if we can support this feature */
1381 # ifdef MPFR_USE_THREAD_SAFE
1382 # error "Enable either `Logging' or `thread-safe', not both"
1383 # endif
1384 # if !__MPFR_GNUC(3,0)
1385 # error "Logging not supported (GCC >= 3.0)"
1386 # endif
1387
1388 #if defined (__cplusplus)
1389 extern "C" {
1390 #endif
1391
1392 __MPFR_DECLSPEC extern FILE *mpfr_log_file;
1393 __MPFR_DECLSPEC extern int mpfr_log_type;
1394 __MPFR_DECLSPEC extern int mpfr_log_level;
1395 __MPFR_DECLSPEC extern int mpfr_log_current;
1396 __MPFR_DECLSPEC extern int mpfr_log_base;
1397 __MPFR_DECLSPEC extern mp_prec_t mpfr_log_prec;
1398
1399 #if defined (__cplusplus)
1400 }
1401 #endif
1402
1403 #define MPFR_LOG_VAR(x) \
1404 if((MPFR_LOG_INTERNAL_F&mpfr_log_type)&&(mpfr_log_current<=mpfr_log_level))\
1405 fprintf (mpfr_log_file, "%s.%d:%s[%#R]=%R\n", __func__,__LINE__, #x, x, x);
1406
1407 #define MPFR_LOG_MSG2(format, ...) \
1408 if ((MPFR_LOG_MSG_F&mpfr_log_type)&&(mpfr_log_current<=mpfr_log_level)) \
1409 fprintf (mpfr_log_file, "%s.%d: "format, __func__, __LINE__, __VA_ARGS__);
1410 #define MPFR_LOG_MSG(x) MPFR_LOG_MSG2 x
1411
1412 #define MPFR_LOG_BEGIN2(format, ...) \
1413 mpfr_log_current ++; \
1414 if ((MPFR_LOG_INPUT_F&mpfr_log_type)&&(mpfr_log_current<=mpfr_log_level)) \
1415 fprintf (mpfr_log_file, "%s:IN "format"\n",__func__,__VA_ARGS__); \
1416 if ((MPFR_LOG_TIME_F&mpfr_log_type)&&(mpfr_log_current<=mpfr_log_level)) \
1417 __gmpfr_log_time = mpfr_get_cputime ();
1418 #define MPFR_LOG_BEGIN(x) \
1419 int __gmpfr_log_time = 0; \
1420 MPFR_LOG_BEGIN2 x
1421
1422 #define MPFR_LOG_END2(format, ...) \
1423 if ((MPFR_LOG_TIME_F&mpfr_log_type)&&(mpfr_log_current<=mpfr_log_level)) \
1424 fprintf (mpfr_log_file, "%s:TIM %dms\n", __mpfr_log_fname, \
1425 mpfr_get_cputime () - __gmpfr_log_time); \
1426 if ((MPFR_LOG_OUTPUT_F&mpfr_log_type)&&(mpfr_log_current<=mpfr_log_level)) \
1427 fprintf (mpfr_log_file, "%s:OUT "format"\n",__mpfr_log_fname,__VA_ARGS__);\
1428 mpfr_log_current --;
1429 #define MPFR_LOG_END(x) \
1430 static const char *__mpfr_log_fname = __func__; \
1431 MPFR_LOG_END2 x
1432
1433 #define MPFR_LOG_FUNC(begin,end) \
1434 static const char *__mpfr_log_fname = __func__; \
1435 auto void __mpfr_log_cleanup (int *time); \
1436 void __mpfr_log_cleanup (int *time) { \
1437 int __gmpfr_log_time = *time; \
1438 MPFR_LOG_END2 end; } \
1439 int __gmpfr_log_time __attribute__ ((cleanup (__mpfr_log_cleanup))); \
1440 __gmpfr_log_time = 0; \
1441 MPFR_LOG_BEGIN2 begin
1442
1443 #else /* MPFR_USE_LOGGING */
1444
1445 /* Define void macro for logging */
1446
1447 #define MPFR_LOG_VAR(x)
1448 #define MPFR_LOG_BEGIN(x)
1449 #define MPFR_LOG_END(x)
1450 #define MPFR_LOG_MSG(x)
1451 #define MPFR_LOG_FUNC(x,y)
1452
1453 #endif /* MPFR_USE_LOGGING */
1454
1455
1456 /**************************************************************
1457 ************ Group Initialize Functions Macros *************
1458 **************************************************************/
1459
1460 #ifndef MPFR_GROUP_STATIC_SIZE
1461 # define MPFR_GROUP_STATIC_SIZE 16
1462 #endif
1463
1464 struct mpfr_group_t {
1465 size_t alloc;
1466 mp_limb_t *mant;
1467 mp_limb_t tab[MPFR_GROUP_STATIC_SIZE];
1468 };
1469
1470 #define MPFR_GROUP_DECL(g) struct mpfr_group_t g
1471 #define MPFR_GROUP_CLEAR(g) do { \
1472 if (MPFR_UNLIKELY ((g).alloc != 0)) { \
1473 MPFR_ASSERTD ((g).mant != (g).tab); \
1474 (*__gmp_free_func) ((g).mant, (g).alloc); \
1475 }} while (0)
1476
1477 #define MPFR_GROUP_INIT_TEMPLATE(g, prec, num, handler) do { \
1478 mp_prec_t _prec = (prec); \
1479 mp_size_t _size; \
1480 MPFR_ASSERTD (_prec >= MPFR_PREC_MIN); \
1481 if (MPFR_UNLIKELY (_prec > MPFR_PREC_MAX)) \
1482 mpfr_abort_prec_max (); \
1483 _size = (mp_prec_t) (_prec + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB; \
1484 if (MPFR_UNLIKELY (_size * (num) > MPFR_GROUP_STATIC_SIZE)) \
1485 { \
1486 (g).alloc = (num) * _size * sizeof (mp_limb_t); \
1487 (g).mant = (mp_limb_t *) (*__gmp_allocate_func) ((g).alloc); \
1488 } \
1489 else \
1490 { \
1491 (g).alloc = 0; \
1492 (g).mant = (g).tab; \
1493 } \
1494 handler; \
1495 } while (0)
1496 #define MPFR_GROUP_TINIT(g, n, x) \
1497 MPFR_TMP_INIT1 ((g).mant + _size * (n), x, _prec)
1498
1499 #define MPFR_GROUP_INIT_1(g, prec, x) \
1500 MPFR_GROUP_INIT_TEMPLATE(g, prec, 1, MPFR_GROUP_TINIT(g, 0, x))
1501 #define MPFR_GROUP_INIT_2(g, prec, x, y) \
1502 MPFR_GROUP_INIT_TEMPLATE(g, prec, 2, \
1503 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y))
1504 #define MPFR_GROUP_INIT_3(g, prec, x, y, z) \
1505 MPFR_GROUP_INIT_TEMPLATE(g, prec, 3, \
1506 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y); \
1507 MPFR_GROUP_TINIT(g, 2, z))
1508 #define MPFR_GROUP_INIT_4(g, prec, x, y, z, t) \
1509 MPFR_GROUP_INIT_TEMPLATE(g, prec, 4, \
1510 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y); \
1511 MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t))
1512 #define MPFR_GROUP_INIT_5(g, prec, x, y, z, t, a) \
1513 MPFR_GROUP_INIT_TEMPLATE(g, prec, 5, \
1514 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y); \
1515 MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t); \
1516 MPFR_GROUP_TINIT(g, 4, a))
1517 #define MPFR_GROUP_INIT_6(g, prec, x, y, z, t, a, b) \
1518 MPFR_GROUP_INIT_TEMPLATE(g, prec, 6, \
1519 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y); \
1520 MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t); \
1521 MPFR_GROUP_TINIT(g, 4, a);MPFR_GROUP_TINIT(g, 5, b))
1522
1523 #define MPFR_GROUP_REPREC_TEMPLATE(g, prec, num, handler) do { \
1524 mp_prec_t _prec = (prec); \
1525 size_t _oalloc = (g).alloc; \
1526 mp_size_t _size; \
1527 MPFR_ASSERTD (_prec >= MPFR_PREC_MIN); \
1528 if (MPFR_UNLIKELY (_prec > MPFR_PREC_MAX)) \
1529 mpfr_abort_prec_max (); \
1530 _size = (mp_prec_t) (_prec + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB; \
1531 (g).alloc = (num) * _size * sizeof (mp_limb_t); \
1532 if (MPFR_LIKELY (_oalloc == 0)) \
1533 (g).mant = (mp_limb_t *) (*__gmp_allocate_func) ((g).alloc); \
1534 else \
1535 (g).mant = (mp_limb_t *) \
1536 (*__gmp_reallocate_func) ((g).mant, _oalloc, (g).alloc); \
1537 handler; \
1538 } while (0)
1539
1540 #define MPFR_GROUP_REPREC_1(g, prec, x) \
1541 MPFR_GROUP_REPREC_TEMPLATE(g, prec, 1, MPFR_GROUP_TINIT(g, 0, x))
1542 #define MPFR_GROUP_REPREC_2(g, prec, x, y) \
1543 MPFR_GROUP_REPREC_TEMPLATE(g, prec, 2, \
1544 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y))
1545 #define MPFR_GROUP_REPREC_3(g, prec, x, y, z) \
1546 MPFR_GROUP_REPREC_TEMPLATE(g, prec, 3, \
1547 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y); \
1548 MPFR_GROUP_TINIT(g, 2, z))
1549 #define MPFR_GROUP_REPREC_4(g, prec, x, y, z, t) \
1550 MPFR_GROUP_REPREC_TEMPLATE(g, prec, 4, \
1551 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y); \
1552 MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t))
1553 #define MPFR_GROUP_REPREC_5(g, prec, x, y, z, t, a) \
1554 MPFR_GROUP_REPREC_TEMPLATE(g, prec, 5, \
1555 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y); \
1556 MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t); \
1557 MPFR_GROUP_TINIT(g, 4, a))
1558 #define MPFR_GROUP_REPREC_6(g, prec, x, y, z, t, a, b) \
1559 MPFR_GROUP_REPREC_TEMPLATE(g, prec, 6, \
1560 MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y); \
1561 MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t); \
1562 MPFR_GROUP_TINIT(g, 4, a);MPFR_GROUP_TINIT(g, 5, b))
1563
1564
1565 /******************************************************
1566 *************** Internal Functions *****************
1567 ******************************************************/
1568
1569 #if defined (__cplusplus)
1570 extern "C" {
1571 #endif
1572
1573 __MPFR_DECLSPEC int mpfr_underflow _MPFR_PROTO ((mpfr_ptr, mp_rnd_t, int));
1574 __MPFR_DECLSPEC int mpfr_overflow _MPFR_PROTO ((mpfr_ptr, mp_rnd_t, int));
1575
1576 __MPFR_DECLSPEC int mpfr_add1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
1577 mpfr_srcptr, mp_rnd_t));
1578 __MPFR_DECLSPEC int mpfr_sub1 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
1579 mpfr_srcptr, mp_rnd_t));
1580 __MPFR_DECLSPEC int mpfr_add1sp _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
1581 mpfr_srcptr, mp_rnd_t));
1582 __MPFR_DECLSPEC int mpfr_sub1sp _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
1583 mpfr_srcptr, mp_rnd_t));
1584 __MPFR_DECLSPEC int mpfr_can_round_raw _MPFR_PROTO ((const mp_limb_t *,
1585 mp_size_t, int, mp_exp_t, mp_rnd_t, mp_rnd_t, mp_prec_t));
1586
1587 __MPFR_DECLSPEC int mpfr_cmp2 _MPFR_PROTO ((mpfr_srcptr, mpfr_srcptr,
1588 mp_prec_t *));
1589
1590 __MPFR_DECLSPEC long __gmpfr_ceil_log2 _MPFR_PROTO ((double));
1591 __MPFR_DECLSPEC long __gmpfr_floor_log2 _MPFR_PROTO ((double));
1592 __MPFR_DECLSPEC double __gmpfr_ceil_exp2 _MPFR_PROTO ((double));
1593 __MPFR_DECLSPEC unsigned long __gmpfr_isqrt _MPFR_PROTO ((unsigned long));
1594 __MPFR_DECLSPEC unsigned long __gmpfr_cuberoot _MPFR_PROTO ((unsigned long));
1595 __MPFR_DECLSPEC int __gmpfr_int_ceil_log2 _MPFR_PROTO ((unsigned long));
1596
1597 __MPFR_DECLSPEC int mpfr_exp_2 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mp_rnd_t));
1598 __MPFR_DECLSPEC int mpfr_exp_3 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mp_rnd_t));
1599 __MPFR_DECLSPEC int mpfr_powerof2_raw _MPFR_PROTO ((mpfr_srcptr));
1600
1601 __MPFR_DECLSPEC int mpfr_pow_general _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
1602 mpfr_srcptr, mp_rnd_t, int, mpfr_save_expo_t *));
1603
1604 __MPFR_DECLSPEC void mpfr_setmax _MPFR_PROTO ((mpfr_ptr, mp_exp_t));
1605 __MPFR_DECLSPEC void mpfr_setmin _MPFR_PROTO ((mpfr_ptr, mp_exp_t));
1606
1607 __MPFR_DECLSPEC long mpfr_mpn_exp _MPFR_PROTO ((mp_limb_t *, mp_exp_t *, int,
1608 mp_exp_t, size_t));
1609
1610 #ifdef _MPFR_H_HAVE_FILE
1611 __MPFR_DECLSPEC void mpfr_fprint_binary _MPFR_PROTO ((FILE *, mpfr_srcptr));
1612 #endif
1613 __MPFR_DECLSPEC void mpfr_print_binary _MPFR_PROTO ((mpfr_srcptr));
1614 __MPFR_DECLSPEC void mpfr_print_mant_binary _MPFR_PROTO ((const char*,
1615 const mp_limb_t*, mp_prec_t));
1616 __MPFR_DECLSPEC void mpfr_set_str_binary _MPFR_PROTO((mpfr_ptr, const char*));
1617
1618 __MPFR_DECLSPEC int mpfr_round_raw _MPFR_PROTO ((mp_limb_t *,
1619 const mp_limb_t *, mp_prec_t, int, mp_prec_t, mp_rnd_t, int *));
1620 __MPFR_DECLSPEC int mpfr_round_raw_2 _MPFR_PROTO ((const mp_limb_t *,
1621 mp_prec_t, int, mp_prec_t, mp_rnd_t));
1622 __MPFR_DECLSPEC int mpfr_round_raw_3 _MPFR_PROTO ((const mp_limb_t *,
1623 mp_prec_t, int, mp_prec_t, mp_rnd_t, int *));
1624 __MPFR_DECLSPEC int mpfr_round_raw_4 _MPFR_PROTO ((mp_limb_t *,
1625 const mp_limb_t *, mp_prec_t, int, mp_prec_t, mp_rnd_t));
1626
1627 #define mpfr_round_raw2(xp, xn, neg, r, prec) \
1628 mpfr_round_raw_2((xp),(xn)*BITS_PER_MP_LIMB,(neg),(prec),(r))
1629
1630 __MPFR_DECLSPEC int mpfr_check _MPFR_PROTO ((mpfr_srcptr));
1631
1632 __MPFR_DECLSPEC int mpfr_sum_sort _MPFR_PROTO ((mpfr_srcptr *const,
1633 unsigned long, mpfr_srcptr *));
1634
1635 __MPFR_DECLSPEC int mpfr_get_cputime _MPFR_PROTO ((void));
1636
1637 __MPFR_DECLSPEC void mpfr_nexttozero _MPFR_PROTO ((mpfr_ptr));
1638 __MPFR_DECLSPEC void mpfr_nexttoinf _MPFR_PROTO ((mpfr_ptr));
1639
1640 __MPFR_DECLSPEC int mpfr_const_pi_internal _MPFR_PROTO ((mpfr_ptr,mp_rnd_t));
1641 __MPFR_DECLSPEC int mpfr_const_log2_internal _MPFR_PROTO((mpfr_ptr,mp_rnd_t));
1642 __MPFR_DECLSPEC int mpfr_const_euler_internal _MPFR_PROTO((mpfr_ptr, mp_rnd_t));
1643 __MPFR_DECLSPEC int mpfr_const_catalan_internal _MPFR_PROTO((mpfr_ptr, mp_rnd_t));
1644
1645 __MPFR_DECLSPEC void mpfr_init_cache _MPFR_PROTO ((mpfr_cache_t,
1646 int(*)(mpfr_ptr,mpfr_rnd_t)));
1647 __MPFR_DECLSPEC void mpfr_clear_cache _MPFR_PROTO ((mpfr_cache_t));
1648 __MPFR_DECLSPEC int mpfr_cache _MPFR_PROTO ((mpfr_ptr, mpfr_cache_t,
1649 mpfr_rnd_t));
1650
1651 __MPFR_DECLSPEC void mpfr_mulhigh_n _MPFR_PROTO ((mp_ptr, mp_srcptr,
1652 mp_srcptr, mp_size_t));
1653 __MPFR_DECLSPEC void mpfr_sqrhigh_n _MPFR_PROTO ((mp_ptr, mp_srcptr, mp_size_t));
1654
1655 __MPFR_DECLSPEC int mpfr_round_p _MPFR_PROTO ((mp_limb_t *, mp_size_t,
1656 mp_exp_t, mp_prec_t));
1657
1658 __MPFR_DECLSPEC void mpfr_dump_mant _MPFR_PROTO ((const mp_limb_t *,
1659 mp_prec_t, mp_prec_t,
1660 mp_prec_t));
1661
1662 __MPFR_DECLSPEC int mpfr_round_near_x _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
1663 mpfr_uexp_t, int,
1664 mp_rnd_t));
1665 __MPFR_DECLSPEC void mpfr_abort_prec_max _MPFR_PROTO ((void))
1666 MPFR_NORETURN_ATTR;
1667
1668 #if defined (__cplusplus)
1669 }
1670 #endif
1671
1672 #endif
DragonFly BSD