sysdeps/i386/fpu/bits/mathinline.h

   1 /* Inline math functions for i387.
   2    Copyright (C) 1995, 1996, 1997 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by John C. Bowman <bowman@ipp-garching.mpg.de>, 1995.
   5
   6    The GNU C Library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Library General Public License as
   8    published by the Free Software Foundation; either version 2 of the
   9    License, or (at your option) any later version.
  10
  11    The GNU C Library is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14    Library General Public License for more details.
  15
  16    You should have received a copy of the GNU Library General Public
  17    License along with the GNU C Library; see the file COPYING.LIB.  If not,
  18    write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19    Boston, MA 02111-1307, USA.  */
  20
  21 #ifndef _BITS_MATHINLINE_H
  22 #define _BITS_MATHINLINE_H      1
  23
  24
  25 #if defined __USE_ISOC9X && defined __GNUC__ && __GNUC__ >= 2
  26 /* ISO C 9X defines some macros to perform unordered comparisons.  The
  27    ix87 FPU supports this with special opcodes and we should use them.
  28    These must not be inline functions since we have to be able to handle
  29    all floating-point types.  */
  30 # define isgreater(x, y) \
  31      ({ int result;                                                           \
  32         __asm__ ("fucompp; fnstsw; andb $0x45, %%ah; setz %%al;"              \
  33                  "andl $0x01, %0"                                             \
  34                  : "=a" (result) : "u" (y), "t" (x) : "cc", "st", "st(1)");   \
  35         result; })
  36
  37 # define isgreaterequal(x, y) \
  38      ({ int result;                                                           \
  39         __asm__ ("fucompp; fnstsw; testb $0x05, %%ah; setz %%al;"             \
  40                  "andl $0x01, %0"                                             \
  41                  : "=a" (result) : "u" (y), "t" (x) : "cc", "st", "st(1)");   \
  42         result; })
  43
  44 # define isless(x, y) \
  45      ({ int result;                                                           \
  46         __asm__ ("fucompp; fnstsw; xorb $0x01, %%ah; testb $0x45, %%ah;"      \
  47                  "setz %%al; andl $0x01, %0"                                  \
  48                  : "=a" (result) : "u" (y), "t" (x) : "cc", "st", "st(1)");   \
  49         result; })
  50
  51 # define islessequal(x, y) \
  52      ({ int result;                                                           \
  53         __asm__ ("fucompp; fnstsw; xorb $0x01, %%ah; testb $0x05, %%ah;"      \
  54                  "setz %%al; andl $0x01, %0"                                  \
  55                  : "=a" (result) : "u" (y), "t" (x) : "cc", "st", "st(1)");   \
  56         result; })
  57
  58 # define islessgreater(x, y) \
  59      ({ int result;                                                           \
  60         __asm__ ("fucompp; fnstsw; testb $0x44, %%ah; setz %%al;"             \
  61                  "andl $0x01, %0"                                             \
  62                  : "=a" (result) : "u" (y), "t" (x) : "cc", "st", "st(1)");   \
  63         result; })
  64
  65 # define isunordered(x, y) \
  66      ({ int result;                                                           \
  67         __asm__ ("fucompp; fnstsw; sahf; setp %%al; andl $0x01, %0"           \
  68                  : "=a" (result) : "u" (y), "t" (x) : "cc", "st", "st(1)");   \
  69         result; })
  70 #endif
  71
  72
  73 #if defined __GNUC__ && \
  74     (__GNUC__ < 2 || (__GNUC__ == 2 && __GNUC_MINOR__ <= 7))
  75 /* gcc 2.7.2 and 2.7.2.1 have problems with inlining `long double'
  76    functions so we disable this now.  */
  77 # undef __NO_MATH_INLINES
  78 # define __NO_MATH_INLINES
  79 #endif
  80
  81
  82 #ifdef  __GNUC__
  83 #if !defined __NO_MATH_INLINES && defined __OPTIMIZE__
  84
  85 #ifdef __cplusplus
  86 # define __MATH_INLINE __inline
  87 #else
  88 # define __MATH_INLINE extern __inline
  89 #endif
  90
  91 __MATH_INLINE double cos (double);
  92 __MATH_INLINE double sin (double);
  93
  94
  95 __MATH_INLINE double __expm1 (double __x);
  96 __MATH_INLINE double
  97 __expm1 (double __x)
  98 {
  99   register double __value, __exponent, __temp;
 100   __asm __volatile__
 101     ("fldl2e                    # e^x - 1 = 2^(x * log2(e)) - 1\n\t"
 102      "fmul      %%st(1)         # x * log2(e)\n\t"
 103      "fstl      %%st(1)\n\t"
 104      "frndint                   # int(x * log2(e))\n\t"
 105      "fxch\n\t"
 106      "fsub      %%st(1)         # fract(x * log2(e))\n\t"
 107      "f2xm1                     # 2^(fract(x * log2(e))) - 1\n\t"
 108      "fscale                    # 2^(x * log2(e)) - 2^(int(x * log2(e)))\n\t"
 109      : "=t" (__value), "=u" (__exponent) : "0" (__x));
 110   __asm __volatile__
 111     ("fscale                    # 2^int(x * log2(e))\n\t"
 112      : "=t" (__temp) : "0" (1.0), "u" (__exponent));
 113   __temp -= 1.0;
 114
 115   return __temp + __value;
 116 }
 117
 118 __MATH_INLINE double __sgn1 (double __x);
 119 __MATH_INLINE double
 120 __sgn1 (double __x)
 121 {
 122   return __x >= 0.0 ? 1.0 : -1.0;
 123 }
 124
 125 __MATH_INLINE double sqrt (double __x);
 126 __MATH_INLINE double
 127 sqrt (double __x)
 128 {
 129   register double __value;
 130   __asm __volatile__
 131     ("fsqrt"
 132      : "=t" (__value) : "0" (__x));
 133
 134   return __value;
 135 }
 136
 137 __MATH_INLINE double fabs (double __x);
 138 __MATH_INLINE double
 139 fabs (double __x)
 140 {
 141   register double __value;
 142   __asm __volatile__
 143     ("fabs"
 144      : "=t" (__value) : "0" (__x));
 145
 146   return __value;
 147 }
 148
 149 /* The argument range of this inline version is limited.  */
 150 __MATH_INLINE double sin (double __x);
 151 __MATH_INLINE double
 152 sin (double __x)
 153 {
 154   register double __value;
 155   __asm __volatile__
 156     ("fsin"
 157      : "=t" (__value) : "0" (__x));
 158
 159   return __value;
 160 }
 161
 162 /* The argument range of this inline version is limited.  */
 163 __MATH_INLINE double cos (double __x);
 164 __MATH_INLINE double
 165 cos (double __x)
 166 {
 167   register double __value;
 168   __asm __volatile__
 169     ("fcos"
 170      : "=t" (__value): "0" (__x));
 171
 172   return __value;
 173 }
 174
 175 __MATH_INLINE double tan (double __x);
 176 __MATH_INLINE double
 177 tan (double __x)
 178 {
 179   register double __value;
 180   register double __value2 __attribute__ ((unused));
 181   __asm __volatile__
 182     ("fptan"
 183      : "=t" (__value2), "=u" (__value) : "0" (__x));
 184
 185   return __value;
 186 }
 187
 188 __MATH_INLINE double atan2 (double __y, double __x);
 189 __MATH_INLINE double
 190 atan2 (double __y, double __x)
 191 {
 192   register double __value;
 193   __asm __volatile__
 194     ("fpatan\n\t"
 195      "fldl %%st(0)"
 196      : "=t" (__value) : "0" (__x), "u" (__y));
 197
 198   return __value;
 199 }
 200
 201 __MATH_INLINE double asin (double __x);
 202 __MATH_INLINE double
 203 asin (double __x)
 204 {
 205   return atan2 (__x, sqrt (1.0 - __x * __x));
 206 }
 207
 208 __MATH_INLINE double acos (double __x);
 209 __MATH_INLINE double
 210 acos (double __x)
 211 {
 212   return atan2 (sqrt (1.0 - __x * __x), __x);
 213 }
 214
 215 __MATH_INLINE double atan (double __x);
 216 __MATH_INLINE double
 217 atan (double __x)
 218 {
 219   register double __value;
 220   __asm __volatile__
 221     ("fld1\n\t"
 222      "fpatan"
 223      : "=t" (__value) : "0" (__x));
 224
 225   return __value;
 226 }
 227
 228 __MATH_INLINE double exp (double __x);
 229 __MATH_INLINE double
 230 exp (double __x)
 231 {
 232   register double __value, __exponent;
 233   __asm __volatile__
 234     ("fldl2e                    # e^x = 2^(x * log2(e))\n\t"
 235      "fmul      %%st(1)         # x * log2(e)\n\t"
 236      "fstl      %%st(1)\n\t"
 237      "frndint                   # int(x * log2(e))\n\t"
 238      "fxch\n\t"
 239      "fsub      %%st(1)         # fract(x * log2(e))\n\t"
 240      "f2xm1                     # 2^(fract(x * log2(e))) - 1\n\t"
 241      : "=t" (__value), "=u" (__exponent) : "0" (__x));
 242   __value += 1.0;
 243   __asm __volatile__
 244     ("fscale"
 245      : "=t" (__value) : "0" (__value), "u" (__exponent));
 246
 247   return __value;
 248 }
 249
 250 __MATH_INLINE double sinh (double __x);
 251 __MATH_INLINE double
 252 sinh (double __x)
 253 {
 254   register double __exm1 = __expm1 (fabs (__x));
 255
 256   return 0.5 * (__exm1 / (__exm1 + 1.0) + __exm1) * __sgn1 (__x);
 257 }
 258
 259 __MATH_INLINE double cosh (double __x);
 260 __MATH_INLINE double
 261 cosh (double __x)
 262 {
 263   register double __ex = exp (__x);
 264
 265   return 0.5 * (__ex + 1.0 / __ex);
 266 }
 267
 268 __MATH_INLINE double tanh (double __x);
 269 __MATH_INLINE double
 270 tanh (double __x)
 271 {
 272   register double __exm1 = __expm1 (-fabs (__x + __x));
 273
 274   return __exm1 / (__exm1 + 2.0) * __sgn1 (-__x);
 275 }
 276
 277 __MATH_INLINE double log (double __x);
 278 __MATH_INLINE double
 279 log (double __x)
 280 {
 281   register double __value;
 282   __asm __volatile__
 283     ("fldln2\n\t"
 284      "fxch\n\t"
 285      "fyl2x"
 286      : "=t" (__value) : "0" (__x));
 287
 288   return __value;
 289 }
 290
 291 __MATH_INLINE double log10 (double __x);
 292 __MATH_INLINE double
 293 log10 (double __x)
 294 {
 295   register double __value;
 296   __asm __volatile__
 297     ("fldlg2\n\t"
 298      "fxch\n\t"
 299      "fyl2x"
 300      : "=t" (__value) : "0" (__x));
 301
 302   return __value;
 303 }
 304
 305 __MATH_INLINE double __log2 (double __x);
 306 __MATH_INLINE double
 307 __log2 (double __x)
 308 {
 309   register double __value;
 310   __asm __volatile__
 311     ("fyl2x"
 312      : "=t" (__value) : "0" (__x), "u" (1.0));
 313
 314   return __value;
 315 }
 316
 317 __MATH_INLINE double fmod (double __x, double __y);
 318 __MATH_INLINE double
 319 fmod (double __x, double __y)
 320 {
 321   register double __value;
 322   __asm __volatile__
 323     ("1:        fprem\n\t"
 324      "fstsw     %%ax\n\t"
 325      "sahf\n\t"
 326      "jp        1b"
 327      : "=t" (__value) : "0" (__x), "u" (__y) : "ax", "cc");
 328
 329   return __value;
 330 }
 331
 332 __MATH_INLINE double ldexp (double __x, int __y);
 333 __MATH_INLINE double
 334 ldexp (double __x, int __y)
 335 {
 336   register double __value;
 337   __asm __volatile__
 338     ("fscale"
 339      : "=t" (__value) : "0" (__x), "u" ((double) __y));
 340
 341   return __value;
 342 }
 343
 344 __MATH_INLINE double pow (double __x, double __y);
 345 __MATH_INLINE double
 346 pow (double __x, double __y)
 347 {
 348   register double __value, __exponent;
 349   long __p = (long) __y;
 350
 351   if (__x == 0.0 && __y > 0.0)
 352     return 0.0;
 353   if (__y == (double) __p)
 354     {
 355       double __r = 1.0;
 356       if (__p == 0)
 357         return 1.0;
 358       if (__p < 0)
 359         {
 360           __p = -__p;
 361           __x = 1.0 / __x;
 362         }
 363       while (1)
 364         {
 365           if (__p & 1)
 366             __r *= __x;
 367           __p >>= 1;
 368           if (__p == 0)
 369             return __r;
 370           __x *= __x;
 371         }
 372       /* NOTREACHED */
 373     }
 374   __asm __volatile__
 375     ("fmul      %%st(1)         # y * log2(x)\n\t"
 376      "fstl      %%st(1)\n\t"
 377      "frndint                   # int(y * log2(x))\n\t"
 378      "fxch\n\t"
 379      "fsub      %%st(1)         # fract(y * log2(x))\n\t"
 380      "f2xm1                     # 2^(fract(y * log2(x))) - 1\n\t"
 381      : "=t" (__value), "=u" (__exponent) :  "0" (__log2 (__x)), "1" (__y));
 382   __value += 1.0;
 383   __asm __volatile__
 384     ("fscale"
 385      : "=t" (__value) : "0" (__value), "u" (__exponent));
 386
 387   return __value;
 388 }
 389
 390 __MATH_INLINE double floor (double __x);
 391 __MATH_INLINE double
 392 floor (double __x)
 393 {
 394   register double __value;
 395   __volatile unsigned short int __cw, __cwtmp;
 396
 397   __asm __volatile ("fnstcw %0" : "=m" (__cw));
 398   __cwtmp = (__cw & 0xf3ff) | 0x0400; /* rounding down */
 399   __asm __volatile ("fldcw %0" : : "m" (__cwtmp));
 400   __asm __volatile ("frndint" : "=t" (__value) : "0" (__x));
 401   __asm __volatile ("fldcw %0" : : "m" (__cw));
 402
 403   return __value;
 404 }
 405
 406 __MATH_INLINE double ceil (double __x);
 407 __MATH_INLINE double
 408 ceil (double __x)
 409 {
 410   register double __value;
 411   __volatile unsigned short int __cw, __cwtmp;
 412
 413   __asm __volatile ("fnstcw %0" : "=m" (__cw));
 414   __cwtmp = (__cw & 0xf3ff) | 0x0800; /* rounding up */
 415   __asm __volatile ("fldcw %0" : : "m" (__cwtmp));
 416   __asm __volatile ("frndint" : "=t" (__value) : "0" (__x));
 417   __asm __volatile ("fldcw %0" : : "m" (__cw));
 418
 419   return __value;
 420 }
 421
 422
 423 /* Optimized versions for some non-standardized functions.  */
 424 #if defined __USE_ISOC9X || defined __USE_MISC
 425
 426 __MATH_INLINE double hypot (double __x, double __y);
 427 __MATH_INLINE double
 428 hypot (double __x, double __y)
 429 {
 430   return sqrt (__x * __x + __y * __y);
 431 }
 432
 433 /* We cannot rely on M_SQRT being defined.  So we do it for ourself
 434    here.  */
 435 # define __M_SQRT2      _Mldbl(1.41421356237309504880)  /* sqrt(2) */
 436
 437 __MATH_INLINE double log1p (double __x);
 438 __MATH_INLINE double
 439 log1p (double __x)
 440 {
 441   register double __value;
 442
 443   if (fabs (__x) >= 1.0 - 0.5 * __M_SQRT2)
 444     __value = log (1.0 + __x);
 445   else
 446     __asm __volatile__
 447       ("fldln2\n\t"
 448        "fxch\n\t"
 449        "fyl2xp1"
 450        : "=t" (__value) : "0" (__x));
 451
 452   return __value;
 453 }
 454
 455 __MATH_INLINE double asinh (double __x);
 456 __MATH_INLINE double
 457 asinh (double __x)
 458 {
 459   register double __y = fabs (__x);
 460
 461   return log1p ((__y * __y / (sqrt (__y * __y + 1.0) + 1.0) + __y)
 462                 * __sgn1 (__x));
 463 }
 464
 465 __MATH_INLINE double acosh (double __x);
 466 __MATH_INLINE double
 467 acosh (double __x)
 468 {
 469   return log (__x + sqrt (__x - 1.0) * sqrt (__x + 1.0));
 470 }
 471
 472 __MATH_INLINE double atanh (double __x);
 473 __MATH_INLINE double
 474 atanh (double __x)
 475 {
 476   register double __y = fabs (__x);
 477
 478   return -0.5 * __log1p (-(__y + __y) / (1.0 + __y)) * __sgn1 (__x);
 479 }
 480
 481 __MATH_INLINE double logb (double __x);
 482 __MATH_INLINE double
 483 logb (double __x)
 484 {
 485   register double __value, __junk;
 486   __asm __volatile__
 487     ("fxtract\n\t"
 488      : "=t" (__junk), "=u" (__value) : "0" (__x));
 489
 490   return __value;
 491 }
 492 #endif
 493
 494
 495 #ifdef __USE_MISC
 496 __MATH_INLINE double drem (double __x, double __y);
 497 __MATH_INLINE double
 498 drem (double __x, double __y)
 499 {
 500   register double __value;
 501   __asm __volatile__
 502     ("1:        fprem1\n\t"
 503      "fstsw     %%ax\n\t"
 504      "sahf\n\t"
 505      "jp        1b"
 506      : "=t" (__value) : "0" (__x), "u" (__y) : "ax", "cc");
 507
 508   return __value;
 509 }
 510
 511 /* This function is used in the `isfinite' macro.  */
 512 __MATH_INLINE int __finite (double __x);
 513 __MATH_INLINE int
 514 __finite (double __x)
 515 {
 516   register int __result;
 517   __asm__ __volatile__
 518     ("orl       $0x800fffff, %0\n\t"
 519      "incl      %0\n\t"
 520      "shrl      $31, %0"
 521      : "=q" (__result) : "0" (((int *) &__x)[1]));
 522   return __result;
 523 }
 524
 525 __MATH_INLINE double coshm1 (double __x);
 526 __MATH_INLINE double
 527 coshm1 (double __x)
 528 {
 529   register double __exm1 = __expm1 (fabs (__x));
 530
 531   return 0.5 * (__exm1 / (__exm1 + 1.0)) * __exm1;
 532 }
 533
 534 __MATH_INLINE double acosh1p (double __x);
 535 __MATH_INLINE double
 536 acosh1p (double __x)
 537 {
 538   return __log1p (__x + sqrt (__x) * sqrt (__x + 2.0));
 539 }
 540
 541 __MATH_INLINE void sincos (double __x, double *__sinx, double *__cosx);
 542 __MATH_INLINE void
 543 sincos (double __x, double *__sinx, double *__cosx)
 544 {
 545   register double __cosr, __sinr;
 546   __asm __volatile__
 547     ("fsincos\n\t"
 548      "fnstsw    %%ax\n\t"
 549      "testl     $0x400, %%eax\n\t"
 550      "jz        1f\n\t"
 551      "fldpi\n\t"
 552      "fadd      %%st(0)\n\t"
 553      "fxch      %%st(1)\n\t"
 554      "2: fprem1\n\t"
 555      "fnstsw    %%ax\n\t"
 556      "testl     $0x400, %%eax\n\t"
 557      "jnz       2b\n\t"
 558      "fstp      %%st(1)\n\t"
 559      "fsincos\n\t"
 560      "1:"
 561      : "=t" (__cosr), "=u" (__sinr) : "0" (__x));
 562
 563   *__sinx = __sinr;
 564   *__cosx = __cosr;
 565 }
 566
 567 __MATH_INLINE double sgn (double __x);
 568 __MATH_INLINE double
 569 sgn (double __x)
 570 {
 571   return __x == 0.0 ? 0.0 : (__x > 0.0 ? 1.0 : -1.0);
 572 }
 573
 574 __MATH_INLINE double pow2 (double __x);
 575 __MATH_INLINE double
 576 pow2 (double __x)
 577 {
 578   register double __value, __exponent;
 579   long __p = (long) __x;
 580
 581   if (__x == (double) __p)
 582     return ldexp (1.0, __p);
 583
 584   __asm __volatile__
 585     ("fldl      %%st(0)\n\t"
 586      "frndint                   # int(x)\n\t"
 587      "fxch\n\t"
 588      "fsub      %%st(1)         # fract(x)\n\t"
 589      "f2xm1                     # 2^(fract(x)) - 1\n\t"
 590      : "=t" (__value), "=u" (__exponent) : "0" (__x));
 591   __value += 1.0;
 592   __asm __volatile__
 593     ("fscale"
 594      : "=t" (__value) : "0" (__value), "u" (__exponent));
 595
 596   return __value;
 597 }
 598
 599 #endif /* __USE_MISC  */
 600
 601 #endif /* __NO_MATH_INLINES  */
 602 #endif /* __GNUC__  */
 603
 604 #endif /* _BITS_MATHINLINE_H  */