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