21 files changed, 298 insertions, 567 deletions

diff --git a/include/tgmath.h b/include/tgmath.h
index 832b052b..e41ccac9 100644
--- a/include/tgmath.h
+++ b/include/tgmath.h
@@ -13,7 +13,7 @@ sizeof(double) == sizeof(long double)
 #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))
 
@@ -27,34 +27,52 @@ sizeof(double) == sizeof(long double)
 /* 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 */
@@ -76,12 +94,12 @@ sizeof(double) == sizeof(long double)
 	__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) ))
@@ -115,7 +133,7 @@ sizeof(double) == sizeof(long double)
 	__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) : \
diff --git a/src/math/acos.c b/src/math/acos.c
index be95d25e..cd5d06a6 100644
--- a/src/math/acos.c
+++ b/src/math/acos.c
@@ -72,7 +72,7 @@ double acos(double 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);
@@ -80,7 +80,7 @@ double acos(double 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 */
diff --git a/src/math/acosl.c b/src/math/acosl.c
index 9e9b01e4..9e7b7fb3 100644
--- a/src/math/acosl.c
+++ b/src/math/acosl.c
@@ -38,14 +38,14 @@ long double acosl(long double x)
 			((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 */
diff --git a/src/math/asin.c b/src/math/asin.c
index a1906b08..d61c04b4 100644
--- a/src/math/asin.c
+++ b/src/math/asin.c
@@ -77,14 +77,14 @@ double asin(double x)
 		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);
diff --git a/src/math/asinh.c b/src/math/asinh.c
index 4152dc39..0829f228 100644
--- a/src/math/asinh.c
+++ b/src/math/asinh.c
@@ -22,7 +22,7 @@ double asinh(double 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;
 }
diff --git a/src/math/asinhl.c b/src/math/asinhl.c
index db966246..3ea88745 100644
--- a/src/math/asinhl.c
+++ b/src/math/asinhl.c
@@ -31,7 +31,7 @@ long double asinhl(long double 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;
 }
diff --git a/src/math/asinl.c b/src/math/asinl.c
index 0ef9853b..8799341d 100644
--- a/src/math/asinl.c
+++ b/src/math/asinl.c
@@ -39,13 +39,13 @@ long double asinl(long double 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);
diff --git a/src/math/atan.c b/src/math/atan.c
index 7fd8a3b8..3c9a59ff 100644
--- a/src/math/atan.c
+++ b/src/math/atan.c
@@ -72,13 +72,13 @@ double atan(double 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;
diff --git a/src/math/atan2l.c b/src/math/atan2l.c
index e86dbffb..7cb42c2f 100644
--- a/src/math/atan2l.c
+++ b/src/math/atan2l.c
@@ -52,39 +52,39 @@ long double atan2l(long double y, long double x)
 		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;
diff --git a/src/math/atanl.c b/src/math/atanl.c
index 6a480a7a..e76693e4 100644
--- a/src/math/atanl.c
+++ b/src/math/atanl.c
@@ -80,7 +80,7 @@ long double atanl(long double x)
 		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. */
@@ -89,7 +89,7 @@ long double atanl(long double x)
 	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;
diff --git a/src/math/cosh.c b/src/math/cosh.c
index 2fcdea8f..100f8231 100644
--- a/src/math/cosh.c
+++ b/src/math/cosh.c
@@ -1,71 +1,40 @@
-/* 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;
 }
diff --git a/src/math/coshf.c b/src/math/coshf.c
index 2bed1258..b09f2ee5 100644
--- a/src/math/coshf.c
+++ b/src/math/coshf.c
@@ -1,54 +1,33 @@
-/* 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;
 }
diff --git a/src/math/coshl.c b/src/math/coshl.c
index 3ea56e00..d09070bb 100644
--- a/src/math/coshl.c
+++ b/src/math/coshl.c
@@ -1,35 +1,3 @@
-/* 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
@@ -45,41 +13,32 @@ long double coshl(long double x)
 		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
diff --git a/src/math/exp2f.c b/src/math/exp2f.c
index 279f32de..ea50db4a 100644
--- a/src/math/exp2f.c
+++ b/src/math/exp2f.c
@@ -97,11 +97,11 @@ float exp2f(float x)
 			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 */
diff --git a/src/math/sinh.c b/src/math/sinh.c
index 0c67ba39..47e36bfa 100644
--- a/src/math/sinh.c
+++ b/src/math/sinh.c
@@ -1,71 +1,39 @@
-/* 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;
 }
diff --git a/src/math/sinhf.c b/src/math/sinhf.c
index b8d88224..6ad19ea2 100644
--- a/src/math/sinhf.c
+++ b/src/math/sinhf.c
@@ -1,57 +1,31 @@
-/* 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;
 }
diff --git a/src/math/sinhl.c b/src/math/sinhl.c
index 8a54677e..14e3371b 100644
--- a/src/math/sinhl.c
+++ b/src/math/sinhl.c
@@ -1,32 +1,3 @@
-/* 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
@@ -35,47 +6,35 @@ long double sinhl(long double x)
 	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
diff --git a/src/math/tanh.c b/src/math/tanh.c
index 21138643..0e766c5c 100644
--- a/src/math/tanh.c
+++ b/src/math/tanh.c
@@ -1,73 +1,41 @@
-/* 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;
 }
diff --git a/src/math/tanhf.c b/src/math/tanhf.c
index 7cb459d0..8099ec30 100644
--- a/src/math/tanhf.c
+++ b/src/math/tanhf.c
@@ -1,55 +1,35 @@
-/* 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;
 }
diff --git a/src/math/tanhl.c b/src/math/tanhl.c
index 92efb20d..66559e9f 100644
--- a/src/math/tanhl.c
+++ b/src/math/tanhl.c
@@ -1,38 +1,3 @@
-/* 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
@@ -41,43 +6,40 @@ long double tanhl(long double x)
 	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
diff --git a/src/math/tgammal.c b/src/math/tgammal.c
index 86782a96..5c1a02a6 100644
--- a/src/math/tgammal.c
+++ b/src/math/tgammal.c
@@ -205,42 +205,36 @@ static long double stirf(long double x)
 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;
@@ -268,8 +262,9 @@ goverf:
 	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));