diff options
author | Danny Smith <dannysmith@users.sourceforge.net> | 2004-10-07 00:31:32 +0400 |
---|---|---|
committer | Danny Smith <dannysmith@users.sourceforge.net> | 2004-10-07 00:31:32 +0400 |
commit | 72db1c11e9887439125ce798bb1b6d571fe7d840 (patch) | |
tree | 9ad7335a58ca086ea29b788d867728ea23fb7cc9 /winsup | |
parent | 74b2bdc351dcbab192faca58b5544bce4f92b973 (diff) |
* include/math.h (ashinh, asinhf, asinhl, acosh, acoshf, acoshl,
atanh, atanhf, atanhl): Add prototypes.
* mingwex/Makefile.in (MATH_OBJS): Add objects for above to list.
(MATH_DISTFILES): Add sources for above and fastmath.h to list.
Specify dependency on fastmath.h for new objects.
* mingwex/math/fastmath.h: New file.
* mingwex/math/ashinh.c: New file.
* mingwex/math/asinhf.c: New file.
* mingwex/math/asinhl.c: New file.
* mingwex/math/acosh.c: New file.
* mingwex/math/acoshf.c: New file.
* mingwex/math/acoshl.c: New file.
* mingwex/math/atanh.c: New file.
* mingwex/math/atanhf.c: New file.
* mingwex/math/atanhl.c: New file.
Diffstat (limited to 'winsup')
-rw-r--r-- | winsup/mingw/ChangeLog | 17 | ||||
-rw-r--r-- | winsup/mingw/include/math.h | 21 | ||||
-rw-r--r-- | winsup/mingw/mingwex/Makefile.in | 12 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/acosh.c | 26 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/acoshf.c | 25 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/acoshl.c | 27 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/asinh.c | 28 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/asinhf.c | 28 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/asinhl.c | 28 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/atanh.c | 31 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/atanhf.c | 30 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/atanhl.c | 29 | ||||
-rwxr-xr-x | winsup/mingw/mingwex/math/fastmath.h | 115 |
13 files changed, 411 insertions, 6 deletions
diff --git a/winsup/mingw/ChangeLog b/winsup/mingw/ChangeLog index 94c6aec72..5e2a6b644 100644 --- a/winsup/mingw/ChangeLog +++ b/winsup/mingw/ChangeLog @@ -1,3 +1,20 @@ +2004-10-07 Danny Smith <dannysmith@users.sourceforge.net> + + * mingwex/math/fastmath.h: New file. + * mingwex/math/asinh.c: New file. + * mingwex/math/asinhf.c: New file. + * mingwex/math/asinhl.c: New file. + * mingwex/math/acosh.c: New file. + * mingwex/math/acoshf.c: New file. + * mingwex/math/acoshl.c: New file. + * mingwex/math/atanh.c: New file. + * mingwex/math/atanhf.c: New file. + * include/math.h (asinh, asinhf, asinhl, acosh, acoshf, acoshl, + atanh, atanhf, atanhl): Add prototypes. + * mingwex/Makefile.in (MATH_OBJS): Add objects for above to list. + (MATH_DISTFILES): Add sources for above and fastmath.h to list. + Specify dependency on fastmath.h for new objects. + 2004-09-08 Earnie Boyd <earnie@users.sf.net> * include/sys/stat.h (_S_IFLNK): Add definition. diff --git a/winsup/mingw/include/math.h b/winsup/mingw/include/math.h index abcd6a4b8..a05fa841b 100644 --- a/winsup/mingw/include/math.h +++ b/winsup/mingw/include/math.h @@ -427,10 +427,23 @@ __CRT_INLINE float __cdecl tanhf (float x) {return (float) tanh (x);} extern long double __cdecl tanhl (long double); -/* - * TODO: asinh, acosh, atanh - */ - +/* Inverse hyperbolic trig functions */ +/* 7.12.5.1 */ +extern double __cdecl acosh (double); +extern float __cdecl acoshf (float); +extern long double __cdecl acoshl (long double); + +/* 7.12.5.2 */ +extern double __cdecl asinh (double); +extern float __cdecl asinhf (float); +extern long double __cdecl asinhl (long double); + +/* 7.12.5.3 */ +extern double __cdecl atanh (double); +extern float __cdecl atanf (float); +extern long double __cdecl atanhl (long double); + +/* Exponentials and logarithms */ /* 7.12.6.1 Double in C89 */ __CRT_INLINE float __cdecl expf (float x) {return (float) exp (x);} diff --git a/winsup/mingw/mingwex/Makefile.in b/winsup/mingw/mingwex/Makefile.in index 38a02b929..0dcf1ed3b 100644 --- a/winsup/mingw/mingwex/Makefile.in +++ b/winsup/mingw/mingwex/Makefile.in @@ -60,7 +60,9 @@ MATH_DISTFILES = \ roundl.c scalbn.S scalbnf.S scalbnl.S s_erf.c sf_erf.c \ signbit.c signbitf.c signbitl.c sinf.S sinhf.c sinhl.c sinl.S \ sqrtf.c sqrtl.c tanf.S tanhf.c tanhl.c tanl.S tgamma.c \ - tgammaf.c tgammal.c trunc.c truncf.c truncl.c + tgammaf.c tgammal.c trunc.c truncf.c truncl.c \ + acosh.c acoshf.c acoshl.c asinh.c asinhf.c asinhl.c \ + atanh.c atanhf.c atanhl.c fastmath.h STDIO_DISTFILES = \ fopen64.c fseeko64.c ftello64.c lseek64.c \ @@ -142,7 +144,9 @@ MATH_OBJS = \ roundl.o scalbn.o scalbnf.o scalbnl.o s_erf.o sf_erf.o \ signbit.o signbitf.o signbitl.o sinf.o sinhf.o sinhl.o sinl.o \ sqrtf.o sqrtl.o tanf.o tanhf.o tanhl.o tanl.o tgamma.o \ - tgammaf.o tgammal.o trunc.o truncf.o truncl.o + tgammaf.o tgammal.o trunc.o truncf.o truncl.o \ + acosh.o acoshf.o acoshl.o asinh.o asinhf.o asinhl.o \ + atanh.o atanhf.o atanhl.o FENV_OBJS = fesetround.o fegetround.o \ fegetenv.o fesetenv.o feupdateenv.o \ feclearexcept.o feholdexcept.o fegetexceptflag.o \ @@ -211,6 +215,10 @@ wdirent.o: $(srcdir)/dirent.c $(srcdir)/wdirent.c strtold.o: $(srcdir)/strtold.c $(srcdir)/math/cephes_emath.h wcstold.o: $(srcdir)/wcstold.c $(srcdir)/math/cephes_emath.h +acosh.o acoshf.o acoshl.o \ +asinh.o asinhf.o asinhl.o \ +atanh.o atanhf.o atanhl.o: fastmath.h + dist: mkdir $(distdir)/mingwex diff --git a/winsup/mingw/mingwex/math/acosh.c b/winsup/mingw/mingwex/math/acosh.c new file mode 100755 index 000000000..1497883cf --- /dev/null +++ b/winsup/mingw/mingwex/math/acosh.c @@ -0,0 +1,26 @@ +#include <math.h> +#include <errno.h> +#include "fastmath.h" + +/* acosh(x) = log (x + sqrt(x * x - 1)) */ +double acosh (double x) +{ + if (isnan (x)) + return x; + + if (x < 1.0) + { + errno = EDOM; + return nan(""); + } + + if (x > 0x1p32) + /* Avoid overflow (and unnecessary calculation when + sqrt (x * x - 1) == x). GCC optimizes by replacing + the long double M_LN2 const with a fldln2 insn. */ + return __fast_log (x) + 6.9314718055994530941723E-1L; + + /* Since x >= 1, the arg to log will always be greater than + the fyl2xp1 limit (approx 0.29) so just use logl. */ + return __fast_log (x + __fast_sqrt((x + 1.0) * (x - 1.0))); +} diff --git a/winsup/mingw/mingwex/math/acoshf.c b/winsup/mingw/mingwex/math/acoshf.c new file mode 100755 index 000000000..08f190fcb --- /dev/null +++ b/winsup/mingw/mingwex/math/acoshf.c @@ -0,0 +1,25 @@ +#include <math.h> +#include <errno.h> +#include "fastmath.h" + +/* acosh(x) = log (x + sqrt(x * x - 1)) */ +float acoshf (float x) +{ + if (isnan (x)) + return x; + if (x < 1.0f) + { + errno = EDOM; + return nan(""); + } + + if (x > 0x1p32f) + /* Avoid overflow (and unnecessary calculation when + sqrt (x * x - 1) == x). GCC optimizes by replacing + the long double M_LN2 const with a fldln2 insn. */ + return __fast_log (x) + 6.9314718055994530941723E-1L; + + /* Since x >= 1, the arg to log will always be greater than + the fyl2xp1 limit (approx 0.29) so just use logl. */ + return __fast_log (x + __fast_sqrt((x + 1.0) * (x - 1.0))); +} diff --git a/winsup/mingw/mingwex/math/acoshl.c b/winsup/mingw/mingwex/math/acoshl.c new file mode 100755 index 000000000..c461176bb --- /dev/null +++ b/winsup/mingw/mingwex/math/acoshl.c @@ -0,0 +1,27 @@ +#include <math.h> +#include <errno.h> +#include "fastmath.h" + +/* acosh(x) = log (x + sqrt(x * x - 1)) */ +long double acoshl (long double x) +{ + if (isnan (x)) + return x; + + if (x < 1.0L) + { + errno = EDOM; + return nanl(""); + } + if (x > 0x1p32L) + /* Avoid overflow (and unnecessary calculation when + sqrt (x * x - 1) == x). + The M_LN2 define doesn't have enough precison for + long double so use this one. GCC optimizes by replacing + the const with a fldln2 insn. */ + return __fast_logl (x) + 6.9314718055994530941723E-1L; + + /* Since x >= 1, the arg to log will always be greater than + the fyl2xp1 limit (approx 0.29) so just use logl. */ + return __fast_logl (x + __fast_sqrtl((x + 1.0L) * (x - 1.0L))); +} diff --git a/winsup/mingw/mingwex/math/asinh.c b/winsup/mingw/mingwex/math/asinh.c new file mode 100755 index 000000000..30404497d --- /dev/null +++ b/winsup/mingw/mingwex/math/asinh.c @@ -0,0 +1,28 @@ +#include <math.h> +#include <errno.h> +#include "fastmath.h" + + /* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */ +double asinh(double x) +{ + double z; + if (!isfinite (x)) + return x; + z = fabs (x); + + /* Avoid setting FPU underflow exception flag in x * x. */ +#if 0 + if ( z < 0x1p-32) + return x; +#endif + + /* Use log1p to avoid cancellation with small x. Put + x * x in denom, so overflow is harmless. + asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0) + = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0)) */ + + z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0)); + + return ( x > 0.0 ? z : -z); +} + diff --git a/winsup/mingw/mingwex/math/asinhf.c b/winsup/mingw/mingwex/math/asinhf.c new file mode 100755 index 000000000..080a9278d --- /dev/null +++ b/winsup/mingw/mingwex/math/asinhf.c @@ -0,0 +1,28 @@ +#include <math.h> +#include <errno.h> +#include "fastmath.h" + + /* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */ +float asinhf(float x) +{ + float z; + if (!isfinite (x)) + return x; + z = fabsf (x); + + /* Avoid setting FPU underflow exception flag in x * x. */ +#if 0 + if ( z < 0x1p-32) + return x; +#endif + + + /* Use log1p to avoid cancellation with small x. Put + x * x in denom, so overflow is harmless. + asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0) + = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0)) */ + + z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0)); + + return ( x > 0.0 ? z : -z); +} diff --git a/winsup/mingw/mingwex/math/asinhl.c b/winsup/mingw/mingwex/math/asinhl.c new file mode 100755 index 000000000..8f027e83d --- /dev/null +++ b/winsup/mingw/mingwex/math/asinhl.c @@ -0,0 +1,28 @@ +#include <math.h> +#include <errno.h> +#include "fastmath.h" + + /* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */ +long double asinhl(long double x) +{ + long double z; + if (!isfinite (x)) + return x; + + z = fabsl (x); + + /* Avoid setting FPU underflow exception flag in x * x. */ +#if 0 + if ( z < 0x1p-32) + return x; +#endif + + /* Use log1p to avoid cancellation with small x. Put + x * x in denom, so overflow is harmless. + asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0) + = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0)) */ + + z = __fast_log1pl (z + z * z / (__fast_sqrtl (z * z + 1.0L) + 1.0L)); + + return ( x > 0.0 ? z : -z); +} diff --git a/winsup/mingw/mingwex/math/atanh.c b/winsup/mingw/mingwex/math/atanh.c new file mode 100755 index 000000000..b5d9ce100 --- /dev/null +++ b/winsup/mingw/mingwex/math/atanh.c @@ -0,0 +1,31 @@ +#include <math.h> +#include <errno.h> +#include "fastmath.h" + +/* atanh (x) = 0.5 * log ((1.0 + x)/(1.0 - x)) */ + +double atanh(double x) +{ + double z; + if isnan (x) + return x; + z = fabs (x); + if (z == 1.0) + { + errno = ERANGE; + return (x > 0 ? INFINITY : -INFINITY); + } + if (z > 1.0) + { + errno = EDOM; + return nan(""); + } + /* Rearrange formula to avoid precision loss for small x. + + atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x)) + = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0) + = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x)) + = 0.5 * log1p ((2.0 * x ) / (1.0 - x)) */ + z = 0.5 * __fast_log1p ((z + z) / (1.0 - z)); + return x >= 0 ? z : -z; +} diff --git a/winsup/mingw/mingwex/math/atanhf.c b/winsup/mingw/mingwex/math/atanhf.c new file mode 100755 index 000000000..b7c30823e --- /dev/null +++ b/winsup/mingw/mingwex/math/atanhf.c @@ -0,0 +1,30 @@ +#include <math.h> +#include <errno.h> +#include "fastmath.h" + +/* atanh (x) = 0.5 * log ((1.0 + x)/(1.0 - x)) */ +float atanhf (float x) +{ + float z; + if isnan (x) + return x; + z = fabsf (x); + if (z == 1.0) + { + errno = ERANGE; + return (x > 0 ? INFINITY : -INFINITY); + } + if ( z > 1.0) + { + errno = EDOM; + return nanf(""); + } + /* Rearrange formula to avoid precision loss for small x. + + atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x)) + = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0) + = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x)) + = 0.5 * log1p ((2.0 * x ) / (1.0 - x)) */ + z = 0.5 * __fast_log1p ((z + z) / (1.0 - z)); + return x >= 0 ? z : -z; +} diff --git a/winsup/mingw/mingwex/math/atanhl.c b/winsup/mingw/mingwex/math/atanhl.c new file mode 100755 index 000000000..2d5fec02a --- /dev/null +++ b/winsup/mingw/mingwex/math/atanhl.c @@ -0,0 +1,29 @@ +#include <math.h> +#include <errno.h> +#include "fastmath.h" + +/* atanh (x) = 0.5 * log ((1.0 + x)/(1.0 - x)) */ +long double atanhl (long double x) +{ + long double z; + if isnan (x) + return x; + z = fabsl (x); + if (z == 1.0L) + { + errno = ERANGE; + return (x > 0 ? INFINITY : -INFINITY); + } + if ( z > 1.0L) + { + errno = EDOM; + return nanl(""); + } + /* Rearrange formula to avoid precision loss for small x. + atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x)) + = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0) + = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x)) + = 0.5 * log1p ((2.0 * x ) / (1.0 - x)) */ + z = 0.5L * __fast_log1pl ((z + z) / (1.0L - z)); + return x >= 0 ? z : -z; +} diff --git a/winsup/mingw/mingwex/math/fastmath.h b/winsup/mingw/mingwex/math/fastmath.h new file mode 100755 index 000000000..01b06b3eb --- /dev/null +++ b/winsup/mingw/mingwex/math/fastmath.h @@ -0,0 +1,115 @@ +#ifndef _MINGWEX_FASTMATH_H_ +#define _MINGWEX_FASTMATH_H_ + +/* Fast math inlines + No range or domain checks. No setting of errno. No tweaks to + protect precision near range limits. */ + +/* For now this is an internal header with just the functions that + are currently used in building libmingwex.a math components */ + +/* FIXME: We really should get rid of the code duplication using euther + C++ templates or tgmath-type macros. */ + +static __inline__ double __fast_sqrt (double x) +{ + double res; + asm __volatile__ ("fsqrt" : "=t" (res) : "0" (x)); + return res; +} + +static __inline__ long double __fast_sqrtl (long double x) +{ + long double res; + asm __volatile__ ("fsqrt" : "=t" (res) : "0" (x)); + return res; +} + +static __inline__ float __fast_sqrtf (float x) +{ + float res; + asm __volatile__ ("fsqrt" : "=t" (res) : "0" (x)); + return res; +} + + +static __inline__ double __fast_log (double x) +{ + double res; + asm __volatile__ + ("fldln2\n\t" + "fxch\n\t" + "fyl2x" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} + +static __inline__ long double __fast_logl (long double x) +{ + long double res; + asm __volatile__ + ("fldln2\n\t" + "fxch\n\t" + "fyl2x" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} + + +static __inline__ float __fast_logf (float x) +{ + float res; + asm __volatile__ + ("fldln2\n\t" + "fxch\n\t" + "fyl2x" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} + +static __inline__ double __fast_log1p (double x) +{ + double res; + /* fyl2xp1 accurate only for |x| <= 1.0 - 0.5 * sqrt (2.0) */ + if (fabs (x) >= 1.0 - 0.5 * 1.41421356237309504880) + res = __fast_log (1.0 + x); + else + asm __volatile__ + ("fldln2\n\t" + "fxch\n\t" + "fyl2xp1" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} + +static __inline__ long double __fast_log1pl (long double x) +{ + long double res; + /* fyl2xp1 accurate only for |x| <= 1.0 - 0.5 * sqrt (2.0) */ + if (fabsl (x) >= 1.0L - 0.5L * 1.41421356237309504880L) + res = __fast_logl (1.0L + x); + else + asm __volatile__ + ("fldln2\n\t" + "fxch\n\t" + "fyl2xp1" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} + +static __inline__ float __fast_log1pf (float x) +{ + float res; + /* fyl2xp1 accurate only for |x| <= 1.0 - 0.5 * sqrt (2.0) */ + if (fabsf (x) >= 1.0 - 0.5 * 1.41421356237309504880) + res = __fast_logf (1.0 + x); + else + asm __volatile__ + ("fldln2\n\t" + "fxch\n\t" + "fyl2xp1" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} + +#endif |