Merge remote-tracking branch 'nsz/math'

38 files changed, 182 insertions, 210 deletions

diff --git a/src/math/__invtrigl.c b/src/math/__invtrigl.c
index a821842c..a51330e2 100644
--- a/src/math/__invtrigl.c
+++ b/src/math/__invtrigl.c
@@ -28,9 +28,8 @@
 #include "__invtrigl.h"
 
 #if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-/*
- * asinl() and acosl()
- */
+
+/* coefficients used in asinl() and acosl() */
 const long double
 pS0 =  1.66666666666666666631e-01L,
 pS1 = -4.16313987993683104320e-01L,
@@ -45,38 +44,9 @@ qS3 = -1.68285799854822427013e+00L,
 qS4 =  3.90699412641738801874e-01L,
 qS5 = -3.14365703596053263322e-02L;
 
-/*
- * atanl()
- */
-const long double atanhi[] = {
-	 4.63647609000806116202e-01L,
-	 7.85398163397448309628e-01L,
-	 9.82793723247329067960e-01L,
-	 1.57079632679489661926e+00L,
-};
-
-const long double atanlo[] = {
-	 1.18469937025062860669e-20L,
-	-1.25413940316708300586e-20L,
-	 2.55232234165405176172e-20L,
-	-2.50827880633416601173e-20L,
-};
-
-const long double aT[] = {
-	 3.33333333333333333017e-01L,
-	-1.99999999999999632011e-01L,
-	 1.42857142857046531280e-01L,
-	-1.11111111100562372733e-01L,
-	 9.09090902935647302252e-02L,
-	-7.69230552476207730353e-02L,
-	 6.66661718042406260546e-02L,
-	-5.88158892835030888692e-02L,
-	 5.25499891539726639379e-02L,
-	-4.70119845393155721494e-02L,
-	 4.03539201366454414072e-02L,
-	-2.91303858419364158725e-02L,
-	 1.24822046299269234080e-02L,
-};
-
+const long double pi_hi = 3.1415926535897932384626433832795L;
 const long double pi_lo = -5.01655761266833202345e-20L;
+const long double pio2_hi = 1.57079632679489661926L;
+const long double pio2_lo = -2.50827880633416601173e-20L;
+
 #endif
diff --git a/src/math/__invtrigl.h b/src/math/__invtrigl.h
index c3ad3c49..22748b68 100644
--- a/src/math/__invtrigl.h
+++ b/src/math/__invtrigl.h
@@ -32,15 +32,6 @@
 #define BIAS            (LDBL_MAX_EXP - 1)
 #define MANH_SIZE       LDBL_MANH_SIZE
 
-/* Approximation thresholds. */
-#define ASIN_LINEAR     (BIAS - 32)     /* 2**-32 */
-#define ACOS_CONST      (BIAS - 65)     /* 2**-65 */
-#define ATAN_CONST      (BIAS + 65)     /* 2**65 */
-#define ATAN_LINEAR     (BIAS - 32)     /* 2**-32 */
-
-/* 0.95 */
-#define THRESH  ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)
-
 /* Constants shared by the long double inverse trig functions. */
 #define pS0     __pS0
 #define pS1     __pS1
@@ -54,56 +45,24 @@
 #define qS3     __qS3
 #define qS4     __qS4
 #define qS5     __qS5
-#define atanhi  __atanhi
-#define atanlo  __atanlo
-#define aT      __aT
+#define pi_hi   __pi_hi
 #define pi_lo   __pi_lo
+#define pio2_hi __pio2_hi
+#define pio2_lo __pio2_lo
 
-#define pio2_hi atanhi[3]
-#define pio2_lo atanlo[3]
-#define pio4_hi atanhi[1]
-
-#ifdef STRUCT_DECLS
-typedef struct longdouble {
-	uint64_t mant;
-	uint16_t expsign;
-} LONGDOUBLE;
-#else
-typedef long double LONGDOUBLE;
-#endif
-
-extern const LONGDOUBLE pS0, pS1, pS2, pS3, pS4, pS5, pS6;
-extern const LONGDOUBLE qS1, qS2, qS3, qS4, qS5;
-extern const LONGDOUBLE atanhi[], atanlo[], aT[];
-extern const LONGDOUBLE pi_lo;
-
-#ifndef STRUCT_DECLS
-static inline long double
-P(long double x)
-{
-	return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \
-		(pS4 + x * (pS5 + x * pS6)))))));
-}
-
-static inline long double
-Q(long double x)
-{
-	return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5)))));
-}
+extern const long double pS0, pS1, pS2, pS3, pS4, pS5, pS6;
+extern const long double qS1, qS2, qS3, qS4, qS5;
+extern const long double pi_hi, pi_lo, pio2_hi, pio2_lo;
 
-static inline long double
-T_even(long double x)
+static long double P(long double x)
 {
-	return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \
-		(aT[8] + x * (aT[10] + x * aT[12]))))));
+	return x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 +
+		x * (pS4 + x * (pS5 + x * pS6))))));
 }
 
-static inline long double
-T_odd(long double x)
+static long double Q(long double x)
 {
-	return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \
-		(aT[9] + x * aT[11])))));
+	return 1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5))));
 }
-#endif
 
 #endif
diff --git a/src/math/acos.c b/src/math/acos.c
index 54d266ee..0eb15bed 100644
--- a/src/math/acos.c
+++ b/src/math/acos.c
@@ -38,6 +38,7 @@
 static const double
 pi      = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
 pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
+// FIXME
 static const volatile double
 pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
 static const double
diff --git a/src/math/acosl.c b/src/math/acosl.c
index cc565336..83857d49 100644
--- a/src/math/acosl.c
+++ b/src/math/acosl.c
@@ -23,9 +23,7 @@ long double acosl(long double x)
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 #include "__invtrigl.h"
-
-static const long double
-pi = 3.14159265358979323846264338327950280e+00L;
+#define ACOS_CONST      (BIAS - 65)     /* 2**-65 */
 
 long double acosl(long double x)
 {
@@ -41,7 +39,8 @@ long double acosl(long double x)
 			if (expsign > 0)
 				return 0.0;  /* acos(1) = 0 */
 			else
-				return pi + 2.0 * pio2_lo;  /* acos(-1)= pi */
+				// FIXME
+				return pi_hi + 2.0 * pio2_lo;  /* acos(-1)= pi */
 		}
 		return (x - x) / (x - x);  /* acos(|x|>1) is NaN */
 	}
@@ -60,7 +59,7 @@ long double acosl(long double x)
 		s = sqrtl(z);
 		r = p / q;
 		w = r * s - pio2_lo;
-		return pi - 2.0 * (s + w);
+		return pi_hi - 2.0 * (s + w);
 	} else {                   /* x > 0.5 */
 		z = (1.0 - x) * 0.5;
 		s = sqrtl(z);
diff --git a/src/math/asinl.c b/src/math/asinl.c
index ddd807e2..7572767c 100644
--- a/src/math/asinl.c
+++ b/src/math/asinl.c
@@ -24,6 +24,10 @@ long double asinl(long double x)
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 #include "__invtrigl.h"
 static const long double huge = 1.000e+300;
+static const long double pio4_hi = 7.85398163397448309628e-01L;
+#define ASIN_LINEAR     (BIAS - 32)     /* 2**-32 */
+/* 0.95 */
+#define THRESH  ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)
 
 long double asinl(long double x)
 {
diff --git a/src/math/atan2.c b/src/math/atan2.c
index 143c3834..825a38fe 100644
--- a/src/math/atan2.c
+++ b/src/math/atan2.c
@@ -39,6 +39,7 @@
 
 #include "libm.h"
 
+// FIXME
 static const volatile double
 tiny  = 1.0e-300;
 static const double
diff --git a/src/math/atan2l.c b/src/math/atan2l.c
index 0fc901c8..45cbfcc1 100644
--- a/src/math/atan2l.c
+++ b/src/math/atan2l.c
@@ -24,10 +24,8 @@ long double atan2l(long double y, long double x)
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 #include "__invtrigl.h"
-static const volatile long double
-tiny = 1.0e-300;
-static const long double
-pi = 3.14159265358979323846264338327950280e+00L;
+// FIXME:
+static const volatile long double tiny = 1.0e-300;
 
 long double atan2l(long double y, long double x)
 {
@@ -55,9 +53,9 @@ long double atan2l(long double y, long double x)
 	if (expty==0 && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)==0) {
 		switch(m) {
 		case 0:
-		case 1: return y;        /* atan(+-0,+anything)=+-0 */
-		case 2: return  pi+tiny; /* atan(+0,-anything) = pi */
-		case 3: return -pi-tiny; /* atan(-0,-anything) =-pi */
+		case 1: return y;           /* atan(+-0,+anything)=+-0 */
+		case 2: return  pi_hi+tiny; /* atan(+0,-anything) = pi */
+		case 3: return -pi_hi-tiny; /* atan(-0,-anything) =-pi */
 		}
 	}
 	/* when x = 0 */
@@ -69,15 +67,15 @@ long double atan2l(long double y, long double x)
 			switch(m) {
 			case 0: return  pio2_hi*0.5+tiny; /* atan(+INF,+INF) */
 			case 1: return -pio2_hi*0.5-tiny; /* atan(-INF,+INF) */
-			case 2: return  1.5*pio2_hi+tiny; /*atan(+INF,-INF)*/
-			case 3: return -1.5*pio2_hi-tiny; /*atan(-INF,-INF)*/
+			case 2: return  1.5*pio2_hi+tiny; /* atan(+INF,-INF) */
+			case 3: return -1.5*pio2_hi-tiny; /* atan(-INF,-INF) */
 			}
 		} else {
 			switch(m) {
-			case 0: return  0.0;     /* atan(+...,+INF) */
-			case 1: return -0.0;     /* atan(-...,+INF) */
-			case 2: return  pi+tiny; /* atan(+...,-INF) */
-			case 3: return -pi-tiny; /* atan(-...,-INF) */
+			case 0: return  0.0;        /* atan(+...,+INF) */
+			case 1: return -0.0;        /* atan(-...,+INF) */
+			case 2: return  pi_hi+tiny; /* atan(+...,-INF) */
+			case 3: return -pi_hi-tiny; /* atan(-...,-INF) */
 			}
 		}
 	}
@@ -95,11 +93,11 @@ long double atan2l(long double y, long double x)
 	else                     /* safe to do y/x */
 		z = atanl(fabsl(y/x));
 	switch (m) {
-	case 0: return z;              /* atan(+,+) */
-	case 1: return -z;             /* atan(-,+) */
-	case 2: return pi - (z-pi_lo); /* atan(+,-) */
+	case 0: return z;               /* atan(+,+) */
+	case 1: return -z;              /* atan(-,+) */
+	case 2: return pi_hi-(z-pi_lo); /* atan(+,-) */
 	default: /* case 3 */
-		return (z-pi_lo) - pi; /* atan(-,-) */
+		return (z-pi_lo)-pi_hi; /* atan(-,-) */
 	}
 }
 #endif
diff --git a/src/math/atanl.c b/src/math/atanl.c
index 36072c17..33ecff0d 100644
--- a/src/math/atanl.c
+++ b/src/math/atanl.c
@@ -23,8 +23,53 @@ long double atanl(long double x)
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 #include "__invtrigl.h"
+
+#define ATAN_CONST      (BIAS + 65)     /* 2**65 */
+#define ATAN_LINEAR     (BIAS - 32)     /* 2**-32 */
 static const long double huge = 1.0e300;
 
+static const long double atanhi[] = {
+	 4.63647609000806116202e-01L,
+	 7.85398163397448309628e-01L,
+	 9.82793723247329067960e-01L,
+	 1.57079632679489661926e+00L,
+};
+
+static const long double atanlo[] = {
+	 1.18469937025062860669e-20L,
+	-1.25413940316708300586e-20L,
+	 2.55232234165405176172e-20L,
+	-2.50827880633416601173e-20L,
+};
+
+static const long double aT[] = {
+	 3.33333333333333333017e-01L,
+	-1.99999999999999632011e-01L,
+	 1.42857142857046531280e-01L,
+	-1.11111111100562372733e-01L,
+	 9.09090902935647302252e-02L,
+	-7.69230552476207730353e-02L,
+	 6.66661718042406260546e-02L,
+	-5.88158892835030888692e-02L,
+	 5.25499891539726639379e-02L,
+	-4.70119845393155721494e-02L,
+	 4.03539201366454414072e-02L,
+	-2.91303858419364158725e-02L,
+	 1.24822046299269234080e-02L,
+};
+
+static long double T_even(long double x)
+{
+	return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] +
+		x * (aT[8] + x * (aT[10] + x * aT[12])))));
+}
+
+static long double T_odd(long double x)
+{
+	return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] +
+		x * (aT[9] + x * aT[11]))));
+}
+
 long double atanl(long double x)
 {
 	union IEEEl2bits u;
diff --git a/src/math/exp10l.c b/src/math/exp10l.c
index 4b2b7ef6..22a4636a 100644
--- a/src/math/exp10l.c
+++ b/src/math/exp10l.c
@@ -5,8 +5,8 @@
 long double exp10l(long double x)
 {
 	static const long double p10[] = {
-		1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10,
-		1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
+		1e-15L, 1e-14L, 1e-13L, 1e-12L, 1e-11L, 1e-10L,
+		1e-9L, 1e-8L, 1e-7L, 1e-6L, 1e-5L, 1e-4L, 1e-3L, 1e-2L, 1e-1L,
 		1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15
 	};
diff --git a/src/math/expm1l.c b/src/math/expm1l.c
index ca0d141f..21a86c00 100644
--- a/src/math/expm1l.c
+++ b/src/math/expm1l.c
@@ -44,13 +44,7 @@
  *
  *                      Relative error:
  * arithmetic   domain     # trials      peak         rms
- *    IEEE    -45,+MAXLOG   200,000     1.2e-19     2.5e-20
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * expm1l overflow   x > MAXLOG         MAXNUM
- *
+ *    IEEE    -45,+maxarg   200,000     1.2e-19     2.5e-20
  */
 
 #include "libm.h"
@@ -61,7 +55,6 @@ long double expm1l(long double x)
 	return expm1(x);
 }
 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double MAXLOGL = 1.1356523406294143949492E4L;
 
 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
    -.5 ln 2  <  x  <  .5 ln 2
@@ -83,19 +76,20 @@ C1 = 6.93145751953125E-1L,
 C2 = 1.428606820309417232121458176568075500134E-6L,
 /* ln 2^-65 */
 minarg = -4.5054566736396445112120088E1L,
-huge = 0x1p10000L;
+/* ln 2^16384 */
+maxarg = 1.1356523406294143949492E4L;
 
 long double expm1l(long double x)
 {
 	long double px, qx, xx;
 	int k;
 
-	/* Overflow.  */
-	if (x > MAXLOGL)
-		return huge*huge;  /* overflow */
+	if (isnan(x))
+		return x;
+	if (x > maxarg)
+		return x*0x1p16383L; /* overflow, unless x==inf */
 	if (x == 0.0)
 		return x;
-	/* Minimum value.*/
 	if (x < minarg)
 		return -1.0;
 
diff --git a/src/math/fma.c b/src/math/fma.c
index 7076d4b9..17f1cdcc 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -89,6 +89,7 @@ static int getexp(long double x)
 
 double fma(double x, double y, double z)
 {
+	#pragma STDC FENV_ACCESS ON
 	long double hi, lo1, lo2, xy;
 	int round, ez, exy;
 
@@ -306,6 +307,7 @@ static inline struct dd dd_mul(double a, double b)
  */
 double fma(double x, double y, double z)
 {
+	#pragma STDC FENV_ACCESS ON
 	double xs, ys, zs, adj;
 	struct dd xy, r;
 	int oround;
diff --git a/src/math/fmaf.c b/src/math/fmaf.c
index cc427fb5..a1c7f4f8 100644
--- a/src/math/fmaf.c
+++ b/src/math/fmaf.c
@@ -37,6 +37,7 @@
  */
 float fmaf(float x, float y, float z)
 {
+	#pragma STDC FENV_ACCESS ON
 	double xy, result;
 	uint32_t hr, lr;
 
diff --git a/src/math/fmal.c b/src/math/fmal.c
index be64f145..ccbe434d 100644
--- a/src/math/fmal.c
+++ b/src/math/fmal.c
@@ -162,6 +162,7 @@ static inline struct dd dd_mul(long double a, long double b)
  */
 long double fmal(long double x, long double y, long double z)
 {
+	#pragma STDC FENV_ACCESS ON
 	long double xs, ys, zs, adj;
 	struct dd xy, r;
 	int oround;
diff --git a/src/math/hypot.c b/src/math/hypot.c
index ba4c7575..9a4cbdb3 100644
--- a/src/math/hypot.c
+++ b/src/math/hypot.c
@@ -117,12 +117,7 @@ double hypot(double x, double y)
 		t2 = a - t1;
 		w  = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 	}
-	if (k != 0) {
-		uint32_t high;
-		t1 = 1.0;
-		GET_HIGH_WORD(high, t1);
-		SET_HIGH_WORD(t1, high+(k<<20));
-		return t1*w;
-	}
+	if (k)
+		w = scalbn(w, k);
 	return w;
 }
diff --git a/src/math/hypotf.c b/src/math/hypotf.c
index 9fd77e6a..4d80178d 100644
--- a/src/math/hypotf.c
+++ b/src/math/hypotf.c
@@ -80,9 +80,7 @@ float hypotf(float x, float y)
 		t2 = a - t1;
 		w  = sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 	}
-	if (k != 0) {
-		SET_FLOAT_WORD(t1, 0x3f800000+(k<<23));
-		return t1*w;
-	}
+	if (k)
+		w = scalbnf(w, k);
 	return w;
 }
diff --git a/src/math/ilogb.c b/src/math/ilogb.c
index 0a3a6a46..5a1819d8 100644
--- a/src/math/ilogb.c
+++ b/src/math/ilogb.c
@@ -8,13 +8,17 @@ int ilogb(double x)
 
 	if (!e) {
 		u.bits <<= 12;
-		if (u.bits == 0)
+		if (u.bits == 0) {
+			FORCE_EVAL(0/0.0f);
 			return FP_ILOGB0;
+		}
 		/* subnormal x */
 		for (e = -0x3ff; u.bits < (uint64_t)1<<63; e--, u.bits<<=1);
 		return e;
 	}
-	if (e == 0x7ff)
+	if (e == 0x7ff) {
+		FORCE_EVAL(0/0.0f);
 		return u.bits<<12 ? FP_ILOGBNAN : INT_MAX;
+	}
 	return e - 0x3ff;
 }
diff --git a/src/math/ilogbf.c b/src/math/ilogbf.c
index 272cbdac..42cd62e2 100644
--- a/src/math/ilogbf.c
+++ b/src/math/ilogbf.c
@@ -8,13 +8,17 @@ int ilogbf(float x)
 
 	if (!e) {
 		u.bits <<= 9;
-		if (u.bits == 0)
+		if (u.bits == 0) {
+			FORCE_EVAL(0/0.0f);
 			return FP_ILOGB0;
+		}
 		/* subnormal x */
 		for (e = -0x7f; u.bits < (uint32_t)1<<31; e--, u.bits<<=1);
 		return e;
 	}
-	if (e == 0xff)
+	if (e == 0xff) {
+		FORCE_EVAL(0/0.0f);
 		return u.bits<<9 ? FP_ILOGBNAN : INT_MAX;
+	}
 	return e - 0x7f;
 }
diff --git a/src/math/ilogbl.c b/src/math/ilogbl.c
index ed9ddcbc..1512934f 100644
--- a/src/math/ilogbl.c
+++ b/src/math/ilogbl.c
@@ -14,15 +14,19 @@ int ilogbl(long double x)
 	int e = u.bits.exp;
 
 	if (!e) {
-		if (m == 0)
+		if (m == 0) {
+			FORCE_EVAL(0/0.0f);
 			return FP_ILOGB0;
+		}
 		/* subnormal x */
 		for (e = -0x3fff+1; m < (uint64_t)1<<63; e--, m<<=1);
 		return e;
 	}
-	if (e == 0x7fff)
+	if (e == 0x7fff) {
+		FORCE_EVAL(0/0.0f);
 		/* in ld80 msb is set in inf */
 		return m & (uint64_t)-1>>1 ? FP_ILOGBNAN : INT_MAX;
+	}
 	return e - 0x3fff;
 }
 #endif
diff --git a/src/math/llrintl.c b/src/math/llrintl.c
index 32bb8b03..3449f6f2 100644
--- a/src/math/llrintl.c
+++ b/src/math/llrintl.c
@@ -18,6 +18,7 @@ raises inexact (with tonearest or upward rounding mode)
 */
 long long llrintl(long double x)
 {
+	#pragma STDC FENV_ACCESS ON
 	int e;
 
 	e = fetestexcept(FE_INEXACT);
diff --git a/src/math/log1pl.c b/src/math/log1pl.c
index 1400d365..edb48df1 100644
--- a/src/math/log1pl.c
+++ b/src/math/log1pl.c
@@ -46,11 +46,6 @@
  *                      Relative error:
  * arithmetic   domain     # trials      peak         rms
  *    IEEE     -1.0, 9.0    100000      8.2e-20    2.5e-20
- *
- * ERROR MESSAGES:
- *
- * log singularity:  x-1 = 0; returns -INFINITY
- * log domain:       x-1 < 0; returns NAN
  */
 
 #include "libm.h"
@@ -123,8 +118,8 @@ long double log1pl(long double xm1)
 	/* Test for domain errors.  */
 	if (x <= 0.0) {
 		if (x == 0.0)
-			return -INFINITY;
-		return NAN;
+			return -1/x; /* -inf with divbyzero */
+		return 0/0.0f; /* nan with invalid */
 	}
 
 	/* Separate mantissa from exponent.
diff --git a/src/math/log2l.c b/src/math/log2l.c
index 8ebce9c4..345b395d 100644
--- a/src/math/log2l.c
+++ b/src/math/log2l.c
@@ -50,11 +50,6 @@
  * In the tests over the interval exp(+-10000), the logarithms
  * of the random arguments were uniformly distributed over
  * [-10000, +10000].
- *
- * ERROR MESSAGES:
- *
- * log singularity:  x = 0; returns -INFINITY
- * log domain:       x < 0; returns NAN
  */
 
 #include "libm.h"
@@ -113,8 +108,7 @@ static const long double S[4] = {
 
 long double log2l(long double x)
 {
-	volatile long double z;
-	long double y;
+	long double y, z;
 	int e;
 
 	if (isnan(x))
@@ -123,8 +117,8 @@ long double log2l(long double x)
 		return x;
 	if (x <= 0.0) {
 		if (x == 0.0)
-			return -INFINITY;
-		return NAN;
+			return -1/(x+0); /* -inf with divbyzero */
+		return 0/0.0f; /* nan with invalid */
 	}
 
 	/* separate mantissa from exponent */
diff --git a/src/math/logb.c b/src/math/logb.c
index f7cd7613..624425a8 100644
--- a/src/math/logb.c
+++ b/src/math/logb.c
@@ -1,20 +1,17 @@
-#include <limits.h>
 #include "libm.h"
 
 /*
 special cases:
-	logb(+-0) = -inf
+	logb(+-0) = -inf, and raise divbyzero
 	logb(+-inf) = +inf
 	logb(nan) = nan
-these are calculated at runtime to raise fp exceptions
 */
 
-double logb(double x) {
-	int i = ilogb(x);
-
-	if (i == FP_ILOGB0)
-		return -1.0/fabs(x);
-	if (i == FP_ILOGBNAN || i == INT_MAX)
+double logb(double x)
+{
+	if (!isfinite(x))
 		return x * x;
-	return i;
+	if (x == 0)
+		return -1/(x+0);
+	return ilogb(x);
 }
diff --git a/src/math/logbf.c b/src/math/logbf.c
index 934827f8..950d3569 100644
--- a/src/math/logbf.c
+++ b/src/math/logbf.c
@@ -1,12 +1,10 @@
-#include <limits.h>
 #include "libm.h"
 
-float logbf(float x) {
-	int i = ilogbf(x);
-
-	if (i == FP_ILOGB0)
-		return -1.0f/fabsf(x);
-	if (i == FP_ILOGBNAN || i == INT_MAX)
+float logbf(float x)
+{
+	if (!isfinite(x))
 		return x * x;
-	return i;
+	if (x == 0)
+		return -1/(x+0);
+	return ilogbf(x);
 }
diff --git a/src/math/logbl.c b/src/math/logbl.c
index 5d04abd7..f3850745 100644
--- a/src/math/logbl.c
+++ b/src/math/logbl.c
@@ -1,4 +1,3 @@
-#include <limits.h>
 #include "libm.h"
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double logbl(long double x)
@@ -8,12 +7,10 @@ long double logbl(long double x)
 #else
 long double logbl(long double x)
 {
-	int i = ilogbl(x);
-
-	if (i == FP_ILOGB0)
-		return -1.0/fabsl(x);
-	if (i == FP_ILOGBNAN || i == INT_MAX)
+	if (!isfinite(x))
 		return x * x;
-	return i;
+	if (x == 0)
+		return -1/(x+0);
+	return ilogbl(x);
 }
 #endif
diff --git a/src/math/logl.c b/src/math/logl.c
index ffd81038..ef2b5515 100644
--- a/src/math/logl.c
+++ b/src/math/logl.c
@@ -50,11 +50,6 @@
  * In the tests over the interval exp(+-10000), the logarithms
  * of the random arguments were uniformly distributed over
  * [-10000, +10000].
- *
- * ERROR MESSAGES:
- *
- * log singularity:  x = 0; returns -INFINITY
- * log domain:       x < 0; returns NAN
  */
 
 #include "libm.h"
@@ -121,8 +116,8 @@ long double logl(long double x)
 		return x;
 	if (x <= 0.0) {
 		if (x == 0.0)
-			return -INFINITY;
-		return NAN;
+			return -1/(x+0); /* -inf with divbyzero */
+		return 0/0.0f; /* nan with invalid */
 	}
 
 	/* separate mantissa from exponent */
diff --git a/src/math/lrint.c b/src/math/lrint.c
index 67091835..bdca8b7c 100644
--- a/src/math/lrint.c
+++ b/src/math/lrint.c
@@ -28,6 +28,7 @@ as a double.
 #if LONG_MAX < 1U<<53 && defined(FE_INEXACT)
 long lrint(double x)
 {
+	#pragma STDC FENV_ACCESS ON
 	int e;
 
 	e = fetestexcept(FE_INEXACT);
diff --git a/src/math/lrintl.c b/src/math/lrintl.c
index 5eb1ba7e..b2a8106d 100644
--- a/src/math/lrintl.c
+++ b/src/math/lrintl.c
@@ -18,6 +18,7 @@ raises inexact (with tonearest or upward rounding mode)
 */
 long lrintl(long double x)
 {
+	#pragma STDC FENV_ACCESS ON
 	int e;
 
 	e = fetestexcept(FE_INEXACT);
diff --git a/src/math/modf.c b/src/math/modf.c
index de45069f..8f337ef0 100644
--- a/src/math/modf.c
+++ b/src/math/modf.c
@@ -1,5 +1,4 @@
-#include <math.h>
-#include <stdint.h>
+#include "libm.h"
 
 double modf(double x, double *iptr)
 {
@@ -33,5 +32,6 @@ double modf(double x, double *iptr)
 	}
 	u.n &= ~mask;
 	*iptr = u.x;
-	return x - *iptr;
+	STRICT_ASSIGN(double, x, x - *iptr);
+	return x;
 }
diff --git a/src/math/modff.c b/src/math/modff.c
index 84d0b82a..90caf316 100644
--- a/src/math/modff.c
+++ b/src/math/modff.c
@@ -1,5 +1,4 @@
-#include <math.h>
-#include <stdint.h>
+#include "libm.h"
 
 float modff(float x, float *iptr)
 {
@@ -33,5 +32,6 @@ float modff(float x, float *iptr)
 	}
 	u.n &= ~mask;
 	*iptr = u.x;
-	return x - *iptr;
+	STRICT_ASSIGN(float, x, x - *iptr);
+	return x;
 }
diff --git a/src/math/nearbyint.c b/src/math/nearbyint.c
index 7a4c58cf..f4e8aac4 100644
--- a/src/math/nearbyint.c
+++ b/src/math/nearbyint.c
@@ -6,6 +6,7 @@
 double nearbyint(double x)
 {
 #ifdef FE_INEXACT
+	#pragma STDC FENV_ACCESS ON
 	int e;
 
 	e = fetestexcept(FE_INEXACT);
diff --git a/src/math/nearbyintf.c b/src/math/nearbyintf.c
index 39c3d73b..092e9ffa 100644
--- a/src/math/nearbyintf.c
+++ b/src/math/nearbyintf.c
@@ -4,6 +4,7 @@
 float nearbyintf(float x)
 {
 #ifdef FE_INEXACT
+	#pragma STDC FENV_ACCESS ON
 	int e;
 
 	e = fetestexcept(FE_INEXACT);
diff --git a/src/math/nearbyintl.c b/src/math/nearbyintl.c
index 0ff4b1f9..82852492 100644
--- a/src/math/nearbyintl.c
+++ b/src/math/nearbyintl.c
@@ -11,6 +11,7 @@ long double nearbyintl(long double x)
 long double nearbyintl(long double x)
 {
 #ifdef FE_INEXACT
+	#pragma STDC FENV_ACCESS ON
 	int e;
 
 	e = fetestexcept(FE_INEXACT);
diff --git a/src/math/nextafter.c b/src/math/nextafter.c
index a3b42c99..9ee82518 100644
--- a/src/math/nextafter.c
+++ b/src/math/nextafter.c
@@ -27,7 +27,7 @@ double nextafter(double x, double y)
 	e = ux.bits >> 52 & 0x7ff;
 	/* raise overflow if ux.value is infinite and x is finite */
 	if (e == 0x7ff)
-		return x + x;
+		FORCE_EVAL(x+x);
 	/* raise underflow if ux.value is subnormal or zero */
 	if (e == 0)
 		FORCE_EVAL(x*x + ux.value*ux.value);
diff --git a/src/math/nextafterf.c b/src/math/nextafterf.c
index b703487b..22b61dce 100644
--- a/src/math/nextafterf.c
+++ b/src/math/nextafterf.c
@@ -26,7 +26,7 @@ float nextafterf(float x, float y)
 	e = ux.bits & 0x7f800000;
 	/* raise overflow if ux.value is infinite and x is finite */
 	if (e == 0x7f800000)
-		return x + x;
+		FORCE_EVAL(x+x);
 	/* raise underflow if ux.value is subnormal or zero */
 	if (e == 0)
 		FORCE_EVAL(x*x + ux.value*ux.value);
diff --git a/src/math/nexttoward.c b/src/math/nexttoward.c
index 7355f2f1..6f32eca4 100644
--- a/src/math/nexttoward.c
+++ b/src/math/nexttoward.c
@@ -36,7 +36,7 @@ double nexttoward(double x, long double y)
 	e = ux.bits>>52 & 0x7ff;
 	/* raise overflow if ux.value is infinite and x is finite */
 	if (e == 0x7ff)
-		return x + x;
+		FORCE_EVAL(x+x);
 	/* raise underflow if ux.value is subnormal or zero */
 	if (e == 0)
 		FORCE_EVAL(x*x + ux.value*ux.value);
diff --git a/src/math/nexttowardf.c b/src/math/nexttowardf.c
index 8648be6a..9a693b1a 100644
--- a/src/math/nexttowardf.c
+++ b/src/math/nexttowardf.c
@@ -28,7 +28,7 @@ float nexttowardf(float x, long double y)
 	e = ux.bits & 0x7f800000;
 	/* raise overflow if ux.value is infinite and x is finite */
 	if (e == 0x7f800000)
-		return x + x;
+		FORCE_EVAL(x+x);
 	/* raise underflow if ux.value is subnormal or zero */
 	if (e == 0)
 		FORCE_EVAL(x*x + ux.value*ux.value);
diff --git a/src/math/scalbn.c b/src/math/scalbn.c
index c9c7af80..003141e3 100644
--- a/src/math/scalbn.c
+++ b/src/math/scalbn.c
@@ -10,8 +10,10 @@ double scalbn(double x, int n)
 		if (n > 1023) {
 			x *= 0x1p1023;
 			n -= 1023;
-			if (n > 1023)
-				return x * 0x1p1023;
+			if (n > 1023) {
+				STRICT_ASSIGN(double, x, x * 0x1p1023);
+				return x;
+			}
 		}
 	} else if (n < -1022) {
 		x *= 0x1p-1022;
@@ -19,10 +21,13 @@ double scalbn(double x, int n)
 		if (n < -1022) {
 			x *= 0x1p-1022;
 			n += 1022;
-			if (n < -1022)
-				return x * 0x1p-1022;
+			if (n < -1022) {
+				STRICT_ASSIGN(double, x, x * 0x1p-1022);
+				return x;
+			}
 		}
 	}
 	INSERT_WORDS(scale, (uint32_t)(0x3ff+n)<<20, 0);
-	return x * scale;
+	STRICT_ASSIGN(double, x, x * scale);
+	return x;
 }
diff --git a/src/math/scalbnf.c b/src/math/scalbnf.c
index 243dafdf..f94b5d59 100644
--- a/src/math/scalbnf.c
+++ b/src/math/scalbnf.c
@@ -10,8 +10,10 @@ float scalbnf(float x, int n)
 		if (n > 127) {
 			x *= 0x1p127f;
 			n -= 127;
-			if (n > 127)
-				return x * 0x1p127f;
+			if (n > 127) {
+				STRICT_ASSIGN(float, x, x * 0x1p127f);
+				return x;
+			}
 		}
 	} else if (n < -126) {
 		x *= 0x1p-126f;
@@ -19,10 +21,13 @@ float scalbnf(float x, int n)
 		if (n < -126) {
 			x *= 0x1p-126f;
 			n += 126;
-			if (n < -126)
-				return x * 0x1p-126f;
+			if (n < -126) {
+				STRICT_ASSIGN(float, x, x * 0x1p-126f);
+				return x;
+			}
 		}
 	}
 	SET_FLOAT_WORD(scale, (uint32_t)(0x7f+n)<<23);
-	return x * scale;
+	STRICT_ASSIGN(float, x, x * scale);
+	return x;
 }