#include <math.h>
#include <complex.h>
-#define __IS_FP(x) !!((1?1:(x))/2)
+#define __IS_FP(x) (sizeof((x)+1ULL) == sizeof((x)+1.0f))
#define __IS_CX(x) (__IS_FP(x) && sizeof(x) == sizeof((x)+I))
#define __IS_REAL(x) (__IS_FP(x) && 2*sizeof(x) == sizeof((x)+I))
/* return type */
#ifdef __GNUC__
+/*
+the result must be casted to the right type
+(otherwise the result type is determined by the conversion
+rules applied to all the function return types so it is long
+double or long double complex except for integral functions)
+
+this cannot be done in c99, so the typeof gcc extension is
+used and that the type of ?: depends on wether an operand is
+a null pointer constant or not
+(in c11 _Generic can be used)
+
+the c arguments below must be integer constant expressions
+so they can be in null pointer constants
+(__IS_FP above was carefully chosen this way)
+*/
+/* if c then t else void */
+#define __type1(c,t) __typeof__(*(0?(t*)0:(void*)!(c)))
+/* if c then t1 else t2 */
+#define __type2(c,t1,t2) __typeof__(*(0?(__type1(c,t1)*)0:(__type1(!(c),t2)*)0))
/* cast to double when x is integral, otherwise use typeof(x) */
-#define __RETCAST(x) (__typeof__(*( \
- 0 ? (__typeof__(0 ? (double *)0 : (void *)__IS_FP(x)))0 : \
- (__typeof__(0 ? (__typeof__(x) *)0 : (void *)!__IS_FP(x)))0 )))
-/* 2 args case, consider complex types (for cpow) */
-#define __RETCAST_2(x, y) (__typeof__(*( \
- 0 ? (__typeof__(0 ? (double *)0 : \
- (void *)!((!__IS_FP(x) || !__IS_FP(y)) && __FLT((x)+(y)+1.0f))))0 : \
- 0 ? (__typeof__(0 ? (double complex *)0 : \
- (void *)!((!__IS_FP(x) || !__IS_FP(y)) && __FLTCX((x)+(y)))))0 : \
- (__typeof__(0 ? (__typeof__((x)+(y)) *)0 : \
- (void *)((!__IS_FP(x) || !__IS_FP(y)) && (__FLT((x)+(y)+1.0f) || __FLTCX((x)+(y))))))0 )))
-/* 3 args case, don't consider complex types (fma only) */
-#define __RETCAST_3(x, y, z) (__typeof__(*( \
- 0 ? (__typeof__(0 ? (double *)0 : \
- (void *)!((!__IS_FP(x) || !__IS_FP(y) || !__IS_FP(z)) && __FLT((x)+(y)+(z)+1.0f))))0 : \
- (__typeof__(0 ? (__typeof__((x)+(y)) *)0 : \
- (void *)((!__IS_FP(x) || !__IS_FP(y) || !__IS_FP(z)) && __FLT((x)+(y)+(z)+1.0f))))0 )))
+#define __RETCAST(x) ( \
+ __type2(__IS_FP(x), __typeof__(x), double))
+/* 2 args case, should work for complex types (cpow) */
+#define __RETCAST_2(x, y) ( \
+ __type2(__IS_FP(x) && __IS_FP(y), \
+ __typeof__((x)+(y)), \
+ __typeof__((x)+(y)+1.0)))
+/* 3 args case (fma only) */
+#define __RETCAST_3(x, y, z) ( \
+ __type2(__IS_FP(x) && __IS_FP(y) && __IS_FP(z), \
+ __typeof__((x)+(y)+(z)), \
+ __typeof__((x)+(y)+(z)+1.0)))
/* drop complex from the type of x */
-#define __TO_REAL(x) *( \
- 0 ? (__typeof__(0 ? (double *)0 : (void *)!__DBLCX(x)))0 : \
- 0 ? (__typeof__(0 ? (float *)0 : (void *)!__FLTCX(x)))0 : \
- 0 ? (__typeof__(0 ? (long double *)0 : (void *)!__LDBLCX(x)))0 : \
- (__typeof__(0 ? (__typeof__(x) *)0 : (void *)__IS_CX(x)))0 )
+/* TODO: wrong when sizeof(long double)==sizeof(double) */
+#define __RETCAST_REAL(x) ( \
+ __type2(__IS_FP(x) && sizeof((x)+I) == sizeof(float complex), float, \
+ __type2(sizeof((x)+1.0+I) == sizeof(double complex), double, \
+ long double)))
+/* add complex to the type of x */
+#define __RETCAST_CX(x) (__typeof__(__RETCAST(x)0+I))
#else
#define __RETCAST(x)
#define __RETCAST_2(x, y)
#define __RETCAST_3(x, y, z)
+#define __RETCAST_REAL(x)
+#define __RETCAST_CX(x)
#endif
/* function selection */
__LDBL((x)+(y)) ? fun ## l (x, y) : \
fun(x, y) ))
-#define __tg_complex(fun, x) (__RETCAST((x)+I)( \
+#define __tg_complex(fun, x) (__RETCAST_CX(x)( \
__FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
__LDBLCX((x)+I) ? fun ## l (x) : \
fun(x) ))
-#define __tg_complex_retreal(fun, x) (__RETCAST(__TO_REAL(x))( \
+#define __tg_complex_retreal(fun, x) (__RETCAST_REAL(x)( \
__FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
__LDBLCX((x)+I) ? fun ## l (x) : \
fun(x) ))
__LDBL((x)+(y)) ? powl(x, y) : \
pow(x, y) ))
-#define __tg_real_complex_fabs(x) (__RETCAST(__TO_REAL(x))( \
+#define __tg_real_complex_fabs(x) (__RETCAST_REAL(x)( \
__FLTCX(x) ? cabsf(x) : \
__DBLCX(x) ? cabs(x) : \
__LDBLCX(x) ? cabsl(x) : \
if ((ix-0x3ff00000 | lx) == 0) {
/* acos(1)=0, acos(-1)=pi */
if (hx >> 31)
- return 2*pio2_hi + 0x1p-1000;
+ return 2*pio2_hi + 0x1p-120f;
return 0;
}
return 0/(x-x);
/* |x| < 0.5 */
if (ix < 0x3fe00000) {
if (ix <= 0x3c600000) /* |x| < 2**-57 */
- return pio2_hi + 0x1p-1000;
+ return pio2_hi + 0x1p-120f;
return pio2_hi - (x - (pio2_lo-x*R(x*x)));
}
/* x < -0.5 */
((u.bits.manh & ~LDBL_NBIT) | u.bits.manl) == 0) {
if (expsign > 0)
return 0; /* acos(1) = 0 */
- return 2*pio2_hi + 0x1p-1000; /* acos(-1)= pi */
+ return 2*pio2_hi + 0x1p-120f; /* acos(-1)= pi */
}
return 0/(x-x); /* acos(|x|>1) is NaN */
}
/* |x| < 0.5 */
if (expt < 0x3fff - 1) {
if (expt < 0x3fff - 65)
- return pio2_hi + 0x1p-1000; /* x < 0x1p-65: acosl(x)=pi/2 */
+ return pio2_hi + 0x1p-120f; /* x < 0x1p-65: acosl(x)=pi/2 */
return pio2_hi - (x - (pio2_lo - x * __invtrigl_R(x*x)));
}
/* x < -0.5 */
GET_LOW_WORD(lx, x);
if ((ix-0x3ff00000 | lx) == 0)
/* asin(1) = +-pi/2 with inexact */
- return x*pio2_hi + 0x1p-1000;
+ return x*pio2_hi + 0x1p-120f;
return 0/(x-x);
}
/* |x| < 0.5 */
if (ix < 0x3fe00000) {
if (ix < 0x3e500000) {
/* |x|<0x1p-26, return x with inexact if x!=0*/
- FORCE_EVAL(x + 0x1p1000);
+ FORCE_EVAL(x + 0x1p120f);
return x;
}
return x + x*R(x*x);
x = log1p(x + x*x/(sqrt(x*x+1)+1));
} else {
/* |x| < 0x1p-26, raise inexact if x != 0 */
- FORCE_EVAL(x + 0x1p1000);
+ FORCE_EVAL(x + 0x1p120f);
}
return s ? -x : x;
}
x = log1pl(x + x*x/(sqrtl(x*x+1)+1));
} else {
/* |x| < 0x1p-32, raise inexact if x!=0 */
- FORCE_EVAL(x + 0x1p1000);
+ FORCE_EVAL(x + 0x1p120f);
}
return s ? -x : x;
}
if (expt == 0x3fff &&
((u.bits.manh&~LDBL_NBIT)|u.bits.manl) == 0)
/* asin(+-1)=+-pi/2 with inexact */
- return x*pio2_hi + 0x1p-1000;
+ return x*pio2_hi + 0x1p-120f;
return 0/(x-x);
}
if (expt < 0x3fff - 1) { /* |x| < 0.5 */
if (expt < 0x3fff - 32) { /* |x|<0x1p-32, asinl(x)=x */
/* return x with inexact if x!=0 */
- FORCE_EVAL(x + 0x1p1000);
+ FORCE_EVAL(x + 0x1p120f);
return x;
}
return x + x*__invtrigl_R(x*x);
if (ix >= 0x44100000) { /* if |x| >= 2^66 */
if (isnan(x))
return x;
- z = atanhi[3] + 0x1p-1000;
+ z = atanhi[3] + 0x1p-120f;
return sign ? -z : z;
}
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e400000) { /* |x| < 2^-27 */
/* raise inexact if x!=0 */
- FORCE_EVAL(x + 0x1p1000);
+ FORCE_EVAL(x + 0x1p120f);
return x;
}
id = -1;
switch(m) {
case 0:
case 1: return y; /* atan(+-0,+anything)=+-0 */
- case 2: return 2*pio2_hi+0x1p-1000; /* atan(+0,-anything) = pi */
- case 3: return -2*pio2_hi-0x1p-1000; /* atan(-0,-anything) =-pi */
+ case 2: return 2*pio2_hi+0x1p-120f; /* atan(+0,-anything) = pi */
+ case 3: return -2*pio2_hi-0x1p-120f; /* atan(-0,-anything) =-pi */
}
}
/* when x = 0 */
if (exptx==0 && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0)
- return expsigny < 0 ? -pio2_hi-0x1p-1000 : pio2_hi+0x1p-1000;
+ return expsigny < 0 ? -pio2_hi-0x1p-120f : pio2_hi+0x1p-120f;
/* when x is INF */
if (exptx == 0x7fff) {
if (expty == 0x7fff) {
switch(m) {
- case 0: return pio2_hi*0.5+0x1p-1000; /* atan(+INF,+INF) */
- case 1: return -pio2_hi*0.5-0x1p-1000; /* atan(-INF,+INF) */
- case 2: return 1.5*pio2_hi+0x1p-1000; /* atan(+INF,-INF) */
- case 3: return -1.5*pio2_hi-0x1p-1000; /* atan(-INF,-INF) */
+ case 0: return pio2_hi*0.5+0x1p-120f; /* atan(+INF,+INF) */
+ case 1: return -pio2_hi*0.5-0x1p-120f; /* atan(-INF,+INF) */
+ case 2: return 1.5*pio2_hi+0x1p-120f; /* atan(+INF,-INF) */
+ case 3: return -1.5*pio2_hi-0x1p-120f; /* atan(-INF,-INF) */
}
} else {
switch(m) {
case 0: return 0.0; /* atan(+...,+INF) */
case 1: return -0.0; /* atan(-...,+INF) */
- case 2: return 2*pio2_hi+0x1p-1000; /* atan(+...,-INF) */
- case 3: return -2*pio2_hi-0x1p-1000; /* atan(-...,-INF) */
+ case 2: return 2*pio2_hi+0x1p-120f; /* atan(+...,-INF) */
+ case 3: return -2*pio2_hi-0x1p-120f; /* atan(-...,-INF) */
}
}
}
/* when y is INF */
if (expty == 0x7fff)
- return expsigny < 0 ? -pio2_hi-0x1p-1000 : pio2_hi+0x1p-1000;
+ return expsigny < 0 ? -pio2_hi-0x1p-120f : pio2_hi+0x1p-120f;
/* compute y/x */
k = expty-exptx;
if(k > LDBL_MANT_DIG+2) { /* |y/x| huge */
- z = pio2_hi+0x1p-1000;
+ z = pio2_hi+0x1p-120f;
m &= 1;
} else if (expsignx < 0 && k < -LDBL_MANT_DIG-2) /* |y/x| tiny, x<0 */
z = 0.0;
if (expt == 0x7fff &&
((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0) /* NaN */
return x+x;
- z = atanhi[3] + 0x1p-1000;
+ z = atanhi[3] + 0x1p-120f;
return expsign < 0 ? -z : z;
}
/* Extract the exponent and the first few bits of the mantissa. */
if (expman < ((0x3fff - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
if (expt < 0x3fff - 32) { /* if |x| is small, atanl(x)~=x */
/* raise inexact if x!=0 */
- FORCE_EVAL(x + 0x1p1000);
+ FORCE_EVAL(x + 0x1p120f);
return x;
}
id = -1;
-/* origin: FreeBSD /usr/src/lib/msun/src/e_cosh.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunSoft, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* cosh(x)
- * Method :
- * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
- * 1. Replace x by |x| (cosh(x) = cosh(-x)).
- * 2.
- * [ exp(x) - 1 ]^2
- * 0 <= x <= ln2/2 : cosh(x) := 1 + -------------------
- * 2*exp(x)
- *
- * exp(x) + 1/exp(x)
- * ln2/2 <= x <= 22 : cosh(x) := -------------------
- * 2
- * 22 <= x <= lnovft : cosh(x) := exp(x)/2
- * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
- * ln2ovft < x : cosh(x) := inf (overflow)
- *
- * Special cases:
- * cosh(x) is |x| if x is +INF, -INF, or NaN.
- * only cosh(0)=1 is exact for finite x.
- */
-
#include "libm.h"
+/* cosh(x) = (exp(x) + 1/exp(x))/2
+ * = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x)
+ * = 1 + x*x/2 + o(x^4)
+ */
double cosh(double x)
{
union {double f; uint64_t i;} u = {.f = x};
- uint32_t ix;
+ uint32_t w;
double t;
/* |x| */
u.i &= (uint64_t)-1/2;
x = u.f;
- ix = u.i >> 32;
+ w = u.i >> 32;
- /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
- if (ix < 0x3fd62e43) {
- t = expm1(x);
- if (ix < 0x3c800000)
+ /* |x| < log(2) */
+ if (w < 0x3fe62e42) {
+ if (w < 0x3ff00000 - (26<<20)) {
+ /* raise inexact if x!=0 */
+ FORCE_EVAL(x + 0x1p120f);
return 1;
+ }
+ t = expm1(x);
return 1 + t*t/(2*(1+t));
}
- /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */
- if (ix < 0x40360000) {
+ /* |x| < log(DBL_MAX) */
+ if (w < 0x40862e42) {
t = exp(x);
- return 0.5*t + 0.5/t;
+ /* note: if x>log(0x1p26) then the 1/t is not needed */
+ return 0.5*(t + 1/t);
}
- /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
- if (ix < 0x40862e42)
- return 0.5*exp(x);
-
- /* |x| in [log(maxdouble), overflowthresold] */
- if (ix <= 0x408633ce)
- return __expo2(x);
-
- /* overflow (or nan) */
- x *= 0x1p1023;
- return x;
+ /* |x| > log(DBL_MAX) or nan */
+ /* note: the result is stored to handle overflow */
+ t = __expo2(x);
+ return t;
}
-/* origin: FreeBSD /usr/src/lib/msun/src/e_coshf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
#include "libm.h"
float coshf(float x)
{
union {float f; uint32_t i;} u = {.f = x};
- uint32_t ix;
+ uint32_t w;
float t;
/* |x| */
u.i &= 0x7fffffff;
x = u.f;
- ix = u.i;
+ w = u.i;
- /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
- if (ix < 0x3eb17218) {
- t = expm1f(x);
- if (ix < 0x39800000)
+ /* |x| < log(2) */
+ if (w < 0x3f317217) {
+ if (w < 0x3f800000 - (12<<23)) {
+ FORCE_EVAL(x + 0x1p120f);
return 1;
+ }
+ t = expm1f(x);
return 1 + t*t/(2*(1+t));
}
- /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
- if (ix < 0x41100000) {
+ /* |x| < log(FLT_MAX) */
+ if (w < 0x42b17217) {
t = expf(x);
- return 0.5f*t + 0.5f/t;
+ return 0.5f*(t + 1/t);
}
- /* |x| in [9, log(maxfloat)] return 0.5f*exp(|x|) */
- if (ix < 0x42b17217)
- return 0.5f*expf(x);
-
- /* |x| in [log(maxfloat), overflowthresold] */
- if (ix <= 0x42b2d4fc)
- return __expo2f(x);
-
- /* |x| > overflowthresold or nan */
- x *= 0x1p127f;
- return x;
+ /* |x| > log(FLT_MAX) or nan */
+ t = __expo2f(x);
+ return t;
}
-/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_coshl.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* coshl(x)
- * Method :
- * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
- * 1. Replace x by |x| (coshl(x) = coshl(-x)).
- * 2.
- * [ exp(x) - 1 ]^2
- * 0 <= x <= ln2/2 : coshl(x) := 1 + -------------------
- * 2*exp(x)
- *
- * exp(x) + 1/exp(x)
- * ln2/2 <= x <= 22 : coshl(x) := -------------------
- * 2
- * 22 <= x <= lnovft : coshl(x) := expl(x)/2
- * lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2)
- * ln2ovft < x : coshl(x) := inf (overflow)
- *
- * Special cases:
- * coshl(x) is |x| if x is +INF, -INF, or NaN.
- * only coshl(0)=1 is exact for finite x.
- */
-
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
struct{uint64_t m; uint16_t se; uint16_t pad;} i;
} u = {.f = x};
unsigned ex = u.i.se & 0x7fff;
+ uint32_t w;
long double t;
- uint32_t mx,lx;
/* |x| */
u.i.se = ex;
x = u.f;
- mx = u.i.m >> 32;
- lx = u.i.m;
+ w = u.i.m >> 32;
- /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
- if (ex < 0x3fff-2 || (ex == 0x3fff-2 && mx < 0xb17217f7)) {
- t = expm1l(x);
- if (ex < 0x3fff-64)
+ /* |x| < log(2) */
+ if (ex < 0x3fff-1 || (ex == 0x3fff-1 && w < 0xb17217f7)) {
+ if (ex < 0x3fff-32) {
+ FORCE_EVAL(x + 0x1p120f);
return 1;
+ }
+ t = expm1l(x);
return 1 + t*t/(2*(1+t));
}
- /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
- if (ex < 0x3fff+4 || (ex == 0x3fff+4 && mx < 0xb0000000)) {
+ /* |x| < log(LDBL_MAX) */
+ if (ex < 0x3fff+13 || (ex == 0x3fff+13 && w < 0xb17217f7)) {
t = expl(x);
- return 0.5*t + 0.5/t;
- }
-
- /* |x| in [22, ln(maxdouble)] return 0.5*exp(|x|) */
- if (ex < 0x3fff+13 || (ex == 0x3fff+13 && mx < 0xb1700000))
- return 0.5*expl(x);
-
- /* |x| in [log(maxdouble), log(2*maxdouble)) */
- if (ex == 0x3fff+13 && (mx < 0xb174ddc0 ||
- (mx == 0xb174ddc0 && lx < 0x31aec0eb))) {
- t = expl(0.5*x);
- return 0.5*t*t;
+ return 0.5*(t + 1/t);
}
- /* |x| >= log(2*maxdouble) or nan */
- return x*0x1p16383L;
+ /* |x| > log(LDBL_MAX) or nan */
+ t = expl(0.5*x);
+ return 0.5*t*t;
}
#endif
return x;
}
if (x >= 128) {
- STRICT_ASSIGN(float, x, x * 0x1p127);
+ STRICT_ASSIGN(float, x, x * 0x1p127f);
return x;
}
if (x <= -150) {
- STRICT_ASSIGN(float, x, 0x1p-100*0x1p-100);
+ STRICT_ASSIGN(float, x, 0x1p-100f*0x1p-100f);
return x;
}
} else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */
-/* origin: FreeBSD /usr/src/lib/msun/src/e_sinh.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunSoft, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* sinh(x)
- * Method :
- * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
- * 1. Replace x by |x| (sinh(-x) = -sinh(x)).
- * 2.
- * E + E/(E+1)
- * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x)
- * 2
- *
- * 22 <= x <= lnovft : sinh(x) := exp(x)/2
- * lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2)
- * ln2ovft < x : sinh(x) := x*shuge (overflow)
- *
- * Special cases:
- * sinh(x) is |x| if x is +INF, -INF, or NaN.
- * only sinh(0)=0 is exact for finite x.
- */
-
#include "libm.h"
-static const double huge = 1.0e307;
-
+/* sinh(x) = (exp(x) - 1/exp(x))/2
+ * = (exp(x)-1 + (exp(x)-1)/exp(x))/2
+ * = x + x^3/6 + o(x^5)
+ */
double sinh(double x)
{
- double t, h;
- int32_t ix, jx;
-
- /* High word of |x|. */
- GET_HIGH_WORD(jx, x);
- ix = jx & 0x7fffffff;
-
- /* x is INF or NaN */
- if (ix >= 0x7ff00000)
- return x + x;
+ union {double f; uint64_t i;} u = {.f = x};
+ uint32_t w;
+ double t, h, absx;
h = 0.5;
- if (jx < 0) h = -h;
- /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
- if (ix < 0x40360000) { /* |x|<22 */
- if (ix < 0x3e300000) /* |x|<2**-28 */
- /* raise inexact, return x */
- if (huge+x > 1.0)
+ if (u.i >> 63)
+ h = -h;
+ /* |x| */
+ u.i &= (uint64_t)-1/2;
+ absx = u.f;
+ w = u.i >> 32;
+
+ /* |x| < log(DBL_MAX) */
+ if (w < 0x40862e42) {
+ t = expm1(absx);
+ if (w < 0x3ff00000) {
+ if (w < 0x3ff00000 - (26<<20))
+ /* note: inexact is raised by expm1 */
+ /* note: this branch avoids underflow */
return x;
- t = expm1(fabs(x));
- if (ix < 0x3ff00000)
- return h*(2.0*t - t*t/(t+1.0));
- return h*(t + t/(t+1.0));
+ return h*(2*t - t*t/(t+1));
+ }
+ /* note: |x|>log(0x1p26)+eps could be just h*exp(x) */
+ return h*(t + t/(t+1));
}
- /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
- if (ix < 0x40862E42)
- return h*exp(fabs(x));
-
- /* |x| in [log(maxdouble), overflowthresold] */
- if (ix <= 0x408633CE)
- return h * 2.0 * __expo2(fabs(x)); /* h is for sign only */
-
- /* |x| > overflowthresold, sinh(x) overflow */
- return x*huge;
+ /* |x| > log(DBL_MAX) or nan */
+ /* note: the result is stored to handle overflow */
+ t = 2*h*__expo2(absx);
+ return t;
}
-/* origin: FreeBSD /usr/src/lib/msun/src/e_sinhf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
#include "libm.h"
-static const float huge = 1.0e37;
-
float sinhf(float x)
{
- float t, h;
- int32_t ix, jx;
-
- GET_FLOAT_WORD(jx, x);
- ix = jx & 0x7fffffff;
-
- /* x is INF or NaN */
- if (ix >= 0x7f800000)
- return x + x;
+ union {float f; uint32_t i;} u = {.f = x};
+ uint32_t w;
+ float t, h, absx;
h = 0.5;
- if (jx < 0)
+ if (u.i >> 31)
h = -h;
- /* |x| in [0,9], return sign(x)*0.5*(E+E/(E+1))) */
- if (ix < 0x41100000) { /* |x|<9 */
- if (ix < 0x39800000) /* |x|<2**-12 */
- /* raise inexact, return x */
- if (huge+x > 1.0f)
+ /* |x| */
+ u.i &= 0x7fffffff;
+ absx = u.f;
+ w = u.i;
+
+ /* |x| < log(FLT_MAX) */
+ if (w < 0x42b17217) {
+ t = expm1f(absx);
+ if (w < 0x3f800000) {
+ if (w < 0x3f800000 - (12<<23))
return x;
- t = expm1f(fabsf(x));
- if (ix < 0x3f800000)
- return h*(2.0f*t - t*t/(t+1.0f));
- return h*(t + t/(t+1.0f));
+ return h*(2*t - t*t/(t+1));
+ }
+ return h*(t + t/(t+1));
}
- /* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */
- if (ix < 0x42b17217)
- return h*expf(fabsf(x));
-
- /* |x| in [logf(maxfloat), overflowthresold] */
- if (ix <= 0x42b2d4fc)
- return h * 2.0f * __expo2f(fabsf(x)); /* h is for sign only */
-
- /* |x| > overflowthresold, sinh(x) overflow */
- return x*huge;
+ /* |x| > logf(FLT_MAX) or nan */
+ t = 2*h*__expo2f(absx);
+ return t;
}
-/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_sinhl.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* sinhl(x)
- * Method :
- * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
- * 1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
- * 2.
- * E + E/(E+1)
- * 0 <= x <= 25 : sinhl(x) := --------------, E=expm1l(x)
- * 2
- *
- * 25 <= x <= lnovft : sinhl(x) := expl(x)/2
- * lnovft <= x <= ln2ovft: sinhl(x) := expl(x/2)/2 * expl(x/2)
- * ln2ovft < x : sinhl(x) := x*huge (overflow)
- *
- * Special cases:
- * sinhl(x) is |x| if x is +INF, -INF, or NaN.
- * only sinhl(0)=0 is exact for finite x.
- */
-
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
return sinh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double huge = 1.0e4931L;
-
long double sinhl(long double x)
{
- long double t,w,h;
- uint32_t jx,ix,i0,i1;
-
- /* Words of |x|. */
- GET_LDOUBLE_WORDS(jx, i0, i1, x);
- ix = jx & 0x7fff;
-
- /* x is INF or NaN */
- if (ix == 0x7fff) return x + x;
+ union {
+ long double f;
+ struct{uint64_t m; uint16_t se; uint16_t pad;} i;
+ } u = {.f = x};
+ unsigned ex = u.i.se & 0x7fff;
+ long double h, t, absx;
h = 0.5;
- if (jx & 0x8000)
+ if (u.i.se & 0x8000)
h = -h;
- /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
- if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x| < 25 */
- if (ix < 0x3fdf) /* |x|<2**-32 */
- if (huge + x > 1.0)
- return x;/* sinh(tiny) = tiny with inexact */
- t = expm1l(fabsl(x));
- if (ix < 0x3fff)
- return h*(2.0*t - t*t/(t + 1.0));
- return h*(t + t/(t + 1.0));
- }
-
- /* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */
- if (ix < 0x400c || (ix == 0x400c && i0 < 0xb17217f7))
- return h*expl(fabsl(x));
-
- /* |x| in [log(maxdouble), overflowthreshold] */
- if (ix < 0x400c || (ix == 0x400c && (i0 < 0xb174ddc0 ||
- (i0 == 0xb174ddc0 && i1 <= 0x31aec0ea)))) {
- w = expl(0.5*fabsl(x));
- t = h*w;
- return t*w;
+ /* |x| */
+ u.i.se = ex;
+ absx = u.f;
+
+ /* |x| < log(LDBL_MAX) */
+ if (ex < 0x3fff+13 || (ex == 0x3fff+13 && u.i.m>>32 < 0xb17217f7)) {
+ t = expm1l(absx);
+ if (ex < 0x3fff) {
+ if (ex < 0x3fff-32)
+ return x;
+ return h*(2*t - t*t/(1+t));
+ }
+ return h*(t + t/(t+1));
}
- /* |x| > overflowthreshold, sinhl(x) overflow */
- return x*huge;
+ /* |x| > log(LDBL_MAX) or nan */
+ t = expl(0.5*absx);
+ return h*t*t;
}
#endif
-/* origin: FreeBSD /usr/src/lib/msun/src/s_tanh.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* Tanh(x)
- * Return the Hyperbolic Tangent of x
- *
- * Method :
- * x -x
- * e - e
- * 0. tanh(x) is defined to be -----------
- * x -x
- * e + e
- * 1. reduce x to non-negative by tanh(-x) = -tanh(x).
- * 2. 0 <= x < 2**-28 : tanh(x) := x with inexact if x != 0
- * -t
- * 2**-28 <= x < 1 : tanh(x) := -----; t = expm1(-2x)
- * t + 2
- * 2
- * 1 <= x < 22 : tanh(x) := 1 - -----; t = expm1(2x)
- * t + 2
- * 22 <= x <= INF : tanh(x) := 1.
- *
- * Special cases:
- * tanh(NaN) is NaN;
- * only tanh(0)=0 is exact for finite argument.
- */
-
#include "libm.h"
-static const double tiny = 1.0e-300, huge = 1.0e300;
-
+/* tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x))
+ * = (exp(2*x) - 1)/(exp(2*x) - 1 + 2)
+ * = (1 - exp(-2*x))/(exp(-2*x) - 1 + 2)
+ */
double tanh(double x)
{
- double t,z;
- int32_t jx,ix;
-
- GET_HIGH_WORD(jx, x);
- ix = jx & 0x7fffffff;
+ union {double f; uint64_t i;} u = {.f = x};
+ uint32_t w;
+ int sign;
+ double t;
- /* x is INF or NaN */
- if (ix >= 0x7ff00000) {
- if (jx >= 0)
- return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */
- else
- return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */
- }
+ /* x = |x| */
+ sign = u.i >> 63;
+ u.i &= (uint64_t)-1/2;
+ x = u.f;
+ w = u.i >> 32;
- if (ix < 0x40360000) { /* |x| < 22 */
- if (ix < 0x3e300000) { /* |x| < 2**-28 */
- /* tanh(tiny) = tiny with inexact */
- if (huge+x > 1.0f)
- return x;
- }
- if (ix >= 0x3ff00000) { /* |x| >= 1 */
- t = expm1(2.0f*fabs(x));
- z = 1.0f - 2.0f/(t+2.0f);
+ if (w > 0x3fe193ea) {
+ /* |x| > log(3)/2 ~= 0.5493 or nan */
+ if (w > 0x40340000) {
+ /* |x| > 20 or nan */
+ /* note: this branch avoids raising overflow */
+ /* raise inexact if x!=+-inf and handle nan */
+ t = 1 + 0/(x + 0x1p-120f);
} else {
- t = expm1(-2.0f*fabs(x));
- z= -t/(t+2.0f);
+ t = expm1(2*x);
+ t = 1 - 2/(t+2);
}
- } else { /* |x| >= 22, return +-1 */
- z = 1.0f - tiny; /* raise inexact */
+ } else if (w > 0x3fd058ae) {
+ /* |x| > log(5/3)/2 ~= 0.2554 */
+ t = expm1(2*x);
+ t = t/(t+2);
+ } else {
+ /* |x| is small, up to 2ulp error in [0.1,0.2554] */
+ t = expm1(-2*x);
+ t = -t/(t+2);
}
- return jx >= 0 ? z : -z;
+ return sign ? -t : t;
}
-/* origin: FreeBSD /usr/src/lib/msun/src/s_tanhf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
#include "libm.h"
-static const float
-tiny = 1.0e-30,
-huge = 1.0e30;
-
float tanhf(float x)
{
- float t,z;
- int32_t jx,ix;
+ union {float f; uint32_t i;} u = {.f = x};
+ uint32_t w;
+ int sign;
+ float t;
- GET_FLOAT_WORD(jx, x);
- ix = jx & 0x7fffffff;
+ /* x = |x| */
+ sign = u.i >> 31;
+ u.i &= 0x7fffffff;
+ x = u.f;
+ w = u.i;
- /* x is INF or NaN */
- if(ix >= 0x7f800000) {
- if (jx >= 0)
- return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */
- else
- return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */
- }
-
- if (ix < 0x41100000) { /* |x| < 9 */
- if (ix < 0x39800000) { /* |x| < 2**-12 */
- /* tanh(tiny) = tiny with inexact */
- if (huge+x > 1.0f)
- return x;
- }
- if (ix >= 0x3f800000) { /* |x|>=1 */
- t = expm1f(2.0f*fabsf(x));
- z = 1.0f - 2.0f/(t+2.0f);
+ if (w > 0x3f0c9f54) {
+ /* |x| > log(3)/2 ~= 0.5493 or nan */
+ if (w > 0x41200000) {
+ /* |x| > 10 */
+ t = 1 + 0/(x + 0x1p-120f);
} else {
- t = expm1f(-2.0f*fabsf(x));
- z = -t/(t+2.0f);
+ t = expm1f(2*x);
+ t = 1 - 2/(t+2);
}
- } else { /* |x| >= 9, return +-1 */
- z = 1.0f - tiny; /* raise inexact */
+ } else if (w > 0x3e82c578) {
+ /* |x| > log(5/3)/2 ~= 0.2554 */
+ t = expm1f(2*x);
+ t = t/(t+2);
+ } else {
+ /* |x| is small */
+ t = expm1f(-2*x);
+ t = -t/(t+2);
}
- return jx >= 0 ? z : -z;
+ return sign ? -t : t;
}
-/* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_tanhl.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* tanhl(x)
- * Return the Hyperbolic Tangent of x
- *
- * Method :
- * x -x
- * e - e
- * 0. tanhl(x) is defined to be -----------
- * x -x
- * e + e
- * 1. reduce x to non-negative by tanhl(-x) = -tanhl(x).
- * 2. 0 <= x <= 2**-55 : tanhl(x) := x*(one+x)
- * -t
- * 2**-55 < x <= 1 : tanhl(x) := -----; t = expm1l(-2x)
- * t + 2
- * 2
- * 1 <= x <= 23.0 : tanhl(x) := 1- ----- ; t=expm1l(2x)
- * t + 2
- * 23.0 < x <= INF : tanhl(x) := 1.
- *
- * Special cases:
- * tanhl(NaN) is NaN;
- * only tanhl(0)=0 is exact for finite argument.
- */
-
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
return tanh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double tiny = 1.0e-4900L;
-
long double tanhl(long double x)
{
- long double t,z;
- int32_t se;
- uint32_t jj0,jj1,ix;
+ union {
+ long double f;
+ struct{uint64_t m; uint16_t se; uint16_t pad;} i;
+ } u = {.f = x};
+ unsigned ex = u.i.se & 0x7fff;
+ unsigned sign = u.i.se & 0x8000;
+ uint32_t w;
+ long double t;
- /* High word of |x|. */
- GET_LDOUBLE_WORDS(se, jj0, jj1, x);
- ix = se & 0x7fff;
-
- /* x is INF or NaN */
- if (ix == 0x7fff) {
- /* for NaN it's not important which branch: tanhl(NaN) = NaN */
- if (se & 0x8000)
- return 1.0/x-1.0; /* tanhl(-inf)= -1; */
- return 1.0/x+1.0; /* tanhl(+inf)= +1 */
- }
+ /* x = |x| */
+ u.i.se = ex;
+ x = u.f;
+ w = u.i.m >> 32;
- /* |x| < 23 */
- if (ix < 0x4003 || (ix == 0x4003 && jj0 < 0xb8000000u)) {
- if ((ix|jj0|jj1) == 0) /* x == +- 0 */
- return x;
- if (ix < 0x3fc8) /* |x| < 2**-55 */
- return x*(1.0+tiny); /* tanh(small) = small */
- if (ix >= 0x3fff) { /* |x| >= 1 */
- t = expm1l(2.0*fabsl(x));
- z = 1.0 - 2.0/(t+2.0);
+ if (ex > 0x3ffe || (ex == 0x3ffe && w > 0x8c9f53d5)) {
+ /* |x| > log(3)/2 ~= 0.5493 or nan */
+ if (ex >= 0x3fff+5) {
+ /* |x| >= 32 */
+ t = 1 + 0/(x + 0x1p-120f);
} else {
- t = expm1l(-2.0*fabsl(x));
- z = -t/(t+2.0);
+ t = expm1l(2*x);
+ t = 1 - 2/(t+2);
}
- /* |x| > 23, return +-1 */
+ } else if (ex > 0x3ffd || (ex == 0x3ffd && w > 0x82c577d4)) {
+ /* |x| > log(5/3)/2 ~= 0.2554 */
+ t = expm1l(2*x);
+ t = t/(t+2);
} else {
- z = 1.0 - tiny; /* raise inexact flag */
+ /* |x| is small */
+ t = expm1l(-2*x);
+ t = -t/(t+2);
}
- return se & 0x8000 ? -z : z;
+ return sign ? -t : t;
}
#endif
long double tgammal(long double x)
{
long double p, q, z;
- int i, sign;
- if (isnan(x))
- return NAN;
- if (x == INFINITY)
- return INFINITY;
- if (x == -INFINITY)
- return x - x;
+ if (!isfinite(x))
+ return x + INFINITY;
+
q = fabsl(x);
if (q > 13.0) {
- sign = 1;
- if (q > MAXGAML)
- goto goverf;
if (x < 0.0) {
p = floorl(q);
- if (p == q)
- return (x - x) / (x - x);
- i = p;
- if ((i & 1) == 0)
- sign = -1;
z = q - p;
- if (z > 0.5L) {
- p += 1.0;
- z = q - p;
- }
- z = q * sinl(PIL * z);
- z = fabsl(z) * stirf(q);
- if (z <= PIL/LDBL_MAX) {
-goverf:
- return sign * INFINITY;
+ if (z == 0)
+ return 0 / z;
+ if (q > MAXGAML) {
+ z = 0;
+ } else {
+ if (z > 0.5) {
+ p += 1.0;
+ z = q - p;
+ }
+ z = q * sinl(PIL * z);
+ z = fabsl(z) * stirf(q);
+ z = PIL/z;
}
- z = PIL/z;
+ if (0.5 * p == floorl(q * 0.5))
+ z = -z;
+ } else if (x > MAXGAML) {
+ z = x * 0x1p16383L;
} else {
z = stirf(x);
}
- return sign * z;
+ return z;
}
z = 1.0;
return z;
small:
- if (x == 0.0)
- return (x - x) / (x - x);
+ /* z==1 if x was originally +-0 */
+ if (x == 0 && z != 1)
+ return x / x;
if (x < 0.0) {
x = -x;
q = z / (x * __polevll(x, SN, 8));
projects / musl / commitdiff