contrib/mpfr/mpfr-impl.h

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