Merge remote branch 'nsz/master'

68 files changed, 440 insertions, 237 deletions

diff --git a/COPYRIGHT b/COPYRIGHT
index eec0066d..730e04a4 100644
--- a/COPYRIGHT
+++ b/COPYRIGHT
@@ -9,11 +9,15 @@ under the terms of the GNU LGPL version 2.1 or later. The included
 version was heavily modified in Spring 2006 by Rich Felker in the
 interests of size, simplicity, and namespace cleanliness.
 
-Most of the math library code (src/math/*) is Copyright © 1993 Sun
-Microsystems, Inc. Some files are Copyright © 2003 Steven G. Kargl and
-labelled as such. All have been licensed under extremely permissive
-terms which are compatible with the GNU LGPL. See the comments in the
-individual files for details.
+Most of the math library code (src/math/* and src/complex/*) is
+Copyright © 1993,2004 Sun Microsystems or
+Copyright © 2003-2011 David Schultz or
+Copyright © 2003-2009 Steven G. Kargl or
+Copyright © 2003-2009 Bruce D. Evans or
+Copyright © 2008 Stephen L. Moshier
+and labelled as such. All have been licensed under extremely
+permissive terms which are compatible with the GNU LGPL. See the
+comments in the individual files for details.
 
 The implementation of DES for crypt (src/misc/crypt.c) is Copyright ©
 1994 David Burren. It is licensed under a BSD license compatible with
diff --git a/include/math.h b/include/math.h
index f320b8e9..9c06c967 100644
--- a/include/math.h
+++ b/include/math.h
@@ -349,41 +349,43 @@ double      gamma(double);
 float       gammaf(float);
 long double gammal(long double);
 
+double      j0(double);
+double      j1(double);
+double      jn(int, double);
+
+double      y0(double);
+double      y1(double);
+double      yn(int, double);
+#endif
+
+#ifdef _GNU_SOURCE
+double      scalb(double, double);
+float       scalbf(float, float);
+long double scalbl(long double, long double);
+
+void        sincos(double, double*, double*);
+void        sincosf(float, float*, float*);
+void        sincosl(long double, long double*, long double*);
+
 double      lgamma_r(double, int*);
 float       lgammaf_r(float, int*);
 long double lgammal_r(long double, int*);
 
-double      j0(double);
 float       j0f(float);
 long double j0l(long double);
-
-double      j1(double);
 float       j1f(float);
 long double j1l(long double);
-
-double      jn(int, double);
 float       jnf(int, float);
 long double jnl(int, long double);
 
-double      y0(double);
 float       y0f(float);
 long double y0l(long double);
-
-double      y1(double);
 float       y1f(float);
 long double y1l(long double);
-
-double      yn(int, double);
 float       ynf(int, float);
 long double ynl(int, long double);
 #endif
 
-#ifdef _GNU_SOURCE
-double      scalb(double, double);
-float       scalbf(float, float);
-long double scalbl(long double, long double);
-#endif
-
 #ifdef __cplusplus
 }
 #endif
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 {