sysdeps/generic/math_private.h

   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
  24 /* The original fdlibm code used statements like:
  25         n0 = ((*(int*)&one)>>29)^1;             * index of high word *
  26         ix0 = *(n0+(int*)&x);                   * high word of x *
  27         ix1 = *((1-n0)+(int*)&x);               * low word of x *
  28    to dig two 32 bit words out of the 64 bit IEEE floating point
  29    value.  That is non-ANSI, and, moreover, the gcc instruction
  30    scheduler gets it wrong.  We instead use the following macros.
  31    Unlike the original code, we determine the endianness at compile
  32    time, not at run time; I don't see much benefit to selecting
  33    endianness at run time.  */
  34
  35 /* A union which permits us to convert between a double and two 32 bit
  36    ints.  */
  37
  38 #if __FLOAT_WORD_ORDER == BIG_ENDIAN
  39
  40 typedef union
  41 {
  42   double value;
  43   struct
  44   {
  45     u_int32_t msw;
  46     u_int32_t lsw;
  47   } parts;
  48   uint64_t word;
  49 } ieee_double_shape_type;
  50
  51 #endif
  52
  53 #if __FLOAT_WORD_ORDER == LITTLE_ENDIAN
  54
  55 typedef union
  56 {
  57   double value;
  58   struct
  59   {
  60     u_int32_t lsw;
  61     u_int32_t msw;
  62   } parts;
  63   uint64_t word;
  64 } ieee_double_shape_type;
  65
  66 #endif
  67
  68 /* Get two 32 bit ints from a double.  */
  69
  70 #define EXTRACT_WORDS(ix0,ix1,d)                                \
  71 do {                                                            \
  72   ieee_double_shape_type ew_u;                                  \
  73   ew_u.value = (d);                                             \
  74   (ix0) = ew_u.parts.msw;                                       \
  75   (ix1) = ew_u.parts.lsw;                                       \
  76 } while (0)
  77
  78 /* Get the more significant 32 bit int from a double.  */
  79
  80 #ifndef GET_HIGH_WORD
  81 # define GET_HIGH_WORD(i,d)                                     \
  82 do {                                                            \
  83   ieee_double_shape_type gh_u;                                  \
  84   gh_u.value = (d);                                             \
  85   (i) = gh_u.parts.msw;                                         \
  86 } while (0)
  87 #endif
  88
  89 /* Get the less significant 32 bit int from a double.  */
  90
  91 #ifndef GET_LOW_WORD
  92 # define GET_LOW_WORD(i,d)                                      \
  93 do {                                                            \
  94   ieee_double_shape_type gl_u;                                  \
  95   gl_u.value = (d);                                             \
  96   (i) = gl_u.parts.lsw;                                         \
  97 } while (0)
  98 #endif
  99
 100 /* Get all in one, efficient on 64-bit machines.  */
 101 #ifndef EXTRACT_WORDS64
 102 # define EXTRACT_WORDS64(i,d)                                   \
 103 do {                                                            \
 104   ieee_double_shape_type gh_u;                                  \
 105   gh_u.value = (d);                                             \
 106   (i) = gh_u.word;                                              \
 107 } while (0)
 108 #endif
 109
 110 /* Set a double from two 32 bit ints.  */
 111 #ifndef INSERT_WORDS
 112 # define INSERT_WORDS(d,ix0,ix1)                                \
 113 do {                                                            \
 114   ieee_double_shape_type iw_u;                                  \
 115   iw_u.parts.msw = (ix0);                                       \
 116   iw_u.parts.lsw = (ix1);                                       \
 117   (d) = iw_u.value;                                             \
 118 } while (0)
 119 #endif
 120
 121 /* Get all in one, efficient on 64-bit machines.  */
 122 #ifndef INSERT_WORDS64
 123 # define INSERT_WORDS64(d,i)                                    \
 124 do {                                                            \
 125   ieee_double_shape_type iw_u;                                  \
 126   iw_u.word = (i);                                              \
 127   (d) = iw_u.value;                                             \
 128 } while (0)
 129 #endif
 130
 131 /* Set the more significant 32 bits of a double from an int.  */
 132 #ifndef SET_HIGH_WORD
 133 #define SET_HIGH_WORD(d,v)                                      \
 134 do {                                                            \
 135   ieee_double_shape_type sh_u;                                  \
 136   sh_u.value = (d);                                             \
 137   sh_u.parts.msw = (v);                                         \
 138   (d) = sh_u.value;                                             \
 139 } while (0)
 140 #endif
 141
 142 /* Set the less significant 32 bits of a double from an int.  */
 143 #ifndef SET_LOW_WORD
 144 # define SET_LOW_WORD(d,v)                                      \
 145 do {                                                            \
 146   ieee_double_shape_type sl_u;                                  \
 147   sl_u.value = (d);                                             \
 148   sl_u.parts.lsw = (v);                                         \
 149   (d) = sl_u.value;                                             \
 150 } while (0)
 151 #endif
 152
 153 /* A union which permits us to convert between a float and a 32 bit
 154    int.  */
 155
 156 typedef union
 157 {
 158   float value;
 159   u_int32_t word;
 160 } ieee_float_shape_type;
 161
 162 /* Get a 32 bit int from a float.  */
 163 #ifndef GET_FLOAT_WORD
 164 # define GET_FLOAT_WORD(i,d)                                    \
 165 do {                                                            \
 166   ieee_float_shape_type gf_u;                                   \
 167   gf_u.value = (d);                                             \
 168   (i) = gf_u.word;                                              \
 169 } while (0)
 170 #endif
 171
 172 /* Set a float from a 32 bit int.  */
 173 #ifndef SET_FLOAT_WORD
 174 # define SET_FLOAT_WORD(d,i)                                    \
 175 do {                                                            \
 176   ieee_float_shape_type sf_u;                                   \
 177   sf_u.word = (i);                                              \
 178   (d) = sf_u.value;                                             \
 179 } while (0)
 180 #endif
 181
 182 /* Get long double macros from a separate header.  */
 183 #include <math_ldbl.h>
 184
 185 /* ieee style elementary functions */
 186 extern double __ieee754_sqrt (double);
 187 extern double __ieee754_acos (double);
 188 extern double __ieee754_acosh (double);
 189 extern double __ieee754_log (double);
 190 extern double __ieee754_atanh (double);
 191 extern double __ieee754_asin (double);
 192 extern double __ieee754_atan2 (double,double);
 193 extern double __ieee754_exp (double);
 194 extern double __ieee754_exp2 (double);
 195 extern double __ieee754_exp10 (double);
 196 extern double __ieee754_cosh (double);
 197 extern double __ieee754_fmod (double,double);
 198 extern double __ieee754_pow (double,double);
 199 extern double __ieee754_lgamma_r (double,int *);
 200 extern double __ieee754_gamma_r (double,int *);
 201 extern double __ieee754_lgamma (double);
 202 extern double __ieee754_gamma (double);
 203 extern double __ieee754_log10 (double);
 204 extern double __ieee754_log2 (double);
 205 extern double __ieee754_sinh (double);
 206 extern double __ieee754_hypot (double,double);
 207 extern double __ieee754_j0 (double);
 208 extern double __ieee754_j1 (double);
 209 extern double __ieee754_y0 (double);
 210 extern double __ieee754_y1 (double);
 211 extern double __ieee754_jn (int,double);
 212 extern double __ieee754_yn (int,double);
 213 extern double __ieee754_remainder (double,double);
 214 extern int32_t __ieee754_rem_pio2 (double,double*);
 215 extern double __ieee754_scalb (double,double);
 216 extern int __ieee754_ilogb (double);
 217
 218 /* fdlibm kernel function */
 219 extern double __kernel_standard (double,double,int);
 220 extern float __kernel_standard_f (float,float,int);
 221 extern long double __kernel_standard_l (long double,long double,int);
 222 extern double __kernel_sin (double,double,int);
 223 extern double __kernel_cos (double,double);
 224 extern double __kernel_tan (double,double,int);
 225 extern int    __kernel_rem_pio2 (double*,double*,int,int,int, const int32_t*);
 226
 227 /* internal functions.  */
 228 extern double __copysign (double x, double __y);
 229
 230 extern inline double __copysign (double x, double y)
 231 { return __builtin_copysign (x, y); }
 232
 233 /* ieee style elementary float functions */
 234 extern float __ieee754_sqrtf (float);
 235 extern float __ieee754_acosf (float);
 236 extern float __ieee754_acoshf (float);
 237 extern float __ieee754_logf (float);
 238 extern float __ieee754_atanhf (float);
 239 extern float __ieee754_asinf (float);
 240 extern float __ieee754_atan2f (float,float);
 241 extern float __ieee754_expf (float);
 242 extern float __ieee754_exp2f (float);
 243 extern float __ieee754_exp10f (float);
 244 extern float __ieee754_coshf (float);
 245 extern float __ieee754_fmodf (float,float);
 246 extern float __ieee754_powf (float,float);
 247 extern float __ieee754_lgammaf_r (float,int *);
 248 extern float __ieee754_gammaf_r (float,int *);
 249 extern float __ieee754_lgammaf (float);
 250 extern float __ieee754_gammaf (float);
 251 extern float __ieee754_log10f (float);
 252 extern float __ieee754_log2f (float);
 253 extern float __ieee754_sinhf (float);
 254 extern float __ieee754_hypotf (float,float);
 255 extern float __ieee754_j0f (float);
 256 extern float __ieee754_j1f (float);
 257 extern float __ieee754_y0f (float);
 258 extern float __ieee754_y1f (float);
 259 extern float __ieee754_jnf (int,float);
 260 extern float __ieee754_ynf (int,float);
 261 extern float __ieee754_remainderf (float,float);
 262 extern int32_t __ieee754_rem_pio2f (float,float*);
 263 extern float __ieee754_scalbf (float,float);
 264 extern int __ieee754_ilogbf (float);
 265
 266
 267 /* float versions of fdlibm kernel functions */
 268 extern float __kernel_sinf (float,float,int);
 269 extern float __kernel_cosf (float,float);
 270 extern float __kernel_tanf (float,float,int);
 271 extern int   __kernel_rem_pio2f (float*,float*,int,int,int, const int32_t*);
 272
 273 /* internal functions.  */
 274 extern float __copysignf (float x, float __y);
 275
 276 extern inline float __copysignf (float x, float y)
 277 { return __builtin_copysignf (x, y); }
 278
 279 /* ieee style elementary long double functions */
 280 extern long double __ieee754_sqrtl (long double);
 281 extern long double __ieee754_acosl (long double);
 282 extern long double __ieee754_acoshl (long double);
 283 extern long double __ieee754_logl (long double);
 284 extern long double __ieee754_atanhl (long double);
 285 extern long double __ieee754_asinl (long double);
 286 extern long double __ieee754_atan2l (long double,long double);
 287 extern long double __ieee754_expl (long double);
 288 extern long double __ieee754_exp2l (long double);
 289 extern long double __ieee754_exp10l (long double);
 290 extern long double __ieee754_coshl (long double);
 291 extern long double __ieee754_fmodl (long double,long double);
 292 extern long double __ieee754_powl (long double,long double);
 293 extern long double __ieee754_lgammal_r (long double,int *);
 294 extern long double __ieee754_gammal_r (long double,int *);
 295 extern long double __ieee754_lgammal (long double);
 296 extern long double __ieee754_gammal (long double);
 297 extern long double __ieee754_log10l (long double);
 298 extern long double __ieee754_log2l (long double);
 299 extern long double __ieee754_sinhl (long double);
 300 extern long double __ieee754_hypotl (long double,long double);
 301 extern long double __ieee754_j0l (long double);
 302 extern long double __ieee754_j1l (long double);
 303 extern long double __ieee754_y0l (long double);
 304 extern long double __ieee754_y1l (long double);
 305 extern long double __ieee754_jnl (int,long double);
 306 extern long double __ieee754_ynl (int,long double);
 307 extern long double __ieee754_remainderl (long double,long double);
 308 extern int   __ieee754_rem_pio2l (long double,long double*);
 309 extern long double __ieee754_scalbl (long double,long double);
 310 extern int   __ieee754_ilogbl (long double);
 311
 312 /* long double versions of fdlibm kernel functions */
 313 extern long double __kernel_sinl (long double,long double,int);
 314 extern long double __kernel_cosl (long double,long double);
 315 extern long double __kernel_tanl (long double,long double,int);
 316 extern void __kernel_sincosl (long double,long double,
 317                               long double *,long double *, int);
 318 extern int   __kernel_rem_pio2l (long double*,long double*,int,int,
 319                                  int,const int*);
 320
 321 #ifndef NO_LONG_DOUBLE
 322 /* prototypes required to compile the ldbl-96 support without warnings */
 323 extern int __finitel (long double);
 324 extern int __ilogbl (long double);
 325 extern int __isinfl (long double);
 326 extern int __isnanl (long double);
 327 extern long double __atanl (long double);
 328 extern long double __copysignl (long double, long double);
 329 extern long double __expm1l (long double);
 330 extern long double __floorl (long double);
 331 extern long double __frexpl (long double, int *);
 332 extern long double __ldexpl (long double, int);
 333 extern long double __log1pl (long double);
 334 extern long double __nanl (const char *);
 335 extern long double __rintl (long double);
 336 extern long double __scalbnl (long double, int);
 337 extern long double __sqrtl (long double x);
 338 extern long double fabsl (long double x);
 339 extern void __sincosl (long double, long double *, long double *);
 340 extern long double __logbl (long double x);
 341 extern long double __significandl (long double x);
 342
 343 extern inline long double __copysignl (long double x, long double y)
 344 { return __builtin_copysignl (x, y); }
 345
 346 #endif
 347
 348 /* Prototypes for functions of the IBM Accurate Mathematical Library.  */
 349 extern double __exp1 (double __x, double __xx, double __error);
 350 extern double __sin (double __x);
 351 extern double __cos (double __x);
 352 extern int __branred (double __x, double *__a, double *__aa);
 353 extern void __doasin (double __x, double __dx, double __v[]);
 354 extern void __dubsin (double __x, double __dx, double __v[]);
 355 extern void __dubcos (double __x, double __dx, double __v[]);
 356 extern double __halfulp (double __x, double __y);
 357 extern double __sin32 (double __x, double __res, double __res1);
 358 extern double __cos32 (double __x, double __res, double __res1);
 359 extern double __mpsin (double __x, double __dx);
 360 extern double __mpcos (double __x, double __dx);
 361 extern double __mpsin1 (double __x);
 362 extern double __mpcos1 (double __x);
 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 either X >=
 369    0.75 or Y >= 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 #ifndef math_opt_barrier
 375 # define math_opt_barrier(x) \
 376 ({ __typeof (x) __x = (x); __asm ("" : "+m" (__x)); __x; })
 377 # define math_force_eval(x) \
 378 ({ __typeof (x) __x = (x); __asm __volatile__ ("" : : "m" (__x)); })
 379 #endif
 380
 381
 382 /* The standards only specify one variant of the fenv.h interfaces.
 383    But at least for some architectures we can be more efficient if we
 384    know what operations are going to be performed.  Therefore we
 385    define additional interfaces.  By default they refer to the normal
 386    interfaces.  */
 387
 388 static __always_inline void
 389 default_libc_feholdexcept (fenv_t *e)
 390 {
 391   (void) feholdexcept (e);
 392 }
 393
 394 #ifndef libc_feholdexcept
 395 # define libc_feholdexcept  default_libc_feholdexcept
 396 #endif
 397 #ifndef libc_feholdexceptf
 398 # define libc_feholdexceptf default_libc_feholdexcept
 399 #endif
 400 #ifndef libc_feholdexceptl
 401 # define libc_feholdexceptl default_libc_feholdexcept
 402 #endif
 403
 404 static __always_inline void
 405 default_libc_fesetround (int r)
 406 {
 407   (void) fesetround (r);
 408 }
 409
 410 #ifndef libc_fesetround
 411 # define libc_fesetround  default_libc_fesetround
 412 #endif
 413 #ifndef libc_fesetroundf
 414 # define libc_fesetroundf default_libc_fesetround
 415 #endif
 416 #ifndef libc_fesetroundl
 417 # define libc_fesetroundl default_libc_fesetround
 418 #endif
 419
 420 static __always_inline void
 421 default_libc_feholdexcept_setround (fenv_t *e, int r)
 422 {
 423   feholdexcept (e);
 424   fesetround (r);
 425 }
 426
 427 #ifndef libc_feholdexcept_setround
 428 # define libc_feholdexcept_setround  default_libc_feholdexcept_setround
 429 #endif
 430 #ifndef libc_feholdexcept_setroundf
 431 # define libc_feholdexcept_setroundf default_libc_feholdexcept_setround
 432 #endif
 433 #ifndef libc_feholdexcept_setroundl
 434 # define libc_feholdexcept_setroundl default_libc_feholdexcept_setround
 435 #endif
 436
 437 #ifndef libc_feholdexcept_setround_53bit
 438 # define libc_feholdexcept_setround_53bit libc_feholdexcept_setround
 439 #endif
 440
 441 #ifndef libc_fetestexcept
 442 # define libc_fetestexcept  fetestexcept
 443 #endif
 444 #ifndef libc_fetestexceptf
 445 # define libc_fetestexceptf fetestexcept
 446 #endif
 447 #ifndef libc_fetestexceptl
 448 # define libc_fetestexceptl fetestexcept
 449 #endif
 450
 451 static __always_inline void
 452 default_libc_fesetenv (fenv_t *e)
 453 {
 454   (void) fesetenv (e);
 455 }
 456
 457 #ifndef libc_fesetenv
 458 # define libc_fesetenv  default_libc_fesetenv
 459 #endif
 460 #ifndef libc_fesetenvf
 461 # define libc_fesetenvf default_libc_fesetenv
 462 #endif
 463 #ifndef libc_fesetenvl
 464 # define libc_fesetenvl default_libc_fesetenv
 465 #endif
 466
 467 static __always_inline void
 468 default_libc_feupdateenv (fenv_t *e)
 469 {
 470   (void) feupdateenv (e);
 471 }
 472
 473 #ifndef libc_feupdateenv
 474 # define libc_feupdateenv  default_libc_feupdateenv
 475 #endif
 476 #ifndef libc_feupdateenvf
 477 # define libc_feupdateenvf default_libc_feupdateenv
 478 #endif
 479 #ifndef libc_feupdateenvl
 480 # define libc_feupdateenvl default_libc_feupdateenv
 481 #endif
 482
 483 #ifndef libc_feupdateenv_53bit
 484 # define libc_feupdateenv_53bit libc_feupdateenv
 485 #endif
 486
 487 static __always_inline int
 488 default_libc_feupdateenv_test (fenv_t *e, int ex)
 489 {
 490   int ret = fetestexcept (ex);
 491   feupdateenv (e);
 492   return ret;
 493 }
 494
 495 #ifndef libc_feupdateenv_test
 496 # define libc_feupdateenv_test  default_libc_feupdateenv_test
 497 #endif
 498 #ifndef libc_feupdateenv_testf
 499 # define libc_feupdateenv_testf default_libc_feupdateenv_test
 500 #endif
 501 #ifndef libc_feupdateenv_testl
 502 # define libc_feupdateenv_testl default_libc_feupdateenv_test
 503 #endif
 504
 505 /* Save and set the rounding mode.  The use of fenv_t to store the old mode
 506    allows a target-specific version of this function to avoid converting the
 507    rounding mode from the fpu format.  By default we have no choice but to
 508    manipulate the entire env.  */
 509
 510 #ifndef libc_feholdsetround
 511 # define libc_feholdsetround  libc_feholdexcept_setround
 512 #endif
 513 #ifndef libc_feholdsetroundf
 514 # define libc_feholdsetroundf libc_feholdexcept_setroundf
 515 #endif
 516 #ifndef libc_feholdsetroundl
 517 # define libc_feholdsetroundl libc_feholdexcept_setroundl
 518 #endif
 519
 520 /* ... and the reverse.  */
 521
 522 #ifndef libc_feresetround
 523 # define libc_feresetround  libc_feupdateenv
 524 #endif
 525 #ifndef libc_feresetroundf
 526 # define libc_feresetroundf libc_feupdateenvf
 527 #endif
 528 #ifndef libc_feresetroundl
 529 # define libc_feresetroundl libc_feupdateenvl
 530 #endif
 531
 532 /* ... and a version that may also discard exceptions.  */
 533
 534 #ifndef libc_feresetround_noex
 535 # define libc_feresetround_noex  libc_fesetenv
 536 #endif
 537 #ifndef libc_feresetround_noexf
 538 # define libc_feresetround_noexf libc_fesetenvf
 539 #endif
 540 #ifndef libc_feresetround_noexl
 541 # define libc_feresetround_noexl libc_fesetenvl
 542 #endif
 543
 544 /* Save and restore the rounding mode within a lexical block.  */
 545
 546 #define SET_RESTORE_ROUND(RM) \
 547   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround)));    \
 548   libc_feholdsetround (&__libc_save_rm, (RM))
 549 #define SET_RESTORE_ROUNDF(RM) \
 550   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetroundf)));   \
 551   libc_feholdsetroundf (&__libc_save_rm, (RM))
 552 #define SET_RESTORE_ROUNDL(RM) \
 553   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetroundl)));   \
 554   libc_feholdsetroundl (&__libc_save_rm, (RM))
 555
 556 /* Save and restore the rounding mode within a lexical block, and also
 557    the set of exceptions raised within the block may be discarded.  */
 558
 559 #define SET_RESTORE_ROUND_NOEX(RM) \
 560   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround_noex))); \
 561   libc_feholdsetround (&__libc_save_rm, (RM))
 562 #define SET_RESTORE_ROUND_NOEXF(RM) \
 563   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround_noexf))); \
 564   libc_feholdsetroundf (&__libc_save_rm, (RM))
 565 #define SET_RESTORE_ROUND_NOEXL(RM) \
 566   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround_noexl))); \
 567   libc_feholdsetroundl (&__libc_save_rm, (RM))
 568
 569 /* Like SET_RESTORE_ROUND, but also set rounding precision to 53 bits.  */
 570 #define SET_RESTORE_ROUND_53BIT(RM) \
 571   fenv_t __libc_save_rm __attribute__((cleanup(libc_feupdateenv_53bit))); \
 572   libc_feholdexcept_setround_53bit (&__libc_save_rm, (RM))
 573
 574 #define __nan(str) \
 575   (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str))
 576 #define __nanf(str) \
 577   (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str))
 578 #define __nanl(str) \
 579   (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str))
 580
 581 #endif /* _MATH_PRIVATE_H_ */