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 __MATH_INLINE double log1p (double __x);
 386 __MATH_INLINE double
 387 log1p (double __x)
 388 {
 389   register double __value;
 390
 391   if (fabs (__x) >= 1.0 - 0.5 * M_SQRT2)
 392     __value = log (1.0 + __x);
 393   else
 394     __asm __volatile__
 395       ("fldln2\n\t"
 396        "fxch\n\t"
 397        "fyl2xp1"
 398        : "=t" (__value) : "0" (__x));
 399
 400   return __value;
 401 }
 402
 403 __MATH_INLINE double asinh (double __x);
 404 __MATH_INLINE double
 405 asinh (double __x)
 406 {
 407   register double __y = fabs (__x);
 408
 409   return log1p ((__y * __y / (sqrt (__y * __y + 1.0) + 1.0) + __y)
 410                 * __sgn1 (__x));
 411 }
 412
 413 __MATH_INLINE double acosh (double __x);
 414 __MATH_INLINE double
 415 acosh (double __x)
 416 {
 417   return log (__x + sqrt (__x - 1.0) * sqrt (__x + 1.0));
 418 }
 419
 420 __MATH_INLINE double atanh (double __x);
 421 __MATH_INLINE double
 422 atanh (double __x)
 423 {
 424   register double __y = fabs (__x);
 425
 426   return -0.5 * __log1p (-(__y + __y) / (1.0 + __y)) * __sgn1 (__x);
 427 }
 428
 429 __MATH_INLINE double logb (double __x);
 430 __MATH_INLINE double
 431 logb (double __x)
 432 {
 433   register double __value;
 434   __asm __volatile__
 435     ("fxtract\n\t"
 436      : "=t" (__value) : "0" (__x));
 437
 438   return __value;
 439 }
 440
 441 __MATH_INLINE double drem (double __x, double __y);
 442 __MATH_INLINE double
 443 drem (double __x, double __y)
 444 {
 445   register double __value;
 446   __asm __volatile__
 447     ("1:        fprem1\n\t"
 448      "fstsw     %%ax\n\t"
 449      "sahf\n\t"
 450      "jp        1b"
 451      : "=t" (__value) : "0" (__x), "u" (__y) : "ax", "cc");
 452
 453   return __value;
 454 }
 455
 456 /* This function is used in the `isfinite' macro.  */
 457 __MATH_INLINE int __finite (double __x);
 458 __MATH_INLINE int
 459 __finite (double __x)
 460 {
 461   register int __result;
 462   __asm__ __volatile__
 463     ("orl       $x0x800fffff, %0\n\t"
 464      "incl      %0\n\t"
 465      "shrl      $31, %0"
 466      : "=q" (__result) : "0" (((int *) &__x)[1]));
 467   return __result;
 468 }
 469 #endif
 470
 471 #ifdef __USE_MISC
 472 __MATH_INLINE double coshm1 (double __x);
 473 __MATH_INLINE double
 474 coshm1 (double __x)
 475 {
 476   register double __exm1 = __expm1 (fabs (__x));
 477
 478   return 0.5 * (__exm1 / (__exm1 + 1.0)) * __exm1;
 479 }
 480
 481 __MATH_INLINE double acosh1p (double __x);
 482 __MATH_INLINE double
 483 acosh1p (double __x)
 484 {
 485   return __log1p (__x + sqrt (__x) * sqrt (__x + 2.0));
 486 }
 487
 488 __MATH_INLINE void sincos (double __x, double *__sinx, double *__cosx);
 489 __MATH_INLINE void
 490 sincos (double __x, double *__sinx, double *__cosx)
 491 {
 492   register double __cosr, __sinr;
 493   __asm __volatile__
 494     ("fsincos\n\t"
 495      "fnstsw    %%ax\n\t"
 496      "testl     $0x400, %%eax\n\t"
 497      "jz        1f\n\t"
 498      "fldpi\n\t"
 499      "fadd      %%st(0)\n\t"
 500      "fxch      %%st(1)\n\t"
 501      "2: fprem1\n\t"
 502      "fnstsw    %%ax\n\t"
 503      "testl     $0x400, %%eax\n\t"
 504      "jnz       2b\n\t"
 505      "fstp      %%st(1)\n\t"
 506      "fsincos\n\t"
 507      "1:"
 508      : "=t" (__cosr), "=u" (__sinr) : "0" (__x));
 509
 510   *__sinx = __sinr;
 511   *__cosx = __cosr;
 512 }
 513
 514 __MATH_INLINE double sgn (double __x);
 515 __MATH_INLINE double
 516 sgn (double __x)
 517 {
 518   return __x == 0.0 ? 0.0 : (__x > 0.0 ? 1.0 : -1.0);
 519 }
 520
 521 __MATH_INLINE double pow2 (double __x);
 522 __MATH_INLINE double
 523 pow2 (double __x)
 524 {
 525   register double __value, __exponent;
 526   long __p = (long) __x;
 527
 528   if (__x == (double) __p)
 529     return ldexp (1.0, __p);
 530
 531   __asm __volatile__
 532     ("fldl      %%st(0)\n\t"
 533      "frndint                   # int(x)\n\t"
 534      "fxch\n\t"
 535      "fsub      %%st(1)         # fract(x)\n\t"
 536      "f2xm1                     # 2^(fract(x)) - 1\n\t"
 537      : "=t" (__value), "=u" (__exponent) : "0" (__x));
 538   __value += 1.0;
 539   __asm __volatile__
 540     ("fscale"
 541      : "=t" (__value) : "0" (__value), "u" (__exponent));
 542
 543   return __value;
 544 }
 545
 546 #endif /* __USE_MISC  */
 547
 548 #endif /* __NO_MATH_INLINES  */
 549 #endif /* __GNUC__  */
 550
 551 #endif /* __MATH_H  */