xref: /netbsd/external/lgpl3/mpfr/dist/src/mpfr-impl.h (revision 606004a0)
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    Note, however, that this is not allowed by the C standard, which just
1659    says "character" and not "multibyte character".
1660    In the mean time, in case of non-single-byte character, revert to the
1661    default value. */
1662 #if MPFR_LCONV_DPTS
1663 #include <locale.h>
1664 /* Warning! In case of signed char, the value of MPFR_DECIMAL_POINT may
1665    be negative (the ISO C99 does not seem to forbid negative values). */
1666 #define MPFR_DECIMAL_POINT                      \
1667   (localeconv()->decimal_point[1] != '\0' ?     \
1668    (char) '.' : localeconv()->decimal_point[0])
1669 #define MPFR_THOUSANDS_SEPARATOR                \
1670   (localeconv()->thousands_sep[0] == '\0' ||    \
1671    localeconv()->thousands_sep[1] != '\0' ?     \
1672    (char) '\0' : localeconv()->thousands_sep[0])
1673 #else
1674 #define MPFR_DECIMAL_POINT ((char) '.')
1675 #define MPFR_THOUSANDS_SEPARATOR ((char) '\0')
1676 #endif
1677 
1678 /* Size of an array, as a constant expression. */
1679 #define numberof_const(x)  (sizeof (x) / sizeof ((x)[0]))
1680 
1681 /* Size of an array, safe version but not a constant expression:
1682    Since an array can silently be converted to a pointer, we check
1683    that this macro is applied on an array, not a pointer.
1684    Also make sure that the type is signed ("long" is sufficient
1685    in practice since the sizes come from the MPFR source), so that
1686    the value can be used in arbitrary expressions without the risk
1687    of silently switching to unsigned arithmetic. */
1688 /* TODO: Make numberof() a constant expression and always use it in
1689    the MPFR code instead of numberof_const(). See the tricks at
1690      https://gcc.gnu.org/pipermail/gcc/2020-September/233763.html
1691      "[PATCH v2] <sys/param.h>: Add nitems() and snitems() macros"
1692      by Alejandro Colomar
1693    but this needs to be fully tested on various platforms and with
1694    various compilers and compilation options.
1695    Moreover, change "long" to "ptrdiff_t", as used at the above URL? */
1696 #undef numberof
1697 #if 0
1698 /* The following should work with GCC as documented in its manual,
1699    but fails: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=38377#c10
1700    Thus disabled for now. */
1701 # define numberof(x)                                                    \
1702   ( __extension__ ({                                                    \
1703       int is_array = (void *) &(x) == (void *) &(x)[0];                 \
1704       MPFR_STAT_STATIC_ASSERT (__builtin_constant_p (is_array) ?        \
1705                                is_array : 1);                           \
1706       MPFR_ASSERTN (is_array);                                          \
1707       (long) numberof_const (x);                                        \
1708     }) )
1709 #else
1710 # define numberof(x)                                    \
1711   (MPFR_ASSERTN ((void *) &(x) == (void *) &(x)[0]),    \
1712    (long) numberof_const (x))
1713 #endif
1714 
1715 /* Addition with carry (detected by GCC and other good compilers). */
1716 #define ADD_LIMB(u,v,c) ((u) += (v), (c) = (u) < (v))
1717 
1718 /* umul_hi(h, x, y) puts in h the high part of x*y */
1719 /* MPFR_NEED_LONGLONG_H needs to be defined to use it. */
1720 #define umul_hi(h, x, y)                        \
1721   do {                                          \
1722     mp_limb_t _l;                               \
1723     umul_ppmm (h, _l, x, y);                    \
1724     (void) _l;  /* unused variable */           \
1725   } while (0)
1726 
1727 
1728 /******************************************************
1729  ************  Save exponent/flags macros  ************
1730  ******************************************************/
1731 
1732 /* See README.dev for details on how to use the macros.
1733    They are used to set the exponent range to the maximum
1734    temporarily */
1735 
1736 typedef struct {
1737   mpfr_flags_t saved_flags;
1738   mpfr_exp_t saved_emin;
1739   mpfr_exp_t saved_emax;
1740 } mpfr_save_expo_t;
1741 
1742 #define MPFR_SAVE_EXPO_DECL(x) mpfr_save_expo_t x
1743 #define MPFR_SAVE_EXPO_MARK(x)     \
1744  ((x).saved_flags = __gmpfr_flags, \
1745   (x).saved_emin = __gmpfr_emin,   \
1746   (x).saved_emax = __gmpfr_emax,   \
1747   __gmpfr_emin = MPFR_EMIN_MIN,    \
1748   __gmpfr_emax = MPFR_EMAX_MAX)
1749 #define MPFR_SAVE_EXPO_FREE(x)     \
1750  (__gmpfr_flags = (x).saved_flags, \
1751   __gmpfr_emin = (x).saved_emin,   \
1752   __gmpfr_emax = (x).saved_emax)
1753 #define MPFR_SAVE_EXPO_UPDATE_FLAGS(x, flags)  \
1754   (x).saved_flags |= (flags)
1755 
1756 /* Speed up final checking */
1757 #define mpfr_check_range(x,t,r) \
1758   (MPFR_LIKELY (MPFR_EXP_IN_RANGE (MPFR_EXP (x)))                \
1759    ? ((t) ? (__gmpfr_flags |= MPFR_FLAGS_INEXACT, (t)) : 0)      \
1760    : mpfr_check_range(x,t,r))
1761 
1762 
1763 /******************************************************
1764  *****************  Inline rounding  ******************
1765  ******************************************************/
1766 
1767 /*
1768  * Note: due to the labels, one cannot use a macro MPFR_RNDRAW* more than
1769  * once in a function (otherwise these labels would not be unique).
1770  */
1771 
1772 /*
1773  * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
1774  * assuming dest's sign is sign.
1775  * In rounding to nearest mode, execute MIDDLE_HANDLER when the value
1776  * is the middle of two consecutive numbers in dest precision.
1777  * Execute OVERFLOW_HANDLER in case of overflow when rounding.
1778  *
1779  * Note: the exponent field of dest is not used, possibly except by the
1780  * handlers. It is the caller (via the handlers) who entirely decides
1781  * how to handle it.
1782  */
1783 #define MPFR_RNDRAW_GEN(inexact, dest, srcp, sprec, rnd, sign,              \
1784                         MIDDLE_HANDLER, OVERFLOW_HANDLER)                   \
1785   do {                                                                      \
1786     mp_size_t _dests, _srcs;                                                \
1787     mp_limb_t *_destp;                                                      \
1788     mpfr_prec_t _destprec, _srcprec;                                        \
1789                                                                             \
1790     /* Check Trivial Case when Dest Mantissa has more bits than source */   \
1791     _srcprec = (sprec);                                                     \
1792     _destprec = MPFR_PREC (dest);                                           \
1793     MPFR_ASSERTD (_srcprec >= MPFR_PREC_MIN);                               \
1794     MPFR_ASSERTD (_destprec >= MPFR_PREC_MIN);                              \
1795     _destp = MPFR_MANT (dest);                                              \
1796     if (MPFR_UNLIKELY (_destprec >= _srcprec))                              \
1797       {                                                                     \
1798         _srcs  = MPFR_PREC2LIMBS (_srcprec);                                \
1799         _dests = MPFR_PREC2LIMBS (_destprec) - _srcs;                       \
1800         MPN_COPY (_destp + _dests, srcp, _srcs);                            \
1801         MPN_ZERO (_destp, _dests);                                          \
1802         inexact = 0;                                                        \
1803       }                                                                     \
1804     else                                                                    \
1805       {                                                                     \
1806         /* Non trivial case: rounding needed */                             \
1807         mpfr_prec_t _sh;                                                    \
1808         mp_limb_t *_sp;                                                     \
1809         mp_limb_t _rb, _sb, _ulp;                                           \
1810                                                                             \
1811         /* Compute Position and shift */                                    \
1812         _srcs  = MPFR_PREC2LIMBS (_srcprec);                                \
1813         _dests = MPFR_PREC2LIMBS (_destprec);                               \
1814         MPFR_UNSIGNED_MINUS_MODULO (_sh, _destprec);                        \
1815         _sp = (srcp) + _srcs - _dests;                                      \
1816                                                                             \
1817         /* General case when prec % GMP_NUMB_BITS != 0 */                   \
1818         if (MPFR_LIKELY (_sh != 0))                                         \
1819           {                                                                 \
1820             mp_limb_t _mask;                                                \
1821             /* Compute Rounding Bit and Sticky Bit */                       \
1822             /* Note: in directed rounding modes, if the rounding bit */     \
1823             /* is 1, the behavior does not depend on the sticky bit; */     \
1824             /* thus we will not try to compute it in this case (this */     \
1825             /* can be much faster and avoids to read uninitialized   */     \
1826             /* data in the current mpfr_mul implementation). We just */     \
1827             /* make sure that _sb is initialized.                    */     \
1828             _mask = MPFR_LIMB_ONE << (_sh - 1);                             \
1829             _rb = _sp[0] & _mask;                                           \
1830             _sb = _sp[0] & (_mask - 1);                                     \
1831             if ((rnd) == MPFR_RNDN || _rb == 0)                             \
1832               {                                                             \
1833                 mp_limb_t *_tmp;                                            \
1834                 mp_size_t _n;                                               \
1835                 for (_tmp = _sp, _n = _srcs - _dests ;                      \
1836                      _n != 0 && _sb == 0 ; _n--)                            \
1837                   _sb = *--_tmp;                                            \
1838               }                                                             \
1839             _ulp = 2 * _mask;                                               \
1840           }                                                                 \
1841         else /* _sh == 0 */                                                 \
1842           {                                                                 \
1843             MPFR_ASSERTD (_dests < _srcs);                                  \
1844             /* Compute Rounding Bit and Sticky Bit - see note above */      \
1845             _rb = _sp[-1] & MPFR_LIMB_HIGHBIT;                              \
1846             _sb = _sp[-1] & (MPFR_LIMB_HIGHBIT-1);                          \
1847             if ((rnd) == MPFR_RNDN || _rb == 0)                             \
1848               {                                                             \
1849                 mp_limb_t *_tmp;                                            \
1850                 mp_size_t _n;                                               \
1851                 for (_tmp = _sp - 1, _n = _srcs - _dests - 1 ;              \
1852                      _n != 0 && _sb == 0 ; _n--)                            \
1853                   _sb = *--_tmp;                                            \
1854               }                                                             \
1855             _ulp = MPFR_LIMB_ONE;                                           \
1856           }                                                                 \
1857         /* Rounding */                                                      \
1858         if (rnd == MPFR_RNDF)                                               \
1859           {                                                                 \
1860             inexact = 0;                                                    \
1861             goto trunc_doit;                                                \
1862           }                                                                 \
1863         else if (rnd == MPFR_RNDN)                                          \
1864           {                                                                 \
1865             if (_rb == 0)                                                   \
1866               {                                                             \
1867               trunc:                                                        \
1868                 inexact = MPFR_LIKELY ((_sb | _rb) != 0) ? -sign : 0;       \
1869               trunc_doit:                                                   \
1870                 MPN_COPY (_destp, _sp, _dests);                             \
1871                 _destp[0] &= ~(_ulp - 1);                                   \
1872               }                                                             \
1873             else if (MPFR_UNLIKELY (_sb == 0))                              \
1874               { /* Middle of two consecutive representable numbers */       \
1875                 MIDDLE_HANDLER;                                             \
1876               }                                                             \
1877             else                                                            \
1878               {                                                             \
1879                 if (0)                                                      \
1880                   goto addoneulp_doit; /* dummy code to avoid warning */    \
1881               addoneulp:                                                    \
1882                 inexact = sign;                                             \
1883               addoneulp_doit:                                               \
1884                 if (MPFR_UNLIKELY (mpn_add_1 (_destp, _sp, _dests, _ulp)))  \
1885                   {                                                         \
1886                     _destp[_dests - 1] = MPFR_LIMB_HIGHBIT;                 \
1887                     OVERFLOW_HANDLER;                                       \
1888                   }                                                         \
1889                 _destp[0] &= ~(_ulp - 1);                                   \
1890               }                                                             \
1891           }                                                                 \
1892         else                                                                \
1893           { /* Directed rounding mode */                                    \
1894             if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))           \
1895               goto trunc;                                                   \
1896              else if (MPFR_UNLIKELY ((_sb | _rb) == 0))                     \
1897                {                                                            \
1898                  inexact = 0;                                               \
1899                  goto trunc_doit;                                           \
1900                }                                                            \
1901              else                                                           \
1902               goto addoneulp;                                               \
1903           }                                                                 \
1904       }                                                                     \
1905   } while (0)
1906 
1907 /*
1908  * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
1909  * assuming dest's sign is sign.
1910  * Execute OVERFLOW_HANDLER in case of overflow when rounding.
1911  */
1912 #define MPFR_RNDRAW(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER) \
1913    MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign,                   \
1914      if ((_sp[0] & _ulp) == 0)                                               \
1915        {                                                                     \
1916          inexact = -sign;                                                    \
1917          goto trunc_doit;                                                    \
1918        }                                                                     \
1919      else                                                                    \
1920        goto addoneulp;                                                       \
1921      , OVERFLOW_HANDLER)
1922 
1923 /*
1924  * Round mantissa (srcp, sprec) to mpfr_t dest using rounding mode rnd
1925  * assuming dest's sign is sign.
1926  * Execute OVERFLOW_HANDLER in case of overflow when rounding.
1927  * Set inexact to +/- MPFR_EVEN_INEX in case of even rounding.
1928  */
1929 #define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, \
1930                          OVERFLOW_HANDLER)                      \
1931    MPFR_RNDRAW_GEN (inexact, dest, srcp, sprec, rnd, sign,      \
1932      if ((_sp[0] & _ulp) == 0)                                  \
1933        {                                                        \
1934          inexact = -MPFR_EVEN_INEX * sign;                      \
1935          goto trunc_doit;                                       \
1936        }                                                        \
1937      else                                                       \
1938        {                                                        \
1939          inexact = MPFR_EVEN_INEX * sign;                       \
1940          goto addoneulp_doit;                                   \
1941        }                                                        \
1942      , OVERFLOW_HANDLER)
1943 
1944 /* Return TRUE if b is non singular and we can round it to precision 'prec'
1945    and determine the ternary value, with rounding mode 'rnd', and with
1946    error at most 'error' */
1947 #define MPFR_CAN_ROUND(b,err,prec,rnd)                                       \
1948  (!MPFR_IS_SINGULAR (b) && mpfr_round_p (MPFR_MANT (b), MPFR_LIMB_SIZE (b),  \
1949                                          (err), (prec) + ((rnd)==MPFR_RNDN)))
1950 
1951 /* Copy the sign and the significand, and handle the exponent in exp. */
1952 #define MPFR_SETRAW(inexact,dest,src,exp,rnd)                           \
1953   if (dest != src)                                                      \
1954     {                                                                   \
1955       MPFR_SET_SIGN (dest, MPFR_SIGN (src));                            \
1956       if (MPFR_PREC (dest) == MPFR_PREC (src))                          \
1957         {                                                               \
1958           MPN_COPY (MPFR_MANT (dest), MPFR_MANT (src),                  \
1959                     MPFR_LIMB_SIZE (src));                              \
1960           inexact = 0;                                                  \
1961         }                                                               \
1962       else                                                              \
1963         {                                                               \
1964           MPFR_RNDRAW (inexact, dest, MPFR_MANT (src), MPFR_PREC (src), \
1965                        rnd, MPFR_SIGN (src), exp++);                    \
1966         }                                                               \
1967     }                                                                   \
1968   else                                                                  \
1969     inexact = 0;
1970 
1971 /* TODO: fix this description (see round_near_x.c). */
1972 /* Assuming that the function has a Taylor expansion which looks like:
1973     y=o(f(x)) = o(v + g(x)) with |g(x)| <= 2^(EXP(v)-err)
1974    we can quickly set y to v if x is small (ie err > prec(y)+1) in most
1975    cases. It assumes that f(x) is not representable exactly as a FP number.
1976    v must not be a singular value (NAN, INF or ZERO); usual values are
1977    v=1 or v=x.
1978 
1979    y is the destination (a mpfr_t), v the value to set (a mpfr_t),
1980    err1+err2 with 0 <= err2 <= 3 the error term (mpfr_exp_t's), dir (an int)
1981    is the direction of the committed error (if dir = 0, it rounds toward 0,
1982    if dir=1, it rounds away from 0), rnd the rounding mode.
1983 
1984    It returns from the function a ternary value in case of success.
1985    If you want to free something, you must fill the "extra" field
1986    in consequences, otherwise put nothing in it.
1987 
1988    The test is less restrictive than necessary, but the function
1989    will finish the check itself.
1990 
1991    Note: err1 + err2 is allowed to overflow as mpfr_exp_t, but it must give
1992    its real value as mpfr_uexp_t.
1993 */
1994 #define MPFR_FAST_COMPUTE_IF_SMALL_INPUT(y,v,err1,err2,dir,rnd,extra)   \
1995   do {                                                                  \
1996     mpfr_ptr _y = (y);                                                  \
1997     mpfr_exp_t _err1 = (err1);                                          \
1998     mpfr_exp_t _err2 = (err2);                                          \
1999     if (_err1 > 0)                                                      \
2000       {                                                                 \
2001         mpfr_uexp_t _err = (mpfr_uexp_t) _err1 + _err2;                 \
2002         if (MPFR_UNLIKELY (_err > MPFR_PREC (_y) + 1))                  \
2003           {                                                             \
2004             int _inexact = mpfr_round_near_x (_y,(v),_err,(dir),(rnd)); \
2005             if (_inexact != 0)                                          \
2006               {                                                         \
2007                 extra;                                                  \
2008                 return _inexact;                                        \
2009               }                                                         \
2010           }                                                             \
2011       }                                                                 \
2012   } while (0)
2013 
2014 /* Variant, to be called somewhere after MPFR_SAVE_EXPO_MARK. This variant
2015    is needed when there are some computations before or when some non-zero
2016    real constant is used, such as __gmpfr_one for mpfr_cos. */
2017 #define MPFR_SMALL_INPUT_AFTER_SAVE_EXPO(y,v,err1,err2,dir,rnd,expo,extra) \
2018   do {                                                                  \
2019     mpfr_ptr _y = (y);                                                  \
2020     mpfr_exp_t _err1 = (err1);                                          \
2021     mpfr_exp_t _err2 = (err2);                                          \
2022     if (_err1 > 0)                                                      \
2023       {                                                                 \
2024         mpfr_uexp_t _err = (mpfr_uexp_t) _err1 + _err2;                 \
2025         if (MPFR_UNLIKELY (_err > MPFR_PREC (_y) + 1))                  \
2026           {                                                             \
2027             int _inexact;                                               \
2028             MPFR_CLEAR_FLAGS ();                                        \
2029             _inexact = mpfr_round_near_x (_y,(v),_err,(dir),(rnd));     \
2030             if (_inexact != 0)                                          \
2031               {                                                         \
2032                 extra;                                                  \
2033                 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);      \
2034                 MPFR_SAVE_EXPO_FREE (expo);                             \
2035                 return mpfr_check_range (_y, _inexact, (rnd));          \
2036               }                                                         \
2037           }                                                             \
2038       }                                                                 \
2039   } while (0)
2040 
2041 
2042 /******************************************************
2043  *****************  Ziv loop macros  ******************
2044  ******************************************************/
2045 
2046 /* To safely increase some precision, detecting integer overflows.
2047    This macro is particularly useful when determining the initial
2048    working precision before Ziv's loop. P is a precision, X is an
2049    arbitrary non-negative integer.
2050    Note: On 2012-02-23, the MPFR_PREC_MAX value has been decreased
2051    by 256 from the maximum value representable in the mpfr_prec_t
2052    type, in order to avoid some integer overflows when this macro
2053    is not used (if the result is larger than MPFR_PREC_MAX, this
2054    should be detected with a later assertion, e.g. in mpfr_init2).
2055    But this change is mainly for existing code that has not been
2056    updated yet. So, it is advised to always use MPFR_ADD_PREC or
2057    MPFR_INC_PREC if the result can be larger than MPFR_PREC_MAX. */
2058 #define MPFR_ADD_PREC(P,X) \
2059   (MPFR_ASSERTN ((X) <= MPFR_PREC_MAX - (P)), (P) + (X))
2060 #define MPFR_INC_PREC(P,X) \
2061   (MPFR_ASSERTN ((X) <= MPFR_PREC_MAX - (P)), (P) += (X))
2062 
2063 #ifndef MPFR_USE_LOGGING
2064 
2065 #define MPFR_ZIV_DECL(_x) mpfr_prec_t _x
2066 #define MPFR_ZIV_INIT(_x, _p) (_x) = GMP_NUMB_BITS
2067 #define MPFR_ZIV_NEXT(_x, _p) (MPFR_INC_PREC (_p, _x), (_x) = (_p)/2)
2068 #define MPFR_ZIV_FREE(x)
2069 
2070 #else
2071 
2072 /* The following test on glibc is there mainly for Darwin (Mac OS X), to
2073    obtain a better error message. The real test should have been a test
2074    concerning nested functions in gcc, which are disabled by default on
2075    Darwin; but it is not possible to do that without a configure test. */
2076 # if defined (__cplusplus) || !(__MPFR_GNUC(3,0) && __MPFR_GLIBC(2,0))
2077 #  error "Logging not supported (needs gcc >= 3.0 and GNU C Library >= 2.0)."
2078 # endif
2079 
2080 /* Use LOGGING */
2081 
2082 /* Note: the mpfr_log_level >= 0 below avoids to take into account
2083    Ziv loops used by the MPFR functions called by the mpfr_fprintf
2084    in LOG_PRINT. */
2085 
2086 #define MPFR_ZIV_DECL(_x)                                               \
2087   mpfr_prec_t _x;                                                       \
2088   int _x ## _cpt = 1;                                                   \
2089   static unsigned long  _x ## _loop = 0, _x ## _bad = 0;                \
2090   static const char *_x ## _fname = __func__;                           \
2091   auto void __attribute__ ((destructor)) x ## _f  (void);               \
2092   void __attribute__ ((destructor)) x ## _f  (void) {                   \
2093     if (_x ## _loop != 0 && (MPFR_LOG_STAT_F & mpfr_log_type)) {        \
2094       fprintf (mpfr_log_file,                                           \
2095                "%s: Ziv failed %2.2f%% (%lu bad cases / %lu calls)\n",  \
2096                _x ## _fname, (double) 100.0 * _x ## _bad / _x ## _loop, \
2097                _x ## _bad, _x ## _loop );                               \
2098       if (mpfr_log_flush)                                               \
2099         fflush (mpfr_log_file);                                         \
2100     }                                                                   \
2101   }
2102 
2103 #define MPFR_ZIV_INIT(_x, _p)                                           \
2104   do                                                                    \
2105     {                                                                   \
2106       (_x) = GMP_NUMB_BITS;                                             \
2107       if (mpfr_log_level >= 0)                                          \
2108         _x ## _loop ++;                                                 \
2109       LOG_PRINT (MPFR_LOG_BADCASE_F, "%s:ZIV 1st prec=%Pd\n",           \
2110                  __func__, (mpfr_prec_t) (_p));                         \
2111     }                                                                   \
2112   while (0)
2113 
2114 #define MPFR_ZIV_NEXT(_x, _p)                                           \
2115   do                                                                    \
2116     {                                                                   \
2117       MPFR_INC_PREC (_p, _x);                                           \
2118       (_x) = (_p) / 2;                                                  \
2119       if (mpfr_log_level >= 0)                                          \
2120         _x ## _bad += (_x ## _cpt == 1);                                \
2121       _x ## _cpt ++;                                                    \
2122       LOG_PRINT (MPFR_LOG_BADCASE_F, "%s:ZIV new prec=%Pd\n",           \
2123                  __func__, (mpfr_prec_t) (_p));                         \
2124     }                                                                   \
2125   while (0)
2126 
2127 #define MPFR_ZIV_FREE(_x)                                               \
2128   do                                                                    \
2129     if (_x ## _cpt > 1)                                                 \
2130       LOG_PRINT (MPFR_LOG_BADCASE_F, "%s:ZIV %d loops\n",               \
2131                  __func__, _x ## _cpt);                                 \
2132   while (0)
2133 
2134 #endif
2135 
2136 
2137 /******************************************************
2138  ******************  Logging macros  ******************
2139  ******************************************************/
2140 
2141 /* The different kind of LOG */
2142 #define MPFR_LOG_INPUT_F    1
2143 #define MPFR_LOG_OUTPUT_F   2
2144 #define MPFR_LOG_INTERNAL_F 4
2145 #define MPFR_LOG_TIME_F     8
2146 #define MPFR_LOG_BADCASE_F  16
2147 #define MPFR_LOG_MSG_F      32
2148 #define MPFR_LOG_STAT_F     64
2149 
2150 #ifdef MPFR_USE_LOGGING
2151 
2152 /* Check if we can support this feature */
2153 # ifdef MPFR_USE_THREAD_SAFE
2154 #  error "Enable either `Logging' or `thread-safe', not both"
2155 # endif
2156 # if !__MPFR_GNUC(3,0)
2157 #  error "Logging not supported (GCC >= 3.0)"
2158 # endif
2159 
2160 #if defined (__cplusplus)
2161 extern "C" {
2162 #endif
2163 
2164 __MPFR_DECLSPEC extern FILE *mpfr_log_file;
2165 __MPFR_DECLSPEC extern int   mpfr_log_flush;
2166 __MPFR_DECLSPEC extern int   mpfr_log_type;
2167 __MPFR_DECLSPEC extern int   mpfr_log_level;
2168 __MPFR_DECLSPEC extern int   mpfr_log_current;
2169 __MPFR_DECLSPEC extern mpfr_prec_t mpfr_log_prec;
2170 
2171 #if defined (__cplusplus)
2172  }
2173 #endif
2174 
2175 /* LOG_PRINT calls mpfr_fprintf on mpfr_log_file with logging disabled
2176    (recursive logging is not wanted and freezes MPFR). */
2177 #define LOG_PRINT(type, format, ...)                                    \
2178   do                                                                    \
2179     if ((mpfr_log_type & (type)) && mpfr_log_current <= mpfr_log_level) \
2180       {                                                                 \
2181         int old_level = mpfr_log_level;                                 \
2182         mpfr_log_level = -1;  /* disable logging in mpfr_fprintf */     \
2183         __gmpfr_cache_const_pi = __gmpfr_logging_pi;                    \
2184         __gmpfr_cache_const_log2 = __gmpfr_logging_log2;                \
2185         mpfr_fprintf (mpfr_log_file, format, __VA_ARGS__);              \
2186         if (mpfr_log_flush)                                             \
2187           fflush (mpfr_log_file);                                       \
2188         mpfr_log_level = old_level;                                     \
2189         __gmpfr_cache_const_pi = __gmpfr_normal_pi;                     \
2190         __gmpfr_cache_const_log2 = __gmpfr_normal_log2;                 \
2191       }                                                                 \
2192   while (0)
2193 
2194 #define MPFR_LOG_VAR(x)                                                 \
2195   LOG_PRINT (MPFR_LOG_INTERNAL_F, "%s.%d:%s[%#Pu]=%.*Rg\n", __func__,   \
2196              (int) __LINE__, #x, mpfr_get_prec (x), mpfr_log_prec, x)
2197 
2198 #define MPFR_LOG_MSG2(format, ...)                                      \
2199   LOG_PRINT (MPFR_LOG_MSG_F, "%s.%d: "format, __func__, (int) __LINE__, \
2200              __VA_ARGS__)
2201 #define MPFR_LOG_MSG(x) MPFR_LOG_MSG2 x
2202 
2203 #define MPFR_LOG_BEGIN2(format, ...)                                    \
2204   mpfr_log_current ++;                                                  \
2205   LOG_PRINT (MPFR_LOG_INPUT_F, "%s:IN  flags=%x "format"\n", __func__,  \
2206              (unsigned int) __gmpfr_flags, __VA_ARGS__);                \
2207   if ((MPFR_LOG_TIME_F & mpfr_log_type) &&                              \
2208       (mpfr_log_current <= mpfr_log_level))                             \
2209     __gmpfr_log_time = mpfr_get_cputime ();
2210 #define MPFR_LOG_BEGIN(x)                                               \
2211   int __gmpfr_log_time = 0;                                             \
2212   MPFR_LOG_BEGIN2 x
2213 
2214 #define MPFR_LOG_END2(format, ...)                                      \
2215   LOG_PRINT (MPFR_LOG_TIME_F, "%s:TIM %dms\n", __mpfr_log_fname,        \
2216              mpfr_get_cputime () - __gmpfr_log_time);                   \
2217   LOG_PRINT (MPFR_LOG_OUTPUT_F, "%s:OUT flags=%x "format"\n",           \
2218              __mpfr_log_fname, (unsigned int) __gmpfr_flags,            \
2219              __VA_ARGS__);                                              \
2220   mpfr_log_current --;
2221 #define MPFR_LOG_END(x)                                                 \
2222   static const char *__mpfr_log_fname = __func__;                       \
2223   MPFR_LOG_END2 x
2224 
2225 #define MPFR_LOG_FUNC(begin,end)                                        \
2226   static const char *__mpfr_log_fname = __func__;                       \
2227   auto void __mpfr_log_cleanup (int *time);                             \
2228   void __mpfr_log_cleanup (int *time) {                                 \
2229     int __gmpfr_log_time = *time;                                       \
2230     MPFR_LOG_END2 end; }                                                \
2231   int __gmpfr_log_time __attribute__ ((cleanup (__mpfr_log_cleanup)));  \
2232   __gmpfr_log_time = 0;                                                 \
2233   MPFR_LOG_BEGIN2 begin
2234 
2235 #else /* MPFR_USE_LOGGING */
2236 
2237 /* Define void macro for logging */
2238 
2239 #define MPFR_LOG_VAR(x)
2240 #define MPFR_LOG_BEGIN(x)
2241 #define MPFR_LOG_END(x)
2242 #define MPFR_LOG_MSG(x)
2243 #define MPFR_LOG_FUNC(x,y)
2244 
2245 #endif /* MPFR_USE_LOGGING */
2246 
2247 
2248 /**************************************************************
2249  ************  Group Initialize Functions Macros  *************
2250  **************************************************************/
2251 
2252 #ifndef MPFR_GROUP_STATIC_SIZE
2253 # define MPFR_GROUP_STATIC_SIZE 16
2254 #endif
2255 
2256 struct mpfr_group_t {
2257   size_t     alloc;
2258   mp_limb_t *mant;
2259 #if MPFR_GROUP_STATIC_SIZE != 0
2260   mp_limb_t  tab[MPFR_GROUP_STATIC_SIZE];
2261 #else
2262   /* In order to detect memory leaks when testing, MPFR_GROUP_STATIC_SIZE
2263      can be set to 0, in which case tab will not be used. ISO C does not
2264      support zero-length arrays[*], thus let's use a flexible array member
2265      (which will be equivalent here). Note: this is new in C99, but this
2266      is just used for testing.
2267      [*] Zero-length arrays are a GNU extension:
2268            https://gcc.gnu.org/onlinedocs/gcc/Zero-Length.html
2269          and as such an extension is forbidden in ISO C, it triggers an
2270          error with -Werror=pedantic.
2271   */
2272   mp_limb_t  tab[];
2273 #endif
2274 };
2275 
2276 #define MPFR_GROUP_DECL(g) struct mpfr_group_t g
2277 #define MPFR_GROUP_CLEAR(g) do {                                 \
2278  MPFR_LOG_MSG (("GROUP_CLEAR: ptr = 0x%lX, size = %lu\n",        \
2279                 (unsigned long) (g).mant,                        \
2280                 (unsigned long) (g).alloc));                     \
2281  if ((g).alloc != 0) {                                           \
2282    MPFR_ASSERTD ((g).mant != (g).tab);                           \
2283    mpfr_free_func ((g).mant, (g).alloc);                         \
2284  }} while (0)
2285 
2286 #define MPFR_GROUP_INIT_TEMPLATE(g, prec, num, handler) do {            \
2287  mpfr_prec_t _prec = (prec);                                            \
2288  mp_size_t _size;                                                       \
2289  MPFR_ASSERTD (_prec >= MPFR_PREC_MIN);                                 \
2290  if (MPFR_UNLIKELY (_prec > MPFR_PREC_MAX))                             \
2291    mpfr_abort_prec_max ();                                              \
2292  _size = MPFR_PREC2LIMBS (_prec);                                       \
2293  if (_size * (num) > MPFR_GROUP_STATIC_SIZE)                            \
2294    {                                                                    \
2295      (g).alloc = (num) * _size * sizeof (mp_limb_t);                    \
2296      (g).mant = (mp_limb_t *) mpfr_allocate_func ((g).alloc);           \
2297    }                                                                    \
2298  else                                                                   \
2299    {                                                                    \
2300      (g).alloc = 0;                                                     \
2301      (g).mant = (g).tab;                                                \
2302    }                                                                    \
2303  MPFR_LOG_MSG (("GROUP_INIT: ptr = 0x%lX, size = %lu\n",                \
2304                 (unsigned long) (g).mant, (unsigned long) (g).alloc));  \
2305  handler;                                                               \
2306  } while (0)
2307 #define MPFR_GROUP_TINIT(g, n, x)                       \
2308   MPFR_TMP_INIT1 ((g).mant + _size * (n), x, _prec)
2309 
2310 #define MPFR_GROUP_INIT_1(g, prec, x)                            \
2311  MPFR_GROUP_INIT_TEMPLATE(g, prec, 1, MPFR_GROUP_TINIT(g, 0, x))
2312 #define MPFR_GROUP_INIT_2(g, prec, x, y)                         \
2313  MPFR_GROUP_INIT_TEMPLATE(g, prec, 2,                            \
2314    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y))
2315 #define MPFR_GROUP_INIT_3(g, prec, x, y, z)                      \
2316  MPFR_GROUP_INIT_TEMPLATE(g, prec, 3,                            \
2317    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
2318    MPFR_GROUP_TINIT(g, 2, z))
2319 #define MPFR_GROUP_INIT_4(g, prec, x, y, z, t)                   \
2320  MPFR_GROUP_INIT_TEMPLATE(g, prec, 4,                            \
2321    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
2322    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t))
2323 #define MPFR_GROUP_INIT_5(g, prec, x, y, z, t, a)                \
2324  MPFR_GROUP_INIT_TEMPLATE(g, prec, 5,                            \
2325    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
2326    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t);          \
2327    MPFR_GROUP_TINIT(g, 4, a))
2328 #define MPFR_GROUP_INIT_6(g, prec, x, y, z, t, a, b)             \
2329  MPFR_GROUP_INIT_TEMPLATE(g, prec, 6,                            \
2330    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
2331    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t);          \
2332    MPFR_GROUP_TINIT(g, 4, a);MPFR_GROUP_TINIT(g, 5, b))
2333 
2334 #define MPFR_GROUP_REPREC_TEMPLATE(g, prec, num, handler) do {          \
2335  mpfr_prec_t _prec = (prec);                                            \
2336  size_t    _oalloc = (g).alloc;                                         \
2337  mp_size_t _size;                                                       \
2338  MPFR_LOG_MSG (("GROUP_REPREC: oldptr = 0x%lX, oldsize = %lu\n",        \
2339                 (unsigned long) (g).mant, (unsigned long) _oalloc));    \
2340  MPFR_ASSERTD (_prec >= MPFR_PREC_MIN);                                 \
2341  if (MPFR_UNLIKELY (_prec > MPFR_PREC_MAX))                             \
2342    mpfr_abort_prec_max ();                                              \
2343  _size = MPFR_PREC2LIMBS (_prec);                                       \
2344  (g).alloc = (num) * _size * sizeof (mp_limb_t);                        \
2345  if (_oalloc == 0)                                                      \
2346    (g).mant = (mp_limb_t *) mpfr_allocate_func ((g).alloc);             \
2347  else                                                                   \
2348    (g).mant = (mp_limb_t *)                                             \
2349      mpfr_reallocate_func ((g).mant, _oalloc, (g).alloc);               \
2350  MPFR_LOG_MSG (("GROUP_REPREC: newptr = 0x%lX, newsize = %lu\n",        \
2351                 (unsigned long) (g).mant, (unsigned long) (g).alloc));  \
2352  handler;                                                               \
2353  } while (0)
2354 
2355 #define MPFR_GROUP_REPREC_1(g, prec, x)                          \
2356  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 1, MPFR_GROUP_TINIT(g, 0, x))
2357 #define MPFR_GROUP_REPREC_2(g, prec, x, y)                       \
2358  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 2,                          \
2359    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y))
2360 #define MPFR_GROUP_REPREC_3(g, prec, x, y, z)                    \
2361  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 3,                          \
2362    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
2363    MPFR_GROUP_TINIT(g, 2, z))
2364 #define MPFR_GROUP_REPREC_4(g, prec, x, y, z, t)                 \
2365  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 4,                          \
2366    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
2367    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t))
2368 #define MPFR_GROUP_REPREC_5(g, prec, x, y, z, t, a)              \
2369  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 5,                          \
2370    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
2371    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t);          \
2372    MPFR_GROUP_TINIT(g, 4, a))
2373 #define MPFR_GROUP_REPREC_6(g, prec, x, y, z, t, a, b)           \
2374  MPFR_GROUP_REPREC_TEMPLATE(g, prec, 6,                          \
2375    MPFR_GROUP_TINIT(g, 0, x);MPFR_GROUP_TINIT(g, 1, y);          \
2376    MPFR_GROUP_TINIT(g, 2, z);MPFR_GROUP_TINIT(g, 3, t);          \
2377    MPFR_GROUP_TINIT(g, 4, a);MPFR_GROUP_TINIT(g, 5, b))
2378 
2379 
2380 /******************************************************
2381  ***************  Internal functions  *****************
2382  ******************************************************/
2383 
2384 #if defined (__cplusplus)
2385 extern "C" {
2386 #endif
2387 
2388 MPFR_COLD_FUNCTION_ATTR __MPFR_DECLSPEC int
2389   mpfr_underflow (mpfr_ptr, mpfr_rnd_t, int);
2390 MPFR_COLD_FUNCTION_ATTR __MPFR_DECLSPEC int
2391   mpfr_overflow (mpfr_ptr, mpfr_rnd_t, int);
2392 
2393 __MPFR_DECLSPEC int mpfr_add1 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
2394 __MPFR_DECLSPEC int mpfr_sub1 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
2395 __MPFR_DECLSPEC int mpfr_add1sp (mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
2396                                  mpfr_rnd_t);
2397 __MPFR_DECLSPEC int mpfr_sub1sp (mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
2398                                  mpfr_rnd_t);
2399 __MPFR_DECLSPEC int mpfr_can_round_raw (const mp_limb_t *,
2400              mp_size_t, int, mpfr_exp_t, mpfr_rnd_t, mpfr_rnd_t, mpfr_prec_t);
2401 
2402 __MPFR_DECLSPEC int mpfr_set_1_2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t, int);
2403 
2404 __MPFR_DECLSPEC int mpfr_cmp2 (mpfr_srcptr, mpfr_srcptr, mpfr_prec_t *);
2405 
2406 __MPFR_DECLSPEC long          __gmpfr_ceil_log2     (double);
2407 __MPFR_DECLSPEC long          __gmpfr_floor_log2    (double);
2408 __MPFR_DECLSPEC double        __gmpfr_ceil_exp2     (double);
2409 __MPFR_DECLSPEC unsigned long __gmpfr_isqrt     (unsigned long);
2410 __MPFR_DECLSPEC unsigned long __gmpfr_cuberoot  (unsigned long);
2411 __MPFR_DECLSPEC int       __gmpfr_int_ceil_log2 (unsigned long);
2412 
2413 __MPFR_DECLSPEC mpfr_exp_t mpfr_ceil_mul (mpfr_exp_t, int, int);
2414 
2415 __MPFR_DECLSPEC int mpfr_exp_2 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
2416 __MPFR_DECLSPEC int mpfr_exp_3 (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t);
2417 __MPFR_DECLSPEC int mpfr_powerof2_raw (mpfr_srcptr);
2418 __MPFR_DECLSPEC int mpfr_powerof2_raw2 (const mp_limb_t *, mp_size_t);
2419 
2420 __MPFR_DECLSPEC int mpfr_pow_general (mpfr_ptr, mpfr_srcptr, mpfr_srcptr,
2421                                       mpfr_rnd_t, int, mpfr_save_expo_t *);
2422 
2423 __MPFR_DECLSPEC void mpfr_setmax (mpfr_ptr, mpfr_exp_t);
2424 __MPFR_DECLSPEC void mpfr_setmin (mpfr_ptr, mpfr_exp_t);
2425 
2426 __MPFR_DECLSPEC int mpfr_mpn_exp (mp_limb_t *, mpfr_exp_t *, int,
2427                                   mpfr_exp_t, size_t);
2428 
2429 #ifdef _MPFR_H_HAVE_FILE
2430 __MPFR_DECLSPEC void mpfr_fdump (FILE *, mpfr_srcptr);
2431 #endif
2432 __MPFR_DECLSPEC void mpfr_print_mant_binary (const char*, const mp_limb_t*,
2433                                              mpfr_prec_t);
2434 __MPFR_DECLSPEC void mpfr_set_str_binary (mpfr_ptr, const char*);
2435 
2436 __MPFR_DECLSPEC int mpfr_round_raw (mp_limb_t *,
2437        const mp_limb_t *, mpfr_prec_t, int, mpfr_prec_t, mpfr_rnd_t, int *);
2438 __MPFR_DECLSPEC int mpfr_round_raw_2 (const mp_limb_t *, mpfr_prec_t, int,
2439                                       mpfr_prec_t, mpfr_rnd_t);
2440 /* No longer defined (see round_prec.c).
2441    Uncomment if it needs to be defined again.
2442 __MPFR_DECLSPEC int mpfr_round_raw_3 (const mp_limb_t *,
2443              mpfr_prec_t, int, mpfr_prec_t, mpfr_rnd_t, int *);
2444 */
2445 __MPFR_DECLSPEC int mpfr_round_raw_4 (mp_limb_t *,
2446        const mp_limb_t *, mpfr_prec_t, int, mpfr_prec_t, mpfr_rnd_t);
2447 
2448 #define mpfr_round_raw2(xp, xn, neg, r, prec) \
2449   mpfr_round_raw_2((xp),(xn)*GMP_NUMB_BITS,(neg),(prec),(r))
2450 
2451 __MPFR_DECLSPEC int mpfr_check (mpfr_srcptr);
2452 
2453 __MPFR_DECLSPEC int mpfr_get_cputime (void);
2454 
2455 __MPFR_DECLSPEC void mpfr_nexttozero (mpfr_ptr);
2456 __MPFR_DECLSPEC void mpfr_nexttoinf (mpfr_ptr);
2457 
2458 __MPFR_DECLSPEC int mpfr_const_pi_internal (mpfr_ptr,mpfr_rnd_t);
2459 __MPFR_DECLSPEC int mpfr_const_log2_internal (mpfr_ptr,mpfr_rnd_t);
2460 __MPFR_DECLSPEC int mpfr_const_euler_internal (mpfr_ptr, mpfr_rnd_t);
2461 __MPFR_DECLSPEC int mpfr_const_catalan_internal (mpfr_ptr, mpfr_rnd_t);
2462 
2463 #if 0
2464 __MPFR_DECLSPEC void mpfr_init_cache (mpfr_cache_t,
2465                                       int(*)(mpfr_ptr,mpfr_rnd_t));
2466 #endif
2467 __MPFR_DECLSPEC void mpfr_clear_cache (mpfr_cache_t);
2468 __MPFR_DECLSPEC int  mpfr_cache (mpfr_ptr, mpfr_cache_t, mpfr_rnd_t);
2469 
2470 __MPFR_DECLSPEC void mpfr_mulhigh_n (mpfr_limb_ptr, mpfr_limb_srcptr,
2471                                      mpfr_limb_srcptr, mp_size_t);
2472 __MPFR_DECLSPEC void mpfr_mullow_n  (mpfr_limb_ptr, mpfr_limb_srcptr,
2473                                      mpfr_limb_srcptr, mp_size_t);
2474 __MPFR_DECLSPEC void mpfr_sqrhigh_n (mpfr_limb_ptr, mpfr_limb_srcptr,
2475                                      mp_size_t);
2476 __MPFR_DECLSPEC mp_limb_t mpfr_divhigh_n (mpfr_limb_ptr, mpfr_limb_ptr,
2477                                           mpfr_limb_ptr, mp_size_t);
2478 
2479 __MPFR_DECLSPEC int mpfr_round_p (mp_limb_t *, mp_size_t, mpfr_exp_t,
2480                                   mpfr_prec_t);
2481 
2482 __MPFR_DECLSPEC int mpfr_round_near_x (mpfr_ptr, mpfr_srcptr, mpfr_uexp_t, int,
2483                                        mpfr_rnd_t);
2484 __MPFR_DECLSPEC MPFR_COLD_FUNCTION_ATTR MPFR_NORETURN void
2485   mpfr_abort_prec_max (void);
2486 
2487 __MPFR_DECLSPEC void mpfr_rand_raw (mpfr_limb_ptr, gmp_randstate_t,
2488                                     mpfr_prec_t);
2489 
2490 __MPFR_DECLSPEC mpz_srcptr mpfr_bernoulli_cache (unsigned long);
2491 __MPFR_DECLSPEC void mpfr_bernoulli_freecache (void);
2492 
2493 __MPFR_DECLSPEC int mpfr_sincos_fast (mpfr_ptr, mpfr_ptr, mpfr_srcptr,
2494                                       mpfr_rnd_t);
2495 
2496 __MPFR_DECLSPEC double mpfr_scale2 (double, int);
2497 
2498 __MPFR_DECLSPEC void mpfr_div_ui2 (mpfr_ptr, mpfr_srcptr, unsigned long,
2499                                    unsigned long, mpfr_rnd_t);
2500 
2501 __MPFR_DECLSPEC void mpfr_gamma_one_and_two_third (mpfr_ptr, mpfr_ptr,
2502                                                    mpfr_prec_t);
2503 
2504 __MPFR_DECLSPEC void mpfr_mpz_init (mpz_ptr);
2505 __MPFR_DECLSPEC void mpfr_mpz_init2 (mpz_ptr, mp_bitcnt_t);
2506 __MPFR_DECLSPEC void mpfr_mpz_clear (mpz_ptr);
2507 
2508 __MPFR_DECLSPEC int mpfr_odd_p (mpfr_srcptr);
2509 
2510 __MPFR_DECLSPEC int mpfr_nbits_ulong (unsigned long);
2511 #ifdef _MPFR_H_HAVE_INTMAX_T
2512 __MPFR_DECLSPEC int mpfr_nbits_uj (uintmax_t);
2513 #endif
2514 
2515 #ifdef _MPFR_H_HAVE_VA_LIST
2516 /* Declared only if <stdarg.h> has been included. */
2517 __MPFR_DECLSPEC int mpfr_vasnprintf_aux (char**, char*, size_t, const char*,
2518                                          va_list);
2519 #endif
2520 
2521 #if MPFR_WANT_ASSERT >= 2
2522 __MPFR_DECLSPEC void flags_fout (FILE *, mpfr_flags_t);
2523 #endif
2524 
2525 #if defined (__cplusplus)
2526 }
2527 #endif
2528 
2529 
2530 /*****************************************************
2531  ***************  Internal mpz_t pool  ***************
2532  *****************************************************/
2533 
2534 /* don't use the mpz_t pool with mini-gmp */
2535 #ifdef MPFR_USE_MINI_GMP
2536 # define MPFR_POOL_NENTRIES 0
2537 #endif
2538 
2539 #ifndef MPFR_POOL_NENTRIES
2540 # define MPFR_POOL_NENTRIES 32  /* default number of entries of the pool */
2541 #endif
2542 
2543 #if MPFR_POOL_NENTRIES && !defined(MPFR_POOL_DONT_REDEFINE)
2544 # undef mpz_init
2545 # undef mpz_init2
2546 # undef mpz_clear
2547 # undef mpz_init_set_ui
2548 # undef mpz_init_set
2549 # define mpz_init mpfr_mpz_init
2550 # define mpz_init2 mpfr_mpz_init2
2551 # define mpz_clear mpfr_mpz_clear
2552 # define mpz_init_set_ui(a,b) do { mpz_init (a); mpz_set_ui (a, b); } while (0)
2553 # define mpz_init_set(a,b) do { mpz_init (a); mpz_set (a, b); } while (0)
2554 #endif
2555 
2556 
2557 /******************************************************
2558  ********  Compute LOG2(LOG2(MPFR_PREC_MAX))  *********
2559  ******************************************************/
2560 
2561 #if   _MPFR_PREC_FORMAT == 1
2562 # define MPFR_PREC_MAX_TEMP USHRT_MAX
2563 #elif _MPFR_PREC_FORMAT == 2
2564 # define MPFR_PREC_MAX_TEMP UINT_MAX
2565 #elif _MPFR_PREC_FORMAT == 3
2566 # define MPFR_PREC_MAX_TEMP ULONG_MAX
2567 #else
2568 # error "Invalid MPFR Prec format"
2569 #endif
2570 
2571 /* Note: In the constants below, it is sufficient to use the suffix U.
2572  * A large enough unsigned type will be chosen automatically, but the
2573  * exact type doesn't matter here.
2574  */
2575 
2576 #if MPFR_PREC_MAX_TEMP == 255U
2577 # define MPFR_PREC_BITS 8
2578 # define MPFR_LOG2_PREC_BITS 3
2579 #elif MPFR_PREC_MAX_TEMP == 65535U
2580 # define MPFR_PREC_BITS 16
2581 # define MPFR_LOG2_PREC_BITS 4
2582 #elif MPFR_PREC_MAX_TEMP == 4294967295U
2583 # define MPFR_PREC_BITS 32
2584 # define MPFR_LOG2_PREC_BITS 5
2585 #elif MPFR_PREC_MAX_TEMP == 18446744073709551615U
2586 # define MPFR_PREC_BITS 64
2587 # define MPFR_LOG2_PREC_BITS 6
2588 #else
2589 # error "Unsupported MPFR_PREC_MAX_TEMP value"
2590 #endif
2591 
2592 
2593 /******************************************************
2594  *************  Value coverage checking  **************
2595  ******************************************************/
2596 
2597 #ifdef MPFR_COV_CHECK
2598 
2599 /* Variable names should start with the __gmpfr_cov_ prefix. */
2600 
2601 #define MPFR_COV_SET(X) (__gmpfr_cov_ ## X = 1)
2602 
2603 #if defined (__cplusplus)
2604 extern "C" {
2605 #endif
2606 
2607 __MPFR_DECLSPEC extern int __gmpfr_cov_div_ui_sb[10][2];
2608 __MPFR_DECLSPEC extern int __gmpfr_cov_sum_tmd[MPFR_RND_MAX][2][2][3][2][2];
2609 
2610 #if defined (__cplusplus)
2611 }
2612 #endif
2613 
2614 #else /* MPFR_COV_CHECK */
2615 
2616 #define MPFR_COV_SET(X) ((void) 0)
2617 
2618 #endif /* MPFR_COV_CHECK */
2619 
2620 
2621 /******************************************************
2622  *****************  Unbounded Floats  *****************
2623  ******************************************************/
2624 
2625 #if defined (__cplusplus)
2626 extern "C" {
2627 #endif
2628 
2629 /* An UBF is like a MPFR number, but with an additional mpz_t member,
2630    which is assumed to be present (with a value in it) when the usual
2631    exponent field has the value MPFR_EXP_UBF. The goal of this compatible
2632    representation is to easily be able to support UBF in "normal" code
2633    using the public API. This is some form of "subtyping".
2634 
2635    Unfortunately this breaks aliasing rules, and C does not provide any way
2636    to avoid that (except with additional syntax ugliness and API breakage,
2637    though there is a workaround -- see the end of this comment):
2638 
2639      https://news.ycombinator.com/item?id=11753236
2640 
2641    The alignment requirement for __mpfr_ubf_struct (UBF) needs to be at least
2642    as strong as the one for __mpfr_struct (MPFR number); this is not required
2643    by the C standard, but this should always be the case in practice, since
2644    __mpfr_ubf_struct starts with the same members as those of __mpfr_struct.
2645    If ever this would not be the case with some particular C implementation,
2646    an _Alignas alignment attribute (C11) could be added for UBF.
2647 
2648    When an input of a public function is an UBF, the semantic remains
2649    internal to MPFR and can change in the future. UBF arguments need
2650    to be explicitly converted to mpfr_ptr (or mpfr_srcptr); be careful
2651    with variadic functions, as the compiler will not be able to check
2652    in general. See fmma.c as an example of usage.
2653 
2654    In general, the type used for values that may be UBF must be either
2655    mpfr_ubf_t or mpfr_ubf_ptr. The type mpfr_ptr or mpfr_srcptr may be
2656    used for UBF only in the case where the pointer has been converted
2657    from mpfr_ubf_ptr, in order to ensure valid alignment. For instance,
2658    in mpfr_fmma_aux, one uses mpfr_ubf_t to generate the exact products
2659    as UBF; then the corresponding pointers are converted to mpfr_srcptr
2660    for mpfr_add (even though they point to UBF).
2661 
2662    Functions that can accept either MPFR arguments (mpfr_ptr type) or
2663    UBF arguments (mpfr_ubf_ptr type) must use a pointer type that can
2664    always be converted from both, typically mpfr_ptr or mpfr_srcptr.
2665    For instance, that's why mpfr_ubf_exp_less_p uses mpfr_srcptr.
2666    Note: "void *" could also be used, but mpfr_ptr is more meaningful
2667    and practical.
2668 
2669    Note that functions used for logging need to support UBF (currently
2670    done by printing that a number is an UBF, as it may be difficult to
2671    do more without significant changes).
2672 
2673    --------
2674 
2675    A workaround to avoid breaking aliasing rules should be to use mpfr_ptr
2676    to access the usual mpfr_t members and mpfr_ubf_ptr to access the
2677    additional member _mpfr_zexp. And never use __mpfr_ubf_struct as a
2678    declared type; otherwise this would force __mpfr_ubf_struct to be the
2679    effective type of the whole object. Thus instead of
2680 
2681      typedef __mpfr_ubf_struct mpfr_ubf_t[1];
2682 
2683    one could use the following definition as a trick to allocate an UBF as
2684    an automatic variable with the required alignment but without forcing
2685    the effective type to __mpfr_ubf_struct.
2686 
2687       typedef union {
2688         __mpfr_ubf_struct u[1];
2689         __mpfr_struct m[1];
2690       } mpfr_ubf_t;
2691 
2692    Then adapt the related code to select to right member, depending on the
2693    context. Unfortunately, this triggers -Wstrict-aliasing false positives
2694    with GCC in the MPFR_UBF_CLEAR_EXP expansion:
2695 
2696      https://gcc.gnu.org/bugzilla/show_bug.cgi?id=94337
2697 
2698    (see changeset r13820 in the ubf2 branch). So, for the time being,
2699    as long as the code does not break, do not change anything.
2700 
2701    Note: The condition "use mpfr_ptr to access the usual mpfr_t members and
2702    mpfr_ubf_ptr to access the additional member _mpfr_zexp" may be ignored
2703    if the union type is visible within the function (see ISO C99 6.5.2.3#5
2704    and 6.5.2.3#8 for the example, this implementation being very similar to
2705    the valid fragment of this example), which must be the case as the union
2706    is declared globally. However, this seems to be buggy in GCC:
2707 
2708      https://gcc.gnu.org/bugzilla/show_bug.cgi?id=14319
2709      https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65892
2710 
2711    Alternatively, GCC's may_alias attribute could conditionally be used
2712    on the __mpfr_ubf_struct and __mpfr_struct types (though it would be
2713    much stronger than needed since only these two types may alias each
2714    other).
2715 */
2716 
2717 typedef struct {
2718   mpfr_prec_t  _mpfr_prec;
2719   mpfr_sign_t  _mpfr_sign;
2720   mpfr_exp_t   _mpfr_exp;
2721   mp_limb_t   *_mpfr_d;
2722   mpz_t        _mpfr_zexp;
2723 } __mpfr_ubf_struct;
2724 
2725 typedef __mpfr_ubf_struct mpfr_ubf_t[1];
2726 typedef __mpfr_ubf_struct *mpfr_ubf_ptr;
2727 
2728 __MPFR_DECLSPEC void mpfr_ubf_mul_exact (mpfr_ubf_ptr,
2729                                          mpfr_srcptr, mpfr_srcptr);
2730 __MPFR_DECLSPEC int mpfr_ubf_exp_less_p (mpfr_srcptr, mpfr_srcptr);
2731 __MPFR_DECLSPEC mpfr_exp_t mpfr_ubf_zexp2exp (mpz_ptr);
2732 __MPFR_DECLSPEC mpfr_exp_t mpfr_ubf_diff_exp (mpfr_srcptr, mpfr_srcptr);
2733 
2734 #if defined (__cplusplus)
2735 }
2736 #endif
2737 
2738 /* Get the _mpfr_zexp field (pointer to a mpz_t) of an UBF object.
2739    For practical reasons, the type of the argument x can be either
2740    mpfr_ubf_ptr or mpfr_ptr, since the latter is used in functions
2741    that accept both MPFR numbers and UBF's; this is checked by the
2742    code "(x)->_mpfr_exp" (the "sizeof" prevents an access, which
2743    could be invalid when MPFR_ZEXP(x) is used for an assignment,
2744    and also avoids breaking the aliasing rules if they are dealt
2745    with in the future).
2746    This macro can be used when building an UBF. So we do not check
2747    that the _mpfr_exp field has the value MPFR_EXP_UBF. */
2748 #define MPFR_ZEXP(x)                            \
2749   ((void) sizeof ((x)->_mpfr_exp),              \
2750    ((mpfr_ubf_ptr) (x))->_mpfr_zexp)
2751 
2752 /* If x is an UBF, clear its mpz_t exponent. */
2753 #define MPFR_UBF_CLEAR_EXP(x) \
2754   ((void) (MPFR_IS_UBF (x) && (mpz_clear (MPFR_ZEXP (x)), 0)))
2755 
2756 /* Like MPFR_GET_EXP, but accepts UBF (with exponent saturated to
2757    the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]). */
2758 #define MPFR_UBF_GET_EXP(x)                                     \
2759   (MPFR_IS_UBF (x) ? mpfr_ubf_zexp2exp (MPFR_ZEXP (x)) :        \
2760    MPFR_GET_EXP ((mpfr_ptr) (x)))
2761 
2762 #endif /* __MPFR_IMPL_H__ */
2763