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