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
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
-/* 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;
}
-/* 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;
}
{
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);
}
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 */
}
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);
q = one+z*qS1;
r = p/q;
w = r*s+c;
- return (float)2.0*(df+w);
+ return 2.0f*(df+w);
}
}
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));
}
}
}
/* 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);
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)));
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) {
/* 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;
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;
}
}
}
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;
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;
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;
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;
}
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;
}
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;
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;
}
}
if (xsb != 0) { /* x < -27*ln2 */
/* raise inexact */
- if (x+tiny < (float)0.0)
+ if (x+tiny < 0.0f)
return tiny-one; /* return -1 */
}
}
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;
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;
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;
-#include "libm.h"
+#include <math.h>
double fdim(double x, double y)
{
-#include "libm.h"
+#include <math.h>
float fdimf(float x, float y)
{
-#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)
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)
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;
-#include "libm.h"
+#include <math.h>
double fmax(double x, double y)
{
-#include "libm.h"
+#include <math.h>
float fmaxf(float x, float y)
{
-#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)
-#include "libm.h"
+#include <math.h>
double fmin(double x, double y)
{
-#include "libm.h"
+#include <math.h>
float fminf(float x, float y)
{
-#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)
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;
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;
}
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);
}
}
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;
}
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] = {
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;
}
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! */
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;
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;
a = temp;
di -= two;
/* scale b to avoid spurious overflow */
- if (b > (float)1e10) {
+ if (b > 1e10f) {
a /= b;
t /= b;
b = one;
-#include "libm.h"
+#include <math.h>
double ldexp(double x, int n)
{
-#include "libm.h"
+#include <math.h>
float ldexpf(float x, int n)
{
-#include "libm.h"
+#include <math.h>
long double ldexpl(long double x, int n)
{
+#define _GNU_SOURCE
#include "libm.h"
double lgamma(double x)
+#define _GNU_SOURCE
#include "libm.h"
float lgammaf(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 */
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;
}
} 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;
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;
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;
}
*
*/
+#define _GNU_SOURCE
#include "libm.h"
long double lgammal(long double x)
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long long llrintl(long double x)
{
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long long llroundl(long double 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
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 */
}
/* 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;
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;
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)
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)
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);
/*
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)
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);
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);
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long lrintl(long double x)
{
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long lroundl(long double x)
{
#include <fenv.h>
-#include "libm.h"
+#include <math.h>
/*
rint may raise inexact (and it should not alter the fenv otherwise)
#include <fenv.h>
-#include "libm.h"
+#include <math.h>
float nearbyintf(float x) {
fenv_t e;
-#include "libm.h"
+#include <math.h>
+#include <float.h>
+
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double nearbyintl(long double x)
{
*/
#include "libm.h"
+
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
double nexttoward(double x, long double y)
{
-#include "libm.h"
+#include <math.h>
long double nexttowardl(long double x, long 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
/* 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
/* 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;
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;
x -= p;
}
} else {
- p_half = (float)0.5*p;
+ p_half = 0.5f*p;
if (x > p_half) {
x -= p;
if (x >= p_half)
-#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)
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
-#include "libm.h"
+#include <math.h>
double round(double x)
{
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
-#include "libm.h"
+#include <math.h>
float roundf(float x)
{
* 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)
{
* should use scalbn() instead.
*/
-#include "libm.h"
+#include <math.h>
double scalb(double x, double fn)
{
* ====================================================
*/
-#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);
}
#include <limits.h>
-#include "libm.h"
+#include <math.h>
double scalbln(double x, long n)
{
#include <limits.h>
-#include "libm.h"
+#include <math.h>
float scalblnf(float x, long n)
{
#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)
--- /dev/null
+/* 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;
+ }
+}
--- /dev/null
+/* 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;
+ }
+}
--- /dev/null
+#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
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));
}
#include "libm.h"
-static const float huge = 1.0e30F;
+static const float huge = 1.0e30f;
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 {
projects / musl / commitdiff