Merge remote branch 'nsz/master'

66 files changed, 413 insertions, 216 deletions

diff --git a/src/math/__expo2.c b/src/math/__expo2.c
index ef14e5f5..740ac680 100644
--- a/src/math/__expo2.c
+++ b/src/math/__expo2.c
@@ -1,51 +1,16 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */
-/*-
- * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
 #include "libm.h"
 
-/*
- * We use exp(x) = exp(x - kln2) * 2**k,
- * k is carefully chosen to minimize |exp(kln2) - 2**k|
- */
-static const uint32_t k = 1799;
-static const double kln2 = 1246.97177782734161156;
+/* k is such that k*ln2 has minimal relative error and x - kln2 > log(DBL_MIN) */
+static const int k = 2043;
+static const double kln2 = 0x1.62066151add8bp+10;
 
-/* exp(x)/2 when x is huge */
+/* exp(x)/2 for x >= log(DBL_MAX), slightly better than 0.5*exp(x/2)*exp(x/2) */
 double __expo2(double x)
 {
 	double scale;
-	int n;
 
-	/*
-	 * efficient scalbn(y, k-1):
-	 * 2**(k-1) cannot be represented
-	 * so we use that k-1 is even and scale in two steps
-	 */
-	n = (k - 1)/2;
-	INSERT_WORDS(scale, (0x3ff + n) << 20, 0);
+	/* note that k is odd and scale*scale overflows */
+	INSERT_WORDS(scale, (uint32_t)(0x3ff + k/2) << 20, 0);
+	/* exp(x - k ln2) * 2**(k-1) */
 	return exp(x - kln2) * scale * scale;
 }
diff --git a/src/math/__expo2f.c b/src/math/__expo2f.c
index 192838f7..5163e418 100644
--- a/src/math/__expo2f.c
+++ b/src/math/__expo2f.c
@@ -1,51 +1,16 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/k_expf.c */
-/*-
- * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
 #include "libm.h"
 
-/*
- * We use exp(x) = exp(x - kln2) * 2**k,
- * k is carefully chosen to minimize |exp(kln2) - 2**k|
- */
-static const uint32_t k = 235;
-static const float kln2 = 162.88958740f;
+/* k is such that k*ln2 has minimal relative error and x - kln2 > log(FLT_MIN) */
+static const int k = 235;
+static const float kln2 = 0x1.45c778p+7f;
 
-/* expf(x)/2 when x is huge */
+/* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */
 float __expo2f(float x)
 {
 	float scale;
-	int n;
 
-	/*
-	 * efficient scalbnf(y, k-1):
-	 * 2**(k-1) cannot be represented
-	 * so we use that k-1 is even and scale in two steps
-	 */
-	n = (k - 1)/2;
-	SET_FLOAT_WORD(scale, (0x7f + n) << 23);
+	/* note that k is odd and scale*scale overflows */
+	SET_FLOAT_WORD(scale, (uint32_t)(0x7f + k/2) << 23);
+	/* exp(x - k ln2) * 2**(k-1) */
 	return expf(x - kln2) * scale * scale;
 }
diff --git a/src/math/__log1pf.h b/src/math/__log1pf.h
index 110acecb..99492c5a 100644
--- a/src/math/__log1pf.h
+++ b/src/math/__log1pf.h
@@ -24,12 +24,12 @@ static inline float __log1pf(float f)
 {
 	float hfsq,s,z,R,w,t1,t2;
 
-	s = f/((float)2.0+f);
+	s = f/(2.0f + f);
 	z = s*s;
 	w = z*z;
 	t1 = w*(Lg2+w*Lg4);
 	t2 = z*(Lg1+w*Lg3);
 	R = t2+t1;
-	hfsq = (float)0.5*f*f;
+	hfsq = 0.5f * f * f;
 	return s*(hfsq+R);
 }
diff --git a/src/math/acosf.c b/src/math/acosf.c
index dd3bba29..b4665d02 100644
--- a/src/math/acosf.c
+++ b/src/math/acosf.c
@@ -36,8 +36,8 @@ float acosf(float x)
 	ix = hx & 0x7fffffff;
 	if (ix >= 0x3f800000) {  /* |x| >= 1 */
 		if (ix == 0x3f800000) {  /* |x| == 1 */
-			if(hx>0) return 0.0;  /* acos(1) = 0 */
-			return pi + (float)2.0*pio2_lo;  /* acos(-1)= pi */
+			if (hx > 0) return 0.0f;  /* acos(1) = 0 */
+			return pi + 2.0f*pio2_lo;  /* acos(-1)= pi */
 		}
 		return (x-x)/(x-x);  /* acos(|x|>1) is NaN */
 	}
@@ -50,17 +50,17 @@ float acosf(float x)
 		r = p/q;
 		return pio2_hi - (x - (pio2_lo-x*r));
 	} else if (hx < 0) {     /* x < -0.5 */
-		z = (one+x)*(float)0.5;
+		z = (one+x)*0.5f;
 		p = z*(pS0+z*(pS1+z*pS2));
 		q = one+z*qS1;
 		s = sqrtf(z);
 		r = p/q;
 		w = r*s-pio2_lo;
-		return pi - (float)2.0*(s+w);
+		return pi - 2.0f*(s+w);
 	} else {                 /* x > 0.5 */
 		int32_t idf;
 
-		z = (one-x)*(float)0.5;
+		z = (one-x)*0.5f;
 		s = sqrtf(z);
 		df = s;
 		GET_FLOAT_WORD(idf,df);
@@ -70,6 +70,6 @@ float acosf(float x)
 		q = one+z*qS1;
 		r = p/q;
 		w = r*s+c;
-		return (float)2.0*(df+w);
+		return 2.0f*(df+w);
 	}
 }
diff --git a/src/math/acoshf.c b/src/math/acoshf.c
index 30a3a943..57ce5cb8 100644
--- a/src/math/acoshf.c
+++ b/src/math/acoshf.c
@@ -35,9 +35,9 @@ float acoshf(float x)
 		return 0.0;  /* acosh(1) = 0 */
 	} else if (hx > 0x40000000) {  /* 2**28 > x > 2 */
 		t = x*x;
-		return logf((float)2.0*x - one/(x+sqrtf(t-one)));
+		return logf(2.0f*x - one/(x+sqrtf(t-one)));
 	} else {                /* 1 < x < 2 */
 		t = x-one;
-		return log1pf(t + sqrtf((float)2.0*t+t*t));
+		return log1pf(t + sqrtf(2.0f*t+t*t));
 	}
 }
diff --git a/src/math/asinf.c b/src/math/asinf.c
index 729dd37f..7865a45b 100644
--- a/src/math/asinf.c
+++ b/src/math/asinf.c
@@ -52,7 +52,7 @@ float asinf(float x)
 	}
 	/* 1 > |x| >= 0.5 */
 	w = one - fabsf(x);
-	t = w*(float)0.5;
+	t = w*0.5f;
 	p = t*(pS0+t*(pS1+t*pS2));
 	q = one+t*qS1;
 	s = sqrt(t);
diff --git a/src/math/asinhf.c b/src/math/asinhf.c
index 5f4bb39c..cb6e5b89 100644
--- a/src/math/asinhf.c
+++ b/src/math/asinhf.c
@@ -38,7 +38,7 @@ float asinhf(float x)
 		w = logf(fabsf(x)) + ln2;
 	} else if (ix > 0x40000000) {  /* 2**28 > |x| > 2.0 */
 		t = fabsf(x);
-		w = logf((float)2.0*t + one/(sqrtf(x*x+one)+t));
+		w = logf(2.0f*t + one/(sqrtf(x*x+one)+t));
 	} else {                /* 2.0 > |x| > 2**-28 */
 		t = x*x;
 		w =log1pf(fabsf(x) + t/(one+sqrtf(one+t)));
diff --git a/src/math/atan2f.c b/src/math/atan2f.c
index 4d78840b..b6a7f92d 100644
--- a/src/math/atan2f.c
+++ b/src/math/atan2f.c
@@ -58,8 +58,8 @@ float atan2f(float y, float x)
 			switch (m) {
 			case 0: return  pi_o_4+tiny; /* atan(+INF,+INF) */
 			case 1: return -pi_o_4-tiny; /* atan(-INF,+INF) */
-			case 2: return (float)3.0*pi_o_4+tiny;  /*atan(+INF,-INF)*/
-			case 3: return (float)-3.0*pi_o_4-tiny; /*atan(-INF,-INF)*/
+			case 2: return 3.0f*pi_o_4+tiny;  /*atan(+INF,-INF)*/
+			case 3: return -3.0f*pi_o_4-tiny; /*atan(-INF,-INF)*/
 			}
 		} else {
 			switch (m) {
@@ -77,7 +77,7 @@ float atan2f(float y, float x)
 	/* compute y/x */
 	k = (iy-ix)>>23;
 	if (k > 26) {                  /* |y/x| >  2**26 */
-		z = pi_o_2+(float)0.5*pi_lo;
+		z = pi_o_2 + 0.5f*pi_lo;
 		m &= 1;
 	} else if (k < -26 && hx < 0)  /* 0 > |y|/x > -2**-26 */
 		z = 0.0;
diff --git a/src/math/atanf.c b/src/math/atanf.c
index 8c2b46b0..73cebb5e 100644
--- a/src/math/atanf.c
+++ b/src/math/atanf.c
@@ -69,18 +69,18 @@ float atanf(float x)
 		if (ix < 0x3f980000) {  /* |x| < 1.1875 */
 			if (ix < 0x3f300000) {  /*  7/16 <= |x| < 11/16 */
 				id = 0;
-				x = ((float)2.0*x-one)/((float)2.0+x);
+				x = (2.0f*x - one)/(2.0f + x);
 			} else {                /* 11/16 <= |x| < 19/16 */
 				id = 1;
-				x = (x-one)/(x+one);
+				x = (x - one)/(x + one);
 			}
 		} else {
 			if (ix < 0x401c0000) {  /* |x| < 2.4375 */
 				id = 2;
-				x = (x-(float)1.5)/(one+(float)1.5*x);
+				x = (x - 1.5f)/(one + 1.5f*x);
 			} else {                /* 2.4375 <= |x| < 2**26 */
 				id = 3;
-				x = -(float)1.0/x;
+				x = -1.0f/x;
 			}
 		}
 	}
diff --git a/src/math/atanhf.c b/src/math/atanhf.c
index 2efbd79c..9c65dc7f 100644
--- a/src/math/atanhf.c
+++ b/src/math/atanhf.c
@@ -34,9 +34,9 @@ float atanhf(float x)
 	SET_FLOAT_WORD(x, ix);
 	if (ix < 0x3f000000) {                 /* x < 0.5 */
 		t = x+x;
-		t = (float)0.5*log1pf(t + t*x/(one-x));
+		t = 0.5f*log1pf(t + t*x/(one-x));
 	} else
-		t = (float)0.5*log1pf((x+x)/(one-x));
+		t = 0.5f*log1pf((x+x)/(one-x));
 	if (hx >= 0)
 		return t;
 	return -t;
diff --git a/src/math/ceilf.c b/src/math/ceilf.c
index d83066a5..d22688a7 100644
--- a/src/math/ceilf.c
+++ b/src/math/ceilf.c
@@ -27,7 +27,7 @@ float ceilf(float x)
 	if (j0 < 23) {
 		if (j0 < 0) {
 			/* raise inexact if x != 0 */
-			if (huge+x > (float)0.0) {
+			if (huge+x > 0.0f) {
 				/* return 0*sign(x) if |x|<1 */
 				if (i0 < 0)
 					i0 = 0x80000000;
@@ -39,7 +39,7 @@ float ceilf(float x)
 			if ((i0&i) == 0)
 				return x; /* x is integral */
 			/* raise inexact flag */
-			if (huge+x > (float)0.0) {
+			if (huge+x > 0.0f) {
 				if (i0 > 0)
 					i0 += 0x00800000>>j0;
 				i0 &= ~i;
diff --git a/src/math/erff.c b/src/math/erff.c
index e4e353d7..eef4851e 100644
--- a/src/math/erff.c
+++ b/src/math/erff.c
@@ -106,7 +106,7 @@ float erff(float x)
 		if (ix < 0x31800000) {  /* |x| < 2**-28 */
 			if (ix < 0x04000000)
 				/*avoid underflow */
-				return (float)0.125*((float)8.0*x+efx8*x);
+				return 0.125f*(8.0f*x + efx8*x);
 			return x + efx*x;
 		}
 		z = x*x;
@@ -143,7 +143,7 @@ float erff(float x)
 	}
 	GET_FLOAT_WORD(ix, x);
 	SET_FLOAT_WORD(z, ix&0xfffff000);
-	r = expf(-z*z - (float)0.5625) * expf((z-x)*(z+x) + R/S);
+	r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
 	if (hx >= 0)
 		return one - r/x;
 	return  r/x - one;
@@ -206,7 +206,7 @@ float erfcf(float x)
 		}
 		GET_FLOAT_WORD(ix, x);
 		SET_FLOAT_WORD(z, ix&0xfffff000);
-		r = expf(-z*z - (float)0.5625) * expf((z-x)*(z+x) + R/S);
+		r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
 		if (hx > 0)
 			return r/x;
 		return two - r/x;
diff --git a/src/math/expf.c b/src/math/expf.c
index a0eaa7a7..8925a6f1 100644
--- a/src/math/expf.c
+++ b/src/math/expf.c
@@ -85,11 +85,11 @@ float expf(float x)
 		SET_FLOAT_WORD(twopk, 0x3f800000+((k+100)<<23));
 	c  = x - t*(P1+t*P2);
 	if (k == 0)
-		return one - ((x*c)/(c-(float)2.0)-x);
-	y = one - ((lo-(x*c)/((float)2.0-c))-hi);
+		return one - ((x*c)/(c - 2.0f) - x);
+	y = one - ((lo - (x*c)/(2.0f - c)) - hi);
 	if (k < -125)
 		return y*twopk*twom100;
 	if (k == 128)
-		return y*2.0F*0x1p127F;
+		return y*2.0f*0x1p127f;
 	return y*twopk;
 }
diff --git a/src/math/expm1f.c b/src/math/expm1f.c
index cfab6975..a8b79e88 100644
--- a/src/math/expm1f.c
+++ b/src/math/expm1f.c
@@ -53,7 +53,7 @@ float expm1f(float x)
 		}
 		if (xsb != 0) {  /* x < -27*ln2 */
 			/* raise inexact */
-			if (x+tiny < (float)0.0)
+			if (x+tiny < 0.0f)
 				return tiny-one;  /* return -1 */
 		}
 	}
@@ -71,7 +71,7 @@ float expm1f(float x)
 				k = -1;
 			}
 		} else {
-			k  = invln2*x+((xsb==0)?(float)0.5:(float)-0.5);
+			k  = invln2*x + (xsb==0 ? 0.5f : -0.5f);
 			t  = k;
 			hi = x - t*ln2_hi;      /* t*ln2_hi is exact here */
 			lo = t*ln2_lo;
@@ -85,27 +85,27 @@ float expm1f(float x)
 		k = 0;
 
 	/* x is now in primary range */
-	hfx = (float)0.5*x;
+	hfx = 0.5f*x;
 	hxs = x*hfx;
 	r1 = one+hxs*(Q1+hxs*Q2);
-	t  = (float)3.0 - r1*hfx;
-	e  = hxs*((r1-t)/((float)6.0 - x*t));
+	t  = 3.0f - r1*hfx;
+	e  = hxs*((r1-t)/(6.0f - x*t));
 	if (k == 0)  /* c is 0 */
 		return x - (x*e-hxs);
 	SET_FLOAT_WORD(twopk, 0x3f800000+(k<<23));   /* 2^k */
 	e  = x*(e-c) - c;
 	e -= hxs;
 	if (k == -1)
-		return (float)0.5*(x-e) - (float)0.5;
+		return 0.5f*(x-e) - 0.5f;
 	if (k == 1) {
-		if (x < (float)-0.25)
-			return -(float)2.0*(e-(x+(float)0.5));
-		return one+(float)2.0*(x-e);
+		if (x < -0.25f)
+			return -2.0f*(e-(x+0.5f));
+		return one + 2.0f*(x-e);
 	}
 	if (k <= -2 || k > 56) {   /* suffice to return exp(x)-1 */
 		y = one - (e - x);
 		if (k == 128)
-			y = y*2.0F*0x1p127F;
+			y = y*2.0f*0x1p127f;
 		else
 			y = y*twopk;
 		return y - one;
@@ -116,7 +116,7 @@ float expm1f(float x)
 		y = t - (e - x);
 		y = y*twopk;
 	} else {
-		SET_FLOAT_WORD(t, ((0x7f-k)<<23));  /* 2^-k */
+		SET_FLOAT_WORD(t, (0x7f-k)<<23);  /* 2^-k */
 		y = x - (e + t);
 		y += one;
 		y = y*twopk;
diff --git a/src/math/fdim.c b/src/math/fdim.c
index fb25521c..95854606 100644
--- a/src/math/fdim.c
+++ b/src/math/fdim.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 double fdim(double x, double y)
 {
diff --git a/src/math/fdimf.c b/src/math/fdimf.c
index 5cfeac6b..543c3648 100644
--- a/src/math/fdimf.c
+++ b/src/math/fdimf.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 float fdimf(float x, float y)
 {
diff --git a/src/math/fdiml.c b/src/math/fdiml.c
index cda3022e..62e29b7d 100644
--- a/src/math/fdiml.c
+++ b/src/math/fdiml.c
@@ -1,4 +1,5 @@
-#include "libm.h"
+#include <math.h>
+#include <float.h>
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double fdiml(long double x, long double y)
diff --git a/src/math/floorf.c b/src/math/floorf.c
index 958abf5b..23fa5e7d 100644
--- a/src/math/floorf.c
+++ b/src/math/floorf.c
@@ -36,7 +36,7 @@ float floorf(float x)
 	if (j0 < 23) {
 		if (j0 < 0) {  /* |x| < 1 */
 			/* raise inexact if x != 0 */
-			if (huge+x > (float)0.0) {
+			if (huge+x > 0.0f) {
 				if (i0 >= 0)  /* x >= 0 */
 					i0 = 0;
 				else if ((i0&0x7fffffff) != 0)
@@ -47,7 +47,7 @@ float floorf(float x)
 			if ((i0&i) == 0)
 				return x; /* x is integral */
 			/* raise inexact flag */
-			if (huge+x > (float)0.0) {
+			if (huge+x > 0.0f) {
 				if (i0 < 0)
 					i0 += 0x00800000>>j0;
 				i0 &= ~i;
diff --git a/src/math/fmax.c b/src/math/fmax.c
index 0b6bf6f3..94f0caa1 100644
--- a/src/math/fmax.c
+++ b/src/math/fmax.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 double fmax(double x, double y)
 {
diff --git a/src/math/fmaxf.c b/src/math/fmaxf.c
index 7767c303..695d8179 100644
--- a/src/math/fmaxf.c
+++ b/src/math/fmaxf.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 float fmaxf(float x, float y)
 {
diff --git a/src/math/fmaxl.c b/src/math/fmaxl.c
index 8a1e3652..4b03158e 100644
--- a/src/math/fmaxl.c
+++ b/src/math/fmaxl.c
@@ -1,4 +1,5 @@
-#include "libm.h"
+#include <math.h>
+#include <float.h>
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double fmaxl(long double x, long double y)
diff --git a/src/math/fmin.c b/src/math/fmin.c
index d1f16454..08a8fd17 100644
--- a/src/math/fmin.c
+++ b/src/math/fmin.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 double fmin(double x, double y)
 {
diff --git a/src/math/fminf.c b/src/math/fminf.c
index 0964cdb3..3573c7de 100644
--- a/src/math/fminf.c
+++ b/src/math/fminf.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 float fminf(float x, float y)
 {
diff --git a/src/math/fminl.c b/src/math/fminl.c
index ae7159a5..69bc24a7 100644
--- a/src/math/fminl.c
+++ b/src/math/fminl.c
@@ -1,4 +1,5 @@
-#include "libm.h"
+#include <math.h>
+#include <float.h>
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double fminl(long double x, long double y)
diff --git a/src/math/hypotf.c b/src/math/hypotf.c
index 40acd917..9fd77e6a 100644
--- a/src/math/hypotf.c
+++ b/src/math/hypotf.c
@@ -40,7 +40,7 @@ float hypotf(float x, float y)
 	if (ha > 0x58800000) {    /* a > 2**50 */
 		if(ha >= 0x7f800000) {  /* Inf or NaN */
 			/* Use original arg order iff result is NaN; quieten sNaNs. */
-			w = fabsf(x+0.0F) - fabsf(y+0.0F);
+			w = fabsf(x+0.0f) - fabsf(y+0.0f);
 			if (ha == 0x7f800000) w = a;
 			if (hb == 0x7f800000) w = b;
 			return w;
diff --git a/src/math/ilogb.c b/src/math/ilogb.c
index c5915a0c..0a3a6a46 100644
--- a/src/math/ilogb.c
+++ b/src/math/ilogb.c
@@ -11,7 +11,6 @@ int ilogb(double x)
 		if (u.bits == 0)
 			return FP_ILOGB0;
 		/* subnormal x */
-		// FIXME: scale up subnormals with a *0x1p53 or find top set bit with a better method
 		for (e = -0x3ff; u.bits < (uint64_t)1<<63; e--, u.bits<<=1);
 		return e;
 	}
diff --git a/src/math/j0f.c b/src/math/j0f.c
index 77a2d734..52cb0ab4 100644
--- a/src/math/j0f.c
+++ b/src/math/j0f.c
@@ -74,16 +74,16 @@ float j0f(float x)
 		if (huge+x > one) {
 			if (ix < 0x32000000)  /* |x| < 2**-27 */
 				return one;
-			return one - (float)0.25*x*x;
+			return one - 0.25f*x*x;
 		}
 	}
 	z = x*x;
 	r =  z*(R02+z*(R03+z*(R04+z*R05)));
 	s =  one+z*(S01+z*(S02+z*(S03+z*S04)));
 	if(ix < 0x3F800000) {   /* |x| < 1.00 */
-		return one + z*((float)-0.25+(r/s));
+		return one + z*(-0.25f + (r/s));
 	} else {
-		u = (float)0.5*x;
+		u = 0.5f*x;
 		return (one+u)*(one-u) + z*(r/s);
 	}
 }
@@ -343,5 +343,5 @@ static float qzerof(float x)
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
-	return (-(float).125 + r/s)/x;
+	return (-.125f + r/s)/x;
 }
diff --git a/src/math/j1f.c b/src/math/j1f.c
index 0323ec78..2d617b67 100644
--- a/src/math/j1f.c
+++ b/src/math/j1f.c
@@ -75,13 +75,13 @@ float j1f(float x)
 	if (ix < 0x32000000) {  /* |x| < 2**-27 */
 		/* raise inexact if x!=0 */
 		if (huge+x > one)
-			return (float)0.5*x;
+			return 0.5f*x;
 	}
 	z = x*x;
 	r = z*(r00+z*(r01+z*(r02+z*r03)));
 	s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
 	r *= x;
-	return x*(float)0.5 + r/s;
+	return 0.5f*x + r/s;
 }
 
 static const float U0[5] = {
@@ -338,5 +338,5 @@ static float qonef(float x)
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
-	return ((float).375 + r/s)/x;
+	return (.375f + r/s)/x;
 }
diff --git a/src/math/jnf.c b/src/math/jnf.c
index 7db93ae7..648db32b 100644
--- a/src/math/jnf.c
+++ b/src/math/jnf.c
@@ -64,7 +64,7 @@ float jnf(int n, float x)
 			if (n > 33)  /* underflow */
 				b = zero;
 			else {
-				temp = x*(float)0.5;
+				temp = 0.5f * x;
 				b = temp;
 				for (a=one,i=2; i<=n; i++) {
 					a *= (float)i;    /* a = n! */
@@ -106,13 +106,13 @@ float jnf(int n, float x)
 			float q0,q1,h,tmp;
 			int32_t k,m;
 
-			w = (n+n)/(float)x;
-			h = (float)2.0/(float)x;
+			w = (n+n)/x;
+			h = 2.0f/x;
 			z = w+h;
 			q0 = w;
-			q1 = w*z - (float)1.0;
+			q1 = w*z - 1.0f;
 			k = 1;
-			while (q1 < (float)1.0e9) {
+			while (q1 < 1.0e9f) {
 				k += 1;
 				z += h;
 				tmp = z*q1 - q0;
@@ -135,7 +135,7 @@ float jnf(int n, float x)
 			tmp = n;
 			v = two/x;
 			tmp = tmp*logf(fabsf(v*tmp));
-			if (tmp < (float)8.8721679688e+01) {
+			if (tmp < 88.721679688f) {
 				for (i=n-1,di=(float)(i+i); i>0; i--) {
 					temp = b;
 					b *= di;
@@ -151,7 +151,7 @@ float jnf(int n, float x)
 					a = temp;
 					di -= two;
 					/* scale b to avoid spurious overflow */
-					if (b > (float)1e10) {
+					if (b > 1e10f) {
 						a /= b;
 						t /= b;
 						b = one;
diff --git a/src/math/ldexp.c b/src/math/ldexp.c
index 36835dba..f4d1cd6a 100644
--- a/src/math/ldexp.c
+++ b/src/math/ldexp.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 double ldexp(double x, int n)
 {
diff --git a/src/math/ldexpf.c b/src/math/ldexpf.c
index f0981ae4..3bad5f39 100644
--- a/src/math/ldexpf.c
+++ b/src/math/ldexpf.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 float ldexpf(float x, int n)
 {
diff --git a/src/math/ldexpl.c b/src/math/ldexpl.c
index 885ff6e9..fd145ccc 100644
--- a/src/math/ldexpl.c
+++ b/src/math/ldexpl.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 long double ldexpl(long double x, int n)
 {
diff --git a/src/math/lgamma.c b/src/math/lgamma.c
index d12462b9..9af7eee4 100644
--- a/src/math/lgamma.c
+++ b/src/math/lgamma.c
@@ -1,3 +1,4 @@
+#define _GNU_SOURCE
 #include "libm.h"
 
 double lgamma(double x)
diff --git a/src/math/lgammaf.c b/src/math/lgammaf.c
index f50f2379..aed98ba4 100644
--- a/src/math/lgammaf.c
+++ b/src/math/lgammaf.c
@@ -1,3 +1,4 @@
+#define _GNU_SOURCE
 #include "libm.h"
 
 float lgammaf(float x)
diff --git a/src/math/lgammaf_r.c b/src/math/lgammaf_r.c
index 9955b2f9..c6280f5b 100644
--- a/src/math/lgammaf_r.c
+++ b/src/math/lgammaf_r.c
@@ -104,9 +104,9 @@ static float sin_pif(float x)
 	 */
 	z = floorf(y);
 	if (z != y) {   /* inexact anyway */
-		y *= (float)0.5;
-		y  = (float)2.0*(y - floorf(y));   /* y = |x| mod 2.0 */
-		n  = (int) (y*(float)4.0);
+		y *= 0.5f;
+		y  = 2.0f*(y - floorf(y));   /* y = |x| mod 2.0 */
+		n  = (int)(y*4.0f);
 	} else {
 		if (ix >= 0x4b800000) {
 			y = zero;  /* y must be even */
@@ -123,12 +123,12 @@ static float sin_pif(float x)
 	switch (n) {
 	case 0:  y =  __sindf(pi*y); break;
 	case 1:
-	case 2:  y =  __cosdf(pi*((float)0.5-y)); break;
+	case 2:  y =  __cosdf(pi*(0.5f - y)); break;
 	case 3:
-	case 4:  y =  __sindf(pi*(one-y)); break;
+	case 4:  y =  __sindf(pi*(one - y)); break;
 	case 5:
-	case 6:  y = -__cosdf(pi*(y-(float)1.5)); break;
-	default: y =  __sindf(pi*(y-(float)2.0)); break;
+	case 6:  y = -__cosdf(pi*(y - 1.5f)); break;
+	default: y =  __sindf(pi*(y - 2.0f)); break;
 	}
 	return -y;
 }
@@ -188,7 +188,7 @@ float lgammaf_r(float x, int *signgamp)
 		} else {
 			r = zero;
 			if (ix >= 0x3fdda618) {  /* [1.7316,2] */
-				y = (float)2.0 - x;
+				y = 2.0f - x;
 				i = 0;
 			} else if (ix >= 0x3F9da620) {  /* [1.23,1.73] */
 				y = x - tc;
@@ -204,7 +204,7 @@ float lgammaf_r(float x, int *signgamp)
 			p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
 			p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
 			p = y*p1+p2;
-			r += (p-(float)0.5*y);
+			r += p - 0.5f*y;
 			break;
 		case 1:
 			z = y*y;
@@ -218,21 +218,21 @@ float lgammaf_r(float x, int *signgamp)
 		case 2:
 			p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
 			p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
-			r += (-(float)0.5*y + p1/p2);
+			r += -0.5f*y + p1/p2;
 		}
 	} else if (ix < 0x41000000) {  /* x < 8.0 */
 		i = (int)x;
-		y = x-(float)i;
+		y = x - (float)i;
 		p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
 		q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
 		r = half*y+p/q;
 		z = one;    /* lgamma(1+s) = log(s) + lgamma(s) */
 		switch (i) {
-		case 7: z *= y + (float)6.0;  /* FALLTHRU */
-		case 6: z *= y + (float)5.0;  /* FALLTHRU */
-		case 5: z *= y + (float)4.0;  /* FALLTHRU */
-		case 4: z *= y + (float)3.0;  /* FALLTHRU */
-		case 3: z *= y + (float)2.0;  /* FALLTHRU */
+		case 7: z *= y + 6.0f;  /* FALLTHRU */
+		case 6: z *= y + 5.0f;  /* FALLTHRU */
+		case 5: z *= y + 4.0f;  /* FALLTHRU */
+		case 4: z *= y + 3.0f;  /* FALLTHRU */
+		case 3: z *= y + 2.0f;  /* FALLTHRU */
 			r += logf(z);
 			break;
 		}
diff --git a/src/math/lgammal.c b/src/math/lgammal.c
index 603477c9..a33707ad 100644
--- a/src/math/lgammal.c
+++ b/src/math/lgammal.c
@@ -85,6 +85,7 @@
  *
  */
 
+#define _GNU_SOURCE
 #include "libm.h"
 
 long double lgammal(long double x)
diff --git a/src/math/llrintl.c b/src/math/llrintl.c
index 6b0838d4..ec2cf86b 100644
--- a/src/math/llrintl.c
+++ b/src/math/llrintl.c
@@ -1,4 +1,6 @@
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long long llrintl(long double x)
 {
diff --git a/src/math/llroundl.c b/src/math/llroundl.c
index 9e2cfdc7..107744a9 100644
--- a/src/math/llroundl.c
+++ b/src/math/llroundl.c
@@ -1,4 +1,6 @@
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long long llroundl(long double x)
 {
diff --git a/src/math/log10f.c b/src/math/log10f.c
index 4175cce2..7e4ac9a8 100644
--- a/src/math/log10f.c
+++ b/src/math/log10f.c
@@ -53,8 +53,8 @@ float log10f(float x)
 	SET_FLOAT_WORD(x, hx|(i^0x3f800000));  /* normalize x or x/2 */
 	k += i>>23;
 	y = (float)k;
-	f = x - (float)1.0;
-	hfsq = (float)0.5*f*f;
+	f = x - 1.0f;
+	hfsq = 0.5f * f * f;
 	r = __log1pf(f);
 
 // FIXME
diff --git a/src/math/log1pf.c b/src/math/log1pf.c
index 5c718152..75eeb371 100644
--- a/src/math/log1pf.c
+++ b/src/math/log1pf.c
@@ -40,7 +40,7 @@ float log1pf(float x)
 	k = 1;
 	if (hx < 0x3ed413d0) {  /* 1+x < sqrt(2)+  */
 		if (ax >= 0x3f800000) {  /* x <= -1.0 */
-			if (x == (float)-1.0)
+			if (x == -1.0f)
 				return -two25/zero; /* log1p(-1)=+inf */
 			return (x-x)/(x-x);         /* log1p(x<-1)=NaN */
 		}
@@ -48,7 +48,7 @@ float log1pf(float x)
 			/* raise inexact */
 			if (two25 + x > zero && ax < 0x33800000)  /* |x| < 2**-24 */
 				return x;
-			return x - x*x*(float)0.5;
+			return x - x*x*0.5f;
 		}
 		if (hx > 0 || hx <= (int32_t)0xbe95f619) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
 			k = 0;
@@ -60,11 +60,11 @@ float log1pf(float x)
 		return x+x;
 	if (k != 0) {
 		if (hx < 0x5a000000) {
-			STRICT_ASSIGN(float, u, (float)1.0 + x);
+			STRICT_ASSIGN(float, u, 1.0f + x);
 			GET_FLOAT_WORD(hu, u);
 			k = (hu>>23) - 127;
 			/* correction term */
-			c = k > 0 ? (float)1.0-(u-x) : x-(u-(float)1.0);
+			c = k > 0 ? 1.0f-(u-x) : x-(u-1.0f);
 			c /= u;
 		} else {
 			u = x;
@@ -87,9 +87,9 @@ float log1pf(float x)
 			SET_FLOAT_WORD(u, hu|0x3f000000);  /* normalize u/2 */
 			hu = (0x00800000-hu)>>2;
 		}
-		f = u - (float)1.0;
+		f = u - 1.0f;
 	}
-	hfsq = (float)0.5*f*f;
+	hfsq = 0.5f * f * f;
 	if (hu == 0) {  /* |f| < 2**-20 */
 		if (f == zero) {
 			if (k == 0)
@@ -97,12 +97,12 @@ float log1pf(float x)
 			c += k*ln2_lo;
 			return k*ln2_hi+c;
 		}
-		R = hfsq*((float)1.0-(float)0.66666666666666666*f);
+		R = hfsq*(1.0f - 0.66666666666666666f * f);
 		if (k == 0)
 			return f - R;
 		return k*ln2_hi - ((R-(k*ln2_lo+c))-f);
 	}
-	s = f/((float)2.0+f);
+	s = f/(2.0f + f);
 	z = s*s;
 	R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
 	if (k == 0)
diff --git a/src/math/log2f.c b/src/math/log2f.c
index a968984d..c9b9282c 100644
--- a/src/math/log2f.c
+++ b/src/math/log2f.c
@@ -51,8 +51,8 @@ float log2f(float x)
 	SET_FLOAT_WORD(x, hx|(i^0x3f800000));  /* normalize x or x/2 */
 	k += i>>23;
 	y = (float)k;
-	f = x - (float)1.0;
-	hfsq = (float)0.5*f*f;
+	f = x - 1.0f;
+	hfsq = 0.5f * f * f;
 	r = __log1pf(f);
 
 	/*
diff --git a/src/math/logf.c b/src/math/logf.c
index 285ee615..a4ed697b 100644
--- a/src/math/logf.c
+++ b/src/math/logf.c
@@ -52,7 +52,7 @@ float logf(float x)
 	i = (ix + (0x95f64<<3)) & 0x800000;
 	SET_FLOAT_WORD(x, ix|(i^0x3f800000));  /* normalize x or x/2 */
 	k += i>>23;
-	f = x - (float)1.0;
+	f = x - 1.0f;
 	if ((0x007fffff & (0x8000 + ix)) < 0xc000) {  /* -2**-9 <= f < 2**-9 */
 		if (f == zero) {
 			if (k == 0)
@@ -60,13 +60,13 @@ float logf(float x)
 			dk = (float)k;
 			return dk*ln2_hi + dk*ln2_lo;
 		}
-		R = f*f*((float)0.5 - (float)0.33333333333333333*f);
+		R = f*f*(0.5f - 0.33333333333333333f*f);
 		if (k == 0)
 			return f-R;
 		dk = (float)k;
 		return dk*ln2_hi - ((R-dk*ln2_lo)-f);
 	}
-	s = f/((float)2.0+f);
+	s = f/(2.0f + f);
 	dk = (float)k;
 	z = s*s;
 	i = ix-(0x6147a<<3);
@@ -77,7 +77,7 @@ float logf(float x)
 	i |= j;
 	R = t2 + t1;
 	if (i > 0) {
-		hfsq = (float)0.5*f*f;
+		hfsq = 0.5f * f * f;
 		if (k == 0)
 			return f - (hfsq-s*(hfsq+R));
 		return dk*ln2_hi - ((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
diff --git a/src/math/lrintl.c b/src/math/lrintl.c
index 7c09653e..c076326f 100644
--- a/src/math/lrintl.c
+++ b/src/math/lrintl.c
@@ -1,4 +1,6 @@
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long lrintl(long double x)
 {
diff --git a/src/math/lroundl.c b/src/math/lroundl.c
index 1469127b..7b593f77 100644
--- a/src/math/lroundl.c
+++ b/src/math/lroundl.c
@@ -1,4 +1,6 @@
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long lroundl(long double x)
 {
diff --git a/src/math/nearbyint.c b/src/math/nearbyint.c
index 781769fb..714c55ca 100644
--- a/src/math/nearbyint.c
+++ b/src/math/nearbyint.c
@@ -1,5 +1,5 @@
 #include <fenv.h>
-#include "libm.h"
+#include <math.h>
 
 /*
 rint may raise inexact (and it should not alter the fenv otherwise)
diff --git a/src/math/nearbyintf.c b/src/math/nearbyintf.c
index e4bdb26c..07df8f54 100644
--- a/src/math/nearbyintf.c
+++ b/src/math/nearbyintf.c
@@ -1,5 +1,5 @@
 #include <fenv.h>
-#include "libm.h"
+#include <math.h>
 
 float nearbyintf(float x) {
 	fenv_t e;
diff --git a/src/math/nearbyintl.c b/src/math/nearbyintl.c
index b58527c8..2906f383 100644
--- a/src/math/nearbyintl.c
+++ b/src/math/nearbyintl.c
@@ -1,4 +1,6 @@
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double nearbyintl(long double x)
 {
diff --git a/src/math/nexttoward.c b/src/math/nexttoward.c
index 5e12c48b..e150eaad 100644
--- a/src/math/nexttoward.c
+++ b/src/math/nexttoward.c
@@ -11,6 +11,7 @@
  */
 
 #include "libm.h"
+
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 double nexttoward(double x, long double y)
 {
diff --git a/src/math/nexttowardl.c b/src/math/nexttowardl.c
index c393ce97..67a63403 100644
--- a/src/math/nexttowardl.c
+++ b/src/math/nexttowardl.c
@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 long double nexttowardl(long double x, long double y)
 {
diff --git a/src/math/pow.c b/src/math/pow.c
index f843645d..5bcedd5b 100644
--- a/src/math/pow.c
+++ b/src/math/pow.c
@@ -112,7 +112,7 @@ double pow(double x, double y)
 	/* y != zero: result is NaN if either arg is NaN */
 	if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0) ||
 	    iy > 0x7ff00000 || (iy == 0x7ff00000 && ly != 0))
-		return (x+0.0)+(y+0.0); // FIXME: x+y ?
+		return (x+0.0) + (y+0.0);
 
 	/* determine if y is an odd int when x < 0
 	 * yisint = 0       ... y is not an integer
diff --git a/src/math/powf.c b/src/math/powf.c
index e322ff28..01aced0e 100644
--- a/src/math/powf.c
+++ b/src/math/powf.c
@@ -70,7 +70,7 @@ float powf(float x, float y)
 
 	/* y != zero: result is NaN if either arg is NaN */
 	if (ix > 0x7f800000 || iy > 0x7f800000)
-		return (x+0.0F) + (y+0.0F);
+		return (x+0.0f) + (y+0.0f);
 
 	/* determine if y is an odd int when x < 0
 	 * yisint = 0       ... y is not an integer
@@ -145,7 +145,7 @@ float powf(float x, float y)
 		/* now |1-x| is tiny <= 2**-20, suffice to compute
 		   log(x) by x-x^2/2+x^3/3-x^4/4 */
 		t = ax - 1;     /* t has 20 trailing zeros */
-		w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25));
+		w = (t*t)*(0.5f - t*(0.333333333333f - t*0.25f));
 		u = ivln2_h*t;  /* ivln2_h has 16 sig. bits */
 		v = t*ivln2_l - w*ivln2;
 		t1 = u + v;
@@ -193,10 +193,10 @@ float powf(float x, float y)
 		r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
 		r += s_l*(s_h+s);
 		s2 = s_h*s_h;
-		t_h = (float)3.0 + s2 + r;
+		t_h = 3.0f + s2 + r;
 		GET_FLOAT_WORD(is, t_h);
 		SET_FLOAT_WORD(t_h, is & 0xfffff000);
-		t_l = r - ((t_h - (float)3.0) - s2);
+		t_l = r - ((t_h - 3.0f) - s2);
 		/* u+v = s*(1+...) */
 		u = s_h*t_h;
 		v = s_l*t_h + t_l*s;
diff --git a/src/math/remainderf.c b/src/math/remainderf.c
index c17bb4f4..61c3c660 100644
--- a/src/math/remainderf.c
+++ b/src/math/remainderf.c
@@ -49,7 +49,7 @@ float remainderf(float x, float p)
 				x -= p;
 		}
 	} else {
-		p_half = (float)0.5*p;
+		p_half = 0.5f*p;
 		if (x > p_half) {
 			x -= p;
 			if (x >= p_half)
diff --git a/src/math/remainderl.c b/src/math/remainderl.c
index b99f9381..2a13c1d5 100644
--- a/src/math/remainderl.c
+++ b/src/math/remainderl.c
@@ -1,4 +1,5 @@
-#include "libm.h"
+#include <math.h>
+#include <float.h>
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double remainderl(long double x, long double y)
diff --git a/src/math/round.c b/src/math/round.c
index 21373847..bf0b453f 100644
--- a/src/math/round.c
+++ b/src/math/round.c
@@ -25,7 +25,7 @@
  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include "libm.h"
+#include <math.h>
 
 double round(double x)
 {
diff --git a/src/math/roundf.c b/src/math/roundf.c
index 3cfd8ae5..662a454b 100644
--- a/src/math/roundf.c
+++ b/src/math/roundf.c
@@ -25,7 +25,7 @@
  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include "libm.h"
+#include <math.h>
 
 float roundf(float x)
 {
diff --git a/src/math/roundl.c b/src/math/roundl.c
index ce56e8a9..99146f07 100644
--- a/src/math/roundl.c
+++ b/src/math/roundl.c
@@ -25,7 +25,9 @@
  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double roundl(long double x)
 {
diff --git a/src/math/scalb.c b/src/math/scalb.c
index 7706e9cb..a54bdf2b 100644
--- a/src/math/scalb.c
+++ b/src/math/scalb.c
@@ -15,7 +15,7 @@
  * should use scalbn() instead.
  */
 
-#include "libm.h"
+#include <math.h>
 
 double scalb(double x, double fn)
 {
diff --git a/src/math/scalbf.c b/src/math/scalbf.c
index 0cc091f1..94364974 100644
--- a/src/math/scalbf.c
+++ b/src/math/scalbf.c
@@ -13,19 +13,19 @@
  * ====================================================
  */
 
-#include "libm.h"
+#include <math.h>
 
 float scalbf(float x, float fn)
 {
 	if (isnan(x) || isnan(fn)) return x*fn;
 	if (!isfinite(fn)) {
-		if (fn > (float)0.0)
+		if (fn > 0.0f)
 			return x*fn;
 		else
 			return x/(-fn);
 	}
 	if (rintf(fn) != fn) return (fn-fn)/(fn-fn);
-	if ( fn > (float)65000.0) return scalbnf(x, 65000);
-	if (-fn > (float)65000.0) return scalbnf(x,-65000);
+	if ( fn > 65000.0f) return scalbnf(x, 65000);
+	if (-fn > 65000.0f) return scalbnf(x,-65000);
 	return scalbnf(x,(int)fn);
 }
diff --git a/src/math/scalbln.c b/src/math/scalbln.c
index 53854fda..e6f3f195 100644
--- a/src/math/scalbln.c
+++ b/src/math/scalbln.c
@@ -1,5 +1,5 @@
 #include <limits.h>
-#include "libm.h"
+#include <math.h>
 
 double scalbln(double x, long n)
 {
diff --git a/src/math/scalblnf.c b/src/math/scalblnf.c
index 61600f18..d8e8166b 100644
--- a/src/math/scalblnf.c
+++ b/src/math/scalblnf.c
@@ -1,5 +1,5 @@
 #include <limits.h>
-#include "libm.h"
+#include <math.h>
 
 float scalblnf(float x, long n)
 {
diff --git a/src/math/scalblnl.c b/src/math/scalblnl.c
index 82ebbed0..854c51c4 100644
--- a/src/math/scalblnl.c
+++ b/src/math/scalblnl.c
@@ -1,5 +1,6 @@
 #include <limits.h>
-#include "libm.h"
+#include <math.h>
+#include <float.h>
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double scalblnl(long double x, long n)
diff --git a/src/math/sincos.c b/src/math/sincos.c
new file mode 100644
index 00000000..442e285e
--- /dev/null
+++ b/src/math/sincos.c
@@ -0,0 +1,68 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_sin.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.
+ * ====================================================
+ */
+
+#include "libm.h"
+
+void sincos(double x, double *sin, double *cos)
+{
+	double y[2], s, c;
+	uint32_t n, ix;
+
+	GET_HIGH_WORD(ix, x);
+	ix &= 0x7fffffff;
+
+	/* |x| ~< pi/4 */
+	if (ix <= 0x3fe921fb) {
+		/* if |x| < 2**-27 * sqrt(2) */
+		if (ix < 0x3e46a09e) {
+			/* raise inexact if x != 0 */
+			if ((int)x == 0) {
+				*sin = x;
+				*cos = 1.0;
+			}
+			return;
+		}
+		*sin = __sin(x, 0.0, 0);
+		*cos = __cos(x, 0.0);
+		return;
+	}
+
+	/* sincos(Inf or NaN) is NaN */
+	if (ix >= 0x7ff00000) {
+		*sin = *cos = x - x;
+		return;
+	}
+
+	/* argument reduction needed */
+	n = __rem_pio2(x, y);
+	s = __sin(y[0], y[1], 1);
+	c = __cos(y[0], y[1]);
+	switch (n&3) {
+	case 0:
+		*sin = s;
+		*cos = c;
+		break;
+	case 1:
+		*sin = c;
+		*cos = -s;
+		break;
+	case 2:
+		*sin = -s;
+		*cos = -c;
+		break;
+	case 3:
+	default:
+		*sin = -c;
+		*cos = s;
+		break;
+	}
+}
diff --git a/src/math/sincosf.c b/src/math/sincosf.c
new file mode 100644
index 00000000..fede5f23
--- /dev/null
+++ b/src/math/sincosf.c
@@ -0,0 +1,113 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_sinf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * 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"
+
+/* Small multiples of pi/2 rounded to double precision. */
+static const double
+s1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
+s2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
+s3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
+s4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */
+
+void sincosf(float x, float *sin, float *cos)
+{
+	double y, s, c;
+	uint32_t n, hx, ix;
+
+	GET_FLOAT_WORD(hx, x);
+	ix = hx & 0x7fffffff;
+
+	/* |x| ~<= pi/4 */
+	if (ix <= 0x3f490fda) {
+		/* |x| < 2**-12 */
+		if (ix < 0x39800000) {
+			/* raise inexact if x != 0 */
+			if((int)x == 0) {
+				*sin = x;
+				*cos = 1.0f;
+			}
+			return;
+		}
+		*sin = __sindf(x);
+		*cos = __cosdf(x);
+		return;
+	}
+
+	/* |x| ~<= 5*pi/4 */
+	if (ix <= 0x407b53d1) {  
+		if (ix <= 0x4016cbe3) {  /* |x| ~<= 3pi/4 */
+			if (hx < 0x80000000) {
+				*sin = __cosdf(x - s1pio2);
+				*cos = __sindf(s1pio2 - x);
+			} else {
+				*sin = -__cosdf(x + s1pio2);
+				*cos = __sindf(x + s1pio2);
+			}
+			return;
+		}
+		*sin = __sindf(hx < 0x80000000 ? s2pio2 - x : -s2pio2 - x);
+		*cos = -__cosdf(hx < 0x80000000 ? x - s2pio2 : x + s2pio2);
+		return;
+	}
+
+	/* |x| ~<= 9*pi/4 */
+	if (ix <= 0x40e231d5) {
+		if (ix <= 0x40afeddf) {  /* |x| ~<= 7*pi/4 */
+			if (hx < 0x80000000) {
+				*sin = -__cosdf(x - s3pio2);
+				*cos = __sindf(x - s3pio2);
+			} else {
+				*sin = __cosdf(x + s3pio2);
+				*cos = __sindf(-s3pio2 - x);
+			}
+			return;
+		}
+		*sin = __sindf(hx < 0x80000000 ? x - s4pio2 : x + s4pio2);
+		*cos = __cosdf(hx < 0x80000000 ? x - s4pio2 : x + s4pio2);
+		return;
+	}
+
+	/* sin(Inf or NaN) is NaN */
+	if (ix >= 0x7f800000) {
+		*sin = *cos = x - x;
+		return;
+	}
+
+	/* general argument reduction needed */
+	n = __rem_pio2f(x, &y);
+	s = __sindf(y);
+	c = __cosdf(y);
+	switch (n&3) {
+	case 0:
+		*sin = s;
+		*cos = c;
+		break;
+	case 1:
+		*sin = c;
+		*cos = -s;
+		break;
+	case 2:
+		*sin = -s;
+		*cos = -c;
+		break;
+	case 3:
+	default:
+		*sin = -c;
+		*cos = s;
+		break;
+	}
+}
diff --git a/src/math/sincosl.c b/src/math/sincosl.c
new file mode 100644
index 00000000..378dc979
--- /dev/null
+++ b/src/math/sincosl.c
@@ -0,0 +1,66 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+void sincosl(long double x, long double *sin, long double *cos)
+{
+	double s, c;
+	sincos(x, &s, &c);
+	*sin = s;
+	*cos = c;
+}
+#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+#include "__rem_pio2l.h"
+
+void sincosl(long double x, long double *sin, long double *cos)
+{
+	union IEEEl2bits u;
+	int n;
+	long double y[2], s, c;
+
+	u.e = x;
+	u.bits.sign = 0;
+
+	/* x = +-0 or subnormal */
+	if (!u.bits.exp) {
+		*sin = x;
+		*cos = 1.0;
+		return;
+	}
+
+	/* x = nan or inf */
+	if (u.bits.exp == 0x7fff) {
+		*sin = *cos = x - x;
+		return;
+	}
+
+	/* |x| < pi/4 */
+	if (u.e < M_PI_4) {
+		*sin = __sinl(x, 0, 0);
+		*cos = __cosl(x, 0);
+		return;
+	}
+
+	n = __rem_pio2l(x, y);
+	s = __sinl(y[0], y[1], 1);
+	c = __cosl(y[0], y[1]);
+	switch (n & 3) {
+	case 0:
+		*sin = s;
+		*cos = c;
+		break;
+	case 1:
+		*sin = c;
+		*cos = -s;
+		break;
+	case 2:
+		*sin = -s;
+		*cos = -c;
+		break;
+	case 3:
+	default:
+		*sin = -c;
+		*cos = s;
+		break;
+	}
+}
+#endif
diff --git a/src/math/sinhf.c b/src/math/sinhf.c
index 056b5f86..fd11b849 100644
--- a/src/math/sinhf.c
+++ b/src/math/sinhf.c
@@ -40,7 +40,7 @@ float sinhf(float x)
 				return x;
 		t = expm1f(fabsf(x));
 		if (ix < 0x3f800000)
-			return h*((float)2.0*t - t*t/(t+one));
+			return h*(2.0f*t - t*t/(t+one));
 		return h*(t + t/(t+one));
 	}
 
diff --git a/src/math/truncf.c b/src/math/truncf.c
index 209586e1..0afcdfbf 100644
--- a/src/math/truncf.c
+++ b/src/math/truncf.c
@@ -20,7 +20,7 @@
 
 #include "libm.h"
 
-static const float huge = 1.0e30F;
+static const float huge = 1.0e30f;
 
 float truncf(float x)
 {
@@ -32,14 +32,14 @@ float truncf(float x)
 	if (j0 < 23) {
 		if (j0 < 0) {  /* |x|<1, return 0*sign(x) */
 			/* raise inexact if x != 0 */
-			if (huge+x > 0.0F)
+			if (huge+x > 0.0f)
 				i0 &= 0x80000000;
 		} else {
 			i = 0x007fffff>>j0;
 			if ((i0&i) == 0)
 				return x; /* x is integral */
 			/* raise inexact */
-			if (huge+x > 0.0F)
+			if (huge+x > 0.0f)
 				i0 &= ~i;
 		}
 	} else {