diff options
Diffstat (limited to 'src')
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 { |