mpfr-impl.h

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