Home | History | Annotate | Line # | Download | only in src
      1 /* Utilities for MPFR developers, not exported.
      2 
      3 Copyright 1999-2023 Free Software Foundation, Inc.
      4 Contributed by the AriC and Caramba 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 https://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 /* Include config.h before using ANY configure macros if needed. */
     27 #ifdef HAVE_CONFIG_H
     28 # include "config.h"
     29 #endif
     30 
     31 
     32 /******************************************************
     33  *****************  Standard headers  *****************
     34  ******************************************************/
     35 
     36 /* Let's include some standard headers unconditionally as they are
     37    already needed by several source files or when some options are
     38    enabled/disabled, and it is easy to forget them (some configure
     39    options may hide the error).
     40    Note: If some source file must not have such a header included
     41    (which is very unlikely and probably means something broken in
     42    this source file), we should do that with some macro (that would
     43    also force to disable incompatible features). */
     44 
     45 #if defined (__cplusplus)
     46 # include <cstdio>
     47 # include <cstring>
     48 #else
     49 # include <stdio.h>
     50 # include <string.h>
     51 #endif
     52 
     53 /* Since <stdio.h> (<cstdio> for C++) is unconditionally included... */
     54 #define MPFR_USE_FILE
     55 
     56 #include <stdlib.h>
     57 #include <limits.h>
     58 #include <float.h>  /* for FLT_RADIX, etc., tested below */
     59 
     60 
     61 /******************************************************
     62  *****************  Include files  ********************
     63  ******************************************************/
     64 
     65 /* The macros defined in mpfr-cvers.h do not depend on anything,
     66    so that it is better to include this header file early: then
     67    it can be used by any other header. */
     68 #include "mpfr-cvers.h"
     69 
     70 #if defined(_MPFR_EXP_FORMAT) && _MPFR_EXP_FORMAT == 4
     71 /* mpfr_exp_t will be defined as intmax_t */
     72 # undef MPFR_NEED_INTMAX_H
     73 # define MPFR_NEED_INTMAX_H 1
     74 #endif
     75 
     76 #ifdef MPFR_NEED_INTMAX_H
     77 # include "mpfr-intmax.h"
     78 #endif
     79 
     80 /* Check if we are inside a build of MPFR or inside the test suite.
     81    This is needed in mpfr.h to export or import the functions.
     82    It matters only for Windows DLL */
     83 #ifndef __MPFR_TEST_H__
     84 # define __MPFR_WITHIN_MPFR 1
     85 #endif
     86 
     87 /* For the definition of MPFR_THREAD_ATTR. GCC/ICC detection macros are
     88    no longer used, as they sometimes gave incorrect information about
     89    the support of thread-local variables. A configure check is now done. */
     90 #if defined(MPFR_WANT_SHARED_CACHE)
     91 # define MPFR_NEED_THREAD_LOCK 1
     92 #endif
     93 #include "mpfr-thread.h"
     94 
     95 #ifndef MPFR_USE_MINI_GMP
     96 #include "gmp.h"
     97 #else
     98 #include "mini-gmp.h"
     99 #endif
    100 
    101 /* With the current code, MPFR_LONG_WITHIN_LIMB must be defined if an
    102    unsigned long fits in a limb. Since one cannot rely on the configure
    103    tests entirely (in particular when GMP is involved) and some platforms
    104    may not use configure, handle this definition here. A limb (mp_limb_t)
    105    is normally defined as an unsigned long, but this may not be the case
    106    with mini-gmp (and we can't rely on __GMP_SHORT_LIMB for this reason).
    107    So, concerning mp_limb_t, we can only test GMP_NUMB_BITS.
    108    Chosen heuristic: define MPFR_LONG_WITHIN_LIMB only when
    109      * mp_limb_t and unsigned long have both 32 bits exactly, or
    110      * mp_limb_t has at least 64 bits.
    111    Since we require that mp_limb_t have a size that is a power of 2, we
    112    can currently be wrong only if mini-gmp is used and unsigned long has
    113    more than 64 bits, which is unlikely to occur. */
    114 #if GMP_NUMB_BITS >= 64 || (GMP_NUMB_BITS == 32 && ULONG_MAX == 0xffffffff)
    115 # undef MPFR_LONG_WITHIN_LIMB
    116 # define MPFR_LONG_WITHIN_LIMB 1
    117 #endif
    118 
    119 #ifdef MPFR_HAVE_GMP_IMPL /* Build with gmp internals */
    120 
    121 # ifdef MPFR_USE_MINI_GMP
    122 #  error "MPFR_HAVE_GMP_IMPL and MPFR_USE_MINI_GMP must not be both defined"
    123 # endif
    124 # include "gmp-impl.h"
    125 # ifdef MPFR_NEED_LONGLONG_H
    126 #  include "longlong.h"
    127 # endif
    128 # include "mpfr.h"
    129 # include "mpfr-gmp.h"
    130 
    131 #else /* Build without gmp internals */
    132 
    133 /* if using mini-gmp, include missing definitions in mini-gmp */
    134 # ifdef MPFR_USE_MINI_GMP
    135 #  include "mpfr-mini-gmp.h"
    136 # endif
    137 # include "mpfr.h"
    138 # include "mpfr-gmp.h"
    139 # ifdef MPFR_LONG_WITHIN_LIMB /* longlong.h is not valid otherwise */
    140 #  ifdef MPFR_NEED_LONGLONG_H
    141 #   define LONGLONG_STANDALONE
    142 #   include "mpfr-longlong.h"
    143 #  endif
    144 # endif
    145 
    146 #endif
    147 
    148 #undef MPFR_NEED_LONGLONG_H
    149 
    150 
    151 /******************************************************
    152  *************  Attribute definitions  ****************
    153  ******************************************************/
    154 
    155 #if defined(MPFR_HAVE_NORETURN)
    156 /* _Noreturn is specified by ISO C11 (Section 6.7.4);
    157    in GCC, it is supported as of version 4.7. */
    158 # define MPFR_NORETURN _Noreturn
    159 #elif !defined(noreturn)
    160 /* A noreturn macro could be defined if <stdnoreturn.h> has been included,
    161    in which case it would make sense to #define MPFR_NORETURN noreturn.
    162    But this is unlikely, as MPFR_HAVE_NORETURN should have been defined
    163    in such a case. So, in doubt, let us avoid any code that would use a
    164    noreturn macro, since it could be invalid. */
    165 # if __MPFR_GNUC(3,0) || __MPFR_ICC(8,1,0)
    166 #  define MPFR_NORETURN __attribute__ ((noreturn))
    167 # elif defined(_MSC_VER) && defined(_WIN32) && (_MSC_VER >= 1200)
    168 #  define MPFR_NORETURN __declspec (noreturn)
    169 # endif
    170 #endif
    171 #ifndef MPFR_NORETURN
    172 # define MPFR_NORETURN
    173 #endif
    174 
    175 #if __MPFR_GNUC(3,0) || __MPFR_ICC(8,1,0)
    176 # define MPFR_CONST_FUNCTION_ATTR   __attribute__ ((const))
    177 #else
    178 # define MPFR_CONST_FUNCTION_ATTR
    179 #endif
    180 
    181 #if __MPFR_GNUC(3,0) || __MPFR_ICC(8,1,0)
    182 # define MPFR_PURE_FUNCTION_ATTR    __attribute__ ((pure))
    183 #else
    184 # define MPFR_PURE_FUNCTION_ATTR
    185 #endif
    186 
    187 /* The hot attribute on a function is used to inform the compiler
    188    that the function is a hot spot of the compiled program. */
    189 #if __MPFR_GNUC(4,3)
    190 # define MPFR_HOT_FUNCTION_ATTR     __attribute__ ((hot))
    191 #else
    192 # define MPFR_HOT_FUNCTION_ATTR
    193 #endif
    194 
    195 /* The cold attribute on functions is used to inform the compiler
    196    that the function is unlikely to be executed. */
    197 #if __MPFR_GNUC(4,3)
    198 # define MPFR_COLD_FUNCTION_ATTR    __attribute__ ((cold))
    199 #else
    200 # define MPFR_COLD_FUNCTION_ATTR
    201 #endif
    202 
    203 /* Add MPFR_MAYBE_UNUSED after a variable declaration to avoid compiler
    204    warnings if it is not used.
    205    TODO: To be replaced by the future maybe_unused attribute (C2x) once
    206    supported. */
    207 #if __MPFR_GNUC(3,4)
    208 #define MPFR_MAYBE_UNUSED __attribute__ ((unused))
    209 #else
    210 #define MPFR_MAYBE_UNUSED
    211 #endif
    212 
    213 /* This MPFR_FALLTHROUGH macro allows one to make fallthrough in switch case
    214    explicit. Use this macro at the end of a switch case if it falls through,
    215    in order to avoid a -Wimplicit-fallthrough warning. */
    216 #if __MPFR_GNUC(7,0)
    217 #define MPFR_FALLTHROUGH __attribute__ ((fallthrough))
    218 #else
    219 #define MPFR_FALLTHROUGH
    220 #endif
    221 
    222 
    223 
    224 /******************************************************
    225  ***  Global internal variables and related macros  ***
    226  ******************************************************/
    227 
    228 #if defined (__cplusplus)
    229 extern "C" {
    230 #endif
    231 
    232 #if defined(MPFR_WANT_SHARED_CACHE)
    233 # define MPFR_CACHE_ATTR
    234 #else
    235 # define MPFR_CACHE_ATTR MPFR_THREAD_ATTR
    236 #endif
    237 
    238 /* Note: The following structure and types depend on the MPFR build options
    239    (including compiler options), due to the various locking methods affecting
    240    MPFR_DEFERRED_INIT_SLAVE_DECL and MPFR_LOCK_DECL. But since this is only
    241    internal, that's OK. */
    242 struct __gmpfr_cache_s {
    243   mpfr_t x;
    244   int inexact;
    245   int (*func)(mpfr_ptr, mpfr_rnd_t);
    246   MPFR_DEFERRED_INIT_SLAVE_DECL()
    247   MPFR_LOCK_DECL(lock)
    248 };
    249 typedef struct __gmpfr_cache_s mpfr_cache_t[1];
    250 typedef struct __gmpfr_cache_s *mpfr_cache_ptr;
    251 
    252 #if __GMP_LIBGMP_DLL
    253 # define MPFR_WIN_THREAD_SAFE_DLL 1
    254 #endif
    255 
    256 #if defined(__MPFR_WITHIN_MPFR) || !defined(MPFR_WIN_THREAD_SAFE_DLL)
    257 extern MPFR_THREAD_ATTR mpfr_flags_t __gmpfr_flags;
    258 extern MPFR_THREAD_ATTR mpfr_exp_t   __gmpfr_emin;
    259 extern MPFR_THREAD_ATTR mpfr_exp_t   __gmpfr_emax;
    260 extern MPFR_THREAD_ATTR mpfr_prec_t  __gmpfr_default_fp_bit_precision;
    261 extern MPFR_THREAD_ATTR mpfr_rnd_t   __gmpfr_default_rounding_mode;
    262 extern MPFR_CACHE_ATTR  mpfr_cache_t __gmpfr_cache_const_euler;
    263 extern MPFR_CACHE_ATTR  mpfr_cache_t __gmpfr_cache_const_catalan;
    264 # ifndef MPFR_USE_LOGGING
    265 extern MPFR_CACHE_ATTR  mpfr_cache_t __gmpfr_cache_const_pi;
    266 extern MPFR_CACHE_ATTR  mpfr_cache_t __gmpfr_cache_const_log2;
    267 # else
    268 /* Two constants are used by the logging functions (via mpfr_fprintf,
    269    then mpfr_log, for the base conversion): pi and log(2). Since the
    270    mpfr_cache function isn't re-entrant when working on the same cache,
    271    we need to define two caches for each constant. */
    272 extern MPFR_CACHE_ATTR  mpfr_cache_t   __gmpfr_normal_pi;
    273 extern MPFR_CACHE_ATTR  mpfr_cache_t   __gmpfr_normal_log2;
    274 extern MPFR_CACHE_ATTR  mpfr_cache_t   __gmpfr_logging_pi;
    275 extern MPFR_CACHE_ATTR  mpfr_cache_t   __gmpfr_logging_log2;
    276 extern MPFR_CACHE_ATTR  mpfr_cache_ptr __gmpfr_cache_const_pi;
    277 extern MPFR_CACHE_ATTR  mpfr_cache_ptr __gmpfr_cache_const_log2;
    278 # endif
    279 #endif
    280 
    281 #ifdef MPFR_WIN_THREAD_SAFE_DLL
    282 # define MPFR_MAKE_VARFCT(T,N) T * N ## _f (void) { return &N; }
    283 __MPFR_DECLSPEC mpfr_flags_t * __gmpfr_flags_f (void);
    284 __MPFR_DECLSPEC mpfr_exp_t *   __gmpfr_emin_f (void);
    285 __MPFR_DECLSPEC mpfr_exp_t *   __gmpfr_emax_f (void);
    286 __MPFR_DECLSPEC mpfr_prec_t *  __gmpfr_default_fp_bit_precision_f (void);
    287 __MPFR_DECLSPEC mpfr_rnd_t *   __gmpfr_default_rounding_mode_f (void);
    288 __MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_euler_f (void);
    289 __MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_catalan_f (void);
    290 # ifndef MPFR_USE_LOGGING
    291 __MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_pi_f (void);
    292 __MPFR_DECLSPEC mpfr_cache_t * __gmpfr_cache_const_log2_f (void);
    293 # else
    294 __MPFR_DECLSPEC mpfr_cache_t *   __gmpfr_normal_pi_f (void);
    295 __MPFR_DECLSPEC mpfr_cache_t *   __gmpfr_normal_log2_f (void);
    296 __MPFR_DECLSPEC mpfr_cache_t *   __gmpfr_logging_pi_f (void);
    297 __MPFR_DECLSPEC mpfr_cache_t *   __gmpfr_logging_log2_f (void);
    298 __MPFR_DECLSPEC mpfr_cache_ptr * __gmpfr_cache_const_pi_f (void);
    299 __MPFR_DECLSPEC mpfr_cache_ptr * __gmpfr_cache_const_log2_f (void);
    300 # endif
    301 # ifndef __MPFR_WITHIN_MPFR
    302 #  define __gmpfr_flags                    (*__gmpfr_flags_f())
    303 #  define __gmpfr_emin                     (*__gmpfr_emin_f())
    304 #  define __gmpfr_emax                     (*__gmpfr_emax_f())
    305 #  define __gmpfr_default_fp_bit_precision (*__gmpfr_default_fp_bit_precision_f())
    306 #  define __gmpfr_default_rounding_mode    (*__gmpfr_default_rounding_mode_f())
    307 #  define __gmpfr_cache_const_euler        (*__gmpfr_cache_const_euler_f())
    308 #  define __gmpfr_cache_const_catalan      (*__gmpfr_cache_const_catalan_f())
    309 #  ifndef MPFR_USE_LOGGING
    310 #   define __gmpfr_cache_const_pi         (*__gmpfr_cache_const_pi_f())
    311 #   define __gmpfr_cache_const_log2       (*__gmpfr_cache_const_log2_f())
    312 #  else
    313 #   define __gmpfr_normal_pi              (*__gmpfr_normal_pi_f())
    314 #   define __gmpfr_logging_pi             (*__gmpfr_logging_pi_f())
    315 #   define __gmpfr_logging_log2           (*__gmpfr_logging_log2_f())
    316 #   define __gmpfr_cache_const_pi         (*__gmpfr_cache_const_pi_f())
    317 #   define __gmpfr_cache_const_log2       (*__gmpfr_cache_const_log2_f())
    318 #  endif
    319 # endif
    320 #else
    321 # define MPFR_MAKE_VARFCT(T,N)
    322 #endif
    323 
    324 # define MPFR_THREAD_VAR(T,N,V)    \
    325   MPFR_THREAD_ATTR T N = (V);      \
    326   MPFR_MAKE_VARFCT (T,N)
    327 
    328 #define BASE_MAX 62
    329 __MPFR_DECLSPEC extern const __mpfr_struct __gmpfr_l2b[BASE_MAX-1][2];
    330 
    331 /* Note: do not use the following values when they can be outside the
    332    current exponent range, e.g. when the exponent range has not been
    333    extended yet; under such a condition, they can be used only in
    334    mpfr_cmpabs. */
    335 __MPFR_DECLSPEC extern const mpfr_t __gmpfr_one;
    336 __MPFR_DECLSPEC extern const mpfr_t __gmpfr_two;
    337 __MPFR_DECLSPEC extern const mpfr_t __gmpfr_four;
    338 __MPFR_DECLSPEC extern const mpfr_t __gmpfr_mone;
    339 __MPFR_DECLSPEC extern const mpfr_t __gmpfr_const_log2_RNDD;
    340 __MPFR_DECLSPEC extern const mpfr_t __gmpfr_const_log2_RNDU;
    341 
    342 #if defined (__cplusplus)
    343  }
    344 #endif
    345 
    346 /* Replace some common functions for direct access to the global vars.
    347    The casts prevent these macros from being used as a lvalue (and this
    348    method makes sure that the expressions have the correct type). */
    349 
    350 #define mpfr_get_emin() ((mpfr_exp_t) __gmpfr_emin)
    351 #define mpfr_get_emax() ((mpfr_exp_t) __gmpfr_emax)
    352 #define mpfr_get_default_rounding_mode() \
    353   ((mpfr_rnd_t) __gmpfr_default_rounding_mode)
    354 #define mpfr_get_default_prec() \
    355   ((mpfr_prec_t) __gmpfr_default_fp_bit_precision)
    356 
    357 /* Flags related macros. */
    358 /* Note: Function-like macros that modify __gmpfr_flags are not defined
    359    because of the risk to break the sequence point rules if two such
    360    macros are used in the same expression (without a sequence point
    361    between). For instance, mpfr_sgn currently uses mpfr_set_erangeflag,
    362    which mustn't be implemented as a macro for this reason. */
    363 
    364 #define mpfr_flags_test(mask) \
    365   (__gmpfr_flags & (mpfr_flags_t) (mask))
    366 
    367 #if MPFR_FLAGS_ALL <= INT_MAX
    368 #define mpfr_underflow_p() \
    369   ((int) (__gmpfr_flags & MPFR_FLAGS_UNDERFLOW))
    370 #define mpfr_overflow_p() \
    371   ((int) (__gmpfr_flags & MPFR_FLAGS_OVERFLOW))
    372 #define mpfr_nanflag_p() \
    373   ((int) (__gmpfr_flags & MPFR_FLAGS_NAN))
    374 #define mpfr_inexflag_p() \
    375   ((int) (__gmpfr_flags & MPFR_FLAGS_INEXACT))
    376 #define mpfr_erangeflag_p() \
    377   ((int) (__gmpfr_flags & MPFR_FLAGS_ERANGE))
    378 #define mpfr_divby0_p() \
    379   ((int) (__gmpfr_flags & MPFR_FLAGS_DIVBY0))
    380 #endif
    381 
    382 /* Use a do-while statement for the following macros in order to prevent
    383    one from using them in an expression, as the sequence point rules could
    384    be broken if __gmpfr_flags is assigned twice in the same expression
    385    (via macro expansions). For instance, the mpfr_sgn macro currently uses
    386    mpfr_set_erangeflag, which mustn't be implemented as a function-like
    387    macro for this reason. It is not clear whether an expression with
    388    sequence points, like (void) (0, __gmpfr_flags = 0), would avoid UB. */
    389 #define MPFR_CLEAR_FLAGS() \
    390   do __gmpfr_flags = 0; while (0)
    391 #define MPFR_CLEAR_UNDERFLOW() \
    392   do __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_UNDERFLOW; while (0)
    393 #define MPFR_CLEAR_OVERFLOW() \
    394   do __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_OVERFLOW; while (0)
    395 #define MPFR_CLEAR_DIVBY0() \
    396   do __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_DIVBY0; while (0)
    397 #define MPFR_CLEAR_NANFLAG() \
    398   do __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_NAN; while (0)
    399 #define MPFR_CLEAR_INEXFLAG() \
    400   do __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_INEXACT; while (0)
    401 #define MPFR_CLEAR_ERANGEFLAG() \
    402   do __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE; while (0)
    403 #define MPFR_SET_UNDERFLOW() \
    404   do __gmpfr_flags |= MPFR_FLAGS_UNDERFLOW; while (0)
    405 #define MPFR_SET_OVERFLOW() \
    406   do __gmpfr_flags |= MPFR_FLAGS_OVERFLOW; while (0)
    407 #define MPFR_SET_DIVBY0() \
    408   do __gmpfr_flags |= MPFR_FLAGS_DIVBY0; while (0)
    409 #define MPFR_SET_NANFLAG() \
    410   do __gmpfr_flags |= MPFR_FLAGS_NAN; while (0)
    411 #define MPFR_SET_INEXFLAG() \
    412   do __gmpfr_flags |= MPFR_FLAGS_INEXACT; while (0)
    413 #define MPFR_SET_ERANGEFLAG() \
    414   do __gmpfr_flags |= MPFR_FLAGS_ERANGE; while (0)
    415 
    416 /* Testing an exception flag correctly is tricky. There are mainly two
    417    pitfalls: First, one needs to remember to clear the corresponding
    418    flag, in case it was set before the function call or during some
    419    intermediate computations (in practice, one can clear all the flags).
    420    Secondly, one needs to test the flag early enough, i.e. before it
    421    can be modified by another function. Moreover, it is quite difficult
    422    (if not impossible) to reliably check problems with "make check". To
    423    avoid these pitfalls, it is recommended to use the following macros.
    424    Other use of the exception-flag predicate functions/macros will be
    425    detected by mpfrlint.
    426    Note: _op can be either a statement or an expression.
    427    MPFR_BLOCK_EXCEP should be used only inside a block; it is useful to
    428    detect some exception in order to exit the block as soon as possible. */
    429 #define MPFR_BLOCK_DECL(_flags) mpfr_flags_t _flags
    430 /* The (void) (_flags) makes sure that _flags is read at least once (it
    431    makes sense to use MPFR_BLOCK while _flags will never be read in the
    432    source, so that we wish to avoid the corresponding warning). */
    433 #define MPFR_BLOCK(_flags,_op)          \
    434   do                                    \
    435     {                                   \
    436       MPFR_CLEAR_FLAGS ();              \
    437       _op;                              \
    438       (_flags) = __gmpfr_flags;         \
    439       (void) (_flags);                  \
    440     }                                   \
    441   while (0)
    442 #define MPFR_BLOCK_TEST(_flags,_f) MPFR_UNLIKELY ((_flags) & (_f))
    443 #define MPFR_BLOCK_EXCEP (__gmpfr_flags & (MPFR_FLAGS_UNDERFLOW | \
    444                                            MPFR_FLAGS_OVERFLOW | \
    445                                            MPFR_FLAGS_DIVBY0 | \
    446                                            MPFR_FLAGS_NAN))
    447 /* Let's use a MPFR_ prefix, because e.g. OVERFLOW is defined by glibc's
    448    math.h, though this is not a reserved identifier! */
    449 #define MPFR_UNDERFLOW(_flags)  MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_UNDERFLOW)
    450 #define MPFR_OVERFLOW(_flags)   MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_OVERFLOW)
    451 #define MPFR_NANFLAG(_flags)    MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_NAN)
    452 #define MPFR_INEXFLAG(_flags)   MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_INEXACT)
    453 #define MPFR_ERANGEFLAG(_flags) MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_ERANGE)
    454 #define MPFR_DIVBY0(_flags)     MPFR_BLOCK_TEST (_flags, MPFR_FLAGS_DIVBY0)
    455 
    456 
    457 /******************************************************
    458  *******************  Assertions  *********************
    459  ******************************************************/
    460 
    461 /* MPFR_WANT_ASSERT can take 4 values (the default value is 0):
    462    -1 (or below): Do not check any assertion. Discouraged, in particular
    463      for a shared library (for time-critical applications, LTO with a
    464      static library should also be used anyway).
    465    0: Check normal assertions.
    466    1: Check debugging assertions too.
    467    2 (or above): Additional checks that may take time. For instance,
    468      some functions may be tested by using two different implementations
    469      and comparing the results.
    470 */
    471 
    472 /* Note: do not use GMP macros ASSERT_ALWAYS and ASSERT as they are not
    473    expressions, and as a consequence, they cannot be used in a for(),
    474    with a comma operator and so on. */
    475 
    476 /* MPFR_ASSERTN(expr): assertions that should normally be checked,
    477      otherwise give a hint to the compiler.
    478    MPFR_ASSERTD(expr): assertions that should be checked when testing,
    479      otherwise give a hint to the compiler.
    480    MPFR_DBGRES(assignment): to be used when the result is tested only
    481      in an MPFR_ASSERTD expression (in order to avoid a warning, e.g.
    482      with GCC's -Wunused-but-set-variable, in non-debug mode).
    483      Note: WG14/N2270 proposed a maybe_unused attribute, which could
    484      be useful to avoid MPFR_DBGRES. See:
    485        https://www.open-std.org/jtc1/sc22/wg14/www/docs/n2270.pdf
    486    Note: Evaluating expr might yield side effects, but such side effects
    487    must not change the results (except by yielding an assertion failure).
    488 */
    489 #ifndef MPFR_WANT_ASSERT
    490 # define MPFR_WANT_ASSERT 0
    491 #endif
    492 
    493 #if MPFR_WANT_ASSERT < 0
    494 # undef MPFR_EXP_CHECK
    495 # define MPFR_ASSERTN(expr)  MPFR_ASSUME (expr)
    496 #else
    497 # define MPFR_ASSERTN(expr)  \
    498   ((void) ((MPFR_LIKELY(expr)) || (ASSERT_FAIL(expr),MPFR_ASSUME(expr),0)))
    499 /* Some explanations: mpfr_assert_fail is not marked as "no return",
    500    so that ((void) ((MPFR_LIKELY(expr)) || (ASSERT_FAIL(expr),0)))
    501    cannot be a way to tell the compiler that after this code, expr is
    502    necessarily true. The MPFR_ASSUME(expr) is a way to tell the compiler
    503    that if expr is false, then ASSERT_FAIL(expr) does not return
    504    (otherwise they would be a contradiction / UB when MPFR_ASSUME(expr)
    505    is reached). Such information is useful to avoid warnings like those
    506    from -Wmaybe-uninitialized, e.g. in tests/turandom.c r11663 on t[0]
    507    from "mpfr_equal_p (y, t[0])".
    508    TODO: Remove the MPFR_ASSUME(expr) once mpfr_assert_fail is marked as
    509    "no return".
    510  */
    511 #endif
    512 
    513 #if MPFR_WANT_ASSERT > 0
    514 # define MPFR_EXP_CHECK 1
    515 # define MPFR_ASSERTD(expr)  MPFR_ASSERTN (expr)
    516 # define MPFR_DBGRES(A)      (A)
    517 #else
    518 # define MPFR_ASSERTD(expr)  MPFR_ASSUME (expr)
    519 # define MPFR_DBGRES(A)      ((void) (A))
    520 #endif
    521 
    522 /* MPFR_ASSUME is like assert(), but it is a hint to a compiler about a
    523    statement of fact in a function call free expression, which allows
    524    the compiler to generate better machine code.
    525    __builtin_unreachable has been introduced in GCC 4.5 but it works
    526    fine since 4.8 only (before it may generate unoptimized code if there
    527    are more than one decision).
    528    Note:
    529      The goal of MPFR_ASSUME() is to allow the compiler to optimize even
    530      more. Thus we need to make sure that its use in MPFR will never yield
    531      code generation. Since MPFR_ASSUME() may be used by MPFR_ASSERTN()
    532      and MPFR_ASSERTD(), whose expression might have side effects, we need
    533      to make sure that the expression x is not evaluated in such a case.
    534      This is done with __builtin_constant_p (!!(x) || !(x)), whose value
    535      is 0 if x has side effects, and normally 1 if the compiler knows that
    536      x has no side effects (since here, it can deduce that !!(x) || !(x)
    537      is equivalent to the constant 1). In the former case, MPFR_ASSUME(x)
    538      will give (void) 0, and in the latter case, it will give:
    539        (x) ? (void) 0 : __builtin_unreachable()
    540    In the development code, it is better to use MPFR_ASSERTD than
    541    MPFR_ASSUME, since it'll check if the condition is true in debug
    542    build.
    543 */
    544 #if defined(MPFR_HAVE_BUILTIN_UNREACHABLE) && __MPFR_GNUC(4, 8)
    545 # define MPFR_ASSUME(x)                                 \
    546     (! __builtin_constant_p (!!(x) || !(x)) || (x) ?    \
    547      (void) 0 : __builtin_unreachable())
    548 #elif defined(_MSC_VER)
    549 # define MPFR_ASSUME(x) __assume(x)
    550 #else
    551 # define MPFR_ASSUME(x) ((void) 0)
    552 #endif
    553 
    554 #include "mpfr-sassert.h"
    555 
    556 /* Code to deal with impossible, for functions returning an int.
    557    The "return 0;" avoids an error with current GCC versions and
    558    "-Werror=return-type".
    559    WARNING: It doesn't use do { } while (0) for Insure++ */
    560 #if defined(HAVE_BUILTIN_UNREACHABLE)
    561 # define MPFR_RET_NEVER_GO_HERE() do { __builtin_unreachable(); } while (0)
    562 #else
    563 # define MPFR_RET_NEVER_GO_HERE() do { MPFR_ASSERTN(0); return 0; } while (0)
    564 #endif
    565 
    566 
    567 /******************************************************
    568  *******************  Warnings  ***********************
    569  ******************************************************/
    570 
    571 /* MPFR_WARNING is no longer useful, but let's keep the macro in case
    572    it needs to be used again in the future. */
    573 
    574 #ifdef MPFR_USE_WARNINGS
    575 # define MPFR_WARNING(W)                    \
    576   do                                        \
    577     {                                       \
    578       char *q = getenv ("MPFR_QUIET");      \
    579       if (q == NULL || *q == 0)             \
    580         fprintf (stderr, "MPFR: %s\n", W);  \
    581     }                                       \
    582   while (0)
    583 #else
    584 # define MPFR_WARNING(W)  ((void) 0)
    585 #endif
    586 
    587 
    588 /******************************************************
    589  *****************  double macros  ********************
    590  ******************************************************/
    591 
    592 /* Precision used for lower precision computations */
    593 #define MPFR_SMALL_PRECISION 32
    594 
    595 /* Definition of constants */
    596 #define LOG2 0.69314718055994528622 /* log(2) rounded to zero on 53 bits */
    597 
    598 /* MPFR_DOUBLE_SPEC = 1 if the C type 'double' corresponds to IEEE-754
    599    double precision, 0 if it doesn't, and undefined if one doesn't know.
    600    On all the tested machines, MPFR_DOUBLE_SPEC = 1. To have this macro
    601    defined here, #include <float.h> is needed. If need be, other values
    602    could be defined for other specs (once they are known). */
    603 #if !defined(MPFR_DOUBLE_SPEC) && defined(FLT_RADIX) && \
    604     defined(DBL_MANT_DIG) && defined(DBL_MIN_EXP) && defined(DBL_MAX_EXP)
    605 # if FLT_RADIX == 2 && DBL_MANT_DIG == 53 && \
    606      DBL_MIN_EXP == -1021 && DBL_MAX_EXP == 1024
    607 #  define MPFR_DOUBLE_SPEC 1
    608 # else
    609 #  define MPFR_DOUBLE_SPEC 0
    610 # endif
    611 #endif
    612 
    613 /* With -DMPFR_DISABLE_IEEE_FLOATS, exercise non IEEE floats */
    614 #ifdef MPFR_DISABLE_IEEE_FLOATS
    615 # ifdef _MPFR_IEEE_FLOATS
    616 #  undef _MPFR_IEEE_FLOATS
    617 # endif
    618 # define _MPFR_IEEE_FLOATS 0
    619 # undef HAVE_LDOUBLE_IS_DOUBLE
    620 # undef HAVE_LDOUBLE_IEEE_EXT_LITTLE
    621 # undef HAVE_LDOUBLE_IEEE_EXT_BIG
    622 # undef HAVE_LDOUBLE_IEEE_QUAD_BIG
    623 # undef HAVE_LDOUBLE_IEEE_QUAD_LITTLE
    624 #endif
    625 
    626 #ifndef IEEE_DBL_MANT_DIG
    627 #define IEEE_DBL_MANT_DIG 53
    628 #endif
    629 #define MPFR_LIMBS_PER_DOUBLE ((IEEE_DBL_MANT_DIG-1)/GMP_NUMB_BITS+1)
    630 
    631 #ifndef IEEE_FLT_MANT_DIG
    632 #define IEEE_FLT_MANT_DIG 24
    633 #endif
    634 #define MPFR_LIMBS_PER_FLT ((IEEE_FLT_MANT_DIG-1)/GMP_NUMB_BITS+1)
    635 
    636 /* Visual C++ doesn't support +1.0/0.0, -1.0/0.0 and 0.0/0.0
    637    at compile time.
    638    Clang with -fsanitize=undefined is a bit similar due to a bug:
    639      https://llvm.org/bugs/show_bug.cgi?id=17381 (fixed on 2015-12-03)
    640    but even without its sanitizer, it may be better to use the
    641    double_zero version until IEEE 754 division by zero is properly
    642    supported:
    643      https://llvm.org/bugs/show_bug.cgi?id=17005
    644    Note: DBL_NAN is 0/0, thus its value is a quiet NaN (qNAN).
    645 */
    646 #if (defined(_MSC_VER) && defined(_WIN32) && (_MSC_VER >= 1200)) || \
    647     defined(__clang__)
    648 static double double_zero = 0.0;
    649 # define DBL_NAN (double_zero/double_zero)
    650 # define DBL_POS_INF ((double) 1.0/double_zero)
    651 # define DBL_NEG_INF ((double)-1.0/double_zero)
    652 # define DBL_NEG_ZERO (-double_zero)
    653 #else
    654 # define DBL_POS_INF ((double) 1.0/0.0)
    655 # define DBL_NEG_INF ((double)-1.0/0.0)
    656 # define DBL_NAN     ((double) 0.0/0.0)
    657 # define DBL_NEG_ZERO (-0.0)
    658 #endif
    659 
    660 /* Note: In the past, there was specific code for _MPFR_IEEE_FLOATS, which
    661    was based on NaN and Inf memory representations. This code was breaking
    662    the aliasing rules (see ISO C99, 6.5#6 and 6.5#7 on the effective type)
    663    and for this reason it did not behave correctly with GCC 4.5.0 20091119.
    664    The code needed a memory transfer and was probably not better than the
    665    macros below with a good compiler (a fix based on the NaN / Inf memory
    666    representation would be even worse due to C limitations), and this code
    667    could be selected only when MPFR was built with --with-gmp-build, thus
    668    introducing a difference (bad for maintaining/testing MPFR); therefore
    669    it has been removed. The old code required that the argument x be an
    670    lvalue of type double. We still require that, in case one would need
    671    to change the macros below, e.g. for some broken compiler. But the
    672    LVALUE(x) condition could be removed if really necessary. */
    673 /* Below, the &(x) == &(x) || &(x) != &(x) allows to make sure that x
    674    is a lvalue without (probably) any warning from the compiler.  The
    675    &(x) != &(x) is needed to avoid a failure under Mac OS X 10.4.11
    676    (with Xcode 2.4.1, i.e. the latest one). */
    677 #define LVALUE(x) (&(x) == &(x) || &(x) != &(x))
    678 #define DOUBLE_ISINF(x) (LVALUE(x) && ((x) > DBL_MAX || (x) < -DBL_MAX))
    679 /* The DOUBLE_ISNAN(x) macro must be valid with any real floating type,
    680    thus constants must be of integer type (e.g. 0). */
    681 #if defined(MPFR_NANISNAN) || (__MPFR_GNUC(1,0) && !defined(__STRICT_ANSI__))
    682 /* Avoid MIPSpro / IRIX64 / GCC (incorrect) optimizations.
    683    The + must not be replaced by a ||. With gcc -ffast-math, NaN is
    684    regarded as a positive number or something like that; the second
    685    test catches this case.
    686    [2016-03-01] Various tests now fail with gcc -ffast-math or just
    687    -ffinite-math-only; such options are not supported, but this makes
    688    difficult to test MPFR assuming x == x optimization to 1. Anyway
    689    support of functions/tests of using native FP and special values for
    690    non-IEEE-754 environment will always be on a case-by-case basis.
    691    [2018-06-02] Let's use this macro instead of the usual (x) != (x) test
    692    with all GCC versions except in ISO C mode[*], as due to
    693      https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
    694    there is no guarantee that (x) != (x) will be true only for NaN.
    695    Testing __STRICT_ANSI__ is suggested in:
    696      https://gcc.gnu.org/bugzilla/show_bug.cgi?id=85995
    697    but this is not safe if the user adds a -f option affecting conformance,
    698    in which case this would be a user error (however, note that the
    699    configure test associated with MPFR_NANISNAN will catch some issues).
    700 */
    701 # define DOUBLE_ISNAN(x) \
    702     (LVALUE(x) && !((((x) >= 0) + ((x) <= 0)) && -(x)*(x) <= 0))
    703 #else
    704 # define DOUBLE_ISNAN(x) (LVALUE(x) && (x) != (x))
    705 #endif
    706 
    707 
    708 /******************************************************
    709  **********  long double macros and typedef  **********
    710  ******************************************************/
    711 
    712 /* We try to get the exact value of the precision of long double
    713    (provided by the implementation) in order to provide correct
    714    rounding in this case (not guaranteed if the C implementation
    715    does not have an adequate long double arithmetic). Note that
    716    it may be lower than the precision of some numbers that can
    717    be represented in a long double; e.g. on FreeBSD/x86, it is
    718    53 because the processor is configured to round in double
    719    precision (even when using the long double type -- this is a
    720    limitation of the x87 arithmetic), and on Mac OS X, it is 106
    721    because the implementation is a double-double arithmetic.
    722    Otherwise (e.g. in base 10), we get an upper bound of the
    723    precision, and correct rounding isn't currently provided.
    724 */
    725 
    726 /* Definitions are enabled only if <float.h> is included. */
    727 #if defined (FLT_RADIX)
    728 
    729 #if defined(LDBL_MANT_DIG) && FLT_RADIX == 2
    730 # define MPFR_LDBL_MANT_DIG LDBL_MANT_DIG
    731 #else
    732 # define MPFR_LDBL_MANT_DIG \
    733   (sizeof(long double)*GMP_NUMB_BITS/sizeof(mp_limb_t))
    734 #endif
    735 
    736 /* LONGDOUBLE_NAN_ACTION executes the code "action" if x is a NaN. */
    737 
    738 /* On hppa2.0n-hp-hpux10 with the unbundled HP cc, the test x!=x on a NaN
    739    has been seen false, meaning NaNs are not detected.  This seemed to
    740    happen only after other comparisons, not sure what's really going on.  In
    741    any case we can pick apart the bytes to identify a NaN.  */
    742 #ifdef HAVE_LDOUBLE_IEEE_QUAD_BIG
    743 # define LONGDOUBLE_NAN_ACTION(x, action)                       \
    744   do {                                                          \
    745     union {                                                     \
    746       long double    ld;                                        \
    747       struct {                                                  \
    748         unsigned int sign : 1;                                  \
    749         unsigned int exp  : 15;                                 \
    750         unsigned int man3 : 16;                                 \
    751         unsigned int man2 : 32;                                 \
    752         unsigned int man1 : 32;                                 \
    753         unsigned int man0 : 32;                                 \
    754       } s;                                                      \
    755     } u;                                                        \
    756     u.ld = (x);                                                 \
    757     if (u.s.exp == 0x7FFFL                                      \
    758         && (u.s.man0 | u.s.man1 | u.s.man2 | u.s.man3) != 0)    \
    759       { action; }                                               \
    760   } while (0)
    761 #endif
    762 
    763 #ifdef HAVE_LDOUBLE_IEEE_QUAD_LITTLE
    764 # define LONGDOUBLE_NAN_ACTION(x, action)                       \
    765   do {                                                          \
    766     union {                                                     \
    767       long double    ld;                                        \
    768       struct {                                                  \
    769         unsigned int man0 : 32;                                 \
    770         unsigned int man1 : 32;                                 \
    771         unsigned int man2 : 32;                                 \
    772         unsigned int man3 : 16;                                 \
    773         unsigned int exp  : 15;                                 \
    774         unsigned int sign : 1;                                  \
    775       } s;                                                      \
    776     } u;                                                        \
    777     u.ld = (x);                                                 \
    778     if (u.s.exp == 0x7FFFL                                      \
    779         && (u.s.man0 | u.s.man1 | u.s.man2 | u.s.man3) != 0)    \
    780       { action; }                                               \
    781   } while (0)
    782 #endif
    783 
    784 /* Under IEEE rules, NaN is not equal to anything, including itself.
    785    "volatile" here stops "cc" on mips64-sgi-irix6.5 from optimizing away
    786    x!=x. */
    787 #ifndef LONGDOUBLE_NAN_ACTION
    788 # define LONGDOUBLE_NAN_ACTION(x, action)               \
    789   do {                                                  \
    790     volatile long double __x = LONGDOUBLE_VOLATILE (x); \
    791     if ((x) != __x)                                     \
    792       { action; }                                       \
    793   } while (0)
    794 
    795 /* Some compilers do not have a proper "volatile" and #define volatile
    796    to empty (to avoid a build failure with programs using "volatile"),
    797    i.e. "volatile" is just ignored and will not prevent optimizations
    798    that could potentially break the IEEE rules. In this case, call an
    799    external function, hoping that the compiler will not optimize. */
    800 # ifdef volatile
    801 __MPFR_DECLSPEC long double
    802   __gmpfr_longdouble_volatile (long double) MPFR_CONST_FUNCTION_ATTR;
    803 #  define LONGDOUBLE_VOLATILE(x)  (__gmpfr_longdouble_volatile (x))
    804 #  define WANT_GMPFR_LONGDOUBLE_VOLATILE 1
    805 # else
    806 #  define LONGDOUBLE_VOLATILE(x)  (x)
    807 # endif
    808 #endif
    809 
    810 /* Some special case for IEEE_EXT Little Endian */
    811 #if HAVE_LDOUBLE_IEEE_EXT_LITTLE
    812 
    813 typedef union {
    814   long double    ld;
    815   struct {
    816     unsigned int manl : 32;
    817     unsigned int manh : 32;
    818     unsigned int expl : 8 ;
    819     unsigned int exph : 7;
    820     unsigned int sign : 1;
    821   } s;
    822 } mpfr_long_double_t;
    823 
    824 #endif /* HAVE_LDOUBLE_IEEE_EXT_LITTLE */
    825 
    826 #endif  /* long double macros and typedef */
    827 
    828 
    829 /******************************************************
    830  *****************  _Float128 support  ****************
    831  ******************************************************/
    832 
    833 /* This is standardized by IEEE 754-2008. */
    834 #define IEEE_FLOAT128_MANT_DIG 113
    835 
    836 
    837 /******************************************************
    838  ******************  Decimal support  *****************
    839  ******************************************************/
    840 
    841 #ifdef MPFR_WANT_DECIMAL_FLOATS
    842 
    843 #if defined(__GNUC__) && \
    844     __FLT64_DECIMAL_DIG__ == 17 && \
    845     __FLT128_DECIMAL_DIG__ == 36
    846 
    847 /* GCC has standard _Decimal64 and _Decimal128 support.
    848    We may be able to detect the encoding here at compile time.
    849 
    850    Note: GCC may define __FLT64_DECIMAL_DIG__ and __FLT128_DECIMAL_DIG__
    851    even when it does not support _Decimal64 and _Decimal128, e.g. on
    852    aarch64 and sparc64. But since MPFR_WANT_DECIMAL_FLOATS has been
    853    defined, we already know that these types should be supported.
    854 
    855    Determining which encoding is used via macros is not documented, and
    856    there is the risk of being wrong. Currently __DECIMAL_BID_FORMAT__ is
    857    defined on x86, where the BID encoding is used. But on PowerPC, where
    858    the DPD encoding is used, nothing indicating the encoding is defined.
    859    A possible reason may be that the decimal support is provided by the
    860    hardware (in this case), so that GCC does not need to care about the
    861    encoding. Thus the absence of a __DECIMAL_BID_FORMAT__ macro would
    862    not necessarily imply DPD, as similarly in the future, GCC could
    863    support an ISA-level implementation using the BID encoding. */
    864 
    865 #ifdef __DECIMAL_BID_FORMAT__
    866 
    867 #if defined(DECIMAL_DPD_FORMAT)
    868 # error "Decimal encoding mismatch (DPD assumed, BID detected)"
    869 #elif !defined(DECIMAL_GENERIC_CODE)
    870 # define DECIMAL_BID_FORMAT 1
    871 #endif
    872 
    873 #endif  /* __DECIMAL_BID_FORMAT__ */
    874 
    875 #endif  /* __GNUC__ */
    876 
    877 #if defined(DECIMAL_GENERIC_CODE)
    878 # if defined(DECIMAL_BID_FORMAT)
    879 #  error "DECIMAL_BID_FORMAT and DECIMAL_GENERIC_CODE both defined"
    880 # endif
    881 # if defined(DECIMAL_DPD_FORMAT)
    882 #  error "DECIMAL_DPD_FORMAT and DECIMAL_GENERIC_CODE both defined"
    883 # endif
    884 #elif defined(DECIMAL_BID_FORMAT) || defined(DECIMAL_DPD_FORMAT)
    885 # if defined(DECIMAL_BID_FORMAT) && defined(DECIMAL_DPD_FORMAT)
    886 #  error "DECIMAL_BID_FORMAT and DECIMAL_DPD_FORMAT both defined"
    887 # endif
    888 #else
    889 # define DECIMAL_GENERIC_CODE 1
    890 #endif
    891 
    892 /* TODO: The following is ugly and possibly wrong on some platforms.
    893    Do something like union ieee_decimal128. */
    894 union ieee_double_decimal64 { double d; _Decimal64 d64; };
    895 
    896 /* FIXME: There's no reason to make the _Decimal128 code depend on
    897    whether _MPFR_IEEE_FLOATS is defined or not, as _MPFR_IEEE_FLOATS
    898    is about binary IEEE-754 floating point only. */
    899 #if _MPFR_IEEE_FLOATS
    900 /* TODO: It would be better to define a different structure for DPD,
    901    where the t* bit-fields correspond to the declets. And to avoid
    902    confusion and detect coding errors, these bit-fields should have
    903    different names for BID and DPD. */
    904 union ieee_decimal128
    905 {
    906   struct
    907     {
    908       /* Assume little-endian double implies little-endian for bit-field
    909          allocation (C99 says: "The order of allocation of bit-fields
    910          within a unit (high-order to low-order or low-order to high-order)
    911          is implementation-defined.") */
    912 #if defined(HAVE_DECIMAL128_IEEE_LITTLE_ENDIAN)
    913 #define HAVE_DECIMAL128_IEEE 1
    914       unsigned int t3:32;
    915       unsigned int t2:32;
    916       unsigned int t1:32;
    917       unsigned int t0:14;
    918       unsigned int comb:17;
    919       unsigned int sig:1;
    920 #elif defined(HAVE_DECIMAL128_IEEE_BIG_ENDIAN)
    921 #define HAVE_DECIMAL128_IEEE 1
    922       unsigned int sig:1;
    923       unsigned int comb:17;
    924       unsigned int t0:14;
    925       unsigned int t1:32;
    926       unsigned int t2:32;
    927       unsigned int t3:32;
    928 #else /* unknown bit-field ordering */
    929       /* This will not be used as HAVE_DECIMAL128_IEEE is not defined. */
    930       unsigned int dummy;
    931 #endif
    932     } s;
    933   _Decimal128 d128;
    934 };
    935 #endif /* _MPFR_IEEE_FLOATS */
    936 
    937 #endif /* MPFR_WANT_DECIMAL_FLOATS */
    938 
    939 
    940 /******************************************************
    941  ****************  mpfr_t properties  *****************
    942  ******************************************************/
    943 
    944 #define MPFR_PREC_COND(p) ((p) >= MPFR_PREC_MIN && (p) <= MPFR_PREC_MAX)
    945 #define MPFR_PREC_IN_RANGE(p) (MPFR_ASSERTD (MPFR_PREC_COND(p)), (p))
    946 
    947 /* In the following macro, p is usually a mpfr_prec_t, but this macro
    948    works with other integer types (without integer overflow). Checking
    949    that p >= 1 in debug mode is useful here because this macro can be
    950    used on a computed precision (in particular, this formula does not
    951    work for a degenerate case p = 0, and could give different results
    952    on different platforms). But let us not use an assertion checking
    953    in the MPFR_LAST_LIMB() and MPFR_LIMB_SIZE() macros below to avoid
    954    too much expansion for assertions (in practice, this should be a
    955    problem just when testing MPFR with the --enable-assert configure
    956    option and the -ansi -pedantic-errors gcc compiler flags). */
    957 #define MPFR_PREC2LIMBS(p) \
    958   (MPFR_ASSERTD ((p) >= 1), ((p) - 1) / GMP_NUMB_BITS + 1)
    959 
    960 #define MPFR_PREC(x)      ((x)->_mpfr_prec)
    961 #define MPFR_EXP(x)       ((x)->_mpfr_exp)
    962 #define MPFR_MANT(x)      ((x)->_mpfr_d)
    963 #define MPFR_GET_PREC(x)  MPFR_PREC_IN_RANGE (MPFR_PREC (x))
    964 #define MPFR_LAST_LIMB(x) ((MPFR_GET_PREC (x) - 1) / GMP_NUMB_BITS)
    965 #define MPFR_LIMB_SIZE(x) (MPFR_LAST_LIMB (x) + 1)
    966 
    967 
    968 /******************************************************
    969  ***************  Exponent properties  ****************
    970  ******************************************************/
    971 
    972 /* Limits of the mpfr_exp_t type (NOT those of valid exponent values).
    973    These macros can be used in preprocessor directives. */
    974 #if   _MPFR_EXP_FORMAT == 1
    975 # define MPFR_EXP_MAX (SHRT_MAX)
    976 # define MPFR_EXP_MIN (SHRT_MIN)
    977 #elif _MPFR_EXP_FORMAT == 2
    978 # define MPFR_EXP_MAX (INT_MAX)
    979 # define MPFR_EXP_MIN (INT_MIN)
    980 #elif _MPFR_EXP_FORMAT == 3
    981 # define MPFR_EXP_MAX (LONG_MAX)
    982 # define MPFR_EXP_MIN (LONG_MIN)
    983 #elif _MPFR_EXP_FORMAT == 4
    984 /* Note: MPFR_EXP_MAX and MPFR_EXP_MIN must not be used in #if directives
    985    if _MPFR_EXP_FORMAT == 4 and MPFR_HAVE_INTMAX_MAX is not defined. */
    986 # define MPFR_EXP_MAX (MPFR_INTMAX_MAX)
    987 # define MPFR_EXP_MIN (MPFR_INTMAX_MIN)
    988 #else
    989 # error "Invalid MPFR Exp format"
    990 #endif
    991 
    992 /* Before doing a cast to mpfr_uexp_t, make sure that the value is
    993    non-negative. */
    994 #define MPFR_UEXP(X) (MPFR_ASSERTD ((X) >= 0), (mpfr_uexp_t) (X))
    995 
    996 /* Define mpfr_eexp_t, mpfr_ueexp_t and MPFR_EXP_FSPEC.
    997    Warning! MPFR_EXP_FSPEC is the length modifier associated with
    998    these types mpfr_eexp_t / mpfr_ueexp_t, not with mpfr_exp_t.
    999    (Should we change that, or is this safer to detect bugs, e.g.
   1000    in the context of an expression with computations with long?)
   1001 */
   1002 #if _MPFR_EXP_FORMAT <= 3
   1003 typedef long mpfr_eexp_t;
   1004 typedef unsigned long mpfr_ueexp_t;
   1005 # define mpfr_get_exp_t(x,r) mpfr_get_si((x),(r))
   1006 # define mpfr_set_exp_t(x,e,r) mpfr_set_si((x),(e),(r))
   1007 # define MPFR_EXP_FSPEC "l"
   1008 #else
   1009 typedef intmax_t mpfr_eexp_t;
   1010 typedef uintmax_t mpfr_ueexp_t;
   1011 # define mpfr_get_exp_t(x,r) mpfr_get_sj((x),(r))
   1012 # define mpfr_set_exp_t(x,e,r) mpfr_set_sj((x),(e),(r))
   1013 # define MPFR_EXP_FSPEC "j"
   1014 #endif
   1015 
   1016 /* Size of mpfr_exp_t in limbs */
   1017 #define MPFR_EXP_LIMB_SIZE \
   1018   ((sizeof (mpfr_exp_t) - 1) / MPFR_BYTES_PER_MP_LIMB + 1)
   1019 
   1020 /* Invalid exponent value (to track bugs...) */
   1021 #define MPFR_EXP_INVALID \
   1022  ((mpfr_exp_t) 1 << (GMP_NUMB_BITS*sizeof(mpfr_exp_t)/sizeof(mp_limb_t)-2))
   1023 
   1024 /* Definition of the exponent limits for MPFR numbers.
   1025  * These limits are chosen so that if e is such an exponent, then 2e-1 and
   1026  * 2e+1 are representable. This is useful for intermediate computations,
   1027  * in particular the multiplication.
   1028  */
   1029 #undef MPFR_EMIN_MIN
   1030 #undef MPFR_EMIN_MAX
   1031 #undef MPFR_EMAX_MIN
   1032 #undef MPFR_EMAX_MAX
   1033 #define MPFR_EMIN_MIN (1-MPFR_EXP_INVALID)
   1034 #define MPFR_EMIN_MAX (MPFR_EXP_INVALID-1)
   1035 #define MPFR_EMAX_MIN (1-MPFR_EXP_INVALID)
   1036 #define MPFR_EMAX_MAX (MPFR_EXP_INVALID-1)
   1037 
   1038 /* Use MPFR_GET_EXP and MPFR_SET_EXP instead of MPFR_EXP directly,
   1039    unless when the exponent may be out-of-range, for instance when
   1040    setting the exponent before calling mpfr_check_range.
   1041    MPFR_EXP_CHECK is defined when MPFR_WANT_ASSERT is defined, but if you
   1042    don't use MPFR_WANT_ASSERT (for speed reasons), you can still define
   1043    MPFR_EXP_CHECK by setting -DMPFR_EXP_CHECK in $CFLAGS.
   1044    Note about MPFR_EXP_IN_RANGE and MPFR_SET_EXP:
   1045      The exp expression is required to have a signed type. To avoid spurious
   1046      failures, we could cast (exp) to mpfr_exp_t, but this wouldn't allow us
   1047      to detect some bugs that can occur on particular platforms. Anyway, an
   1048      unsigned type for exp is suspicious and should be regarded as a bug.
   1049 */
   1050 
   1051 #define MPFR_EXP_IN_RANGE(e)                                          \
   1052   (MPFR_ASSERTD (IS_SIGNED (e)), (e) >= __gmpfr_emin && (e) <= __gmpfr_emax)
   1053 
   1054 #ifdef MPFR_EXP_CHECK
   1055 # define MPFR_GET_EXP(x)          (mpfr_get_exp) (x)
   1056 # define MPFR_SET_EXP(x,e)        (MPFR_ASSERTN (MPFR_EXP_IN_RANGE (e)), \
   1057                                    (void) (MPFR_EXP (x) = (e)))
   1058 # define MPFR_SET_INVALID_EXP(x)  ((void) (MPFR_EXP (x) = MPFR_EXP_INVALID))
   1059 #else
   1060 # define MPFR_GET_EXP(x)          MPFR_EXP (x)
   1061 # define MPFR_SET_EXP(x,e)        ((void) (MPFR_EXP (x) = (e)))
   1062 # define MPFR_SET_INVALID_EXP(x)  ((void) 0)
   1063 #endif
   1064 
   1065 /* Compare the exponents of two numbers, which can be either MPFR numbers
   1066    or UBF numbers. */
   1067 #define MPFR_UBF_EXP_LESS_P(x,y) \
   1068   (MPFR_UNLIKELY (MPFR_IS_UBF (x) || MPFR_IS_UBF (y)) ? \
   1069    mpfr_ubf_exp_less_p (x, y) : MPFR_GET_EXP (x) < MPFR_GET_EXP (y))
   1070 
   1071 
   1072 /******************************************************
   1073  *********  Singular values (NAN, INF, ZERO)  *********
   1074  ******************************************************/
   1075 
   1076 /* Enum special value of exponent. */
   1077 # define MPFR_EXP_ZERO (MPFR_EXP_MIN+1)
   1078 # define MPFR_EXP_NAN  (MPFR_EXP_MIN+2)
   1079 # define MPFR_EXP_INF  (MPFR_EXP_MIN+3)
   1080 # define MPFR_EXP_UBF  (MPFR_EXP_MIN+4)
   1081 
   1082 #define MPFR_IS_NAN(x)   (MPFR_EXP(x) == MPFR_EXP_NAN)
   1083 #define MPFR_SET_NAN(x)  (MPFR_EXP(x) =  MPFR_EXP_NAN)
   1084 #define MPFR_IS_INF(x)   (MPFR_EXP(x) == MPFR_EXP_INF)
   1085 #define MPFR_SET_INF(x)  (MPFR_EXP(x) =  MPFR_EXP_INF)
   1086 #define MPFR_IS_ZERO(x)  (MPFR_EXP(x) == MPFR_EXP_ZERO)
   1087 #define MPFR_SET_ZERO(x) (MPFR_EXP(x) =  MPFR_EXP_ZERO)
   1088 #define MPFR_NOTZERO(x)  (MPFR_EXP(x) != MPFR_EXP_ZERO)
   1089 #define MPFR_IS_UBF(x)   (MPFR_EXP(x) == MPFR_EXP_UBF)
   1090 #define MPFR_SET_UBF(x)  (MPFR_EXP(x) =  MPFR_EXP_UBF)
   1091 
   1092 #define MPFR_IS_NORMALIZED(x) \
   1093   (MPFR_LIMB_MSB (MPFR_MANT(x)[MPFR_LAST_LIMB(x)]) != 0)
   1094 
   1095 #define MPFR_IS_FP(x)       (!MPFR_IS_NAN(x) && !MPFR_IS_INF(x))
   1096 
   1097 /* Note: contrary to the MPFR_IS_PURE_*(x) macros, the MPFR_IS_SINGULAR*(x)
   1098    macros may be used even when x is being constructed, i.e. its exponent
   1099    field is already set (possibly out-of-range), but its significand field
   1100    may still contain arbitrary data. Thus MPFR_IS_PURE_FP(x) is not always
   1101    equivalent to !MPFR_IS_SINGULAR(x); see the code below. */
   1102 #define MPFR_IS_SINGULAR(x) (MPFR_EXP(x) <= MPFR_EXP_INF)
   1103 #define MPFR_IS_SINGULAR_OR_UBF(x) (MPFR_EXP(x) <= MPFR_EXP_UBF)
   1104 
   1105 /* The following two macros return true iff the value is a regular number,
   1106    i.e. it is not a singular number. In debug mode, the format is also
   1107    checked: valid exponent, but potentially out of range; normalized value.
   1108    In contexts where UBF's are not accepted or not possible, MPFR_IS_PURE_FP
   1109    is preferable. If UBF's are accepted, MPFR_IS_PURE_UBF must be used. */
   1110 #define MPFR_IS_PURE_FP(x)                          \
   1111   (!MPFR_IS_SINGULAR(x) &&                          \
   1112    (MPFR_ASSERTD (MPFR_EXP (x) >= MPFR_EMIN_MIN &&  \
   1113                   MPFR_EXP (x) <= MPFR_EMAX_MAX &&  \
   1114                   MPFR_IS_NORMALIZED (x)), 1))
   1115 #define MPFR_IS_PURE_UBF(x)                             \
   1116   (!MPFR_IS_SINGULAR(x) &&                              \
   1117    (MPFR_ASSERTD ((MPFR_IS_UBF (x) ||                   \
   1118                    (MPFR_EXP (x) >= MPFR_EMIN_MIN &&    \
   1119                     MPFR_EXP (x) <= MPFR_EMAX_MAX)) &&  \
   1120                   MPFR_IS_NORMALIZED (x)), 1))
   1121 
   1122 /* Ditto for 2 numbers. */
   1123 #define MPFR_ARE_SINGULAR(x,y) \
   1124   (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)) || MPFR_UNLIKELY(MPFR_IS_SINGULAR(y)))
   1125 #define MPFR_ARE_SINGULAR_OR_UBF(x,y)           \
   1126   (MPFR_UNLIKELY(MPFR_IS_SINGULAR_OR_UBF(x)) || \
   1127    MPFR_UNLIKELY(MPFR_IS_SINGULAR_OR_UBF(y)))
   1128 
   1129 
   1130 /******************************************************
   1131  ********************  Sign macros  *******************
   1132  ******************************************************/
   1133 
   1134 /* These are sign macros for MPFR numbers only. */
   1135 
   1136 #define MPFR_SIGN_POS (1)
   1137 #define MPFR_SIGN_NEG (-1)
   1138 
   1139 #define MPFR_IS_STRICTPOS(x)  (MPFR_NOTZERO (x) && MPFR_IS_POS (x))
   1140 #define MPFR_IS_STRICTNEG(x)  (MPFR_NOTZERO (x) && MPFR_IS_NEG (x))
   1141 
   1142 #define MPFR_IS_NEG(x) (MPFR_SIGN(x) < 0)
   1143 #define MPFR_IS_POS(x) (MPFR_SIGN(x) > 0)
   1144 
   1145 #define MPFR_SET_POS(x) (MPFR_SIGN(x) = MPFR_SIGN_POS)
   1146 #define MPFR_SET_NEG(x) (MPFR_SIGN(x) = MPFR_SIGN_NEG)
   1147 
   1148 #define MPFR_CHANGE_SIGN(x) (MPFR_SIGN(x) = -MPFR_SIGN(x))
   1149 #define MPFR_SET_SAME_SIGN(x, y) (MPFR_SIGN(x) = MPFR_SIGN(y))
   1150 #define MPFR_SET_OPPOSITE_SIGN(x, y) (MPFR_SIGN(x) = -MPFR_SIGN(y))
   1151 #define MPFR_ASSERT_SIGN(s) \
   1152  (MPFR_ASSERTD((s) == MPFR_SIGN_POS || (s) == MPFR_SIGN_NEG))
   1153 #define MPFR_SET_SIGN(x, s) \
   1154   (MPFR_ASSERT_SIGN(s), MPFR_SIGN(x) = s)
   1155 #define MPFR_IS_POS_SIGN(s1) ((s1) > 0)
   1156 #define MPFR_IS_NEG_SIGN(s1) ((s1) < 0)
   1157 #define MPFR_MULT_SIGN(s1, s2) ((s1) * (s2))
   1158 /* Transform a sign to 1 or -1 */
   1159 #define MPFR_FROM_SIGN_TO_INT(s) (s)
   1160 #define MPFR_INT_SIGN(x) MPFR_FROM_SIGN_TO_INT(MPFR_SIGN(x))
   1161 
   1162 
   1163 /******************************************************
   1164  ***************  Ternary value macros  ***************
   1165  ******************************************************/
   1166 
   1167 /* Special inexact value */
   1168 #define MPFR_EVEN_INEX 2
   1169 
   1170 /* Note: the addition/subtraction of 2 comparisons below instead of the
   1171    use of the ?: operator allows GCC and Clang to generate better code
   1172    in general; this is the case at least with GCC on x86 (32 & 64 bits),
   1173    PowerPC and Aarch64 (64-bit ARM), and with Clang on x86_64.
   1174    VSIGN code based on mini-gmp's GMP_CMP macro; adapted for INEXPOS. */
   1175 
   1176 /* Macros for functions returning two inexact values in an 'int'
   1177    (exact = 0, positive = 1, negative = 2) */
   1178 #define INEXPOS(y) (((y) != 0) + ((y) < 0))
   1179 #define INEX(y,z) (INEXPOS(y) | (INEXPOS(z) << 2))
   1180 
   1181 /* When returning the ternary inexact value, ALWAYS use one of the
   1182    following two macros, unless the flag comes from another function
   1183    returning the ternary inexact value */
   1184 #define MPFR_RET(I) return \
   1185   (I) != 0 ? ((__gmpfr_flags |= MPFR_FLAGS_INEXACT), (I)) : 0
   1186 #define MPFR_RET_NAN return (__gmpfr_flags |= MPFR_FLAGS_NAN), 0
   1187 
   1188 /* Sign of a native value. */
   1189 #define VSIGN(I) (((I) > 0) - ((I) < 0))
   1190 #define SAME_SIGN(I1,I2) (VSIGN (I1) == VSIGN (I2))
   1191 
   1192 
   1193 /******************************************************
   1194  ***************  Rounding mode macros  ***************
   1195  ******************************************************/
   1196 
   1197 /* MPFR_RND_MAX gives the number of supported rounding modes by all functions.
   1198  */
   1199 #define MPFR_RND_MAX ((mpfr_rnd_t)((MPFR_RNDF)+1))
   1200 
   1201 /* We want to test this :
   1202  *  (rnd == MPFR_RNDU && test) || (rnd == RNDD && !test)
   1203  * i.e. it transforms RNDU or RNDD to away or zero according to the sign.
   1204  * The argument test must be 0 or 1. */
   1205 #define MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd, test) \
   1206   (MPFR_ASSERTD ((test) == 0 || (test) == 1),      \
   1207    ((rnd) + (test)) == MPFR_RNDD)
   1208 
   1209 /* We want to test if rnd rounds toward zero or away from zero.
   1210    'neg' is 1 if negative, and 0 if positive. */
   1211 #define MPFR_IS_LIKE_RNDZ(rnd, neg) \
   1212   ((rnd) == MPFR_RNDZ || MPFR_IS_RNDUTEST_OR_RNDDNOTTEST (rnd, neg))
   1213 #define MPFR_IS_LIKE_RNDA(rnd, neg) \
   1214   ((rnd) == MPFR_RNDA || MPFR_IS_RNDUTEST_OR_RNDDNOTTEST (rnd, (neg) == 0))
   1215 
   1216 #define MPFR_IS_LIKE_RNDU(rnd, sign)                    \
   1217   (((rnd) == MPFR_RNDU) ||                              \
   1218    ((rnd) == MPFR_RNDZ && MPFR_IS_NEG_SIGN (sign)) ||   \
   1219    ((rnd) == MPFR_RNDA && MPFR_IS_POS_SIGN (sign)))
   1220 
   1221 #define MPFR_IS_LIKE_RNDD(rnd, sign)                    \
   1222   (((rnd) == MPFR_RNDD) ||                              \
   1223    ((rnd) == MPFR_RNDZ && MPFR_IS_POS_SIGN (sign)) ||   \
   1224    ((rnd) == MPFR_RNDA && MPFR_IS_NEG_SIGN (sign)))
   1225 
   1226 /* Invert RNDU and RNDD; the other rounding modes are unchanged. */
   1227 #define MPFR_INVERT_RND(rnd) ((rnd) == MPFR_RNDU ? MPFR_RNDD :          \
   1228                               (rnd) == MPFR_RNDD ? MPFR_RNDU : (rnd))
   1229 
   1230 /* Transform RNDU and RNDD to RNDZ according to test */
   1231 #define MPFR_UPDATE_RND_MODE(rnd, test)                             \
   1232   do {                                                              \
   1233     if (MPFR_UNLIKELY(MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd, test)))  \
   1234       rnd = MPFR_RNDZ;                                              \
   1235   } while (0)
   1236 
   1237 /* Transform RNDU and RNDD to RNDZ or RNDA according to sign;
   1238    leave the other modes unchanged.
   1239    Usage: MPFR_UPDATE2_RND_MODE (rnd_mode, MPFR_SIGN (x)) */
   1240 #define MPFR_UPDATE2_RND_MODE(rnd, sign)                       \
   1241   do {                                                         \
   1242     if (rnd == MPFR_RNDU)                                      \
   1243       rnd = MPFR_IS_POS_SIGN (sign) ? MPFR_RNDA : MPFR_RNDZ;   \
   1244     else if (rnd == MPFR_RNDD)                                 \
   1245       rnd = MPFR_IS_NEG_SIGN (sign) ? MPFR_RNDA : MPFR_RNDZ;   \
   1246   } while (0)
   1247 
   1248 
   1249 /******************************************************
   1250  ******************  Limb macros  *********************
   1251  ******************************************************/
   1252 
   1253 /* MPFR_LIMB: Cast to mp_limb_t, assuming that x is based on mp_limb_t
   1254    variables (needed when mp_limb_t is defined as an integer type shorter
   1255    than int, due to the integer promotion rules, which is possible only
   1256    if MPFR_LONG_WITHIN_LIMB is not defined). Warning! This will work
   1257    only when the computed value x is congruent to the expected value
   1258    modulo MPFR_LIMB_MAX + 1. Be aware that this macro may not solve all
   1259    the problems related to the integer promotion rules, because it won't
   1260    have an influence on the evaluation of x itself. Hence the need for...
   1261 
   1262    MPFR_LIMB_LSHIFT: Left shift by making sure that the shifted argument
   1263    is unsigned (use unsigned long due to the MPFR_LONG_WITHIN_LIMB test).
   1264    For instance, due to the integer promotion rules, if mp_limb_t is
   1265    defined as a 16-bit unsigned short and an int has 32 bits, then a
   1266    mp_limb_t will be converted to an int, which is signed.
   1267 */
   1268 #ifdef MPFR_LONG_WITHIN_LIMB
   1269 #define MPFR_LIMB(x) (x)
   1270 #define MPFR_LIMB_LSHIFT(x,c) ((x) << (c))
   1271 #else
   1272 #define MPFR_LIMB(x) ((mp_limb_t) (x))
   1273 #define MPFR_LIMB_LSHIFT(x,c) MPFR_LIMB((unsigned long) (x) << (c))
   1274 #endif
   1275 
   1276 /* Definition of simple mp_limb_t constants */
   1277 #define MPFR_LIMB_ZERO    ((mp_limb_t) 0)
   1278 #define MPFR_LIMB_ONE     ((mp_limb_t) 1)
   1279 #define MPFR_LIMB_HIGHBIT MPFR_LIMB_LSHIFT (MPFR_LIMB_ONE, GMP_NUMB_BITS - 1)
   1280 #define MPFR_LIMB_MAX     ((mp_limb_t) -1)
   1281 
   1282 /* Mask to get the Most Significant Bit of a limb */
   1283 #define MPFR_LIMB_MSB(l) ((mp_limb_t) ((l) & MPFR_LIMB_HIGHBIT))
   1284 
   1285 /* Mask for the low 's' bits of a limb */
   1286 #define MPFR_LIMB_MASK(s)                                               \
   1287   (MPFR_ASSERTD (s >= 0 && s <= GMP_NUMB_BITS),                         \
   1288    s == GMP_NUMB_BITS ? MPFR_LIMB_MAX :                                 \
   1289    (mp_limb_t) (MPFR_LIMB_LSHIFT (MPFR_LIMB_ONE, s) - MPFR_LIMB_ONE))
   1290 
   1291 /******************************************************
   1292  **********************  Memory  **********************
   1293  ******************************************************/
   1294 
   1295 #define MPFR_BYTES_PER_MP_LIMB (GMP_NUMB_BITS/CHAR_BIT)
   1296 
   1297 /* Heap memory handling
   1298    --------------------
   1299    Memory allocated for a significand (mantissa) has the following
   1300    format:
   1301      * A mp_size_t in a mpfr_size_limb_t union (see below).
   1302      * An array of mp_limb_t (not all of them are necessarily used,
   1303        as the precision can change without a reallocation).
   1304    The goal of the mpfr_size_limb_t union is to make sure that
   1305    size and alignment requirements are satisfied if mp_size_t and
   1306    mp_limb_t have different sizes and/or alignment requirements.
   1307    And the casts to void * prevents the compiler from emitting a
   1308    warning (or error), such as:
   1309      cast increases required alignment of target type
   1310    with the -Wcast-align GCC option. Correct alignment is checked
   1311    by MPFR_SET_MANT_PTR (when setting MPFR_MANT(x), the MPFR code
   1312    should use this macro or guarantee a correct alignment at this
   1313    time).
   1314    Moreover, pointer conversions are not fully specified by the
   1315    C standard, and the use of a union (and the double casts below)
   1316    might help even if mp_size_t and mp_limb_t have the same size
   1317    and the same alignment requirements. Still, there is currently
   1318    no guarantee that this code is portable. Note that union members
   1319    are not used at all.
   1320 */
   1321 typedef union { mp_size_t s; mp_limb_t l; } mpfr_size_limb_t;
   1322 #define MPFR_GET_ALLOC_SIZE(x) \
   1323   (((mp_size_t *) (void *) MPFR_MANT(x))[-1] + 0)
   1324 #define MPFR_SET_ALLOC_SIZE(x, n) \
   1325   (((mp_size_t *) (void *) MPFR_MANT(x))[-1] = (n))
   1326 #define MPFR_MALLOC_SIZE(s) \
   1327   (sizeof(mpfr_size_limb_t) + MPFR_BYTES_PER_MP_LIMB * (size_t) (s))
   1328 #define MPFR_SET_MANT_PTR(x,p) \
   1329   (MPFR_MANT(x) = (mp_limb_t *) ((mpfr_size_limb_t *) (p) + 1))
   1330 #define MPFR_GET_REAL_PTR(x) \
   1331   ((void *) ((mpfr_size_limb_t *) (void *) MPFR_MANT(x) - 1))
   1332 
   1333 /* Temporary memory handling */
   1334 #ifndef TMP_SALLOC
   1335 /* GMP 4.1.x or below or internals */
   1336 #define MPFR_TMP_DECL TMP_DECL
   1337 #define MPFR_TMP_MARK TMP_MARK
   1338 #define MPFR_TMP_ALLOC TMP_ALLOC
   1339 #define MPFR_TMP_FREE TMP_FREE
   1340 #else
   1341 #define MPFR_TMP_DECL(x) TMP_DECL
   1342 #define MPFR_TMP_MARK(x) TMP_MARK
   1343 #define MPFR_TMP_ALLOC(s) TMP_ALLOC(s)
   1344 #define MPFR_TMP_FREE(x) TMP_FREE
   1345 #endif
   1346 
   1347 #define MPFR_TMP_LIMBS_ALLOC(N) \
   1348   ((mp_limb_t *) MPFR_TMP_ALLOC ((size_t) (N) * MPFR_BYTES_PER_MP_LIMB))
   1349 
   1350 /* The temporary var doesn't have any size field, but it doesn't matter
   1351  * since only functions dealing with the Heap care about it */
   1352 #define MPFR_TMP_INIT1(xp, x, p)                                     \
   1353  ( MPFR_PREC(x) = (p),                                               \
   1354    MPFR_MANT(x) = (xp),                                              \
   1355    MPFR_SET_POS(x),                                                  \
   1356    MPFR_SET_INVALID_EXP(x))
   1357 
   1358 #define MPFR_TMP_INIT(xp, x, p, s)                                   \
   1359   (xp = MPFR_TMP_LIMBS_ALLOC(s),                                     \
   1360    MPFR_TMP_INIT1(xp, x, p))
   1361 
   1362 /* Set y to s*significand(x)*2^e, for example MPFR_ALIAS(y,x,1,MPFR_EXP(x))
   1363    sets y to |x|, and MPFR_ALIAS(y,x,MPFR_SIGN(x),0) sets y to x*2^f such
   1364    that 1/2 <= |y| < 1. Does not check y is in the valid exponent range.
   1365    WARNING! x and y share the same mantissa. So, some operations are
   1366    not valid if x has been provided via an argument, e.g., trying to
   1367    modify the mantissa of y, even temporarily, or calling mpfr_clear on y.
   1368 */
   1369 #define MPFR_ALIAS(y,x,s,e)                     \
   1370   (MPFR_PREC(y) = MPFR_PREC(x),                 \
   1371    MPFR_SIGN(y) = (s),                          \
   1372    MPFR_EXP(y) = (e),                           \
   1373    MPFR_MANT(y) = MPFR_MANT(x))
   1374 
   1375 #define MPFR_TMP_INIT_ABS(y,x) \
   1376   MPFR_ALIAS (y, x, MPFR_SIGN_POS, MPFR_EXP (x))
   1377 
   1378 #define MPFR_TMP_INIT_NEG(y,x) \
   1379   MPFR_ALIAS (y, x, - MPFR_SIGN (x), MPFR_EXP (x))
   1380 
   1381 
   1382 /******************************************************
   1383  *******************  Cache macros  *******************
   1384  ******************************************************/
   1385 
   1386 /* Cache struct */
   1387 #define mpfr_const_pi(_d,_r)    mpfr_cache(_d, __gmpfr_cache_const_pi,_r)
   1388 #define mpfr_const_log2(_d,_r)  mpfr_cache(_d, __gmpfr_cache_const_log2, _r)
   1389 #define mpfr_const_euler(_d,_r) mpfr_cache(_d, __gmpfr_cache_const_euler, _r)
   1390 #define mpfr_const_catalan(_d,_r) mpfr_cache(_d,__gmpfr_cache_const_catalan,_r)
   1391 
   1392 /* Declare a global cache for a MPFR constant.
   1393    If the shared cache is enabled, and if the constructor/destructor
   1394    attributes are available, we need to initialize the shared lock of
   1395    the cache with a constructor. It is the goal of the macro
   1396    MPFR_DEFERRED_INIT_MASTER_DECL.
   1397    FIXME: When MPFR is built with shared cache, the field "lock" is
   1398    not explicitly initialized, which can yield a warning, e.g. with
   1399    GCC's -Wmissing-field-initializers (and an error with -Werror).
   1400    Since one does not know what is behind the associated typedef name,
   1401    one cannot provide an explicit initialization for such a type. Two
   1402    possible solutions:
   1403      1. Encapsulate the type in a structure or a union and use the
   1404         universal zero initializer: { 0 }
   1405         But: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80454
   1406      2. Use designated initializers when supported. But this needs a
   1407         configure test.
   1408    Using a diagnostic pragma to ignore the warning in this particular case
   1409    is not really possible, because the warning occurs when the macro is
   1410    expanded and one cannot put a pragma in the contents of a #define.
   1411 */
   1412 #define MPFR_DECL_INIT_CACHE(_cache,_func)                           \
   1413   MPFR_DEFERRED_INIT_MASTER_DECL(_func,                              \
   1414                                  MPFR_LOCK_INIT( (_cache)->lock),    \
   1415                                  MPFR_LOCK_CLEAR((_cache)->lock))    \
   1416   MPFR_CACHE_ATTR mpfr_cache_t _cache = {{                           \
   1417       {{ 0, MPFR_SIGN_POS, 0, (mp_limb_t *) 0 }}, 0, _func           \
   1418       MPFR_DEFERRED_INIT_SLAVE_VALUE(_func)                          \
   1419     }};                                                              \
   1420   MPFR_MAKE_VARFCT (mpfr_cache_t,_cache)
   1421 
   1422 /******************************************************
   1423  ***************  Threshold parameters  ***************
   1424  ******************************************************/
   1425 
   1426 #include "mparam.h"
   1427 
   1428 
   1429 /******************************************************
   1430  ******************  Useful macros  *******************
   1431  ******************************************************/
   1432 
   1433 /* The MAX, MIN and ABS macros may already be defined if gmp-impl.h has
   1434    been included. They have the same semantics as in gmp-impl.h, but the
   1435    expressions may be slightly different. So, it's better to undefine
   1436    them first, as required by the ISO C standard. */
   1437 #undef MAX
   1438 #undef MIN
   1439 #undef ABS
   1440 #define MAX(a, b) (((a) > (b)) ? (a) : (b))
   1441 #define MIN(a, b) (((a) < (b)) ? (a) : (b))
   1442 #define ABS(x) (((x)>0) ? (x) : -(x))
   1443 
   1444 /* These macros help the compiler to determine if a test is
   1445    likely or unlikely. The !! is necessary in case x is larger
   1446    than a long. */
   1447 #if defined MPFR_DEBUG_PREDICTION && __MPFR_GNUC(3,0)
   1448 
   1449 /* Code to debug branch prediction, based on Ulrich Drepper's paper
   1450  * "What Every Programmer Should Know About Memory":
   1451  *   https://people.freebsd.org/~lstewart/articles/cpumemory.pdf
   1452  */
   1453 asm (".section predict_data, \"aw\"; .previous\n"
   1454      ".section predict_line, \"a\"; .previous\n"
   1455      ".section predict_file, \"a\"; .previous");
   1456 # if defined __x86_64__
   1457 #  define MPFR_DEBUGPRED__(e,E)                                         \
   1458   ({ long _e = !!(e);                                                   \
   1459     asm volatile (".pushsection predict_data\n"                         \
   1460                   "..predictcnt%=: .quad 0; .quad 0\n"                  \
   1461                   ".section predict_line; .quad %c1\n"                  \
   1462                   ".section predict_file; .quad %c2; .popsection\n"     \
   1463                   "addq $1,..predictcnt%=(,%0,8)"                       \
   1464                   : : "r" (_e == E), "i" (__LINE__), "i" (__FILE__));   \
   1465     __builtin_expect (_e, E);                                           \
   1466   })
   1467 # elif defined __i386__
   1468 #  define MPFR_DEBUGPRED__(e,E)                                         \
   1469   ({ long _e = !!(e);                                                   \
   1470     asm volatile (".pushsection predict_data\n"                         \
   1471                   "..predictcnt%=: .long 0; .long 0\n"                  \
   1472                   ".section predict_line; .long %c1\n"                  \
   1473                   ".section predict_file; .long %c2; .popsection\n"     \
   1474                   "incl ..predictcnt%=(,%0,4)"                          \
   1475                   : : "r" (_e == E), "i" (__LINE__), "i" (__FILE__));   \
   1476     __builtin_expect (_e, E);                                           \
   1477   })
   1478 # else
   1479 #  error "MPFR_DEBUGPRED__ definition missing"
   1480 # endif
   1481 
   1482 # define MPFR_LIKELY(x) MPFR_DEBUGPRED__ ((x), 1)
   1483 # define MPFR_UNLIKELY(x) MPFR_DEBUGPRED__ ((x), 0)
   1484 
   1485 #elif __MPFR_GNUC(3,0) || __MPFR_ICC(8,1,0)
   1486 
   1487 # define MPFR_LIKELY(x) (__builtin_expect(!!(x), 1))
   1488 # define MPFR_UNLIKELY(x) (__builtin_expect(!!(x), 0))
   1489 
   1490 #else
   1491 
   1492 # define MPFR_LIKELY(x) (x)
   1493 # define MPFR_UNLIKELY(x) (x)
   1494 
   1495 #endif
   1496 
   1497 /* Declare that some variable is initialized before being used (without a
   1498    dummy initialization) in order to avoid some compiler warnings. Use the
   1499    VAR = VAR trick (see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=36296#c3)
   1500    only with gcc as this is undefined behavior, and we don't know what other
   1501    compilers do (they may also be smarter). This self-initialization trick
   1502    could be disabled with future gcc versions. */
   1503 #if defined(__GNUC__)
   1504 # define INITIALIZED(VAR) VAR = VAR
   1505 #else
   1506 # define INITIALIZED(VAR) VAR
   1507 #endif
   1508 
   1509 /* Ceil log 2: If GCC, uses a GCC extension, otherwise calls a function */
   1510 /* Warning:
   1511  *   Needs to define MPFR_NEED_LONGLONG.
   1512  *   Computes ceil(log2(x)) only for x integer (unsigned long)
   1513  *   Undefined if x is 0 */
   1514 #if __MPFR_GNUC(2,95) || __MPFR_ICC(8,1,0)
   1515 /* Note: This macro MPFR_INT_CEIL_LOG2 shouldn't be used in an MPFR_ASSERT*
   1516    macro, either directly or indirectly via other macros, otherwise it can
   1517    yield an error due to a too large stringized expression in ASSERT_FAIL.
   1518    A static inline function could be a better solution than this macro. */
   1519 /* FIXME: The current code assumes that x fits in an unsigned long
   1520    (used by __gmpfr_int_ceil_log2) while MPFR_INT_CEIL_LOG2 is used on
   1521    values that might be larger than ULONG_MAX on some platforms and/or
   1522    with some build options; a loop could be used if x > ULONG_MAX. If
   1523    the type of x is <= unsigned long, then no additional code will be
   1524    generated thanks to obvious compiler optimization. */
   1525 #ifdef MPFR_LONG_WITHIN_LIMB
   1526 # define MPFR_INT_CEIL_LOG2(x)                            \
   1527     (MPFR_UNLIKELY ((x) == 1) ? 0 :                       \
   1528      __extension__ ({ int _b; mp_limb_t _limb;            \
   1529       MPFR_ASSERTN ((x) > 1);                             \
   1530       _limb = (x) - 1;                                    \
   1531       MPFR_ASSERTN (_limb == (x) - 1);                    \
   1532       count_leading_zeros (_b, _limb);                    \
   1533       _b = GMP_NUMB_BITS - _b;                            \
   1534       MPFR_ASSERTD (_b >= 0);                             \
   1535       _b; }))
   1536 #else
   1537 # define MPFR_INT_CEIL_LOG2(x)                              \
   1538   (MPFR_UNLIKELY ((x) == 1) ? 0 :                           \
   1539    __extension__ ({ int _c = 0; unsigned long _x = (x) - 1; \
   1540        MPFR_ASSERTN ((x) > 1);                              \
   1541        while (_x != 0)                                      \
   1542          {                                                  \
   1543            _x = _x >> 1;                                    \
   1544            _c ++;                                           \
   1545          };                                                 \
   1546        MPFR_ASSERTD (_c >= 0);                              \
   1547        _c; }))
   1548 #endif /* MPFR_LONG_WITHIN_LIMB */
   1549 #else
   1550 # define MPFR_INT_CEIL_LOG2(x) \
   1551   (MPFR_ASSERTN (x <= ULONG_MAX), __gmpfr_int_ceil_log2(x))
   1552 #endif /* __MPFR_GNUC(2,95) || __MPFR_ICC(8,1,0) */
   1553 
   1554 /* Add two integers with overflow handling */
   1555 /* Example: MPFR_SADD_OVERFLOW (c, a, b, long, unsigned long,
   1556  *                              LONG_MIN, LONG_MAX,
   1557  *                              goto overflow, goto underflow); */
   1558 #define MPFR_UADD_OVERFLOW(c,a,b,ACTION_IF_OVERFLOW)                  \
   1559  do {                                                                 \
   1560   (c) = (a) + (b);                                                    \
   1561   if ((c) < (a)) ACTION_IF_OVERFLOW;                                  \
   1562  } while (0)
   1563 
   1564 #define MPFR_SADD_OVERFLOW(c,a,b,STYPE,UTYPE,MIN,MAX,ACTION_IF_POS_OVERFLOW,ACTION_IF_NEG_OVERFLOW) \
   1565   do {                                                                \
   1566   if ((a) >= 0 && (b) >= 0) {                                         \
   1567          UTYPE uc,ua,ub;                                              \
   1568          ua = (UTYPE) (a); ub = (UTYPE) (b);                          \
   1569          MPFR_UADD_OVERFLOW (uc, ua, ub, ACTION_IF_POS_OVERFLOW);     \
   1570          if (uc > (UTYPE)(MAX)) ACTION_IF_POS_OVERFLOW;               \
   1571          else (c) = (STYPE) uc;                                       \
   1572   } else if ((a) < 0 && (b) < 0) {                                    \
   1573          UTYPE uc,ua,ub;                                              \
   1574          ua = -(UTYPE) (a); ub = -(UTYPE) (b);                        \
   1575          MPFR_UADD_OVERFLOW (uc, ua, ub, ACTION_IF_NEG_OVERFLOW);     \
   1576          if (uc >= -(UTYPE)(MIN) || uc > (UTYPE)(MAX)) {              \
   1577            if (uc ==  -(UTYPE)(MIN)) (c) = (MIN);                     \
   1578            else  ACTION_IF_NEG_OVERFLOW; }                            \
   1579          else (c) = -(STYPE) uc;                                      \
   1580   } else (c) = (a) + (b);                                             \
   1581  } while (0)
   1582 
   1583 
   1584 /* Set a number to 1 (Fast) - It doesn't check if 1 is in the exponent range */
   1585 #define MPFR_SET_ONE(x)                                               \
   1586 do {                                                                  \
   1587   mp_size_t _size = MPFR_LAST_LIMB(x);                                \
   1588   MPFR_SET_POS(x);                                                    \
   1589   MPFR_EXP(x) = 1;                                                    \
   1590   MPN_ZERO ( MPFR_MANT(x), _size);                                    \
   1591   MPFR_MANT(x)[_size] = MPFR_LIMB_HIGHBIT;                            \
   1592 } while (0)
   1593 
   1594 /* Compute s = (-a) % GMP_NUMB_BITS as unsigned */
   1595 #define MPFR_UNSIGNED_MINUS_MODULO(s, a)                              \
   1596   do                                                                  \
   1597     {                                                                 \
   1598       if (IS_POW2 (GMP_NUMB_BITS))                                    \
   1599         (s) = (- (unsigned int) (a)) % GMP_NUMB_BITS;                 \
   1600       else                                                            \
   1601         {                                                             \
   1602           (s) = (a) % GMP_NUMB_BITS;                                  \
   1603           if ((s) != 0)                                               \
   1604             (s) = GMP_NUMB_BITS - (s);                                \
   1605         }                                                             \
   1606       MPFR_ASSERTD ((s) >= 0 && (s) < GMP_NUMB_BITS);                 \
   1607     }                                                                 \
   1608   while (0)
   1609 
   1610 /* Test if X (positive) is a power of 2 */
   1611 #define IS_POW2(X) (((X) & ((X) - 1)) == 0)
   1612 #define NOT_POW2(X) (((X) & ((X) - 1)) != 0)
   1613 
   1614 /* Safe absolute value and difference (to avoid possible integer overflow) */
   1615 /* type is the target (unsigned) type */
   1616 #define SAFE_ABS(type,x) ((x) >= 0 ? (type)(x) : -(type)(x))
   1617 #define SAFE_DIFF(type,x,y) (MPFR_ASSERTD((x) >= (y)), (type)(x) - (type)(y))
   1618 
   1619 #define ULONG2LONG(U) ((U) > LONG_MAX ? -1 - (long) ~(U) : (long) (U))
   1620 
   1621 /* Check whether an integer type (after integer promotion) is signed.
   1622    This can be determined at compilation time, but unfortunately,
   1623    when used in practice, this is not a constant expression (because
   1624    the argument X is not a constant expression, even though the result
   1625    does not depend on its value), so that this cannot be used for a
   1626    static assertion. */
   1627 #define IS_SIGNED(X) ((X) * 0 - 1 < 0)
   1628 
   1629 #define mpfr_get_d1(x) mpfr_get_d(x,__gmpfr_default_rounding_mode)
   1630 
   1631 /* Store in r the size in bits of the mpz_t z */
   1632 #define MPFR_MPZ_SIZEINBASE2(r, z)                      \
   1633   do {                                                  \
   1634     int _cnt;                                           \
   1635     mp_size_t _size;                                    \
   1636     MPFR_ASSERTD (mpz_sgn (z) != 0);                    \
   1637     _size = ABSIZ(z);                                   \
   1638     MPFR_ASSERTD (_size >= 1);                          \
   1639     count_leading_zeros (_cnt, PTR(z)[_size-1]);        \
   1640     (r) = (mp_bitcnt_t) _size * GMP_NUMB_BITS - _cnt;   \
   1641   } while (0)
   1642 
   1643 /* MPFR_LCONV_DPTS can also be forced to 0 or 1 by the user. */
   1644 #ifndef MPFR_LCONV_DPTS
   1645 # if defined(HAVE_LOCALE_H) && \
   1646      defined(HAVE_STRUCT_LCONV_DECIMAL_POINT) && \
   1647      defined(HAVE_STRUCT_LCONV_THOUSANDS_SEP)
   1648 #  define MPFR_LCONV_DPTS 1
   1649 # else
   1650 #  define MPFR_LCONV_DPTS 0
   1651 # endif
   1652 #endif
   1653 
   1654 /* FIXME: Add support for multibyte decimal_point and thousands_sep since
   1655    this can be found in practice: https://reviews.llvm.org/D27167 says:
   1656    "I found this problem on FreeBSD 11, where thousands_sep in fr_FR.UTF-8
   1657    is a no-break space (U+00A0)."
   1658    Under Linux, this is U+202F NARROW NO-BREAK SPACE (e2 80 af).
   1659    And in the ps_AF locale,
   1660      decimal_point = U+066B ARABIC DECIMAL SEPARATOR (d9 ab)
   1661      thousands_sep = U+066C ARABIC THOUSANDS SEPARATOR (d9 ac)
   1662    In the mean time, in case of non-single-byte character, revert to the
   1663    default value. */
   1664 #if MPFR_LCONV_DPTS
   1665 #include <locale.h>
   1666 /* Warning! In case of signed char, the value of MPFR_DECIMAL_POINT may
   1667    be negative (the ISO C99 does not seem to forbid negative values). */
   1668 #define MPFR_DECIMAL_POINT                      \
   1669   (localeconv()->decimal_point[1] != '\0' ?     \
   1670    (char) '.' : localeconv()->decimal_point[0])
   1671 #define MPFR_THOUSANDS_SEPARATOR                \
   1672   (localeconv()->thousands_sep[0] == '\0' ||    \
   1673    localeconv()->thousands_sep[1] != '\0' ?     \
   1674    (char) '\0' : localeconv()->thousands_sep[0])
   1675 #else
   1676 #define MPFR_DECIMAL_POINT ((char) '.')
   1677 #define MPFR_THOUSANDS_SEPARATOR ((char) '\0')
   1678 #endif
   1679 
   1680 /* Size of an array, as a constant expression. */
   1681 #define numberof_const(x)  (sizeof (x) / sizeof ((x)[0]))
   1682 
   1683 /* Size of an array, safe version but not a constant expression:
   1684    Since an array can silently be converted to a pointer, we check
   1685    that this macro is applied on an array, not a pointer.
   1686    Also make sure that the type is signed ("long" is sufficient
   1687    in practice since the sizes come from the MPFR source), so that
   1688    the value can be used in arbitrary expressions without the risk
   1689    of silently switching to unsigned arithmetic. */
   1690 /* TODO: Make numberof() a constant expression and always use it in
   1691    the MPFR code instead of numberof_const(). See the tricks at
   1692      https://gcc.gnu.org/pipermail/gcc/2020-September/233763.html
   1693      "[PATCH v2] <sys/param.h>: Add nitems() and snitems() macros"
   1694      by Alejandro Colomar
   1695    but this needs to be fully tested on various platforms and with
   1696    various compilers and compilation options.
   1697    Moreover, change "long" to "ptrdiff_t", as used at the above URL? */
   1698 #undef numberof
   1699 #if 0
   1700 /* The following should work with GCC as documented in its manual,
   1701    but fails: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=38377#c10
   1702    Thus disabled for now. */
   1703 # define numberof(x)                                                    \
   1704   ( __extension__ ({                                                    \
   1705       int is_array = (void *) &(x) == (void *) &(x)[0];                 \
   1706       MPFR_STAT_STATIC_ASSERT (__builtin_constant_p (is_array) ?        \
   1707                                is_array : 1);                           \
   1708       MPFR_ASSERTN (is_array);                                          \
   1709       (long) numberof_const (x);                                        \
   1710     }) )
   1711 #else
   1712 # define numberof(x)                                    \
   1713   (MPFR_ASSERTN ((void *) &(x) == (void *) &(x)[0]),    \
   1714    (long) numberof_const (x))
   1715 #endif
   1716 
   1717 /* Addition with carry (detected by GCC and other good compilers). */
   1718 #define ADD_LIMB(u,v,c) ((u) += (v), (c) = (u) < (v))
   1719 
   1720 /* umul_hi(h, x, y) puts in h the high part of x*y */
   1721 /* MPFR_NEED_LONGLONG_H needs to be defined to use it. */
   1722 #define umul_hi(h, x, y)                        \
   1723   do {                                          \
   1724     mp_limb_t _l;                               \
   1725     umul_ppmm (h, _l, x, y);                    \
   1726     (void) _l;  /* unused variable */           \
   1727   } while (0)
   1728 
   1729 
   1730 /******************************************************
   1731  ************  Save exponent/flags macros  ************
   1732  ******************************************************/
   1733 
   1734 /* See README.dev for details on how to use the macros.
   1735    They are used to set the exponent range to the maximum
   1736    temporarily */
   1737 
   1738 typedef struct {
   1739   mpfr_flags_t saved_flags;
   1740   mpfr_exp_t saved_emin;
   1741   mpfr_exp_t saved_emax;
   1742 } mpfr_save_expo_t;
   1743 
   1744 #define MPFR_SAVE_EXPO_DECL(x) mpfr_save_expo_t x
   1745 #define MPFR_SAVE_EXPO_MARK(x)     \
   1746  ((x).saved_flags = __gmpfr_flags, \
   1747   (x).saved_emin = __gmpfr_emin,   \
   1748   (x).saved_emax = __gmpfr_emax,   \
   1749   __gmpfr_emin = MPFR_EMIN_MIN,    \
   1750   __gmpfr_emax = MPFR_EMAX_MAX)
   1751 #define MPFR_SAVE_EXPO_FREE(x)     \
   1752  (__gmpfr_flags = (x).saved_flags, \
   1753   __gmpfr_emin = (x).saved_emin,   \
   1754   __gmpfr_emax = (x).saved_emax)
   1755 #define MPFR_SAVE_EXPO_UPDATE_FLAGS(x, flags)  \
   1756   (x).saved_flags |= (flags)
   1757 
   1758 /* Speed up final checking */
   1759 #define mpfr_check_range(x,t,r) \
   1760   (MPFR_LIKELY (MPFR_EXP_IN_RANGE (MPFR_EXP (x)))                \
   1761    ? ((t) ? (__gmpfr_flags |= MPFR_FLAGS_INEXACT, (t)) : 0)      \
   1762    : mpfr_check_range(x,t,r))
   1763 
   1764 
   1765 /******************************************************
   1766  *****************  Inline rounding  ******************
   1767  ******************************************************/
   1768 
   1769 /*
   1770  * Note: due to the labels, one cannot use a macro MPFR_RNDRAW* more than
   1771  * once in a function (otherwise these labels would not be unique).
   1772  */
   1773 
   1774 /*
   1775  * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
   1776  * assuming dest's sign is sign.
   1777  * In rounding to nearest mode, execute MIDDLE_HANDLER when the value
   1778  * is the middle of two consecutive numbers in dest precision.
   1779  * Execute OVERFLOW_HANDLER in case of overflow when rounding.
   1780  *
   1781  * Note: the exponent field of dest is not used, possibly except by the
   1782  * handlers. It is the caller (via the handlers) who entirely decides
   1783  * how to handle it.
   1784  */
   1785 #define MPFR_RNDRAW_GEN(inexact, dest, srcp, sprec, rnd, sign,              \
   1786                         MIDDLE_HANDLER, OVERFLOW_HANDLER)                   \
   1787   do {                                                                      \
   1788     mp_size_t _dests, _srcs;                                                \
   1789     mp_limb_t *_destp;                                                      \
   1790     mpfr_prec_t _destprec, _srcprec;                                        \
   1791                                                                             \
   1792     /* Check Trivial Case when Dest Mantissa has more bits than source */   \
   1793     _srcprec = (sprec);                                                     \
   1794     _destprec = MPFR_PREC (dest);                                           \
   1795     MPFR_ASSERTD (_srcprec >= MPFR_PREC_MIN);                               \
   1796     MPFR_ASSERTD (_destprec >= MPFR_PREC_MIN);                              \
   1797     _destp = MPFR_MANT (dest);                                              \
   1798     if (MPFR_UNLIKELY (_destprec >= _srcprec))                              \
   1799       {                                                                     \
   1800         _srcs  = MPFR_PREC2LIMBS (_srcprec);                                \
   1801         _dests = MPFR_PREC2LIMBS (_destprec) - _srcs;                       \
   1802         MPN_COPY (_destp + _dests, srcp, _srcs);                            \
   1803         MPN_ZERO (_destp, _dests);                                          \
   1804         inexact = 0;                                                        \
   1805       }                                                                     \
   1806     else                                                                    \
   1807       {                                                                     \
   1808         /* Non trivial case: rounding needed */                             \
   1809         mpfr_prec_t _sh;                                                    \
   1810         mp_limb_t *_sp;                                                     \
   1811         mp_limb_t _rb, _sb, _ulp;                                           \
   1812                                                                             \
   1813         /* Compute Position and shift */                                    \
   1814         _srcs  = MPFR_PREC2LIMBS (_srcprec);                                \
   1815         _dests = MPFR_PREC2LIMBS (_destprec);                               \
   1816         MPFR_UNSIGNED_MINUS_MODULO (_sh, _destprec);                        \
   1817         _sp = (srcp) + _srcs - _dests;                                      \
   1818                                                                             \
   1819         /* General case when prec % GMP_NUMB_BITS != 0 */                   \
   1820         if (MPFR_LIKELY (_sh != 0))                                         \
   1821           {                                                                 \
   1822             mp_limb_t _mask;                                                \
   1823             /* Compute Rounding Bit and Sticky Bit */                       \
   1824             /* Note: in directed rounding modes, if the rounding bit */     \
   1825             /* is 1, the behavior does not depend on the sticky bit; */     \
   1826             /* thus we will not try to compute it in this case (this */     \
   1827             /* can be much faster and avoids to read uninitialized   */     \
   1828             /* data in the current mpfr_mul implementation). We just */     \
   1829             /* make sure that _sb is initialized.                    */     \
   1830             _mask = MPFR_LIMB_ONE << (_sh - 1);                             \
   1831             _rb = _sp[0] & _mask;                                           \
   1832             _sb = _sp[0] & (_mask - 1);                                     \
   1833             if ((rnd) == MPFR_RNDN || _rb == 0)                             \
   1834               {                                                             \
   1835                 mp_limb_t *_tmp;                                            \
   1836                 mp_size_t _n;                                               \
   1837                 for (_tmp = _sp, _n = _srcs - _dests ;                      \
   1838                      _n != 0 && _sb == 0 ; _n--)                            \
   1839                   _sb = *--_tmp;                                            \
   1840               }                                                             \
   1841             _ulp = 2 * _mask;                                               \
   1842           }                                                                 \
   1843         else /* _sh == 0 */                                                 \
   1844           {                                                                 \
   1845             MPFR_ASSERTD (_dests < _srcs);                                  \
   1846             /* Compute Rounding Bit and Sticky Bit - see note above */      \
   1847             _rb = _sp[-1] & MPFR_LIMB_HIGHBIT;                              \
   1848             _sb = _sp[-1] & (MPFR_LIMB_HIGHBIT-1);                          \
   1849             if ((rnd) == MPFR_RNDN || _rb == 0)                             \
   1850               {                                                             \
   1851                 mp_limb_t *_tmp;                                            \
   1852                 mp_size_t _n;                                               \
   1853                 for (_tmp = _sp - 1, _n = _srcs - _dests - 1 ;              \
   1854                      _n != 0 && _sb == 0 ; _n--)                            \
   1855                   _sb = *--_tmp;                                            \
   1856               }                                                             \
   1857             _ulp = MPFR_LIMB_ONE;                                           \
   1858           }                                                                 \
   1859         /* Rounding */                                                      \
   1860         if (rnd == MPFR_RNDF)                                               \
   1861           {                                                                 \
   1862             inexact = 0;                                                    \
   1863             goto trunc_doit;                                                \
   1864           }                                                                 \
   1865         else if (rnd == MPFR_RNDN)                                          \
   1866           {                                                                 \
   1867             if (_rb == 0)                                                   \
   1868               {                                                             \
   1869               trunc:                                                        \
   1870                 inexact = MPFR_LIKELY ((_sb | _rb) != 0) ? -sign : 0;       \
   1871               trunc_doit:                                                   \
   1872                 MPN_COPY (_destp, _sp, _dests);                             \
   1873                 _destp[0] &= ~(_ulp - 1);                                   \
   1874               }                                                             \
   1875             else if (MPFR_UNLIKELY (_sb == 0))                              \
   1876               { /* Middle of two consecutive representable numbers */       \
   1877                 MIDDLE_HANDLER;                                             \
   1878               }                                                             \
   1879             else                                                            \
   1880               {                                                             \
   1881                 if (0)                                                      \
   1882                   goto addoneulp_doit; /* dummy code to avoid warning */    \
   1883               addoneulp:                                                    \
   1884                 inexact = sign;                                             \
   1885               addoneulp_doit:                                               \
   1886                 if (MPFR_UNLIKELY (mpn_add_1 (_destp, _sp, _dests, _ulp)))  \
   1887                   {                                                         \
   1888                     _destp[_dests - 1] = MPFR_LIMB_HIGHBIT;                 \
   1889                     OVERFLOW_HANDLER;                                       \
   1890                   }                                                         \
   1891                 _destp[0] &= ~(_ulp - 1);                                   \
   1892               }                                                             \
   1893           }                                                                 \
   1894         else                                                                \
   1895           { /* Directed rounding mode */                                    \
   1896             if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))           \
   1897               goto trunc;                                                   \
   1898              else if (MPFR_UNLIKELY ((_sb | _rb) == 0))                     \
   1899                {                                                            \
   1900                  inexact = 0;                                               \
   1901                  goto trunc_doit;                                           \
   1902                }                                                            \
   1903              else                                                           \
   1904               goto addoneulp;                                               \
   1905           }                                                                 \
   1906       }                                                                     \
   1907   } while (0)
   1908 
   1909 /*
   1910  * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
   1911  * assuming dest's sign is sign.
   1912  * Execute OVERFLOW_HANDLER in case of overflow when rounding.
   1913  */
   1914 #define MPFR_RNDRAW(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER) \
   1915    MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign,                   \
   1916      if ((_sp[0] & _ulp) == 0)                                               \
   1917        {                                                                     \
   1918          inexact = -sign;                                                    \
   1919          goto trunc_doit;                                                    \
   1920        }                                                                     \
   1921      else                                                                    \
   1922        goto addoneulp;                                                       \
   1923      , OVERFLOW_HANDLER)
   1924 
   1925 /*
   1926  * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
   1927  * assuming dest's sign is sign.
   1928  * Execute OVERFLOW_HANDLER in case of overflow when rounding.
   1929  * Set inexact to +/- MPFR_EVEN_INEX in case of even rounding.
   1930  */
   1931 #define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, \
   1932                          OVERFLOW_HANDLER)                      \
   1933    MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign,      \
   1934      if ((_sp[0] & _ulp) == 0)                                  \
   1935        {                                                        \
   1936          inexact = -MPFR_EVEN_INEX * sign;                      \
   1937          goto trunc_doit;                                       \
   1938        }                                                        \
   1939      else                                                       \
   1940        {                                                        \
   1941          inexact = MPFR_EVEN_INEX * sign;                       \
   1942          goto addoneulp_doit;                                   \
   1943        }                                                        \
   1944      , OVERFLOW_HANDLER)
   1945 
   1946 /* Return TRUE if b is non singular and we can round it to precision 'prec'
   1947    and determine the ternary value, with rounding mode 'rnd', and with
   1948    error at most 'error' */
   1949 #define MPFR_CAN_ROUND(b,err,prec,rnd)                                       \
   1950  (!MPFR_IS_SINGULAR (b) && mpfr_round_p (MPFR_MANT (b), MPFR_LIMB_SIZE (b),  \
   1951                                          (err), (prec) + ((rnd)==MPFR_RNDN)))
   1952 
   1953 /* Copy the sign and the significand, and handle the exponent in exp. */
   1954 #define MPFR_SETRAW(inexact,dest,src,exp,rnd)                           \
   1955   if (dest != src)                                                      \
   1956     {                                                                   \
   1957       MPFR_SET_SIGN (dest, MPFR_SIGN (src));                            \
   1958       if (MPFR_PREC (dest) == MPFR_PREC (src))                          \
   1959         {                                                               \
   1960           MPN_COPY (MPFR_MANT (dest), MPFR_MANT (src),                  \
   1961                     MPFR_LIMB_SIZE (src));                              \
   1962           inexact = 0;                                                  \
   1963         }                                                               \
   1964       else                                                              \
   1965         {                                                               \
   1966           MPFR_RNDRAW (inexact, dest, MPFR_MANT (src), MPFR_PREC (src), \
   1967                        rnd, MPFR_SIGN (src), exp++);                    \
   1968         }                                                               \
   1969     }                                                                   \
   1970   else                                                                  \
   1971     inexact = 0;
   1972 
   1973 /* TODO: fix this description (see round_near_x.c). */
   1974 /* Assuming that the function has a Taylor expansion which looks like:
   1975     y=o(f(x)) = o(v + g(x)) with |g(x)| <= 2^(EXP(v)-err)
   1976    we can quickly set y to v if x is small (ie err > prec(y)+1) in most
   1977    cases. It assumes that f(x) is not representable exactly as a FP number.
   1978    v must not be a singular value (NAN, INF or ZERO); usual values are
   1979    v=1 or v=x.
   1980 
   1981    y is the destination (a mpfr_t), v the value to set (a mpfr_t),
   1982    err1+err2 with 0 <= err2 <= 3 the error term (mpfr_exp_t's), dir (an int)
   1983    is the direction of the committed error (if dir = 0, it rounds toward 0,
   1984    if dir=1, it rounds away from 0), rnd the rounding mode.
   1985 
   1986    It returns from the function a ternary value in case of success.
   1987    If you want to free something, you must fill the "extra" field
   1988    in consequences, otherwise put nothing in it.
   1989 
   1990    The test is less restrictive than necessary, but the function
   1991    will finish the check itself.
   1992 
   1993    Note: err1 + err2 is allowed to overflow as mpfr_exp_t, but it must give
   1994    its real value as mpfr_uexp_t.
   1995 */
   1996 #define MPFR_FAST_COMPUTE_IF_SMALL_INPUT(y,v,err1,err2,dir,rnd,extra)   \
   1997   do {                                                                  \
   1998     mpfr_ptr _y = (y);                                                  \
   1999     mpfr_exp_t _err1 = (err1);                                          \
   2000     mpfr_exp_t _err2 = (err2);                                          \
   2001     if (_err1 > 0)                                                      \
   2002       {                                                                 \
   2003         mpfr_uexp_t _err = (mpfr_uexp_t) _err1 + _err2;                 \
   2004         if (MPFR_UNLIKELY (_err > MPFR_PREC (_y) + 1))                  \
   2005           {                                                             \
   2006             int _inexact = mpfr_round_near_x (_y,(v),_err,(dir),(rnd)); \
   2007             if (_inexact != 0)                                          \
   2008               {                                                         \
   2009                 extra;                                                  \
   2010                 return _inexact;                                        \
   2011               }                                                         \
   2012           }                                                             \
   2013       }                                                                 \
   2014   } while (0)
   2015 
   2016 /* Variant, to be called somewhere after MPFR_SAVE_EXPO_MARK. This variant
   2017    is needed when there are some computations before or when some non-zero
   2018    real constant is used, such as __gmpfr_one for mpfr_cos. */
   2019 #define MPFR_SMALL_INPUT_AFTER_SAVE_EXPO(y,v,err1,err2,dir,rnd,expo,extra) \
   2020   do {                                                                  \
   2021     mpfr_ptr _y = (y);                                                  \
   2022     mpfr_exp_t _err1 = (err1);                                          \
   2023     mpfr_exp_t _err2 = (err2);                                          \
   2024     if (_err1 > 0)                                                      \
   2025       {                                                                 \
   2026         mpfr_uexp_t _err = (mpfr_uexp_t) _err1 + _err2;                 \
   2027         if (MPFR_UNLIKELY (_err > MPFR_PREC (_y) + 1))                  \
   2028           {                                                             \
   2029             int _inexact;                                               \
   2030             MPFR_CLEAR_FLAGS ();                                        \
   2031             _inexact = mpfr_round_near_x (_y,(v),_err,(dir),(rnd));     \
   2032             if (_inexact != 0)                                          \
   2033               {                                                         \
   2034                 extra;                                                  \
   2035                 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);      \
   2036                 MPFR_SAVE_EXPO_FREE (expo);                             \
   2037                 return mpfr_check_range (_y, _inexact, (rnd));          \
   2038               }                                                         \
   2039           }                                                             \
   2040       }                                                                 \
   2041   } while (0)
   2042 
   2043 
   2044 /******************************************************
   2045  *****************  Ziv loop macros  ******************
   2046  ******************************************************/
   2047 
   2048 /* To safely increase some precision, detecting integer overflows.
   2049    This macro is particularly useful when determining the initial
   2050    working precision before Ziv's loop. P is a precision, X is an
   2051    arbitrary non-negative integer.
   2052    Note: On 2012-02-23, the MPFR_PREC_MAX value has been decreased
   2053    by 256 from the maximum value representable in the mpfr_prec_t
   2054    type, in order to avoid some integer overflows when this macro
   2055    is not used (if the result is larger than MPFR_PREC_MAX, this
   2056    should be detected with a later assertion, e.g. in mpfr_init2).
   2057    But this change is mainly for existing code that has not been
   2058    updated yet. So, it is advised to always use MPFR_ADD_PREC or
   2059    MPFR_INC_PREC if the result can be larger than MPFR_PREC_MAX. */
   2060 #define MPFR_ADD_PREC(P,X) \
   2061   (MPFR_ASSERTN ((X) <= MPFR_PREC_MAX - (P)), (P) + (X))
   2062 #define MPFR_INC_PREC(P,X) \
   2063   (MPFR_ASSERTN ((X) <= MPFR_PREC_MAX - (P)), (P) += (X))
   2064 
   2065 #ifndef MPFR_USE_LOGGING
   2066 
   2067 #define MPFR_ZIV_DECL(_x) mpfr_prec_t _x
   2068 #define MPFR_ZIV_INIT(_x, _p) (_x) = GMP_NUMB_BITS
   2069 #define MPFR_ZIV_NEXT(_x, _p) (MPFR_INC_PREC (_p, _x), (_x) = (_p)/2)
   2070 #define MPFR_ZIV_FREE(x)
   2071 
   2072 #else
   2073 
   2074 /* The following test on glibc is there mainly for Darwin (Mac OS X), to
   2075    obtain a better error message. The real test should have been a test
   2076    concerning nested functions in gcc, which are disabled by default on
   2077    Darwin; but it is not possible to do that without a configure test. */
   2078 # if defined (__cplusplus) || !(__MPFR_GNUC(3,0) && __MPFR_GLIBC(2,0))
   2079 #  error "Logging not supported (needs gcc >= 3.0 and GNU C Library >= 2.0)."
   2080 # endif
   2081 
   2082 /* Use LOGGING */
   2083 
   2084 /* Note: the mpfr_log_level >= 0 below avoids to take into account
   2085    Ziv loops used by the MPFR functions called by the mpfr_fprintf
   2086    in LOG_PRINT. */
   2087 
   2088 #define MPFR_ZIV_DECL(_x)                                               \
   2089   mpfr_prec_t _x;                                                       \
   2090   int _x ## _cpt = 1;                                                   \
   2091   static unsigned long  _x ## _loop = 0, _x ## _bad = 0;                \
   2092   static const char *_x ## _fname = __func__;                           \
   2093   auto void __attribute__ ((destructor)) x ## _f  (void);               \
   2094   void __attribute__ ((destructor)) x ## _f  (void) {                   \
   2095     if (_x ## _loop != 0 && (MPFR_LOG_STAT_F & mpfr_log_type)) {        \
   2096       fprintf (mpfr_log_file,                                           \
   2097                "%s: Ziv failed %2.2f%% (%lu bad cases / %lu calls)\n",  \
   2098                _x ## _fname, (double) 100.0 * _x ## _bad / _x ## _loop, \
   2099                _x ## _bad, _x ## _loop );                               \
   2100       if (mpfr_log_flush)                                               \
   2101         fflush (mpfr_log_file);                                         \
   2102     }                                                                   \
   2103   }
   2104 
   2105 #define MPFR_ZIV_INIT(_x, _p)                                           \
   2106   do                                                                    \
   2107     {                                                                   \
   2108       (_x) = GMP_NUMB_BITS;                                             \
   2109       if (mpfr_log_level >= 0)                                          \
   2110         _x ## _loop ++;                                                 \
   2111       LOG_PRINT (MPFR_LOG_BADCASE_F, "%s:ZIV 1st prec=%Pd\n",           \
   2112                  __func__, (mpfr_prec_t) (_p));                         \
   2113     }                                                                   \
   2114   while (0)
   2115 
   2116 #define MPFR_ZIV_NEXT(_x, _p)                                           \
   2117   do                                                                    \
   2118     {                                                                   \
   2119       MPFR_INC_PREC (_p, _x);                                           \
   2120       (_x) = (_p) / 2;                                                  \
   2121       if (mpfr_log_level >= 0)                                          \
   2122         _x ## _bad += (_x ## _cpt == 1);                                \
   2123       _x ## _cpt ++;                                                    \
   2124       LOG_PRINT (MPFR_LOG_BADCASE_F, "%s:ZIV new prec=%Pd\n",           \
   2125                  __func__, (mpfr_prec_t) (_p));                         \
   2126     }                                                                   \
   2127   while (0)
   2128 
   2129 #define MPFR_ZIV_FREE(_x)                                               \
   2130   do                                                                    \
   2131     if (_x ## _cpt > 1)                                                 \
   2132       LOG_PRINT (MPFR_LOG_BADCASE_F, "%s:ZIV %d loops\n",               \
   2133                  __func__, _x ## _cpt);                                 \
   2134   while (0)
   2135 
   2136 #endif
   2137 
   2138 
   2139 /******************************************************
   2140  ******************  Logging macros  ******************
   2141  ******************************************************/
   2142 
   2143 /* The different kind of LOG */
   2144 #define MPFR_LOG_INPUT_F    1
   2145 #define MPFR_LOG_OUTPUT_F   2
   2146 #define MPFR_LOG_INTERNAL_F 4
   2147 #define MPFR_LOG_TIME_F     8
   2148 #define MPFR_LOG_BADCASE_F  16
   2149 #define MPFR_LOG_MSG_F      32
   2150 #define MPFR_LOG_STAT_F     64
   2151 
   2152 #ifdef MPFR_USE_LOGGING
   2153 
   2154 /* Check if we can support this feature */
   2155 # ifdef MPFR_USE_THREAD_SAFE
   2156 #  error "Enable either `Logging' or `thread-safe', not both"
   2157 # endif
   2158 # if !__MPFR_GNUC(3,0)
   2159 #  error "Logging not supported (GCC >= 3.0)"
   2160 # endif
   2161 
   2162 #if defined (__cplusplus)
   2163 extern "C" {
   2164 #endif
   2165 
   2166 __MPFR_DECLSPEC extern FILE *mpfr_log_file;
   2167 __MPFR_DECLSPEC extern int   mpfr_log_flush;
   2168 __MPFR_DECLSPEC extern int   mpfr_log_type;
   2169 __MPFR_DECLSPEC extern int   mpfr_log_level;
   2170 __MPFR_DECLSPEC extern int   mpfr_log_current;
   2171 __MPFR_DECLSPEC extern mpfr_prec_t mpfr_log_prec;
   2172 
   2173 #if defined (__cplusplus)
   2174  }
   2175 #endif
   2176 
   2177 /* LOG_PRINT calls mpfr_fprintf on mpfr_log_file with logging disabled
   2178    (recursive logging is not wanted and freezes MPFR). */
   2179 #define LOG_PRINT(type, format, ...)                                    \
   2180   do                                                                    \
   2181     if ((mpfr_log_type & (type)) && mpfr_log_current <= mpfr_log_level) \
   2182       {                                                                 \
   2183         int old_level = mpfr_log_level;                                 \
   2184         mpfr_log_level = -1;  /* disable logging in mpfr_fprintf */     \
   2185         __gmpfr_cache_const_pi = __gmpfr_logging_pi;                    \
   2186         __gmpfr_cache_const_log2 = __gmpfr_logging_log2;                \
   2187         mpfr_fprintf (mpfr_log_file, format, __VA_ARGS__);              \
   2188         if (mpfr_log_flush)                                             \
   2189           fflush (mpfr_log_file);                                       \
   2190         mpfr_log_level = old_level;                                     \
   2191         __gmpfr_cache_const_pi = __gmpfr_normal_pi;                     \
   2192         __gmpfr_cache_const_log2 = __gmpfr_normal_log2;                 \
   2193       }                                                                 \
   2194   while (0)
   2195 
   2196 #define MPFR_LOG_VAR(x)                                                 \
   2197   LOG_PRINT (MPFR_LOG_INTERNAL_F, "%s.%d:%s[%#Pd]=%.*Rg\n", __func__,   \
   2198              (int) __LINE__, #x, mpfr_get_prec (x), mpfr_log_prec, x)
   2199 
   2200 #define MPFR_LOG_MSG2(format, ...)                                      \
   2201   LOG_PRINT (MPFR_LOG_MSG_F, "%s.%d: "format, __func__, (int) __LINE__, \
   2202              __VA_ARGS__)
   2203 #define MPFR_LOG_MSG(x) MPFR_LOG_MSG2 x
   2204 
   2205 #define MPFR_LOG_BEGIN2(format, ...)                                    \
   2206   mpfr_log_current ++;                                                  \
   2207   LOG_PRINT (MPFR_LOG_INPUT_F, "%s:IN  flags=%x "format"\n", __func__,  \
   2208              (unsigned int) __gmpfr_flags, __VA_ARGS__);                \
   2209   if ((MPFR_LOG_TIME_F & mpfr_log_type) &&                              \
   2210       (mpfr_log_current <= mpfr_log_level))                             \
   2211     __gmpfr_log_time = mpfr_get_cputime ();
   2212 #define MPFR_LOG_BEGIN(x)                                               \
   2213   int __gmpfr_log_time = 0;                                             \
   2214   MPFR_LOG_BEGIN2 x
   2215 
   2216 #define MPFR_LOG_END2(format, ...)                                      \
   2217   LOG_PRINT (MPFR_LOG_TIME_F, "%s:TIM %dms\n", __mpfr_log_fname,        \
   2218              mpfr_get_cputime () - __gmpfr_log_time);                   \
   2219   LOG_PRINT (MPFR_LOG_OUTPUT_F, "%s:OUT flags=%x "format"\n",           \
   2220              __mpfr_log_fname, (unsigned int) __gmpfr_flags,            \
   2221              __VA_ARGS__);                                              \
   2222   mpfr_log_current --;
   2223 #define MPFR_LOG_END(x)                                                 \
   2224   static const char *__mpfr_log_fname = __func__;                       \
   2225   MPFR_LOG_END2 x
   2226 
   2227 #define MPFR_LOG_FUNC(begin,end)                                        \
   2228   static const char *__mpfr_log_fname = __func__;                       \
   2229   auto void __mpfr_log_cleanup (int *time);                             \
   2230   void __mpfr_log_cleanup (int *time) {                                 \
   2231     int __gmpfr_log_time = *time;                                       \
   2232     MPFR_LOG_END2 end; }                                                \
   2233   int __gmpfr_log_time __attribute__ ((cleanup (__mpfr_log_cleanup)));  \
   2234   __gmpfr_log_time = 0;                                                 \
   2235   MPFR_LOG_BEGIN2 begin
   2236 
   2237 #else /* MPFR_USE_LOGGING */
   2238 
   2239 /* Define void macro for logging */
   2240 
   2241 #define MPFR_LOG_VAR(x)
   2242 #define MPFR_LOG_BEGIN(x)
   2243 #define MPFR_LOG_END(x)
   2244 #define MPFR_LOG_MSG(x)
   2245 #define MPFR_LOG_FUNC(x,y)
   2246 
   2247 #endif /* MPFR_USE_LOGGING */
   2248 
   2249 
   2250 /**************************************************************
   2251  ************  Group Initialize Functions Macros  *************
   2252  **************************************************************/
   2253 
   2254 #ifndef MPFR_GROUP_STATIC_SIZE
   2255 # define MPFR_GROUP_STATIC_SIZE 16
   2256 #endif
   2257 
   2258 struct mpfr_group_t {
   2259   size_t     alloc;
   2260   mp_limb_t *mant;
   2261 #if MPFR_GROUP_STATIC_SIZE != 0
   2262   mp_limb_t  tab[MPFR_GROUP_STATIC_SIZE];
   2263 #else
   2264   /* In order to detect memory leaks when testing, MPFR_GROUP_STATIC_SIZE
   2265      can be set to 0, in which case tab will not be used. ISO C does not
   2266      support zero-length arrays[*], thus let's use a flexible array member
   2267      (which will be equivalent here). Note: this is new in C99, but this
   2268      is just used for testing.
   2269      [*] Zero-length arrays are a GNU extension:
   2270            https://gcc.gnu.org/onlinedocs/gcc/Zero-Length.html
   2271          and as such an extension is forbidden in ISO C, it triggers an
   2272          error with -Werror=pedantic.
   2273   */
   2274   mp_limb_t  tab[];
   2275 #endif
   2276 };
   2277 
   2278 #define MPFR_GROUP_DECL(g) struct mpfr_group_t g
   2279 #define MPFR_GROUP_CLEAR(g) do {                                 \
   2280  MPFR_LOG_MSG (("GROUP_CLEAR: ptr = 0x%lX, size = %lu\n",        \
   2281                 (unsigned long) (g).mant,                        \
   2282                 (unsigned long) (g).alloc));                     \
   2283  if ((g).alloc != 0) {                                           \
   2284    MPFR_ASSERTD ((g).mant != (g).tab);                           \
   2285    mpfr_free_func ((g).mant, (g).alloc);                         \
   2286  }} while (0)
   2287 
   2288 #define MPFR_GROUP_INIT_TEMPLATE(g, prec, num, handler) do {            \
   2289  mpfr_prec_t _prec = (prec);                                            \
   2290  mp_size_t _size;                                                       \
   2291  MPFR_ASSERTD (_prec >= MPFR_PREC_MIN);                                 \
   2292  if (MPFR_UNLIKELY (_prec > MPFR_PREC_MAX))                             \
   2293    mpfr_abort_prec_max ();                                              \
   2294  _size = MPFR_PREC2LIMBS (_prec);                                       \
   2295  if (_size * (num) > MPFR_GROUP_STATIC_SIZE)                            \
   2296    {                                                                    \
   2297      (g).alloc = (num) * _size * sizeof (mp_limb_t);                    \
   2298      (g).mant = (mp_limb_t *) mpfr_allocate_func ((g).alloc);           \
   2299    }                                                                    \
   2300  else                                                                   \
   2301    {                                                                    \
   2302      (g).alloc = 0;                                                     \
   2303      (g).mant = (g).tab;                                                \
   2304    }                                                                    \
   2305  MPFR_LOG_MSG (("GROUP_INIT: ptr = 0x%lX, size = %lu\n",                \
   2306                 (unsigned long) (g).mant, (unsigned long) (g).alloc));  \
   2307  handler;                                                               \
   2308  } while (0)
   2309 #define MPFR_GROUP_TINIT(g, n, x)                       \
   2310   MPFR_TMP_INIT1 ((g).mant + _size * (n), x, _prec)
   2311 
   2312 #define MPFR_GROUP_INIT_1(g, prec, x)                            \
   2313  MPFR_GROUP_INIT_TEMPLATE(g, prec, 1, MPFR_GROUP_TINIT(g, 0, x))
   2314 #define MPFR_GROUP_INIT_2(g, prec, x, y)                         \
   2315  MPFR_GROUP_INIT_TEMPLATE(g, prec, 2,                            \
   2316    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y))
   2317 #define MPFR_GROUP_INIT_3(g, prec, x, y, z)                      \
   2318  MPFR_GROUP_INIT_TEMPLATE(g, prec, 3,                            \
   2319    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
   2320    MPFR_GROUP_TINIT(g, 2, z))
   2321 #define MPFR_GROUP_INIT_4(g, prec, x, y, z, t)                   \
   2322  MPFR_GROUP_INIT_TEMPLATE(g, prec, 4,                            \
   2323    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
   2324    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t))
   2325 #define MPFR_GROUP_INIT_5(g, prec, x, y, z, t, a)                \
   2326  MPFR_GROUP_INIT_TEMPLATE(g, prec, 5,                            \
   2327    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
   2328    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t);          \
   2329    MPFR_GROUP_TINIT(g, 4, a))
   2330 #define MPFR_GROUP_INIT_6(g, prec, x, y, z, t, a, b)             \
   2331  MPFR_GROUP_INIT_TEMPLATE(g, prec, 6,                            \
   2332    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
   2333    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t);          \
   2334    MPFR_GROUP_TINIT(g, 4, a);MPFR_GROUP_TINIT(g, 5, b))
   2335 
   2336 #define MPFR_GROUP_REPREC_TEMPLATE(g, prec, num, handler) do {          \
   2337  mpfr_prec_t _prec = (prec);                                            \
   2338  size_t    _oalloc = (g).alloc;                                         \
   2339  mp_size_t _size;                                                       \
   2340  MPFR_LOG_MSG (("GROUP_REPREC: oldptr = 0x%lX, oldsize = %lu\n",        \
   2341                 (unsigned long) (g).mant, (unsigned long) _oalloc));    \
   2342  MPFR_ASSERTD (_prec >= MPFR_PREC_MIN);                                 \
   2343  if (MPFR_UNLIKELY (_prec > MPFR_PREC_MAX))                             \
   2344    mpfr_abort_prec_max ();                                              \
   2345  _size = MPFR_PREC2LIMBS (_prec);                                       \
   2346  (g).alloc = (num) * _size * sizeof (mp_limb_t);                        \
   2347  if (_oalloc == 0)                                                      \
   2348    (g).mant = (mp_limb_t *) mpfr_allocate_func ((g).alloc);             \
   2349  else                                                                   \
   2350    (g).mant = (mp_limb_t *)                                             \
   2351      mpfr_reallocate_func ((g).mant, _oalloc, (g).alloc);               \
   2352  MPFR_LOG_MSG (("GROUP_REPREC: newptr = 0x%lX, newsize = %lu\n",        \
   2353                 (unsigned long) (g).mant, (unsigned long) (g).alloc));  \
   2354  handler;                                                               \
   2355  } while (0)
   2356 
   2357 #define MPFR_GROUP_REPREC_1(g, prec, x)                          \
   2358  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 1, MPFR_GROUP_TINIT(g, 0, x))
   2359 #define MPFR_GROUP_REPREC_2(g, prec, x, y)                       \
   2360  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 2,                          \
   2361    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y))
   2362 #define MPFR_GROUP_REPREC_3(g, prec, x, y, z)                    \
   2363  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 3,                          \
   2364    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
   2365    MPFR_GROUP_TINIT(g, 2, z))
   2366 #define MPFR_GROUP_REPREC_4(g, prec, x, y, z, t)                 \
   2367  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 4,                          \
   2368    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
   2369    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t))
   2370 #define MPFR_GROUP_REPREC_5(g, prec, x, y, z, t, a)              \
   2371  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 5,                          \
   2372    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
   2373    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t);          \
   2374    MPFR_GROUP_TINIT(g, 4, a))
   2375 #define MPFR_GROUP_REPREC_6(g, prec, x, y, z, t, a, b)           \
   2376  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 6,                          \
   2377    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
   2378    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t);          \
   2379    MPFR_GROUP_TINIT(g, 4, a);MPFR_GROUP_TINIT(g, 5, b))
   2380 
   2381 
   2382 /******************************************************
   2383  ***************  Internal functions  *****************
   2384  ******************************************************/
   2385 
   2386 #if defined (__cplusplus)
   2387 extern "C" {
   2388 #endif
   2389 
   2390 MPFR_COLD_FUNCTION_ATTR __MPFR_DECLSPEC int
   2391   mpfr_underflow (mpfr_ptr, mpfr_rnd_t, int);
   2392 MPFR_COLD_FUNCTION_ATTR __MPFR_DECLSPEC int
   2393   mpfr_overflow (mpfr_ptr, mpfr_rnd_t, int);
   2394 
   2395 __MPFR_DECLSPEC int mpfr_add1 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
   2396 __MPFR_DECLSPEC int mpfr_sub1 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
   2397 __MPFR_DECLSPEC int mpfr_add1sp (mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
   2398                                  mpfr_rnd_t);
   2399 __MPFR_DECLSPEC int mpfr_sub1sp (mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
   2400                                  mpfr_rnd_t);
   2401 __MPFR_DECLSPEC int mpfr_can_round_raw (const mp_limb_t *,
   2402              mp_size_t, int, mpfr_exp_t, mpfr_rnd_t, mpfr_rnd_t, mpfr_prec_t);
   2403 
   2404 __MPFR_DECLSPEC int mpfr_set_1_2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t, int);
   2405 
   2406 __MPFR_DECLSPEC int mpfr_cmp2 (mpfr_srcptr, mpfr_srcptr, mpfr_prec_t *);
   2407 
   2408 __MPFR_DECLSPEC long          __gmpfr_ceil_log2     (double);
   2409 __MPFR_DECLSPEC long          __gmpfr_floor_log2    (double);
   2410 __MPFR_DECLSPEC double        __gmpfr_ceil_exp2     (double);
   2411 __MPFR_DECLSPEC unsigned long __gmpfr_isqrt     (unsigned long);
   2412 __MPFR_DECLSPEC unsigned long __gmpfr_cuberoot  (unsigned long);
   2413 __MPFR_DECLSPEC int       __gmpfr_int_ceil_log2 (unsigned long);
   2414 
   2415 __MPFR_DECLSPEC mpfr_exp_t mpfr_ceil_mul (mpfr_exp_t, int, int);
   2416 
   2417 __MPFR_DECLSPEC int mpfr_exp_2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
   2418 __MPFR_DECLSPEC int mpfr_exp_3 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
   2419 __MPFR_DECLSPEC int mpfr_powerof2_raw (mpfr_srcptr);
   2420 __MPFR_DECLSPEC int mpfr_powerof2_raw2 (const mp_limb_t *, mp_size_t);
   2421 
   2422 __MPFR_DECLSPEC int mpfr_pow_general (mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
   2423                                       mpfr_rnd_t, int, mpfr_save_expo_t *);
   2424 
   2425 __MPFR_DECLSPEC void mpfr_setmax (mpfr_ptr, mpfr_exp_t);
   2426 __MPFR_DECLSPEC void mpfr_setmin (mpfr_ptr, mpfr_exp_t);
   2427 
   2428 __MPFR_DECLSPEC int mpfr_mpn_exp (mp_limb_t *, mpfr_exp_t *, int,
   2429                                   mpfr_exp_t, size_t);
   2430 
   2431 #ifdef _MPFR_H_HAVE_FILE
   2432 __MPFR_DECLSPEC void mpfr_fdump (FILE *, mpfr_srcptr);
   2433 #endif
   2434 __MPFR_DECLSPEC void mpfr_print_mant_binary (const char*, const mp_limb_t*,
   2435                                              mpfr_prec_t);
   2436 __MPFR_DECLSPEC void mpfr_set_str_binary (mpfr_ptr, const char*);
   2437 
   2438 __MPFR_DECLSPEC int mpfr_round_raw (mp_limb_t *,
   2439        const mp_limb_t *, mpfr_prec_t, int, mpfr_prec_t, mpfr_rnd_t, int *);
   2440 __MPFR_DECLSPEC int mpfr_round_raw_2 (const mp_limb_t *, mpfr_prec_t, int,
   2441                                       mpfr_prec_t, mpfr_rnd_t);
   2442 /* No longer defined (see round_prec.c).
   2443    Uncomment if it needs to be defined again.
   2444 __MPFR_DECLSPEC int mpfr_round_raw_3 (const mp_limb_t *,
   2445              mpfr_prec_t, int, mpfr_prec_t, mpfr_rnd_t, int *);
   2446 */
   2447 __MPFR_DECLSPEC int mpfr_round_raw_4 (mp_limb_t *,
   2448        const mp_limb_t *, mpfr_prec_t, int, mpfr_prec_t, mpfr_rnd_t);
   2449 
   2450 #define mpfr_round_raw2(xp, xn, neg, r, prec) \
   2451   mpfr_round_raw_2((xp),(xn)*GMP_NUMB_BITS,(neg),(prec),(r))
   2452 
   2453 __MPFR_DECLSPEC int mpfr_check (mpfr_srcptr);
   2454 
   2455 __MPFR_DECLSPEC int mpfr_get_cputime (void);
   2456 
   2457 __MPFR_DECLSPEC void mpfr_nexttozero (mpfr_ptr);
   2458 __MPFR_DECLSPEC void mpfr_nexttoinf (mpfr_ptr);
   2459 
   2460 __MPFR_DECLSPEC int mpfr_const_pi_internal (mpfr_ptr,mpfr_rnd_t);
   2461 __MPFR_DECLSPEC int mpfr_const_log2_internal (mpfr_ptr,mpfr_rnd_t);
   2462 __MPFR_DECLSPEC int mpfr_const_euler_internal (mpfr_ptr, mpfr_rnd_t);
   2463 __MPFR_DECLSPEC int mpfr_const_catalan_internal (mpfr_ptr, mpfr_rnd_t);
   2464 
   2465 #if 0
   2466 __MPFR_DECLSPEC void mpfr_init_cache (mpfr_cache_t,
   2467                                       int(*)(mpfr_ptr,mpfr_rnd_t));
   2468 #endif
   2469 __MPFR_DECLSPEC void mpfr_clear_cache (mpfr_cache_t);
   2470 __MPFR_DECLSPEC int  mpfr_cache (mpfr_ptr, mpfr_cache_t, mpfr_rnd_t);
   2471 
   2472 __MPFR_DECLSPEC void mpfr_mulhigh_n (mpfr_limb_ptr, mpfr_limb_srcptr,
   2473                                      mpfr_limb_srcptr, mp_size_t);
   2474 __MPFR_DECLSPEC void mpfr_mullow_n  (mpfr_limb_ptr, mpfr_limb_srcptr,
   2475                                      mpfr_limb_srcptr, mp_size_t);
   2476 __MPFR_DECLSPEC void mpfr_sqrhigh_n (mpfr_limb_ptr, mpfr_limb_srcptr,
   2477                                      mp_size_t);
   2478 __MPFR_DECLSPEC mp_limb_t mpfr_divhigh_n (mpfr_limb_ptr, mpfr_limb_ptr,
   2479                                           mpfr_limb_ptr, mp_size_t);
   2480 
   2481 __MPFR_DECLSPEC int mpfr_round_p (mp_limb_t *, mp_size_t, mpfr_exp_t,
   2482                                   mpfr_prec_t);
   2483 
   2484 __MPFR_DECLSPEC int mpfr_round_near_x (mpfr_ptr, mpfr_srcptr, mpfr_uexp_t, int,
   2485                                        mpfr_rnd_t);
   2486 __MPFR_DECLSPEC MPFR_COLD_FUNCTION_ATTR MPFR_NORETURN void
   2487   mpfr_abort_prec_max (void);
   2488 
   2489 __MPFR_DECLSPEC void mpfr_rand_raw (mpfr_limb_ptr, gmp_randstate_t,
   2490                                     mpfr_prec_t);
   2491 
   2492 __MPFR_DECLSPEC mpz_srcptr mpfr_bernoulli_cache (unsigned long);
   2493 __MPFR_DECLSPEC void mpfr_bernoulli_freecache (void);
   2494 
   2495 __MPFR_DECLSPEC int mpfr_sincos_fast (mpfr_ptr, mpfr_ptr, mpfr_srcptr,
   2496                                       mpfr_rnd_t);
   2497 
   2498 __MPFR_DECLSPEC double mpfr_scale2 (double, int);
   2499 
   2500 __MPFR_DECLSPEC void mpfr_div_ui2 (mpfr_ptr, mpfr_srcptr, unsigned long,
   2501                                    unsigned long, mpfr_rnd_t);
   2502 
   2503 __MPFR_DECLSPEC void mpfr_gamma_one_and_two_third (mpfr_ptr, mpfr_ptr,
   2504                                                    mpfr_prec_t);
   2505 
   2506 __MPFR_DECLSPEC void mpfr_mpz_init (mpz_ptr);
   2507 __MPFR_DECLSPEC void mpfr_mpz_init2 (mpz_ptr, mp_bitcnt_t);
   2508 __MPFR_DECLSPEC void mpfr_mpz_clear (mpz_ptr);
   2509 
   2510 __MPFR_DECLSPEC int mpfr_odd_p (mpfr_srcptr);
   2511 
   2512 __MPFR_DECLSPEC int mpfr_nbits_ulong (unsigned long);
   2513 #ifdef _MPFR_H_HAVE_INTMAX_T
   2514 __MPFR_DECLSPEC int mpfr_nbits_uj (uintmax_t);
   2515 #endif
   2516 
   2517 #ifdef _MPFR_H_HAVE_VA_LIST
   2518 /* Declared only if <stdarg.h> has been included. */
   2519 __MPFR_DECLSPEC int mpfr_vasnprintf_aux (char**, char*, size_t, const char*,
   2520                                          va_list);
   2521 #endif
   2522 
   2523 #if MPFR_WANT_ASSERT >= 2
   2524 __MPFR_DECLSPEC void flags_fout (FILE *, mpfr_flags_t);
   2525 #endif
   2526 
   2527 #if defined (__cplusplus)
   2528 }
   2529 #endif
   2530 
   2531 
   2532 /*****************************************************
   2533  ***************  Internal mpz_t pool  ***************
   2534  *****************************************************/
   2535 
   2536 /* don't use the mpz_t pool with mini-gmp */
   2537 #ifdef MPFR_USE_MINI_GMP
   2538 # define MPFR_POOL_NENTRIES 0
   2539 #endif
   2540 
   2541 #ifndef MPFR_POOL_NENTRIES
   2542 # define MPFR_POOL_NENTRIES 32  /* default number of entries of the pool */
   2543 #endif
   2544 
   2545 #if MPFR_POOL_NENTRIES && !defined(MPFR_POOL_DONT_REDEFINE)
   2546 # undef mpz_init
   2547 # undef mpz_init2
   2548 # undef mpz_clear
   2549 # undef mpz_init_set_ui
   2550 # undef mpz_init_set
   2551 # define mpz_init mpfr_mpz_init
   2552 # define mpz_init2 mpfr_mpz_init2
   2553 # define mpz_clear mpfr_mpz_clear
   2554 # define mpz_init_set_ui(a,b) do { mpz_init (a); mpz_set_ui (a, b); } while (0)
   2555 # define mpz_init_set(a,b) do { mpz_init (a); mpz_set (a, b); } while (0)
   2556 #endif
   2557 
   2558 
   2559 /******************************************************
   2560  ********  Compute LOG2(LOG2(MPFR_PREC_MAX))  *********
   2561  ******************************************************/
   2562 
   2563 #if   _MPFR_PREC_FORMAT == 1
   2564 # define MPFR_PREC_MAX_TEMP USHRT_MAX
   2565 #elif _MPFR_PREC_FORMAT == 2
   2566 # define MPFR_PREC_MAX_TEMP UINT_MAX
   2567 #elif _MPFR_PREC_FORMAT == 3
   2568 # define MPFR_PREC_MAX_TEMP ULONG_MAX
   2569 #else
   2570 # error "Invalid MPFR Prec format"
   2571 #endif
   2572 
   2573 /* Note: In the constants below, it is sufficient to use the suffix U.
   2574  * A large enough unsigned type will be chosen automatically, but the
   2575  * exact type doesn't matter here.
   2576  */
   2577 
   2578 #if MPFR_PREC_MAX_TEMP == 255U
   2579 # define MPFR_PREC_BITS 8
   2580 # define MPFR_LOG2_PREC_BITS 3
   2581 #elif MPFR_PREC_MAX_TEMP == 65535U
   2582 # define MPFR_PREC_BITS 16
   2583 # define MPFR_LOG2_PREC_BITS 4
   2584 #elif MPFR_PREC_MAX_TEMP == 4294967295U
   2585 # define MPFR_PREC_BITS 32
   2586 # define MPFR_LOG2_PREC_BITS 5
   2587 #elif MPFR_PREC_MAX_TEMP == 18446744073709551615U
   2588 # define MPFR_PREC_BITS 64
   2589 # define MPFR_LOG2_PREC_BITS 6
   2590 #else
   2591 # error "Unsupported MPFR_PREC_MAX_TEMP value"
   2592 #endif
   2593 
   2594 
   2595 /******************************************************
   2596  *************  Value coverage checking  **************
   2597  ******************************************************/
   2598 
   2599 #ifdef MPFR_COV_CHECK
   2600 
   2601 /* Variable names should start with the __gmpfr_cov_ prefix. */
   2602 
   2603 #define MPFR_COV_SET(X) (__gmpfr_cov_ ## X = 1)
   2604 
   2605 #if defined (__cplusplus)
   2606 extern "C" {
   2607 #endif
   2608 
   2609 __MPFR_DECLSPEC extern int __gmpfr_cov_div_ui_sb[10][2];
   2610 __MPFR_DECLSPEC extern int __gmpfr_cov_sum_tmd[MPFR_RND_MAX][2][2][3][2][2];
   2611 
   2612 #if defined (__cplusplus)
   2613 }
   2614 #endif
   2615 
   2616 #else /* MPFR_COV_CHECK */
   2617 
   2618 #define MPFR_COV_SET(X) ((void) 0)
   2619 
   2620 #endif /* MPFR_COV_CHECK */
   2621 
   2622 
   2623 /******************************************************
   2624  *****************  Unbounded Floats  *****************
   2625  ******************************************************/
   2626 
   2627 #if defined (__cplusplus)
   2628 extern "C" {
   2629 #endif
   2630 
   2631 /* An UBF is like a MPFR number, but with an additional mpz_t member,
   2632    which is assumed to be present (with a value in it) when the usual
   2633    exponent field has the value MPFR_EXP_UBF. The goal of this compatible
   2634    representation is to easily be able to support UBF in "normal" code
   2635    using the public API. This is some form of "subtyping".
   2636 
   2637    Unfortunately this breaks aliasing rules, and C does not provide any way
   2638    to avoid that (except with additional syntax ugliness and API breakage,
   2639    though there is a workaround -- see the end of this comment):
   2640 
   2641      https://news.ycombinator.com/item?id=11753236
   2642 
   2643    The alignment requirement for __mpfr_ubf_struct (UBF) needs to be at least
   2644    as strong as the one for __mpfr_struct (MPFR number); this is not required
   2645    by the C standard, but this should always be the case in practice, since
   2646    __mpfr_ubf_struct starts with the same members as those of __mpfr_struct.
   2647    If ever this would not be the case with some particular C implementation,
   2648    an _Alignas alignment attribute (C11) could be added for UBF.
   2649 
   2650    When an input of a public function is an UBF, the semantic remains
   2651    internal to MPFR and can change in the future. UBF arguments need
   2652    to be explicitly converted to mpfr_ptr (or mpfr_srcptr); be careful
   2653    with variadic functions, as the compiler will not be able to check
   2654    in general. See fmma.c as an example of usage.
   2655 
   2656    In general, the type used for values that may be UBF must be either
   2657    mpfr_ubf_t or mpfr_ubf_ptr. The type mpfr_ptr or mpfr_srcptr may be
   2658    used for UBF only in the case where the pointer has been converted
   2659    from mpfr_ubf_ptr, in order to ensure valid alignment. For instance,
   2660    in mpfr_fmma_aux, one uses mpfr_ubf_t to generate the exact products
   2661    as UBF; then the corresponding pointers are converted to mpfr_srcptr
   2662    for mpfr_add (even though they point to UBF).
   2663 
   2664    Functions that can accept either MPFR arguments (mpfr_ptr type) or
   2665    UBF arguments (mpfr_ubf_ptr type) must use a pointer type that can
   2666    always be converted from both, typically mpfr_ptr or mpfr_srcptr.
   2667    For instance, that's why mpfr_ubf_exp_less_p uses mpfr_srcptr.
   2668    Note: "void *" could also be used, but mpfr_ptr is more meaningful
   2669    and practical.
   2670 
   2671    Note that functions used for logging need to support UBF (currently
   2672    done by printing that a number is an UBF, as it may be difficult to
   2673    do more without significant changes).
   2674 
   2675    --------
   2676 
   2677    A workaround to avoid breaking aliasing rules should be to use mpfr_ptr
   2678    to access the usual mpfr_t members and mpfr_ubf_ptr to access the
   2679    additional member _mpfr_zexp. And never use __mpfr_ubf_struct as a
   2680    declared type; otherwise this would force __mpfr_ubf_struct to be the
   2681    effective type of the whole object. Thus instead of
   2682 
   2683      typedef __mpfr_ubf_struct mpfr_ubf_t[1];
   2684 
   2685    one could use the following definition as a trick to allocate an UBF as
   2686    an automatic variable with the required alignment but without forcing
   2687    the effective type to __mpfr_ubf_struct.
   2688 
   2689       typedef union {
   2690         __mpfr_ubf_struct u[1];
   2691         __mpfr_struct m[1];
   2692       } mpfr_ubf_t;
   2693 
   2694    Then adapt the related code to select to right member, depending on the
   2695    context. Unfortunately, this triggers -Wstrict-aliasing false positives
   2696    with GCC in the MPFR_UBF_CLEAR_EXP expansion:
   2697 
   2698      https://gcc.gnu.org/bugzilla/show_bug.cgi?id=94337
   2699 
   2700    (see changeset r13820 in the ubf2 branch). So, for the time being,
   2701    as long as the code does not break, do not change anything.
   2702 
   2703    Note: The condition "use mpfr_ptr to access the usual mpfr_t members and
   2704    mpfr_ubf_ptr to access the additional member _mpfr_zexp" may be ignored
   2705    if the union type is visible within the function (see ISO C99 6.5.2.3#5
   2706    and 6.5.2.3#8 for the example, this implementation being very similar to
   2707    the valid fragment of this example), which must be the case as the union
   2708    is declared globally. However, this seems to be buggy in GCC:
   2709 
   2710      https://gcc.gnu.org/bugzilla/show_bug.cgi?id=14319
   2711      https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65892
   2712 
   2713    Alternatively, GCC's may_alias attribute could conditionally be used
   2714    on the __mpfr_ubf_struct and __mpfr_struct types (though it would be
   2715    much stronger than needed since only these two types may alias each
   2716    other).
   2717 */
   2718 
   2719 typedef struct {
   2720   mpfr_prec_t  _mpfr_prec;
   2721   mpfr_sign_t  _mpfr_sign;
   2722   mpfr_exp_t   _mpfr_exp;
   2723   mp_limb_t   *_mpfr_d;
   2724   mpz_t        _mpfr_zexp;
   2725 } __mpfr_ubf_struct;
   2726 
   2727 typedef __mpfr_ubf_struct mpfr_ubf_t[1];
   2728 typedef __mpfr_ubf_struct *mpfr_ubf_ptr;
   2729 
   2730 __MPFR_DECLSPEC void mpfr_ubf_mul_exact (mpfr_ubf_ptr,
   2731                                          mpfr_srcptr, mpfr_srcptr);
   2732 __MPFR_DECLSPEC int mpfr_ubf_exp_less_p (mpfr_srcptr, mpfr_srcptr);
   2733 __MPFR_DECLSPEC mpfr_exp_t mpfr_ubf_zexp2exp (mpz_ptr);
   2734 __MPFR_DECLSPEC mpfr_exp_t mpfr_ubf_diff_exp (mpfr_srcptr, mpfr_srcptr);
   2735 
   2736 #if defined (__cplusplus)
   2737 }
   2738 #endif
   2739 
   2740 /* Get the _mpfr_zexp field (pointer to a mpz_t) of an UBF object.
   2741    For practical reasons, the type of the argument x can be either
   2742    mpfr_ubf_ptr or mpfr_ptr, since the latter is used in functions
   2743    that accept both MPFR numbers and UBF's; this is checked by the
   2744    code "(x)->_mpfr_exp" (the "sizeof" prevents an access, which
   2745    could be invalid when MPFR_ZEXP(x) is used for an assignment,
   2746    and also avoids breaking the aliasing rules if they are dealt
   2747    with in the future).
   2748    This macro can be used when building an UBF. So we do not check
   2749    that the _mpfr_exp field has the value MPFR_EXP_UBF. */
   2750 #define MPFR_ZEXP(x)                            \
   2751   ((void) sizeof ((x)->_mpfr_exp),              \
   2752    ((mpfr_ubf_ptr) (x))->_mpfr_zexp)
   2753 
   2754 /* If x is an UBF, clear its mpz_t exponent. */
   2755 #define MPFR_UBF_CLEAR_EXP(x) \
   2756   ((void) (MPFR_IS_UBF (x) && (mpz_clear (MPFR_ZEXP (x)), 0)))
   2757 
   2758 /* Like MPFR_GET_EXP, but accepts UBF (with exponent saturated to
   2759    the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]). */
   2760 #define MPFR_UBF_GET_EXP(x)                                     \
   2761   (MPFR_IS_UBF (x) ? mpfr_ubf_zexp2exp (MPFR_ZEXP (x)) :        \
   2762    MPFR_GET_EXP ((mpfr_ptr) (x)))
   2763 
   2764 #endif /* __MPFR_IMPL_H__ */
   2765