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