Merge remote branch 'nsz/master'

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

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

src/math/__expo2.c

@@ -1,51 +1,16 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */
-/*-
- * Copyright (c) 2011 David Schultz <[email protected]>
- * 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;
 }

src/math/__expo2f.c

@@ -1,51 +1,16 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/k_expf.c */
-/*-
- * Copyright (c) 2011 David Schultz <[email protected]>
- * 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;
 }

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);
 }

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);
 	}
 }

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));
 	}
 }

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);

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)));

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;

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;
 			}
 		}
 	}

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;

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;

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;

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;
 }

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;

src/math/fdim.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 double fdim(double x, double y)
 {

src/math/fdimf.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 float fdimf(float x, float y)
 {

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)

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;

src/math/fmax.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 double fmax(double x, double y)
 {

src/math/fmaxf.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 float fmaxf(float x, float y)
 {

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)

src/math/fmin.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 double fmin(double x, double y)
 {

src/math/fminf.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 float fminf(float x, float y)
 {

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)

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;

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;
 	}

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;
 }

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;
 }

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;

src/math/ldexp.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 double ldexp(double x, int n)
 {

src/math/ldexpf.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 float ldexpf(float x, int n)
 {

src/math/ldexpl.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 long double ldexpl(long double x, int n)
 {

src/math/lgamma.c

@@ -1,3 +1,4 @@
+#define _GNU_SOURCE
 #include "libm.h"
 
 double lgamma(double x)

src/math/lgammaf.c

@@ -1,3 +1,4 @@
+#define _GNU_SOURCE
 #include "libm.h"
 
 float lgammaf(float x)

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;
 		}

src/math/lgammal.c

@@ -85,6 +85,7 @@
  *
  */
 
+#define _GNU_SOURCE
 #include "libm.h"
 
 long double lgammal(long double x)

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)
 {

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)
 {

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

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)

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);
 
 	/*

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);

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)
 {

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)
 {

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)

src/math/nearbyintf.c

@@ -1,5 +1,5 @@
 #include <fenv.h>
-#include "libm.h"
+#include <math.h>
 
 float nearbyintf(float x) {
 	fenv_t e;

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)
 {

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)
 {

src/math/nexttowardl.c

@@ -1,4 +1,4 @@
-#include "libm.h"
+#include <math.h>
 
 long double nexttowardl(long double x, long double y)
 {

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

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;

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)

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)

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)
 {

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)
 {

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)
 {

src/math/scalb.c

@@ -15,7 +15,7 @@
  * should use scalbn() instead.
  */
 
-#include "libm.h"
+#include <math.h>
 
 double scalb(double x, double fn)
 {

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);
 }

src/math/scalbln.c

@@ -1,5 +1,5 @@
 #include <limits.h>
-#include "libm.h"
+#include <math.h>
 
 double scalbln(double x, long n)
 {

src/math/scalblnf.c

@@ -1,5 +1,5 @@
 #include <limits.h>
-#include "libm.h"
+#include <math.h>
 
 float scalblnf(float x, long n)
 {

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)

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;
+	}
+}

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, [email protected].
+ * 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;
+	}
+}

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

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));
 	}
 

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 {