1/* Utilities for MPFR developers, not exported. 2 3Copyright 1999-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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) 229extern "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. */ 242struct __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}; 249typedef struct __gmpfr_cache_s mpfr_cache_t[1]; 250typedef 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) 257extern MPFR_THREAD_ATTR mpfr_flags_t __gmpfr_flags; 258extern MPFR_THREAD_ATTR mpfr_exp_t __gmpfr_emin; 259extern MPFR_THREAD_ATTR mpfr_exp_t __gmpfr_emax; 260extern MPFR_THREAD_ATTR mpfr_prec_t __gmpfr_default_fp_bit_precision; 261extern MPFR_THREAD_ATTR mpfr_rnd_t __gmpfr_default_rounding_mode; 262extern MPFR_CACHE_ATTR mpfr_cache_t __gmpfr_cache_const_euler; 263extern MPFR_CACHE_ATTR mpfr_cache_t __gmpfr_cache_const_catalan; 264# ifndef MPFR_USE_LOGGING 265extern MPFR_CACHE_ATTR mpfr_cache_t __gmpfr_cache_const_pi; 266extern 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. */ 272extern MPFR_CACHE_ATTR mpfr_cache_t __gmpfr_normal_pi; 273extern MPFR_CACHE_ATTR mpfr_cache_t __gmpfr_normal_log2; 274extern MPFR_CACHE_ATTR mpfr_cache_t __gmpfr_logging_pi; 275extern MPFR_CACHE_ATTR mpfr_cache_t __gmpfr_logging_log2; 276extern MPFR_CACHE_ATTR mpfr_cache_ptr __gmpfr_cache_const_pi; 277extern 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__) 648static 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 813typedef 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. */ 894union 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. */ 904union 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 1003typedef long mpfr_eexp_t; 1004typedef 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 1009typedef intmax_t mpfr_eexp_t; 1010typedef 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*/ 1321typedef 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 */ 1453asm (".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) \ 1586do { \ 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 1736typedef 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) 2161extern "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 2256struct 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) 2385extern "C" { 2386#endif 2387 2388MPFR_COLD_FUNCTION_ATTR __MPFR_DECLSPEC int 2389 mpfr_underflow (mpfr_ptr, mpfr_rnd_t, int); 2390MPFR_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) 2604extern "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) 2626extern "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 2717typedef 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 2725typedef __mpfr_ubf_struct mpfr_ubf_t[1]; 2726typedef __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