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_feholdexcept_setround (fenv_t *e, int r)
 406 {
 407   feholdexcept (e);
 408   fesetround (r);
 409 }
 410
 411 #ifndef libc_feholdexcept_setround
 412 # define libc_feholdexcept_setround  default_libc_feholdexcept_setround
 413 #endif
 414 #ifndef libc_feholdexcept_setroundf
 415 # define libc_feholdexcept_setroundf default_libc_feholdexcept_setround
 416 #endif
 417 #ifndef libc_feholdexcept_setroundl
 418 # define libc_feholdexcept_setroundl default_libc_feholdexcept_setround
 419 #endif
 420
 421 #ifndef libc_feholdexcept_setround_53bit
 422 # define libc_feholdexcept_setround_53bit libc_feholdexcept_setround
 423 #endif
 424
 425 #ifndef libc_fetestexcept
 426 # define libc_fetestexcept  fetestexcept
 427 #endif
 428 #ifndef libc_fetestexceptf
 429 # define libc_fetestexceptf fetestexcept
 430 #endif
 431 #ifndef libc_fetestexceptl
 432 # define libc_fetestexceptl fetestexcept
 433 #endif
 434
 435 static __always_inline void
 436 default_libc_fesetenv (fenv_t *e)
 437 {
 438   (void) fesetenv (e);
 439 }
 440
 441 #ifndef libc_fesetenv
 442 # define libc_fesetenv  default_libc_fesetenv
 443 #endif
 444 #ifndef libc_fesetenvf
 445 # define libc_fesetenvf default_libc_fesetenv
 446 #endif
 447 #ifndef libc_fesetenvl
 448 # define libc_fesetenvl default_libc_fesetenv
 449 #endif
 450
 451 static __always_inline void
 452 default_libc_feupdateenv (fenv_t *e)
 453 {
 454   (void) feupdateenv (e);
 455 }
 456
 457 #ifndef libc_feupdateenv
 458 # define libc_feupdateenv  default_libc_feupdateenv
 459 #endif
 460 #ifndef libc_feupdateenvf
 461 # define libc_feupdateenvf default_libc_feupdateenv
 462 #endif
 463 #ifndef libc_feupdateenvl
 464 # define libc_feupdateenvl default_libc_feupdateenv
 465 #endif
 466
 467 #ifndef libc_feupdateenv_53bit
 468 # define libc_feupdateenv_53bit libc_feupdateenv
 469 #endif
 470
 471 static __always_inline int
 472 default_libc_feupdateenv_test (fenv_t *e, int ex)
 473 {
 474   int ret = fetestexcept (ex);
 475   feupdateenv (e);
 476   return ret;
 477 }
 478
 479 #ifndef libc_feupdateenv_test
 480 # define libc_feupdateenv_test  default_libc_feupdateenv_test
 481 #endif
 482 #ifndef libc_feupdateenv_testf
 483 # define libc_feupdateenv_testf default_libc_feupdateenv_test
 484 #endif
 485 #ifndef libc_feupdateenv_testl
 486 # define libc_feupdateenv_testl default_libc_feupdateenv_test
 487 #endif
 488
 489 /* Save and set the rounding mode.  The use of fenv_t to store the old mode
 490    allows a target-specific version of this function to avoid converting the
 491    rounding mode from the fpu format.  By default we have no choice but to
 492    manipulate the entire env.  */
 493
 494 #ifndef libc_feholdsetround
 495 # define libc_feholdsetround  libc_feholdexcept_setround
 496 #endif
 497 #ifndef libc_feholdsetroundf
 498 # define libc_feholdsetroundf libc_feholdexcept_setroundf
 499 #endif
 500 #ifndef libc_feholdsetroundl
 501 # define libc_feholdsetroundl libc_feholdexcept_setroundl
 502 #endif
 503
 504 /* ... and the reverse.  */
 505
 506 #ifndef libc_feresetround
 507 # define libc_feresetround  libc_feupdateenv
 508 #endif
 509 #ifndef libc_feresetroundf
 510 # define libc_feresetroundf libc_feupdateenvf
 511 #endif
 512 #ifndef libc_feresetroundl
 513 # define libc_feresetroundl libc_feupdateenvl
 514 #endif
 515
 516 /* ... and a version that may also discard exceptions.  */
 517
 518 #ifndef libc_feresetround_noex
 519 # define libc_feresetround_noex  libc_fesetenv
 520 #endif
 521 #ifndef libc_feresetround_noexf
 522 # define libc_feresetround_noexf libc_fesetenvf
 523 #endif
 524 #ifndef libc_feresetround_noexl
 525 # define libc_feresetround_noexl libc_fesetenvl
 526 #endif
 527
 528 /* Save and restore the rounding mode within a lexical block.  */
 529
 530 #define SET_RESTORE_ROUND(RM) \
 531   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround)));    \
 532   libc_feholdsetround (&__libc_save_rm, (RM))
 533 #define SET_RESTORE_ROUNDF(RM) \
 534   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetroundf)));   \
 535   libc_feholdsetroundf (&__libc_save_rm, (RM))
 536 #define SET_RESTORE_ROUNDL(RM) \
 537   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetroundl)));   \
 538   libc_feholdsetroundl (&__libc_save_rm, (RM))
 539
 540 /* Save and restore the rounding mode within a lexical block, and also
 541    the set of exceptions raised within the block may be discarded.  */
 542
 543 #define SET_RESTORE_ROUND_NOEX(RM) \
 544   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround_noex))); \
 545   libc_feholdsetround (&__libc_save_rm, (RM))
 546 #define SET_RESTORE_ROUND_NOEXF(RM) \
 547   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround_noexf))); \
 548   libc_feholdsetroundf (&__libc_save_rm, (RM))
 549 #define SET_RESTORE_ROUND_NOEXL(RM) \
 550   fenv_t __libc_save_rm __attribute__((cleanup(libc_feresetround_noexl))); \
 551   libc_feholdsetroundl (&__libc_save_rm, (RM))
 552
 553 /* Like SET_RESTORE_ROUND, but also set rounding precision to 53 bits.  */
 554 #define SET_RESTORE_ROUND_53BIT(RM) \
 555   fenv_t __libc_save_rm __attribute__((cleanup(libc_feupdateenv_53bit))); \
 556   libc_feholdexcept_setround_53bit (&__libc_save_rm, (RM))
 557
 558 #define __nan(str) \
 559   (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str))
 560 #define __nanf(str) \
 561   (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str))
 562 #define __nanl(str) \
 563   (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str))
 564
 565 #endif /* _MATH_PRIVATE_H_ */