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 /*
  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