* 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.
This commit is contained in:
parent
74b2bdc351
commit
72db1c11e9
|
@ -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.
|
||||
|
|
|
@ -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);}
|
||||
|
|
|
@ -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
|
||||
|
|
|
@ -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)));
|
||||
}
|
|
@ -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)));
|
||||
}
|
|
@ -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)));
|
||||
}
|
|
@ -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);
|
||||
}
|
||||
|
|
@ -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);
|
||||
}
|
|
@ -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);
|
||||
}
|
|
@ -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;
|
||||
}
|
|
@ -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;
|
||||
}
|
|
@ -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;
|
||||
}
|
|
@ -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
|
Loading…
Reference in New Issue