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