xf.li | bfc6e71 | 2025-02-07 01:54:34 -0800 | [diff] [blame^] | 1 | /* |
| 2 | * ==================================================== |
| 3 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 4 | * |
| 5 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
| 6 | * Permission to use, copy, modify, and distribute this |
| 7 | * software is freely granted, provided that this notice |
| 8 | * is preserved. |
| 9 | * ==================================================== |
| 10 | */ |
| 11 | |
| 12 | /* |
| 13 | * from: @(#)fdlibm.h 5.1 93/09/24 |
| 14 | */ |
| 15 | |
| 16 | #ifndef _MATH_PRIVATE_H_ |
| 17 | #define _MATH_PRIVATE_H_ |
| 18 | |
| 19 | #include <endian.h> |
| 20 | #include <stdint.h> |
| 21 | #include <sys/types.h> |
| 22 | #include <fenv.h> |
| 23 | #include <float.h> |
| 24 | #include <get-rounding-mode.h> |
| 25 | |
| 26 | /* The original fdlibm code used statements like: |
| 27 | n0 = ((*(int*)&one)>>29)^1; * index of high word * |
| 28 | ix0 = *(n0+(int*)&x); * high word of x * |
| 29 | ix1 = *((1-n0)+(int*)&x); * low word of x * |
| 30 | to dig two 32 bit words out of the 64 bit IEEE floating point |
| 31 | value. That is non-ANSI, and, moreover, the gcc instruction |
| 32 | scheduler gets it wrong. We instead use the following macros. |
| 33 | Unlike the original code, we determine the endianness at compile |
| 34 | time, not at run time; I don't see much benefit to selecting |
| 35 | endianness at run time. */ |
| 36 | |
| 37 | /* A union which permits us to convert between a double and two 32 bit |
| 38 | ints. */ |
| 39 | |
| 40 | #if __FLOAT_WORD_ORDER == BIG_ENDIAN |
| 41 | |
| 42 | typedef union |
| 43 | { |
| 44 | double value; |
| 45 | struct |
| 46 | { |
| 47 | u_int32_t msw; |
| 48 | u_int32_t lsw; |
| 49 | } parts; |
| 50 | uint64_t word; |
| 51 | } ieee_double_shape_type; |
| 52 | |
| 53 | #endif |
| 54 | |
| 55 | #if __FLOAT_WORD_ORDER == LITTLE_ENDIAN |
| 56 | |
| 57 | typedef union |
| 58 | { |
| 59 | double value; |
| 60 | struct |
| 61 | { |
| 62 | u_int32_t lsw; |
| 63 | u_int32_t msw; |
| 64 | } parts; |
| 65 | uint64_t word; |
| 66 | } ieee_double_shape_type; |
| 67 | |
| 68 | #endif |
| 69 | |
| 70 | /* Get two 32 bit ints from a double. */ |
| 71 | |
| 72 | #define EXTRACT_WORDS(ix0,ix1,d) \ |
| 73 | do { \ |
| 74 | ieee_double_shape_type ew_u; \ |
| 75 | ew_u.value = (d); \ |
| 76 | (ix0) = ew_u.parts.msw; \ |
| 77 | (ix1) = ew_u.parts.lsw; \ |
| 78 | } while (0) |
| 79 | |
| 80 | /* Get the more significant 32 bit int from a double. */ |
| 81 | |
| 82 | #ifndef GET_HIGH_WORD |
| 83 | # define GET_HIGH_WORD(i,d) \ |
| 84 | do { \ |
| 85 | ieee_double_shape_type gh_u; \ |
| 86 | gh_u.value = (d); \ |
| 87 | (i) = gh_u.parts.msw; \ |
| 88 | } while (0) |
| 89 | #endif |
| 90 | |
| 91 | /* Get the less significant 32 bit int from a double. */ |
| 92 | |
| 93 | #ifndef GET_LOW_WORD |
| 94 | # define GET_LOW_WORD(i,d) \ |
| 95 | do { \ |
| 96 | ieee_double_shape_type gl_u; \ |
| 97 | gl_u.value = (d); \ |
| 98 | (i) = gl_u.parts.lsw; \ |
| 99 | } while (0) |
| 100 | #endif |
| 101 | |
| 102 | /* Get all in one, efficient on 64-bit machines. */ |
| 103 | #ifndef EXTRACT_WORDS64 |
| 104 | # define EXTRACT_WORDS64(i,d) \ |
| 105 | do { \ |
| 106 | ieee_double_shape_type gh_u; \ |
| 107 | gh_u.value = (d); \ |
| 108 | (i) = gh_u.word; \ |
| 109 | } while (0) |
| 110 | #endif |
| 111 | |
| 112 | /* Set a double from two 32 bit ints. */ |
| 113 | #ifndef INSERT_WORDS |
| 114 | # define INSERT_WORDS(d,ix0,ix1) \ |
| 115 | do { \ |
| 116 | ieee_double_shape_type iw_u; \ |
| 117 | iw_u.parts.msw = (ix0); \ |
| 118 | iw_u.parts.lsw = (ix1); \ |
| 119 | (d) = iw_u.value; \ |
| 120 | } while (0) |
| 121 | #endif |
| 122 | |
| 123 | /* Get all in one, efficient on 64-bit machines. */ |
| 124 | #ifndef INSERT_WORDS64 |
| 125 | # define INSERT_WORDS64(d,i) \ |
| 126 | do { \ |
| 127 | ieee_double_shape_type iw_u; \ |
| 128 | iw_u.word = (i); \ |
| 129 | (d) = iw_u.value; \ |
| 130 | } while (0) |
| 131 | #endif |
| 132 | |
| 133 | /* Set the more significant 32 bits of a double from an int. */ |
| 134 | #ifndef SET_HIGH_WORD |
| 135 | #define SET_HIGH_WORD(d,v) \ |
| 136 | do { \ |
| 137 | ieee_double_shape_type sh_u; \ |
| 138 | sh_u.value = (d); \ |
| 139 | sh_u.parts.msw = (v); \ |
| 140 | (d) = sh_u.value; \ |
| 141 | } while (0) |
| 142 | #endif |
| 143 | |
| 144 | /* Set the less significant 32 bits of a double from an int. */ |
| 145 | #ifndef SET_LOW_WORD |
| 146 | # define SET_LOW_WORD(d,v) \ |
| 147 | do { \ |
| 148 | ieee_double_shape_type sl_u; \ |
| 149 | sl_u.value = (d); \ |
| 150 | sl_u.parts.lsw = (v); \ |
| 151 | (d) = sl_u.value; \ |
| 152 | } while (0) |
| 153 | #endif |
| 154 | |
| 155 | /* A union which permits us to convert between a float and a 32 bit |
| 156 | int. */ |
| 157 | |
| 158 | typedef union |
| 159 | { |
| 160 | float value; |
| 161 | u_int32_t word; |
| 162 | } ieee_float_shape_type; |
| 163 | |
| 164 | /* Get a 32 bit int from a float. */ |
| 165 | #ifndef GET_FLOAT_WORD |
| 166 | # define GET_FLOAT_WORD(i,d) \ |
| 167 | do { \ |
| 168 | ieee_float_shape_type gf_u; \ |
| 169 | gf_u.value = (d); \ |
| 170 | (i) = gf_u.word; \ |
| 171 | } while (0) |
| 172 | #endif |
| 173 | |
| 174 | /* Set a float from a 32 bit int. */ |
| 175 | #ifndef SET_FLOAT_WORD |
| 176 | # define SET_FLOAT_WORD(d,i) \ |
| 177 | do { \ |
| 178 | ieee_float_shape_type sf_u; \ |
| 179 | sf_u.word = (i); \ |
| 180 | (d) = sf_u.value; \ |
| 181 | } while (0) |
| 182 | #endif |
| 183 | |
| 184 | /* Get long double macros from a separate header. */ |
| 185 | #include <math_ldbl.h> |
| 186 | |
| 187 | /* ieee style elementary functions */ |
| 188 | extern double __ieee754_sqrt (double); |
| 189 | extern double __ieee754_acos (double); |
| 190 | extern double __ieee754_acosh (double); |
| 191 | extern double __ieee754_log (double); |
| 192 | extern double __ieee754_atanh (double); |
| 193 | extern double __ieee754_asin (double); |
| 194 | extern double __ieee754_atan2 (double,double); |
| 195 | extern double __ieee754_exp (double); |
| 196 | extern double __ieee754_exp2 (double); |
| 197 | extern double __ieee754_exp10 (double); |
| 198 | extern double __ieee754_cosh (double); |
| 199 | extern double __ieee754_fmod (double,double); |
| 200 | extern double __ieee754_pow (double,double); |
| 201 | extern double __ieee754_lgamma_r (double,int *); |
| 202 | extern double __ieee754_gamma_r (double,int *); |
| 203 | extern double __ieee754_lgamma (double); |
| 204 | extern double __ieee754_gamma (double); |
| 205 | extern double __ieee754_log10 (double); |
| 206 | extern double __ieee754_log2 (double); |
| 207 | extern double __ieee754_sinh (double); |
| 208 | extern double __ieee754_hypot (double,double); |
| 209 | extern double __ieee754_j0 (double); |
| 210 | extern double __ieee754_j1 (double); |
| 211 | extern double __ieee754_y0 (double); |
| 212 | extern double __ieee754_y1 (double); |
| 213 | extern double __ieee754_jn (int,double); |
| 214 | extern double __ieee754_yn (int,double); |
| 215 | extern double __ieee754_remainder (double,double); |
| 216 | extern int32_t __ieee754_rem_pio2 (double,double*); |
| 217 | extern double __ieee754_scalb (double,double); |
| 218 | extern int __ieee754_ilogb (double); |
| 219 | |
| 220 | /* fdlibm kernel function */ |
| 221 | extern double __kernel_standard (double,double,int); |
| 222 | extern float __kernel_standard_f (float,float,int); |
| 223 | extern long double __kernel_standard_l (long double,long double,int); |
| 224 | extern double __kernel_sin (double,double,int); |
| 225 | extern double __kernel_cos (double,double); |
| 226 | extern double __kernel_tan (double,double,int); |
| 227 | extern int __kernel_rem_pio2 (double*,double*,int,int,int, const int32_t*); |
| 228 | |
| 229 | /* internal functions. */ |
| 230 | extern double __copysign (double x, double __y); |
| 231 | |
| 232 | extern inline double __copysign (double x, double y) |
| 233 | { return __builtin_copysign (x, y); } |
| 234 | |
| 235 | /* ieee style elementary float functions */ |
| 236 | extern float __ieee754_sqrtf (float); |
| 237 | extern float __ieee754_acosf (float); |
| 238 | extern float __ieee754_acoshf (float); |
| 239 | extern float __ieee754_logf (float); |
| 240 | extern float __ieee754_atanhf (float); |
| 241 | extern float __ieee754_asinf (float); |
| 242 | extern float __ieee754_atan2f (float,float); |
| 243 | extern float __ieee754_expf (float); |
| 244 | extern float __ieee754_exp2f (float); |
| 245 | extern float __ieee754_exp10f (float); |
| 246 | extern float __ieee754_coshf (float); |
| 247 | extern float __ieee754_fmodf (float,float); |
| 248 | extern float __ieee754_powf (float,float); |
| 249 | extern float __ieee754_lgammaf_r (float,int *); |
| 250 | extern float __ieee754_gammaf_r (float,int *); |
| 251 | extern float __ieee754_lgammaf (float); |
| 252 | extern float __ieee754_gammaf (float); |
| 253 | extern float __ieee754_log10f (float); |
| 254 | extern float __ieee754_log2f (float); |
| 255 | extern float __ieee754_sinhf (float); |
| 256 | extern float __ieee754_hypotf (float,float); |
| 257 | extern float __ieee754_j0f (float); |
| 258 | extern float __ieee754_j1f (float); |
| 259 | extern float __ieee754_y0f (float); |
| 260 | extern float __ieee754_y1f (float); |
| 261 | extern float __ieee754_jnf (int,float); |
| 262 | extern float __ieee754_ynf (int,float); |
| 263 | extern float __ieee754_remainderf (float,float); |
| 264 | extern int32_t __ieee754_rem_pio2f (float,float*); |
| 265 | extern float __ieee754_scalbf (float,float); |
| 266 | extern int __ieee754_ilogbf (float); |
| 267 | |
| 268 | |
| 269 | /* float versions of fdlibm kernel functions */ |
| 270 | extern float __kernel_sinf (float,float,int); |
| 271 | extern float __kernel_cosf (float,float); |
| 272 | extern float __kernel_tanf (float,float,int); |
| 273 | extern int __kernel_rem_pio2f (float*,float*,int,int,int, const int32_t*); |
| 274 | |
| 275 | /* internal functions. */ |
| 276 | extern float __copysignf (float x, float __y); |
| 277 | |
| 278 | extern inline float __copysignf (float x, float y) |
| 279 | { return __builtin_copysignf (x, y); } |
| 280 | |
| 281 | /* ieee style elementary long double functions */ |
| 282 | extern long double __ieee754_sqrtl (long double); |
| 283 | extern long double __ieee754_acosl (long double); |
| 284 | extern long double __ieee754_acoshl (long double); |
| 285 | extern long double __ieee754_logl (long double); |
| 286 | extern long double __ieee754_atanhl (long double); |
| 287 | extern long double __ieee754_asinl (long double); |
| 288 | extern long double __ieee754_atan2l (long double,long double); |
| 289 | extern long double __ieee754_expl (long double); |
| 290 | extern long double __ieee754_exp2l (long double); |
| 291 | extern long double __ieee754_exp10l (long double); |
| 292 | extern long double __ieee754_coshl (long double); |
| 293 | extern long double __ieee754_fmodl (long double,long double); |
| 294 | extern long double __ieee754_powl (long double,long double); |
| 295 | extern long double __ieee754_lgammal_r (long double,int *); |
| 296 | extern long double __ieee754_gammal_r (long double,int *); |
| 297 | extern long double __ieee754_lgammal (long double); |
| 298 | extern long double __ieee754_gammal (long double); |
| 299 | extern long double __ieee754_log10l (long double); |
| 300 | extern long double __ieee754_log2l (long double); |
| 301 | extern long double __ieee754_sinhl (long double); |
| 302 | extern long double __ieee754_hypotl (long double,long double); |
| 303 | extern long double __ieee754_j0l (long double); |
| 304 | extern long double __ieee754_j1l (long double); |
| 305 | extern long double __ieee754_y0l (long double); |
| 306 | extern long double __ieee754_y1l (long double); |
| 307 | extern long double __ieee754_jnl (int,long double); |
| 308 | extern long double __ieee754_ynl (int,long double); |
| 309 | extern long double __ieee754_remainderl (long double,long double); |
| 310 | extern int __ieee754_rem_pio2l (long double,long double*); |
| 311 | extern long double __ieee754_scalbl (long double,long double); |
| 312 | extern int __ieee754_ilogbl (long double); |
| 313 | |
| 314 | /* long double versions of fdlibm kernel functions */ |
| 315 | extern long double __kernel_sinl (long double,long double,int); |
| 316 | extern long double __kernel_cosl (long double,long double); |
| 317 | extern long double __kernel_tanl (long double,long double,int); |
| 318 | extern void __kernel_sincosl (long double,long double, |
| 319 | long double *,long double *, int); |
| 320 | extern int __kernel_rem_pio2l (long double*,long double*,int,int, |
| 321 | int,const int*); |
| 322 | |
| 323 | #ifndef NO_LONG_DOUBLE |
| 324 | /* prototypes required to compile the ldbl-96 support without warnings */ |
| 325 | extern int __finitel (long double); |
| 326 | extern int __ilogbl (long double); |
| 327 | extern int __isinfl (long double); |
| 328 | extern int __isnanl (long double); |
| 329 | extern long double __atanl (long double); |
| 330 | extern long double __copysignl (long double, long double); |
| 331 | extern long double __expm1l (long double); |
| 332 | extern long double __floorl (long double); |
| 333 | extern long double __frexpl (long double, int *); |
| 334 | extern long double __ldexpl (long double, int); |
| 335 | extern long double __log1pl (long double); |
| 336 | extern long double __nanl (const char *); |
| 337 | extern long double __rintl (long double); |
| 338 | extern long double __scalbnl (long double, int); |
| 339 | extern long double __sqrtl (long double x); |
| 340 | extern long double fabsl (long double x); |
| 341 | extern void __sincosl (long double, long double *, long double *); |
| 342 | extern long double __logbl (long double x); |
| 343 | extern long double __significandl (long double x); |
| 344 | |
| 345 | extern inline long double __copysignl (long double x, long double y) |
| 346 | { return __builtin_copysignl (x, y); } |
| 347 | |
| 348 | #endif |
| 349 | |
| 350 | /* Prototypes for functions of the IBM Accurate Mathematical Library. */ |
| 351 | extern double __exp1 (double __x, double __xx, double __error); |
| 352 | extern double __sin (double __x); |
| 353 | extern double __cos (double __x); |
| 354 | extern int __branred (double __x, double *__a, double *__aa); |
| 355 | extern void __doasin (double __x, double __dx, double __v[]); |
| 356 | extern void __dubsin (double __x, double __dx, double __v[]); |
| 357 | extern void __dubcos (double __x, double __dx, double __v[]); |
| 358 | extern double __halfulp (double __x, double __y); |
| 359 | extern double __sin32 (double __x, double __res, double __res1); |
| 360 | extern double __cos32 (double __x, double __res, double __res1); |
| 361 | extern double __mpsin (double __x, double __dx, bool __range_reduce); |
| 362 | extern double __mpcos (double __x, double __dx, bool __range_reduce); |
| 363 | extern double __slowexp (double __x); |
| 364 | extern double __slowpow (double __x, double __y, double __z); |
| 365 | extern void __docos (double __x, double __dx, double __v[]); |
| 366 | |
| 367 | /* Return X^2 + Y^2 - 1, computed without large cancellation error. |
| 368 | It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >= |
| 369 | 0.5. */ |
| 370 | extern float __x2y2m1f (float x, float y); |
| 371 | extern double __x2y2m1 (double x, double y); |
| 372 | extern long double __x2y2m1l (long double x, long double y); |
| 373 | |
| 374 | /* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N |
| 375 | - 1, in the form R * (1 + *EPS) where the return value R is an |
| 376 | approximation to the product and *EPS is set to indicate the |
| 377 | approximate error in the return value. X is such that all the |
| 378 | values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / |
| 379 | X is small enough that factors quadratic in it can be |
| 380 | neglected. */ |
| 381 | extern float __gamma_productf (float x, float x_eps, int n, float *eps); |
| 382 | extern double __gamma_product (double x, double x_eps, int n, double *eps); |
| 383 | extern long double __gamma_productl (long double x, long double x_eps, |
| 384 | int n, long double *eps); |
| 385 | |
| 386 | /* Compute lgamma of a negative argument X, if it is in a range |
| 387 | (depending on the floating-point format) for which expansion around |
| 388 | zeros is used, setting *SIGNGAMP accordingly. */ |
| 389 | extern float __lgamma_negf (float x, int *signgamp); |
| 390 | extern double __lgamma_neg (double x, int *signgamp); |
| 391 | extern long double __lgamma_negl (long double x, int *signgamp); |
| 392 | |
| 393 | /* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS + |
| 394 | 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that |
| 395 | all the values X + 1, ..., X + N - 1 are exactly representable, and |
| 396 | X_EPS / X is small enough that factors quadratic in it can be |
| 397 | neglected. */ |
| 398 | extern double __lgamma_product (double t, double x, double x_eps, int n); |
| 399 | extern long double __lgamma_productl (long double t, long double x, |
| 400 | long double x_eps, int n); |
| 401 | |
| 402 | #ifndef math_opt_barrier |
| 403 | # define math_opt_barrier(x) \ |
| 404 | ({ __typeof (x) __x = (x); __asm ("" : "+m" (__x)); __x; }) |
| 405 | # define math_force_eval(x) \ |
| 406 | ({ __typeof (x) __x = (x); __asm __volatile__ ("" : : "m" (__x)); }) |
| 407 | #endif |
| 408 | |
| 409 | /* math_narrow_eval reduces its floating-point argument to the range |
| 410 | and precision of its semantic type. (The original evaluation may |
| 411 | still occur with excess range and precision, so the result may be |
| 412 | affected by double rounding.) */ |
| 413 | #if FLT_EVAL_METHOD == 0 |
| 414 | # define math_narrow_eval(x) (x) |
| 415 | #else |
| 416 | # if FLT_EVAL_METHOD == 1 |
| 417 | # define excess_precision(type) __builtin_types_compatible_p (type, float) |
| 418 | # else |
| 419 | # define excess_precision(type) (__builtin_types_compatible_p (type, float) \ |
| 420 | || __builtin_types_compatible_p (type, \ |
| 421 | double)) |
| 422 | # endif |
| 423 | # define math_narrow_eval(x) \ |
| 424 | ({ \ |
| 425 | __typeof (x) math_narrow_eval_tmp = (x); \ |
| 426 | if (excess_precision (__typeof (math_narrow_eval_tmp))) \ |
| 427 | __asm__ ("" : "+m" (math_narrow_eval_tmp)); \ |
| 428 | math_narrow_eval_tmp; \ |
| 429 | }) |
| 430 | #endif |
| 431 | |
| 432 | #define fabs_tg(x) __builtin_choose_expr \ |
| 433 | (__builtin_types_compatible_p (__typeof (x), float), \ |
| 434 | __builtin_fabsf (x), \ |
| 435 | __builtin_choose_expr \ |
| 436 | (__builtin_types_compatible_p (__typeof (x), double), \ |
| 437 | __builtin_fabs (x), __builtin_fabsl (x))) |
| 438 | #define min_of_type(type) __builtin_choose_expr \ |
| 439 | (__builtin_types_compatible_p (type, float), \ |
| 440 | FLT_MIN, \ |
| 441 | __builtin_choose_expr \ |
| 442 | (__builtin_types_compatible_p (type, double), \ |
| 443 | DBL_MIN, LDBL_MIN)) |
| 444 | |
| 445 | /* If X (which is not a NaN) is subnormal, force an underflow |
| 446 | exception. */ |
| 447 | #define math_check_force_underflow(x) \ |
| 448 | do \ |
| 449 | { \ |
| 450 | __typeof (x) force_underflow_tmp = (x); \ |
| 451 | if (fabs_tg (force_underflow_tmp) \ |
| 452 | < min_of_type (__typeof (force_underflow_tmp))) \ |
| 453 | { \ |
| 454 | __typeof (force_underflow_tmp) force_underflow_tmp2 \ |
| 455 | = force_underflow_tmp * force_underflow_tmp; \ |
| 456 | math_force_eval (force_underflow_tmp2); \ |
| 457 | } \ |
| 458 | } \ |
| 459 | while (0) |
| 460 | /* Likewise, but X is also known to be nonnegative. */ |
| 461 | #define math_check_force_underflow_nonneg(x) \ |
| 462 | do \ |
| 463 | { \ |
| 464 | __typeof (x) force_underflow_tmp = (x); \ |
| 465 | if (force_underflow_tmp \ |
| 466 | < min_of_type (__typeof (force_underflow_tmp))) \ |
| 467 | { \ |
| 468 | __typeof (force_underflow_tmp) force_underflow_tmp2 \ |
| 469 | = force_underflow_tmp * force_underflow_tmp; \ |
| 470 | math_force_eval (force_underflow_tmp2); \ |
| 471 | } \ |
| 472 | } \ |
| 473 | while (0) |
| 474 | /* Likewise, for both real and imaginary parts of a complex |
| 475 | result. */ |
| 476 | #define math_check_force_underflow_complex(x) \ |
| 477 | do \ |
| 478 | { \ |
| 479 | __typeof (x) force_underflow_complex_tmp = (x); \ |
| 480 | math_check_force_underflow (__real__ force_underflow_complex_tmp); \ |
| 481 | math_check_force_underflow (__imag__ force_underflow_complex_tmp); \ |
| 482 | } \ |
| 483 | while (0) |
| 484 | |
| 485 | /* The standards only specify one variant of the fenv.h interfaces. |
| 486 | But at least for some architectures we can be more efficient if we |
| 487 | know what operations are going to be performed. Therefore we |
| 488 | define additional interfaces. By default they refer to the normal |
| 489 | interfaces. */ |
| 490 | |
| 491 | static __always_inline void |
| 492 | default_libc_feholdexcept (fenv_t *e) |
| 493 | { |
| 494 | (void) __feholdexcept (e); |
| 495 | } |
| 496 | |
| 497 | #ifndef libc_feholdexcept |
| 498 | # define libc_feholdexcept default_libc_feholdexcept |
| 499 | #endif |
| 500 | #ifndef libc_feholdexceptf |
| 501 | # define libc_feholdexceptf default_libc_feholdexcept |
| 502 | #endif |
| 503 | #ifndef libc_feholdexceptl |
| 504 | # define libc_feholdexceptl default_libc_feholdexcept |
| 505 | #endif |
| 506 | |
| 507 | static __always_inline void |
| 508 | default_libc_fesetround (int r) |
| 509 | { |
| 510 | (void) __fesetround (r); |
| 511 | } |
| 512 | |
| 513 | #ifndef libc_fesetround |
| 514 | # define libc_fesetround default_libc_fesetround |
| 515 | #endif |
| 516 | #ifndef libc_fesetroundf |
| 517 | # define libc_fesetroundf default_libc_fesetround |
| 518 | #endif |
| 519 | #ifndef libc_fesetroundl |
| 520 | # define libc_fesetroundl default_libc_fesetround |
| 521 | #endif |
| 522 | |
| 523 | static __always_inline void |
| 524 | default_libc_feholdexcept_setround (fenv_t *e, int r) |
| 525 | { |
| 526 | __feholdexcept (e); |
| 527 | __fesetround (r); |
| 528 | } |
| 529 | |
| 530 | #ifndef libc_feholdexcept_setround |
| 531 | # define libc_feholdexcept_setround default_libc_feholdexcept_setround |
| 532 | #endif |
| 533 | #ifndef libc_feholdexcept_setroundf |
| 534 | # define libc_feholdexcept_setroundf default_libc_feholdexcept_setround |
| 535 | #endif |
| 536 | #ifndef libc_feholdexcept_setroundl |
| 537 | # define libc_feholdexcept_setroundl default_libc_feholdexcept_setround |
| 538 | #endif |
| 539 | |
| 540 | #ifndef libc_feholdsetround_53bit |
| 541 | # define libc_feholdsetround_53bit libc_feholdsetround |
| 542 | #endif |
| 543 | |
| 544 | #ifndef libc_fetestexcept |
| 545 | # define libc_fetestexcept fetestexcept |
| 546 | #endif |
| 547 | #ifndef libc_fetestexceptf |
| 548 | # define libc_fetestexceptf fetestexcept |
| 549 | #endif |
| 550 | #ifndef libc_fetestexceptl |
| 551 | # define libc_fetestexceptl fetestexcept |
| 552 | #endif |
| 553 | |
| 554 | static __always_inline void |
| 555 | default_libc_fesetenv (fenv_t *e) |
| 556 | { |
| 557 | (void) __fesetenv (e); |
| 558 | } |
| 559 | |
| 560 | #ifndef libc_fesetenv |
| 561 | # define libc_fesetenv default_libc_fesetenv |
| 562 | #endif |
| 563 | #ifndef libc_fesetenvf |
| 564 | # define libc_fesetenvf default_libc_fesetenv |
| 565 | #endif |
| 566 | #ifndef libc_fesetenvl |
| 567 | # define libc_fesetenvl default_libc_fesetenv |
| 568 | #endif |
| 569 | |
| 570 | static __always_inline void |
| 571 | default_libc_feupdateenv (fenv_t *e) |
| 572 | { |
| 573 | (void) __feupdateenv (e); |
| 574 | } |
| 575 | |
| 576 | #ifndef libc_feupdateenv |
| 577 | # define libc_feupdateenv default_libc_feupdateenv |
| 578 | #endif |
| 579 | #ifndef libc_feupdateenvf |
| 580 | # define libc_feupdateenvf default_libc_feupdateenv |
| 581 | #endif |
| 582 | #ifndef libc_feupdateenvl |
| 583 | # define libc_feupdateenvl default_libc_feupdateenv |
| 584 | #endif |
| 585 | |
| 586 | #ifndef libc_feresetround_53bit |
| 587 | # define libc_feresetround_53bit libc_feresetround |
| 588 | #endif |
| 589 | |
| 590 | static __always_inline int |
| 591 | default_libc_feupdateenv_test (fenv_t *e, int ex) |
| 592 | { |
| 593 | int ret = fetestexcept (ex); |
| 594 | __feupdateenv (e); |
| 595 | return ret; |
| 596 | } |
| 597 | |
| 598 | #ifndef libc_feupdateenv_test |
| 599 | # define libc_feupdateenv_test default_libc_feupdateenv_test |
| 600 | #endif |
| 601 | #ifndef libc_feupdateenv_testf |
| 602 | # define libc_feupdateenv_testf default_libc_feupdateenv_test |
| 603 | #endif |
| 604 | #ifndef libc_feupdateenv_testl |
| 605 | # define libc_feupdateenv_testl default_libc_feupdateenv_test |
| 606 | #endif |
| 607 | |
| 608 | /* Save and set the rounding mode. The use of fenv_t to store the old mode |
| 609 | allows a target-specific version of this function to avoid converting the |
| 610 | rounding mode from the fpu format. By default we have no choice but to |
| 611 | manipulate the entire env. */ |
| 612 | |
| 613 | #ifndef libc_feholdsetround |
| 614 | # define libc_feholdsetround libc_feholdexcept_setround |
| 615 | #endif |
| 616 | #ifndef libc_feholdsetroundf |
| 617 | # define libc_feholdsetroundf libc_feholdexcept_setroundf |
| 618 | #endif |
| 619 | #ifndef libc_feholdsetroundl |
| 620 | # define libc_feholdsetroundl libc_feholdexcept_setroundl |
| 621 | #endif |
| 622 | |
| 623 | /* ... and the reverse. */ |
| 624 | |
| 625 | #ifndef libc_feresetround |
| 626 | # define libc_feresetround libc_feupdateenv |
| 627 | #endif |
| 628 | #ifndef libc_feresetroundf |
| 629 | # define libc_feresetroundf libc_feupdateenvf |
| 630 | #endif |
| 631 | #ifndef libc_feresetroundl |
| 632 | # define libc_feresetroundl libc_feupdateenvl |
| 633 | #endif |
| 634 | |
| 635 | /* ... and a version that may also discard exceptions. */ |
| 636 | |
| 637 | #ifndef libc_feresetround_noex |
| 638 | # define libc_feresetround_noex libc_fesetenv |
| 639 | #endif |
| 640 | #ifndef libc_feresetround_noexf |
| 641 | # define libc_feresetround_noexf libc_fesetenvf |
| 642 | #endif |
| 643 | #ifndef libc_feresetround_noexl |
| 644 | # define libc_feresetround_noexl libc_fesetenvl |
| 645 | #endif |
| 646 | |
| 647 | #ifndef HAVE_RM_CTX |
| 648 | # define HAVE_RM_CTX 0 |
| 649 | #endif |
| 650 | |
| 651 | #if HAVE_RM_CTX |
| 652 | /* Set/Restore Rounding Modes only when necessary. If defined, these functions |
| 653 | set/restore floating point state only if the state needed within the lexical |
| 654 | block is different from the current state. This saves a lot of time when |
| 655 | the floating point unit is much slower than the fixed point units. */ |
| 656 | |
| 657 | # ifndef libc_feholdsetround_noex_ctx |
| 658 | # define libc_feholdsetround_noex_ctx libc_feholdsetround_ctx |
| 659 | # endif |
| 660 | # ifndef libc_feholdsetround_noexf_ctx |
| 661 | # define libc_feholdsetround_noexf_ctx libc_feholdsetroundf_ctx |
| 662 | # endif |
| 663 | # ifndef libc_feholdsetround_noexl_ctx |
| 664 | # define libc_feholdsetround_noexl_ctx libc_feholdsetroundl_ctx |
| 665 | # endif |
| 666 | |
| 667 | # ifndef libc_feresetround_noex_ctx |
| 668 | # define libc_feresetround_noex_ctx libc_fesetenv_ctx |
| 669 | # endif |
| 670 | # ifndef libc_feresetround_noexf_ctx |
| 671 | # define libc_feresetround_noexf_ctx libc_fesetenvf_ctx |
| 672 | # endif |
| 673 | # ifndef libc_feresetround_noexl_ctx |
| 674 | # define libc_feresetround_noexl_ctx libc_fesetenvl_ctx |
| 675 | # endif |
| 676 | |
| 677 | #else |
| 678 | |
| 679 | /* Default implementation using standard fenv functions. |
| 680 | Avoid unnecessary rounding mode changes by first checking the |
| 681 | current rounding mode. Note the use of __glibc_unlikely is |
| 682 | important for performance. */ |
| 683 | |
| 684 | static __always_inline void |
| 685 | libc_feholdsetround_ctx (struct rm_ctx *ctx, int round) |
| 686 | { |
| 687 | ctx->updated_status = false; |
| 688 | |
| 689 | /* Update rounding mode only if different. */ |
| 690 | if (__glibc_unlikely (round != get_rounding_mode ())) |
| 691 | { |
| 692 | ctx->updated_status = true; |
| 693 | __fegetenv (&ctx->env); |
| 694 | __fesetround (round); |
| 695 | } |
| 696 | } |
| 697 | |
| 698 | static __always_inline void |
| 699 | libc_feresetround_ctx (struct rm_ctx *ctx) |
| 700 | { |
| 701 | /* Restore the rounding mode if updated. */ |
| 702 | if (__glibc_unlikely (ctx->updated_status)) |
| 703 | __feupdateenv (&ctx->env); |
| 704 | } |
| 705 | |
| 706 | static __always_inline void |
| 707 | libc_feholdsetround_noex_ctx (struct rm_ctx *ctx, int round) |
| 708 | { |
| 709 | /* Save exception flags and rounding mode. */ |
| 710 | __fegetenv (&ctx->env); |
| 711 | |
| 712 | /* Update rounding mode only if different. */ |
| 713 | if (__glibc_unlikely (round != get_rounding_mode ())) |
| 714 | __fesetround (round); |
| 715 | } |
| 716 | |
| 717 | static __always_inline void |
| 718 | libc_feresetround_noex_ctx (struct rm_ctx *ctx) |
| 719 | { |
| 720 | /* Restore exception flags and rounding mode. */ |
| 721 | __fesetenv (&ctx->env); |
| 722 | } |
| 723 | |
| 724 | # define libc_feholdsetroundf_ctx libc_feholdsetround_ctx |
| 725 | # define libc_feholdsetroundl_ctx libc_feholdsetround_ctx |
| 726 | # define libc_feresetroundf_ctx libc_feresetround_ctx |
| 727 | # define libc_feresetroundl_ctx libc_feresetround_ctx |
| 728 | |
| 729 | # define libc_feholdsetround_noexf_ctx libc_feholdsetround_noex_ctx |
| 730 | # define libc_feholdsetround_noexl_ctx libc_feholdsetround_noex_ctx |
| 731 | # define libc_feresetround_noexf_ctx libc_feresetround_noex_ctx |
| 732 | # define libc_feresetround_noexl_ctx libc_feresetround_noex_ctx |
| 733 | |
| 734 | #endif |
| 735 | |
| 736 | #ifndef libc_feholdsetround_53bit_ctx |
| 737 | # define libc_feholdsetround_53bit_ctx libc_feholdsetround_ctx |
| 738 | #endif |
| 739 | #ifndef libc_feresetround_53bit_ctx |
| 740 | # define libc_feresetround_53bit_ctx libc_feresetround_ctx |
| 741 | #endif |
| 742 | |
| 743 | #define SET_RESTORE_ROUND_GENERIC(RM,ROUNDFUNC,CLEANUPFUNC) \ |
| 744 | struct rm_ctx ctx __attribute__((cleanup (CLEANUPFUNC ## _ctx))); \ |
| 745 | ROUNDFUNC ## _ctx (&ctx, (RM)) |
| 746 | |
| 747 | /* Set the rounding mode within a lexical block. Restore the rounding mode to |
| 748 | the value at the start of the block. The exception mode must be preserved. |
| 749 | Exceptions raised within the block must be set in the exception flags. |
| 750 | Non-stop mode may be enabled inside the block. */ |
| 751 | |
| 752 | #define SET_RESTORE_ROUND(RM) \ |
| 753 | SET_RESTORE_ROUND_GENERIC (RM, libc_feholdsetround, libc_feresetround) |
| 754 | #define SET_RESTORE_ROUNDF(RM) \ |
| 755 | SET_RESTORE_ROUND_GENERIC (RM, libc_feholdsetroundf, libc_feresetroundf) |
| 756 | #define SET_RESTORE_ROUNDL(RM) \ |
| 757 | SET_RESTORE_ROUND_GENERIC (RM, libc_feholdsetroundl, libc_feresetroundl) |
| 758 | |
| 759 | /* Set the rounding mode within a lexical block. Restore the rounding mode to |
| 760 | the value at the start of the block. The exception mode must be preserved. |
| 761 | Exceptions raised within the block must be discarded, and exception flags |
| 762 | are restored to the value at the start of the block. |
| 763 | Non-stop mode may be enabled inside the block. */ |
| 764 | |
| 765 | #define SET_RESTORE_ROUND_NOEX(RM) \ |
| 766 | SET_RESTORE_ROUND_GENERIC (RM, libc_feholdsetround_noex, \ |
| 767 | libc_feresetround_noex) |
| 768 | #define SET_RESTORE_ROUND_NOEXF(RM) \ |
| 769 | SET_RESTORE_ROUND_GENERIC (RM, libc_feholdsetround_noexf, \ |
| 770 | libc_feresetround_noexf) |
| 771 | #define SET_RESTORE_ROUND_NOEXL(RM) \ |
| 772 | SET_RESTORE_ROUND_GENERIC (RM, libc_feholdsetround_noexl, \ |
| 773 | libc_feresetround_noexl) |
| 774 | |
| 775 | /* Like SET_RESTORE_ROUND, but also set rounding precision to 53 bits. */ |
| 776 | #define SET_RESTORE_ROUND_53BIT(RM) \ |
| 777 | SET_RESTORE_ROUND_GENERIC (RM, libc_feholdsetround_53bit, \ |
| 778 | libc_feresetround_53bit) |
| 779 | |
| 780 | #define __nan(str) \ |
| 781 | (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str)) |
| 782 | #define __nanf(str) \ |
| 783 | (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str)) |
| 784 | #define __nanl(str) \ |
| 785 | (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str)) |
| 786 | |
| 787 | #endif /* _MATH_PRIVATE_H_ */ |