diff options
Diffstat (limited to 'sysdeps/i386/fpu/bits/mathinline.h')
| -rw-r--r-- | sysdeps/i386/fpu/bits/mathinline.h | 608 |
1 files changed, 608 insertions, 0 deletions
diff --git a/sysdeps/i386/fpu/bits/mathinline.h b/sysdeps/i386/fpu/bits/mathinline.h new file mode 100644 index 0000000000..d5c480a51a --- /dev/null +++ b/sysdeps/i386/fpu/bits/mathinline.h @@ -0,0 +1,608 @@ +/* Inline math functions for i387. + Copyright (C) 1995, 1996, 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by John C. Bowman <bowman@ipp-garching.mpg.de>, 1995. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#ifndef _BITS_MATHINLINE_H +#define _BITS_MATHINLINE_H 1 + +#if defined __GNUG__ && \ + (__GNUC__ < 2 || (__GNUC__ == 2 && __GNUC_MINOR__ <= 7)) +/* gcc 2.7.2 and 2.7.2.1 have problems with inlining `long double' + functions so we disable this now. */ +#undef __NO_MATH_INLINES +#define __NO_MATH_INLINES +#endif + +#ifdef __GNUC__ +#ifndef __NO_MATH_INLINES + +#ifdef __cplusplus +#define __MATH_INLINE __inline +#else +#define __MATH_INLINE extern __inline +#endif + +__MATH_INLINE double cos (double); +__MATH_INLINE double sin (double); + + +__MATH_INLINE double __expm1 (double __x); +__MATH_INLINE double +__expm1 (double __x) +{ + register double __value, __exponent, __temp; + __asm __volatile__ + ("fldl2e # e^x - 1 = 2^(x * log2(e)) - 1\n\t" + "fmul %%st(1) # x * log2(e)\n\t" + "fstl %%st(1)\n\t" + "frndint # int(x * log2(e))\n\t" + "fxch\n\t" + "fsub %%st(1) # fract(x * log2(e))\n\t" + "f2xm1 # 2^(fract(x * log2(e))) - 1\n\t" + "fscale # 2^(x * log2(e)) - 2^(int(x * log2(e)))\n\t" + : "=t" (__value), "=u" (__exponent) : "0" (__x)); + __asm __volatile__ + ("fscale # 2^int(x * log2(e))\n\t" + : "=t" (__temp) : "0" (1.0), "u" (__exponent)); + __temp -= 1.0; + + return __temp + __value; +} + +__MATH_INLINE double __sgn1 (double __x); +__MATH_INLINE double +__sgn1 (double __x) +{ + return __x >= 0.0 ? 1.0 : -1.0; +} + +__MATH_INLINE double sqrt (double __x); +__MATH_INLINE double +sqrt (double __x) +{ + register double __value; + __asm __volatile__ + ("fsqrt" + : "=t" (__value) : "0" (__x)); + + return __value; +} + +__MATH_INLINE double fabs (double __x); +__MATH_INLINE double +fabs (double __x) +{ + register double __value; + __asm __volatile__ + ("fabs" + : "=t" (__value) : "0" (__x)); + + return __value; +} + +/* The argument range of this inline version is limited. */ +__MATH_INLINE double sin (double __x); +__MATH_INLINE double +sin (double __x) +{ + register double __value; + __asm __volatile__ + ("fsin" + : "=t" (__value) : "0" (__x)); + + return __value; +} + +/* The argument range of this inline version is limited. */ +__MATH_INLINE double cos (double __x); +__MATH_INLINE double +cos (double __x) +{ + register double __value; + __asm __volatile__ + ("fcos" + : "=t" (__value): "0" (__x)); + + return __value; +} + +__MATH_INLINE double tan (double __x); +__MATH_INLINE double +tan (double __x) +{ + register double __value; + register double __value2 __attribute__ ((unused)); + __asm __volatile__ + ("fptan" + : "=t" (__value2), "=u" (__value) : "0" (__x)); + + return __value; +} + +__MATH_INLINE double atan2 (double __y, double __x); +__MATH_INLINE double +atan2 (double __y, double __x) +{ + register double __value; + __asm __volatile__ + ("fpatan\n\t" + "fldl %%st(0)" + : "=t" (__value) : "0" (__x), "u" (__y)); + + return __value; +} + +__MATH_INLINE double asin (double __x); +__MATH_INLINE double +asin (double __x) +{ + return atan2 (__x, sqrt (1.0 - __x * __x)); +} + +__MATH_INLINE double acos (double __x); +__MATH_INLINE double +acos (double __x) +{ + return atan2 (sqrt (1.0 - __x * __x), __x); +} + +__MATH_INLINE double atan (double __x); +__MATH_INLINE double +atan (double __x) +{ + register double __value; + __asm __volatile__ + ("fld1\n\t" + "fpatan" + : "=t" (__value) : "0" (__x)); + + return __value; +} + +__MATH_INLINE double exp (double __x); +__MATH_INLINE double +exp (double __x) +{ + register double __value, __exponent; + __asm __volatile__ + ("fldl2e # e^x = 2^(x * log2(e))\n\t" + "fmul %%st(1) # x * log2(e)\n\t" + "fstl %%st(1)\n\t" + "frndint # int(x * log2(e))\n\t" + "fxch\n\t" + "fsub %%st(1) # fract(x * log2(e))\n\t" + "f2xm1 # 2^(fract(x * log2(e))) - 1\n\t" + : "=t" (__value), "=u" (__exponent) : "0" (__x)); + __value += 1.0; + __asm __volatile__ + ("fscale" + : "=t" (__value) : "0" (__value), "u" (__exponent)); + + return __value; +} + +__MATH_INLINE double sinh (double __x); +__MATH_INLINE double +sinh (double __x) +{ + register double __exm1 = __expm1 (fabs (__x)); + + return 0.5 * (__exm1 / (__exm1 + 1.0) + __exm1) * __sgn1 (__x); +} + +__MATH_INLINE double cosh (double __x); +__MATH_INLINE double +cosh (double __x) +{ + register double __ex = exp (__x); + + return 0.5 * (__ex + 1.0 / __ex); +} + +__MATH_INLINE double tanh (double __x); +__MATH_INLINE double +tanh (double __x) +{ + register double __exm1 = __expm1 (-fabs (__x + __x)); + + return __exm1 / (__exm1 + 2.0) * __sgn1 (-__x); +} + +__MATH_INLINE double log (double __x); +__MATH_INLINE double +log (double __x) +{ + register double __value; + __asm __volatile__ + ("fldln2\n\t" + "fxch\n\t" + "fyl2x" + : "=t" (__value) : "0" (__x)); + + return __value; +} + +__MATH_INLINE double log10 (double __x); +__MATH_INLINE double +log10 (double __x) +{ + register double __value; + __asm __volatile__ + ("fldlg2\n\t" + "fxch\n\t" + "fyl2x" + : "=t" (__value) : "0" (__x)); + + return __value; +} + +__MATH_INLINE double __log2 (double __x); +__MATH_INLINE double +__log2 (double __x) +{ + register double __value; + __asm __volatile__ + ("fld1\n\t" + "fxch\n\t" + "fyl2x" + : "=t" (__value) : "0" (__x)); + + return __value; +} + +__MATH_INLINE double fmod (double __x, double __y); +__MATH_INLINE double +fmod (double __x, double __y) +{ + register double __value; + __asm __volatile__ + ("1: fprem\n\t" + "fstsw %%ax\n\t" + "sahf\n\t" + "jp 1b" + : "=t" (__value) : "0" (__x), "u" (__y) : "ax", "cc"); + + return __value; +} + +__MATH_INLINE double ldexp (double __x, int __y); +__MATH_INLINE double +ldexp (double __x, int __y) +{ + register double __value; + __asm __volatile__ + ("fscale" + : "=t" (__value) : "0" (__x), "u" ((double) __y)); + + return __value; +} + +__MATH_INLINE double pow (double __x, double __y); +__MATH_INLINE double +pow (double __x, double __y) +{ + register double __value, __exponent; + long __p = (long) __y; + + if (__x == 0.0 && __y > 0.0) + return 0.0; + if (__y == (double) __p) + { + double __r = 1.0; + if (__p == 0) + return 1.0; + if (__p < 0) + { + __p = -__p; + __x = 1.0 / __x; + } + while (1) + { + if (__p & 1) + __r *= __x; + __p >>= 1; + if (__p == 0) + return __r; + __x *= __x; + } + /* NOTREACHED */ + } + __asm __volatile__ + ("fmul %%st(1) # y * log2(x)\n\t" + "fstl %%st(1)\n\t" + "frndint # int(y * log2(x))\n\t" + "fxch\n\t" + "fsub %%st(1) # fract(y * log2(x))\n\t" + "f2xm1 # 2^(fract(y * log2(x))) - 1\n\t" + : "=t" (__value), "=u" (__exponent) : "0" (__log2 (__x)), "1" (__y)); + __value += 1.0; + __asm __volatile__ + ("fscale" + : "=t" (__value) : "0" (__value), "u" (__exponent)); + + return __value; +} + +__MATH_INLINE double floor (double __x); +__MATH_INLINE double +floor (double __x) +{ + register double __value; + __volatile unsigned short int __cw, __cwtmp; + + __asm __volatile ("fnstcw %0" : "=m" (__cw)); + __cwtmp = (__cw & 0xf3ff) | 0x0400; /* rounding down */ + __asm __volatile ("fldcw %0" : : "m" (__cwtmp)); + __asm __volatile ("frndint" : "=t" (__value) : "0" (__x)); + __asm __volatile ("fldcw %0" : : "m" (__cw)); + + return __value; +} + +__MATH_INLINE double ceil (double __x); +__MATH_INLINE double +ceil (double __x) +{ + register double __value; + __volatile unsigned short int __cw, __cwtmp; + + __asm __volatile ("fnstcw %0" : "=m" (__cw)); + __cwtmp = (__cw & 0xf3ff) | 0x0800; /* rounding up */ + __asm __volatile ("fldcw %0" : : "m" (__cwtmp)); + __asm __volatile ("frndint" : "=t" (__value) : "0" (__x)); + __asm __volatile ("fldcw %0" : : "m" (__cw)); + + return __value; +} + + +/* Optimized versions for some non-standardized functions. */ +#if defined __USE_ISOC9X || defined __USE_MISC + +__MATH_INLINE double hypot (double __x, double __y); +__MATH_INLINE double +hypot (double __x, double __y) +{ + return sqrt (__x * __x + __y * __y); +} + +/* We cannot rely on M_SQRT being defined. So we do it for ourself + here. */ +#define __M_SQRT2 _Mldbl(1.41421356237309504880) /* sqrt(2) */ + +__MATH_INLINE double log1p (double __x); +__MATH_INLINE double +log1p (double __x) +{ + register double __value; + + if (fabs (__x) >= 1.0 - 0.5 * __M_SQRT2) + __value = log (1.0 + __x); + else + __asm __volatile__ + ("fldln2\n\t" + "fxch\n\t" + "fyl2xp1" + : "=t" (__value) : "0" (__x)); + + return __value; +} + +__MATH_INLINE double asinh (double __x); +__MATH_INLINE double +asinh (double __x) +{ + register double __y = fabs (__x); + + return log1p ((__y * __y / (sqrt (__y * __y + 1.0) + 1.0) + __y) + * __sgn1 (__x)); +} + +__MATH_INLINE double acosh (double __x); +__MATH_INLINE double +acosh (double __x) +{ + return log (__x + sqrt (__x - 1.0) * sqrt (__x + 1.0)); +} + +__MATH_INLINE double atanh (double __x); +__MATH_INLINE double +atanh (double __x) +{ + register double __y = fabs (__x); + + return -0.5 * __log1p (-(__y + __y) / (1.0 + __y)) * __sgn1 (__x); +} + +__MATH_INLINE double logb (double __x); +__MATH_INLINE double +logb (double __x) +{ + register double __value, __junk; + __asm __volatile__ + ("fxtract\n\t" + : "=t" (__junk), "=u" (__value) : "0" (__x)); + + return __value; +} + +__MATH_INLINE double drem (double __x, double __y); +__MATH_INLINE double +drem (double __x, double __y) +{ + register double __value; + __asm __volatile__ + ("1: fprem1\n\t" + "fstsw %%ax\n\t" + "sahf\n\t" + "jp 1b" + : "=t" (__value) : "0" (__x), "u" (__y) : "ax", "cc"); + + return __value; +} + +/* This function is used in the `isfinite' macro. */ +__MATH_INLINE int __finite (double __x); +__MATH_INLINE int +__finite (double __x) +{ + register int __result; + __asm__ __volatile__ + ("orl $0x800fffff, %0\n\t" + "incl %0\n\t" + "shrl $31, %0" + : "=q" (__result) : "0" (((int *) &__x)[1])); + return __result; +} + + +/* ISO C 9X defines some macros to perform unordered comparisons. The + ix87 FPU supports this with special opcodes and we should use them. + These must not be inline functions since we have to be able to handle + all floating-point types. */ +#undef isgreater +#define isgreater(x, y) \ + ({ int result; \ + __asm__ ("fucompp; fnstsw; andb $0x45, %%ah; setz %%al;" \ + "andl $0xff, %0" \ + : "=a" (result) : "t" (x), "u" (y) : "cc"); \ + result; }) + +#undef isgreaterequal +#define isgreaterequal(x, y) \ + ({ int result; \ + __asm__ ("fucompp; fnstsw; testb $0x05, %%ah; setz %%al;" \ + "andl $0xff, %0" \ + : "=a" (result) : "t" (x), "u" (y) : "cc"); \ + result; }) + +#undef isless +#define isless(x, y) \ + ({ int result; \ + __asm__ ("fucompp; fnstsw; xorb $0x01, %%ah; testb $0x45, %%ah;" \ + "setz %%al; andl $0xff, %0" \ + : "=a" (result) : "t" (x), "u" (y) : "cc"); \ + result; }) + +#undef islessequal +#define islessequal(x, y) \ + ({ int result; \ + __asm__ ("fucompp; fnstsw; xorb $0x01, %%ah; testb $0x05, %%ah;" \ + "setz %%al; andl $0xff, %0" \ + : "=a" (result) : "t" (x), "u" (y) : "cc"); \ + result; }) + +#undef islessgreater +#define islessgreater(x, y) \ + ({ int result; \ + __asm__ ("fucompp; fnstsw; testb $0x44, %%ah; setz %%al;" \ + "andl $0xff, %0" \ + : "=a" (result) : "t" (x), "u" (y) : "cc"); \ + result; }) + +#undef isunordered +#define isunordered(x, y) \ + ({ int result; \ + __asm__ ("fucompp; fnstsw; sahf; setp %%al; andl $0xff, %0" \ + : "=a" (result) : "t" (x), "u" (y) : "cc"); \ + result; }) +#endif + + +#ifdef __USE_MISC +__MATH_INLINE double coshm1 (double __x); +__MATH_INLINE double +coshm1 (double __x) +{ + register double __exm1 = __expm1 (fabs (__x)); + + return 0.5 * (__exm1 / (__exm1 + 1.0)) * __exm1; +} + +__MATH_INLINE double acosh1p (double __x); +__MATH_INLINE double +acosh1p (double __x) +{ + return __log1p (__x + sqrt (__x) * sqrt (__x + 2.0)); +} + +__MATH_INLINE void sincos (double __x, double *__sinx, double *__cosx); +__MATH_INLINE void +sincos (double __x, double *__sinx, double *__cosx) +{ + register double __cosr, __sinr; + __asm __volatile__ + ("fsincos\n\t" + "fnstsw %%ax\n\t" + "testl $0x400, %%eax\n\t" + "jz 1f\n\t" + "fldpi\n\t" + "fadd %%st(0)\n\t" + "fxch %%st(1)\n\t" + "2: fprem1\n\t" + "fnstsw %%ax\n\t" + "testl $0x400, %%eax\n\t" + "jnz 2b\n\t" + "fstp %%st(1)\n\t" + "fsincos\n\t" + "1:" + : "=t" (__cosr), "=u" (__sinr) : "0" (__x)); + + *__sinx = __sinr; + *__cosx = __cosr; +} + +__MATH_INLINE double sgn (double __x); +__MATH_INLINE double +sgn (double __x) +{ + return __x == 0.0 ? 0.0 : (__x > 0.0 ? 1.0 : -1.0); +} + +__MATH_INLINE double pow2 (double __x); +__MATH_INLINE double +pow2 (double __x) +{ + register double __value, __exponent; + long __p = (long) __x; + + if (__x == (double) __p) + return ldexp (1.0, __p); + + __asm __volatile__ + ("fldl %%st(0)\n\t" + "frndint # int(x)\n\t" + "fxch\n\t" + "fsub %%st(1) # fract(x)\n\t" + "f2xm1 # 2^(fract(x)) - 1\n\t" + : "=t" (__value), "=u" (__exponent) : "0" (__x)); + __value += 1.0; + __asm __volatile__ + ("fscale" + : "=t" (__value) : "0" (__value), "u" (__exponent)); + + return __value; +} + +#endif /* __USE_MISC */ + +#endif /* __NO_MATH_INLINES */ +#endif /* __GNUC__ */ + +#endif /* _BITS_MATHINLINE_H */ |