RepoMirrors
/
musl
mirror of git://git.musl-libc.org/musl
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge remote branch 'nsz/master'

This commit is contained in:
Rich Felker 2012-03-19 22:07:43 -04:00
parent acb744921b 0cbb654791
commit 97721a5508
95 changed files with 655 additions and 806 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/internal/libm.h
View File

@ -135,6 +135,7 @@ float __tandf(double,int);
float __expo2f(float);
float complex __ldexp_cexpf(float complex,int);
int __rem_pio2l(long double, long double *);
long double __sinl(long double, long double, int);
long double __cosl(long double, long double);
long double __tanl(long double, long double, int);

5
src/math/__cos.c
View File

@ -51,7 +51,6 @@
#include "libm.h"
static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
@ -67,6 +66,6 @@ double __cos(double x, double y)
w = z*z;
r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
hz = 0.5*z;
w = one-hz;
return w + (((one-w)-hz) + (z*r-x*y));
w = 1.0-hz;
return w + (((1.0-w)-hz) + (z*r-x*y));
}

3
src/math/__cosdf.c
View File

@ -18,7 +18,6 @@
/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
static const double
one = 1.0,
C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */
C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
@ -32,5 +31,5 @@ float __cosdf(double x)
z = x*x;
w = z*z;
r = C2+z*C3;
return ((one+z*C0) + w*C1) + (w*z)*r;
return ((1.0+z*C0) + w*C1) + (w*z)*r;
}

6
src/math/__cosl.c
View File

@ -41,8 +41,6 @@
* almost for free from the complications needed to search for the best
* higher coefficients.
*/
static const double one = 1.0;
static const long double
C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */
@ -61,7 +59,7 @@ long double __cosl(long double x, long double y)
z = x*x;
r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))));
hz = 0.5*z;
w = one-hz;
return w + (((one-w)-hz) + (z*r-x*y));
w = 1.0-hz;
return w + (((1.0-w)-hz) + (z*r-x*y));
}
#endif

3
src/math/__rem_pio2.c
View File

@ -29,7 +29,6 @@
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const double
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
@ -163,7 +162,7 @@ medium:
}
tx[2] = z;
nx = 3;
while (tx[nx-1] == zero) nx--; /* skip zero term */
while (tx[nx-1] == 0.0) nx--; /* skip zero term */
n = __rem_pio2_large(tx,ty,e0,nx,1);
if (hx < 0) {
y[0] = -ty[0];

12
src/math/__rem_pio2_large.c
View File

@ -271,8 +271,6 @@ static const double PIo2[] = {
};
static const double
zero = 0.0,
one = 1.0,
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
@ -293,7 +291,7 @@ int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec)
/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
j = jv-jx; m = jx+jk;
for (i=0; i<=m; i++,j++)
f[i] = j<0 ? zero : (double)ipio2[j];
f[i] = j<0 ? 0.0 : (double)ipio2[j];
/* compute q[0],q[1],...q[jk] */
for (i=0; i<=jk; i++) {
@ -346,14 +344,14 @@ recompute:
}
}
if (ih == 2) {
z = one - z;
z = 1.0 - z;
if (carry != 0)
z -= scalbn(one,q0);
z -= scalbn(1.0,q0);
}
}
/* check if recomputation is needed */
if (z == zero) {
if (z == 0.0) {
j = 0;
for (i=jz-1; i>=jk; i--) j |= iq[i];
if (j == 0) { /* need recomputation */
@ -391,7 +389,7 @@ recompute:
}
/* convert integer "bit" chunk to floating-point value */
fw = scalbn(one,q0);
fw = scalbn(1.0,q0);
for (i=jz; i>=0; i--) {
q[i] = fw*(double)iq[i];
fw *= twon24;

5
src/math/__rem_pio2l.h → src/math/__rem_pio2l.c
View File

@ -32,7 +32,6 @@
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const double
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */
pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */
@ -44,7 +43,7 @@ pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */
pio2_2t = 6.36831716351095013979e-25L, /* 0xc51701b839a25205.0p-144 */
pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */
static inline int __rem_pio2l(long double x, long double *y)
int __rem_pio2l(long double x, long double *y)
{
union IEEEl2bits u,u1;
long double z,w,t,r,fn;
@ -114,7 +113,7 @@ static inline int __rem_pio2l(long double x, long double *y)
}
tx[2] = z;
nx = 3;
while (tx[nx-1] == zero)
while (tx[nx-1] == 0.0)
nx--; /* skip zero term */
n = __rem_pio2_large(tx,ty,e0,nx,2);
r = (long double)ty[0] + ty[1];

3
src/math/__sin.c
View File

@ -42,7 +42,6 @@
#include "libm.h"
static const double
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
@ -61,5 +60,5 @@ double __sin(double x, double y, int iy)
if (iy == 0)
return x + v*(S1 + z*r);
else
return x - ((z*(half*y - v*r) - y) - v*S1);
return x - ((z*(0.5*y - v*r) - y) - v*S1);
}

4
src/math/__sinl.c
View File

@ -24,8 +24,6 @@
* See __cosl.c for more details about the polynomial.
*/
static const double half = 0.5;
static const long double
S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */
@ -47,6 +45,6 @@ long double __sinl(long double x, long double y, int iy)
r = S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8)))));
if (iy == 0)
return x+v*(S1+z*r);
return x-((z*(half*y-v*r)-y)-v*S1);
return x-((z*(0.5*y-v*r)-y)-v*S1);
}
#endif

10
src/math/__tan.c
View File

@ -59,13 +59,9 @@ static const double T[] = {
7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */
-1.85586374855275456654e-05, /* BEF375CB, DB605373 */
2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */
/* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */
/* pio4 */ 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
/* pio4lo */ 3.06161699786838301793e-17 /* 3C81A626, 33145C07 */
};
#define one T[13]
#define pio4 T[14]
#define pio4lo T[15]
},
pio4 = 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
pio4lo = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */
double __tan(double x, double y, int iy)
{

3
src/math/__tanl.c
View File

@ -51,8 +51,7 @@ long double __tanl(long double x, long double y, int iy) {
int i;
iy = iy == 1 ? -1 : 1; /* XXX recover original interface */
// FIXME: this is wrong, use copysign, signbit or union bithack
osign = x >= 0 ? 1.0 : -1.0; /* XXX slow, probably wrong for -0 */
osign = copysignl(1.0, x);
if (fabsl(x) >= 0.67434) {
if (x < 0) {
x = -x;

13
src/math/acos.c
View File

@ -36,8 +36,7 @@
#include "libm.h"
static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
static const volatile double
pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
@ -75,25 +74,25 @@ double acos(double x)
return pio2_hi + pio2_lo;
z = x*x;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q;
return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx < 0) { /* x < -0.5 */
z = (one+x)*0.5;
z = (1.0+x)*0.5;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
s = sqrt(z);
r = p/q;
w = r*s-pio2_lo;
return pi - 2.0*(s+w);
} else { /* x > 0.5 */
z = (one-x)*0.5;
z = (1.0-x)*0.5;
s = sqrt(z);
df = s;
SET_LOW_WORD(df,0);
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q;
w = r*s+c;
return 2.0*(df+w);

13
src/math/acosf.c
View File

@ -16,8 +16,7 @@
#include "libm.h"
static const float
one = 1.0000000000e+00, /* 0x3F800000 */
pi = 3.1415925026e+00, /* 0x40490fda */
pi = 3.1415925026e+00, /* 0x40490fda */
pio2_hi = 1.5707962513e+00; /* 0x3fc90fda */
static const volatile float
pio2_lo = 7.5497894159e-08; /* 0x33a22168 */
@ -46,13 +45,13 @@ float acosf(float x)
return pio2_hi + pio2_lo;
z = x*x;
p = z*(pS0+z*(pS1+z*pS2));
q = one+z*qS1;
q = 1.0f+z*qS1;
r = p/q;
return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx < 0) { /* x < -0.5 */
z = (one+x)*0.5f;
z = (1.0f+x)*0.5f;
p = z*(pS0+z*(pS1+z*pS2));
q = one+z*qS1;
q = 1.0f+z*qS1;
s = sqrtf(z);
r = p/q;
w = r*s-pio2_lo;
@ -60,14 +59,14 @@ float acosf(float x)
} else { /* x > 0.5 */
int32_t idf;
z = (one-x)*0.5f;
z = (1.0f-x)*0.5f;
s = sqrtf(z);
df = s;
GET_FLOAT_WORD(idf,df);
SET_FLOAT_WORD(df,idf&0xfffff000);
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*pS2));
q = one+z*qS1;
q = 1.0f+z*qS1;
r = p/q;
w = r*s+c;
return 2.0f*(df+w);

5
src/math/acosh.c
View File

@ -27,7 +27,6 @@
#include "libm.h"
static const double
one = 1.0,
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
double acosh(double x)
@ -47,9 +46,9 @@ double acosh(double x)
return 0.0; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t = x*x;
return log(2.0*x - one/(x+sqrt(t-one)));
return log(2.0*x - 1.0/(x+sqrt(t-1.0)));
} else { /* 1 < x < 2 */
t = x-one;
t = x-1.0;
return log1p(t + sqrt(2.0*t+t*t));
}
}

7
src/math/acoshf.c
View File

@ -16,7 +16,6 @@
#include "libm.h"
static const float
one = 1.0,
ln2 = 6.9314718246e-01; /* 0x3f317218 */
float acoshf(float x)
@ -32,12 +31,12 @@ float acoshf(float x)
return x + x;
return logf(x) + ln2; /* acosh(huge)=log(2x) */
} else if (hx == 0x3f800000) {
return 0.0; /* acosh(1) = 0 */
return 0.0f; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t = x*x;
return logf(2.0f*x - one/(x+sqrtf(t-one)));
return logf(2.0f*x - 1.0f/(x+sqrtf(t-1.0f)));
} else { /* 1 < x < 2 */
t = x-one;
t = x-1.0f;
return log1pf(t + sqrtf(2.0f*t+t*t));
}
}

5
src/math/acoshl.c
View File

@ -32,7 +32,6 @@ long double acoshl(long double x)
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
one = 1.0,
ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
long double acoshl(long double x)
@ -51,10 +50,10 @@ long double acoshl(long double x)
return 0.0; /* acosh(1) = 0 */
} else if (se > 0x4000) { /* x > 2 */
t = x*x;
return logl(2.0*x - one/(x + sqrtl(t - one)));
return logl(2.0*x - 1.0/(x + sqrtl(t - 1.0)));
}
/* 1 < x <= 2 */
t = x - one;
t = x - 1.0;
return log1pl(t + sqrtl(2.0*t + t*t));
}
#endif

5
src/math/acosl.c
View File

@ -25,7 +25,6 @@ long double acosl(long double x)
#include "__invtrigl.h"
static const long double
one = 1.00000000000000000000e+00,
pi = 3.14159265358979323846264338327950280e+00L;
long double acosl(long double x)
@ -55,7 +54,7 @@ long double acosl(long double x)
r = p / q;
return pio2_hi - (x - (pio2_lo - x * r));
} else if (expsign < 0) { /* x < -0.5 */
z = (one + x) * 0.5;
z = (1.0 + x) * 0.5;
p = P(z);
q = Q(z);
s = sqrtl(z);
@ -63,7 +62,7 @@ long double acosl(long double x)
w = r * s - pio2_lo;
return pi - 2.0 * (s + w);
} else { /* x > 0.5 */
z = (one - x) * 0.5;
z = (1.0 - x) * 0.5;
s = sqrtl(z);
u.e = s;
u.bits.manl = 0;

9
src/math/asin.c
View File

@ -42,7 +42,6 @@
#include "libm.h"
static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
huge = 1.000e+300,
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
@ -76,20 +75,20 @@ double asin(double x)
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3fe00000) { /* |x|<0.5 */
if (ix < 0x3e500000) { /* if |x| < 2**-26 */
if (huge+x > one)
if (huge+x > 1.0)
return x; /* return x with inexact if x!=0*/
}
t = x*x;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
w = p/q;
return x + x*w;
}
/* 1 > |x| >= 0.5 */
w = one - fabs(x);
w = 1.0 - fabs(x);
t = w*0.5;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
s = sqrt(t);
if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */
w = p/q;

9
src/math/asinf.c
View File

@ -16,7 +16,6 @@
#include "libm.h"
static const float
one = 1.0000000000e+00, /* 0x3F800000 */
huge = 1.000e+30,
/* coefficients for R(x^2) */
pS0 = 1.6666586697e-01,
@ -41,20 +40,20 @@ float asinf(float x)
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3f000000) { /* |x|<0.5 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
if (huge+x > one)
if (huge+x > 1.0f)
return x; /* return x with inexact if x!=0 */
}
t = x*x;
p = t*(pS0+t*(pS1+t*pS2));
q = one+t*qS1;
q = 1.0f+t*qS1;
w = p/q;
return x + x*w;
}
/* 1 > |x| >= 0.5 */
w = one - fabsf(x);
w = 1.0f - fabsf(x);
t = w*0.5f;
p = t*(pS0+t*(pS1+t*pS2));
q = one+t*qS1;
q = 1.0f+t*qS1;
s = sqrt(t);
w = p/q;
t = pio2-2.0*(s+s*w);

7
src/math/asinh.c
View File

@ -23,7 +23,6 @@
#include "libm.h"
static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
huge= 1.00000000000000000000e+300;
@ -38,17 +37,17 @@ double asinh(double x)
return x+x;
if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* return x inexact except 0 */
if (huge+x > one)
if (huge+x > 1.0)
return x;
}
if (ix > 0x41b00000) { /* |x| > 2**28 */
w = log(fabs(x)) + ln2;
} else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabs(x);
w = log(2.0*t + one/(sqrt(x*x+one)+t));
w = log(2.0*t + 1.0/(sqrt(x*x+1.0)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1p(fabs(x) + t/(one+sqrt(one+t)));
w =log1p(fabs(x) + t/(1.0+sqrt(1.0+t)));
}
if (hx > 0)
return w;

7
src/math/asinhf.c
View File

@ -16,7 +16,6 @@
#include "libm.h"
static const float
one = 1.0000000000e+00, /* 0x3F800000 */
ln2 = 6.9314718246e-01, /* 0x3f317218 */
huge= 1.0000000000e+30;
@ -31,17 +30,17 @@ float asinhf(float x)
return x+x;
if (ix < 0x31800000) { /* |x| < 2**-28 */
/* return x inexact except 0 */
if (huge+x > one)
if (huge+x > 1.0f)
return x;
}
if (ix > 0x4d800000) { /* |x| > 2**28 */
w = logf(fabsf(x)) + ln2;
} else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabsf(x);
w = logf(2.0f*t + one/(sqrtf(x*x+one)+t));
w = logf(2.0f*t + 1.0f/(sqrtf(x*x+1.0f)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1pf(fabsf(x) + t/(one+sqrtf(one+t)));
w =log1pf(fabsf(x) + t/(1.0f+sqrtf(1.0f+t)));
}
if (hx > 0)
return w;

7
src/math/asinhl.c
View File

@ -29,7 +29,6 @@ long double asinhl(long double x)
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
one = 1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */
ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
huge = 1.000000000000000000e+4900L;
@ -44,17 +43,17 @@ long double asinhl(long double x)
return x + x; /* x is inf or NaN */
if (ix < 0x3fde) { /* |x| < 2**-34 */
/* return x, raise inexact if x != 0 */
if (huge+x > one)
if (huge+x > 1.0)
return x;
}
if (ix > 0x4020) { /* |x| > 2**34 */
w = logl(fabsl(x)) + ln2;
} else if (ix > 0x4000) { /* 2**34 > |x| > 2.0 */
t = fabsl(x);
w = logl(2.0*t + one/(sqrtl(x*x + one) + t));
w = logl(2.0*t + 1.0/(sqrtl(x*x + 1.0) + t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1pl(fabsl(x) + t/(one + sqrtl(one + t)));
w =log1pl(fabsl(x) + t/(1.0 + sqrtl(1.0 + t)));
}
if (hx & 0x8000)
return -w;

8
src/math/asinl.c
View File

@ -23,9 +23,7 @@ long double asinl(long double x)
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
static const long double
one = 1.00000000000000000000e+00,
huge = 1.000e+300;
static const long double huge = 1.000e+300;
long double asinl(long double x)
{
@ -45,7 +43,7 @@ long double asinl(long double x)
} else if (expt < BIAS-1) { /* |x|<0.5 */
if (expt < ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */
/* return x with inexact if x!=0 */
if (huge+x > one)
if (huge+x > 1.0)
return x;
}
t = x*x;
@ -55,7 +53,7 @@ long double asinl(long double x)
return x + x*w;
}
/* 1 > |x| >= 0.5 */
w = one - fabsl(x);
w = 1.0 - fabsl(x);
t = w*0.5;
p = P(t);
q = Q(t);

12
src/math/atan.c
View File

@ -60,9 +60,7 @@ static const double aT[] = {
1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
};
static const double
one = 1.0,
huge = 1.0e300;
static const double huge = 1.0e300;
double atan(double x)
{
@ -86,7 +84,7 @@ double atan(double x)
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e400000) { /* |x| < 2^-27 */
/* raise inexact */
if (huge+x > one)
if (huge+x > 1.0)
return x;
}
id = -1;
@ -95,15 +93,15 @@ double atan(double x)
if (ix < 0x3ff30000) { /* |x| < 1.1875 */
if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */
id = 0;
x = (2.0*x-one)/(2.0+x);
x = (2.0*x-1.0)/(2.0+x);
} else { /* 11/16 <= |x| < 19/16 */
id = 1;
x = (x-one)/(x+one);
x = (x-1.0)/(x+1.0);
}
} else {
if (ix < 0x40038000) { /* |x| < 2.4375 */
id = 2;
x = (x-1.5)/(one+1.5*x);
x = (x-1.5)/(1.0+1.5*x);
} else { /* 2.4375 <= |x| < 2^66 */
id = 3;
x = -1.0/x;

5
src/math/atan2.c
View File

@ -42,7 +42,6 @@
static const volatile double
tiny = 1.0e-300;
static const double
zero = 0.0,
pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */
@ -89,8 +88,8 @@ double atan2(double y, double x)
}
} else {
switch(m) {
case 0: return zero; /* atan(+...,+INF) */
case 1: return -zero; /* atan(-...,+INF) */
case 0: return 0.0; /* atan(+...,+INF) */
case 1: return -0.0; /* atan(-...,+INF) */
case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */
}

5
src/math/atan2f.c
View File

@ -18,7 +18,6 @@
static const volatile float
tiny = 1.0e-30;
static const float
zero = 0.0,
pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
pi = 3.1415927410e+00; /* 0x40490fdb */
@ -63,8 +62,8 @@ float atan2f(float y, float x)
}
} else {
switch (m) {
case 0: return zero; /* atan(+...,+INF) */
case 1: return -zero; /* atan(-...,+INF) */
case 0: return 0.0f; /* atan(+...,+INF) */
case 1: return -0.0f; /* atan(-...,+INF) */
case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */
}

5
src/math/atan2l.c
View File

@ -27,7 +27,6 @@ long double atan2l(long double y, long double x)
static const volatile long double
tiny = 1.0e-300;
static const long double
zero = 0.0,
pi = 3.14159265358979323846264338327950280e+00L;
long double atan2l(long double y, long double x)
@ -75,8 +74,8 @@ long double atan2l(long double y, long double x)
}
} else {
switch(m) {
case 0: return zero; /* atan(+...,+INF) */
case 1: return -zero; /* atan(-...,+INF) */
case 0: return 0.0; /* atan(+...,+INF) */
case 1: return -0.0; /* atan(-...,+INF) */
case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */
}

12
src/math/atanf.c
View File

@ -38,9 +38,7 @@ static const float aT[] = {
6.1687607318e-02,
};
static const float
one = 1.0,
huge = 1.0e30;
static const float huge = 1.0e30;
float atanf(float x)
{
@ -60,7 +58,7 @@ float atanf(float x)
if (ix < 0x3ee00000) { /* |x| < 0.4375 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
/* raise inexact */
if(huge+x>one)
if(huge+x>1.0f)
return x;
}
id = -1;
@ -69,15 +67,15 @@ float atanf(float x)
if (ix < 0x3f980000) { /* |x| < 1.1875 */
if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */
id = 0;
x = (2.0f*x - one)/(2.0f + x);
x = (2.0f*x - 1.0f)/(2.0f + x);
} else { /* 11/16 <= |x| < 19/16 */
id = 1;
x = (x - one)/(x + one);
x = (x - 1.0f)/(x + 1.0f);
}
} else {
if (ix < 0x401c0000) { /* |x| < 2.4375 */
id = 2;
x = (x - 1.5f)/(one + 1.5f*x);
x = (x - 1.5f)/(1.0f + 1.5f*x);
} else { /* 2.4375 <= |x| < 2**26 */
id = 3;
x = -1.0f/x;

11
src/math/atanh.c
View File

@ -30,8 +30,7 @@
#include "libm.h"
static const double one = 1.0, huge = 1e300;
static const double zero = 0.0;
static const double huge = 1e300;
double atanh(double x)
{
@ -44,15 +43,15 @@ double atanh(double x)
if ((ix | ((lx|-lx)>>31)) > 0x3ff00000) /* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3ff00000)
return x/zero;
if (ix < 0x3e300000 && (huge+x) > zero) /* x < 2**-28 */
return x/0.0;
if (ix < 0x3e300000 && (huge+x) > 0.0) /* x < 2**-28 */
return x;
SET_HIGH_WORD(x, ix);
if (ix < 0x3fe00000) { /* x < 0.5 */
t = x+x;
t = 0.5*log1p(t + t*x/(one-x));
t = 0.5*log1p(t + t*x/(1.0-x));
} else
t = 0.5*log1p((x+x)/(one-x));
t = 0.5*log1p((x+x)/(1.0-x));
if (hx >= 0)
return t;
return -t;

11
src/math/atanhf.c
View File

@ -15,8 +15,7 @@
#include "libm.h"
static const float one = 1.0, huge = 1e30;
static const float zero = 0.0;
static const float huge = 1e30;
float atanhf(float x)
{
@ -28,15 +27,15 @@ float atanhf(float x)
if (ix > 0x3f800000) /* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3f800000)
return x/zero;
if (ix < 0x31800000 && huge+x > zero) /* x < 2**-28 */
return x/0.0f;
if (ix < 0x31800000 && huge+x > 0.0f) /* x < 2**-28 */
return x;
SET_FLOAT_WORD(x, ix);
if (ix < 0x3f000000) { /* x < 0.5 */
t = x+x;
t = 0.5f*log1pf(t + t*x/(one-x));
t = 0.5f*log1pf(t + t*x/(1.0f-x));
} else
t = 0.5f*log1pf((x+x)/(one-x));
t = 0.5f*log1pf((x+x)/(1.0f-x));
if (hx >= 0)
return t;
return -t;

10
src/math/atanhl.c
View File

@ -34,7 +34,7 @@ long double atanhl(long double x)
return atanh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double zero = 0.0, one = 1.0, huge = 1e4900L;
static const long double huge = 1e4900L;
long double atanhl(long double x)
{
@ -48,15 +48,15 @@ long double atanhl(long double x)
/* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3fff)
return x/zero;
if (ix < 0x3fe3 && huge+x > zero) /* x < 2**-28 */
return x/0.0;
if (ix < 0x3fe3 && huge+x > 0.0) /* x < 2**-28 */
return x;
SET_LDOUBLE_EXP(x, ix);
if (ix < 0x3ffe) { /* x < 0.5 */
t = x + x;
t = 0.5*log1pl(t + t*x/(one-x));
t = 0.5*log1pl(t + t*x/(1.0 - x));
} else
t = 0.5*log1pl((x + x)/(one - x));
t = 0.5*log1pl((x + x)/(1.0 - x));
if (se <= 0x7fff)
return t;
return -t;

12
src/math/atanl.c
View File

@ -23,9 +23,7 @@ long double atanl(long double x)
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
static const long double
one = 1.0,
huge = 1.0e300;
static const long double huge = 1.0e300;
long double atanl(long double x)
{
@ -53,7 +51,7 @@ long double atanl(long double x)
if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */
/* raise inexact */
if (huge+x > one)
if (huge+x > 1.0)
return x;
}
id = -1;
@ -62,15 +60,15 @@ long double atanl(long double x)
if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */
if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */
id = 0;
x = (2.0*x-one)/(2.0+x);
x = (2.0*x-1.0)/(2.0+x);
} else { /* 11/16 <= |x| < 19/16 */
id = 1;
x = (x-one)/(x+one);
x = (x-1.0)/(x+1.0);
}
} else {
if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */
id = 2;
x = (x-1.5)/(one+1.5*x);
x = (x-1.5)/(1.0+1.5*x);
} else { /* 2.4375 <= |x| < 2^ATAN_CONST */
id = 3;
x = -1.0/x;

6
src/math/cbrtl.c
View File

@ -118,11 +118,7 @@ long double cbrtl(long double x)
* Round it away from zero to 32 bits (32 so that t*t is exact, and
* away from zero for technical reasons).
*/
volatile double vd2 = 0x1.0p32;
volatile double vd1 = 0x1.0p-31;
#define vd ((long double)vd2 + vd1)
t = dt + vd - 0x1.0p32;
t = dt + (0x1.0p32L + 0x1.0p-32L) - 0x1.0p32;
#elif LDBL_MANT_DIG == 113
/*
* Round dt away from zero to 47 bits. Since we don't trust the 47,

12
src/math/cosh.c
View File

@ -32,7 +32,7 @@
#include "libm.h"
static const double one = 1.0, half = 0.5, huge = 1.0e300;
static const double huge = 1.0e300;
double cosh(double x)
{
@ -49,21 +49,21 @@ double cosh(double x)
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3fd62e43) {
t = expm1(fabs(x));
w = one+t;
w = 1.0+t;
if (ix < 0x3c800000)
return w; /* cosh(tiny) = 1 */
return one + (t*t)/(w+w);
return 1.0 + (t*t)/(w+w);
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x40360000) {
t = exp(fabs(x));
return half*t + half/t;
return 0.5*t + 0.5/t;
}
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862E42)
return half*exp(fabs(x));
return 0.5*exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */
if (ix <= 0x408633CE)

14
src/math/coshf.c
View File

@ -15,7 +15,7 @@
#include "libm.h"
static const float one = 1.0, half = 0.5, huge = 1.0e30;
static const float huge = 1.0e30;
float coshf(float x)
{
@ -32,21 +32,21 @@ float coshf(float x)
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3eb17218) {
t = expm1f(fabsf(x));
w = one+t;
w = 1.0f+t;
if (ix<0x39800000)
return one; /* cosh(tiny) = 1 */
return one + (t*t)/(w+w);
return 1.0f; /* cosh(tiny) = 1 */
return 1.0f + (t*t)/(w+w);
}
/* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x41100000) {
t = expf(fabsf(x));
return half*t + half/t;
return 0.5f*t + 0.5f/t;
}
/* |x| in [9, log(maxfloat)] return half*exp(|x|) */
/* |x| in [9, log(maxfloat)] return 0.5f*exp(|x|) */
if (ix < 0x42b17217)
return half*expf(fabsf(x));
return 0.5f*expf(fabsf(x));
/* |x| in [log(maxfloat), overflowthresold] */
if (ix <= 0x42b2d4fc)

16
src/math/coshl.c
View File

@ -38,7 +38,7 @@ long double coshl(long double x)
return cosh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double one = 1.0, half = 0.5, huge = 1.0e4900L;
static const long double huge = 1.0e4900L;
long double coshl(long double x)
{
@ -56,27 +56,27 @@ long double coshl(long double x)
/* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
if (ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
t = expm1l(fabsl(x));
w = one + t;
w = 1.0 + t;
if (ex < 0x3fbc) return w; /* cosh(tiny) = 1 */
return one+(t*t)/(w+w);
return 1.0+(t*t)/(w+w);
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
t = expl(fabsl(x));
return half*t + half/t;
return 0.5*t + 0.5/t;
}
/* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
/* |x| in [22, ln(maxdouble)] return 0.5*exp(|x|) */
if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
return half*expl(fabsl(x));
return 0.5*expl(fabsl(x));
/* |x| in [log(maxdouble), log(2*maxdouble)) */
if (ex == 0x400c && (mx < 0xb174ddc0u ||
(mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
{
w = expl(half*fabsl(x));
t = half*w;
w = expl(0.5*fabsl(x));
t = 0.5*w;
return t*w;
}

2
src/math/cosl.c
View File

@ -36,8 +36,6 @@ long double cosl(long double x) {
return cos(x);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__rem_pio2l.h"
long double cosl(long double x)
{
union IEEEl2bits z;

57
src/math/erf.c
View File

@ -107,9 +107,6 @@
static const double
tiny = 1e-300,
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
/* c = (float)0.84506291151 */
erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
/*
@ -190,7 +187,7 @@ double erf(double x)
if (ix >= 0x7ff00000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1;
return (double)(1-i) + one/x;
return (double)(1-i) + 1.0/x;
}
if (ix < 0x3feb0000) { /* |x|<0.84375 */
if (ix < 0x3e300000) { /* |x|<2**-28 */
@ -201,42 +198,42 @@ double erf(double x)
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one;
s = fabs(x)-1.0;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
Q = 1.0+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0)
return erx + P/Q;
return -erx - P/Q;
}
if (ix >= 0x40180000) { /* inf > |x| >= 6 */
if (hx >= 0)
return one-tiny;
return tiny-one;
return 1.0 - tiny;
return tiny - 1.0;
}
x = fabs(x);
s = one/(x*x);
s = 1.0/(x*x);
if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
SET_LOW_WORD(z,0);
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx >= 0)
return one-r/x;
return r/x-one;
return 1.0 - r/x;
return r/x - 1.0;
}
double erfc(double x)
@ -248,49 +245,49 @@ double erfc(double x)
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
return (double)(((uint32_t)hx>>31)<<1) + one/x;
return (double)(((uint32_t)hx>>31)<<1) + 1.0/x;
}
if (ix < 0x3feb0000) { /* |x| < 0.84375 */
if (ix < 0x3c700000) /* |x| < 2**-56 */
return one - x;
return 1.0 - x;
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
if (hx < 0x3fd00000) { /* x < 1/4 */
return one - (x+x*y);
return 1.0 - (x+x*y);
} else {
r = x*y;
r += x-half;
return half - r ;
r += x - 0.5;
return 0.5 - r ;
}
}
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one;
s = fabs(x)-1.0;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
Q = 1.0+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0) {
z = one-erx;
z = 1.0-erx;
return z - P/Q;
} else {
z = erx+P/Q;
return one+z;
return 1.0 + z;
}
}
if (ix < 0x403c0000) { /* |x| < 28 */
x = fabs(x);
s = one/(x*x);
s = 1.0/(x*x);
if (ix < 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40180000) /* x < -6 */
return two-tiny;
return 2.0 - tiny;
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
@ -298,9 +295,9 @@ double erfc(double x)
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx > 0)
return r/x;
return two-r/x;
return 2.0 - r/x;
}
if (hx > 0)
return tiny*tiny;
return two-tiny;
return 2.0 - tiny;
}

57
src/math/erff.c
View File

@ -17,9 +17,6 @@
static const float
tiny = 1e-30,
half = 5.0000000000e-01, /* 0x3F000000 */
one = 1.0000000000e+00, /* 0x3F800000 */
two = 2.0000000000e+00, /* 0x40000000 */
/* c = (subfloat)0.84506291151 */
erx = 8.4506291151e-01, /* 0x3f58560b */
/*
@ -100,7 +97,7 @@ float erff(float x)
if (ix >= 0x7f800000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1;
return (float)(1-i)+one/x;
return (float)(1-i)+1.0f/x;
}
if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x31800000) { /* |x| < 2**-28 */
@ -111,42 +108,42 @@ float erff(float x)
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-one;
s = fabsf(x)-1.0f;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
Q = 1.0f+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0)
return erx + P/Q;
return -erx - P/Q;
}
if (ix >= 0x40c00000) { /* inf > |x| >= 6 */
if (hx >= 0)
return one - tiny;
return tiny - one;
return 1.0f - tiny;
return tiny - 1.0f;
}
x = fabsf(x);
s = one/(x*x);
s = 1.0f/(x*x);
if (ix < 0x4036DB6E) { /* |x| < 1/0.35 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
GET_FLOAT_WORD(ix, x);
SET_FLOAT_WORD(z, ix&0xfffff000);
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx >= 0)
return one - r/x;
return r/x - one;
return 1.0f - r/x;
return r/x - 1.0f;
}
float erfcf(float x)
@ -158,50 +155,50 @@ float erfcf(float x)
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
return (float)(((uint32_t)hx>>31)<<1) + one/x;
return (float)(((uint32_t)hx>>31)<<1) + 1.0f/x;
}
if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x23800000) /* |x| < 2**-56 */
return one - x;
return 1.0f - x;
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
if (hx < 0x3e800000) { /* x<1/4 */
return one - (x+x*y);
return 1.0f - (x+x*y);
} else {
r = x*y;
r += (x-half);
return half - r ;
r += (x-0.5f);
return 0.5f - r ;
}
}
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-one;
s = fabsf(x)-1.0f;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
Q = 1.0f+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if(hx >= 0) {
z = one - erx;
z = 1.0f - erx;
return z - P/Q;
} else {
z = erx + P/Q;
return one + z;
return 1.0f + z;
}
}
if (ix < 0x41e00000) { /* |x| < 28 */
x = fabsf(x);
s = one/(x*x);
s = 1.0f/(x*x);
if (ix < 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40c00000) /* x < -6 */
return two-tiny;
return 2.0f-tiny;
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
GET_FLOAT_WORD(ix, x);
@ -209,9 +206,9 @@ float erfcf(float x)
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx > 0)
return r/x;
return two - r/x;
return 2.0f - r/x;
}
if (hx > 0)
return tiny*tiny;
return two - tiny;
return 2.0f - tiny;
}

54
src/math/erfl.c
View File

@ -108,9 +108,6 @@ long double erfl(long double x)
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
tiny = 1e-4931L,
half = 0.5L,
one = 1.0L,
two = 2.0L,
/* c = (float)0.84506291151 */
erx = 0.845062911510467529296875L,
@ -248,12 +245,12 @@ long double erfl(long double x)
int32_t ix, i;
uint32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
GET_LDOUBLE_WORDS(se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x7fff) { /* erf(nan)=nan */
i = ((se & 0xffff) >> 15) << 1;
return (long double)(1 - i) + one / x; /* erf(+-inf)=+-1 */
return (long double)(1 - i) + 1.0 / x; /* erf(+-inf)=+-1 */
}
ix = (ix << 16) | (i0 >> 16);
@ -272,7 +269,7 @@ long double erfl(long double x)
return x + x * y;
}
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
s = fabsl (x) - one;
s = fabsl(x) - 1.0;
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
@ -283,11 +280,11 @@ long double erfl(long double x)
}
if (ix >= 0x4001d555) { /* inf > |x| >= 6.6666259765625 */
if ((se & 0x8000) == 0)
return one - tiny;
return tiny - one;
return 1.0 - tiny;
return tiny - 1.0;
}
x = fabsl (x);
s = one / (x * x);
s = 1.0 / (x * x);
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
@ -305,8 +302,8 @@ long double erfl(long double x)
SET_LDOUBLE_WORDS(z, i, i0, i1);
r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0)
return one - r / x;
return r / x - one;
return 1.0 - r / x;
return r / x - 1.0;
}
long double erfcl(long double x)
@ -315,16 +312,16 @@ long double erfcl(long double x)
long double R, S, P, Q, s, y, z, r;
uint32_t se, i0, i1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
GET_LDOUBLE_WORDS(se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x7fff) { /* erfc(nan) = nan, erfc(+-inf) = 0,2 */
return (long double)(((se & 0xffff) >> 15) << 1) + one / x;
return (long double)(((se & 0xffff) >> 15) << 1) + 1.0 / x;
}
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffed800) { /* |x| < 0.84375 */
if (ix < 0x3fbe0000) /* |x| < 2**-65 */
return one - x;
return 1.0 - x;
z = x * x;
r = pp[0] + z * (pp[1] +
z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
@ -332,27 +329,27 @@ long double erfcl(long double x)
z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s;
if (ix < 0x3ffd8000) /* x < 1/4 */
return one - (x + x * y);
return 1.0 - (x + x * y);
r = x * y;
r += x - half;
return half - r;
r += x - 0.5L;
return 0.5L - r;
}
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
s = fabsl (x) - one;
s = fabsl(x) - 1.0;
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
if ((se & 0x8000) == 0) {
z = one - erx;
z = 1.0 - erx;
return z - P / Q;
}
z = erx + P / Q;
return one + z;
return 1.0 + z;
}
if (ix < 0x4005d600) { /* |x| < 107 */
x = fabsl (x);
s = one / (x * x);
x = fabsl(x);
s = 1.0 / (x * x);
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
@ -365,26 +362,25 @@ long double erfcl(long double x)
s * (sb[5] + s * (sb[6] + s))))));
} else { /* 107 > |x| >= 6.666 */
if (se & 0x8000)
return two - tiny;/* x < -6.666 */
return 2.0 - tiny;/* x < -6.666 */
R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
s * (rc[4] + s * rc[5]))));
S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
s * (sc[4] + s))));
}
z = x;
GET_LDOUBLE_WORDS (hx, i0, i1, z);
GET_LDOUBLE_WORDS(hx, i0, i1, z);
i1 = 0;
i0 &= 0xffffff00;
SET_LDOUBLE_WORDS (z, hx, i0, i1);
r = expl (-z * z - 0.5625) *
expl ((z - x) * (z + x) + R / S);
SET_LDOUBLE_WORDS(z, hx, i0, i1);
r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0)
return r / x;
return two - r / x;
return 2.0 - r / x;
}
if ((se & 0x8000) == 0)
return tiny * tiny;
return two - tiny;
return 2.0 - tiny;
}
#endif

9
src/math/exp.c
View File

@ -74,7 +74,6 @@
#include "libm.h"
static const double
one = 1.0,
halF[2] = {0.5,-0.5,},
huge = 1.0e+300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
@ -134,8 +133,8 @@ double exp(double x)
STRICT_ASSIGN(double, x, hi - lo);
} else if(hx < 0x3e300000) { /* |x| < 2**-28 */
/* raise inexact */
if (huge+x > one)
return one+x;
if (huge+x > 1.0)
return 1.0+x;
} else
k = 0;
@ -147,8 +146,8 @@ double exp(double x)
INSERT_WORDS(twopk, 0x3ff00000+((k+1000)<<20), 0);
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
if (k == 0)
return one - ((x*c)/(c-2.0) - x);
y = one-((lo-(x*c)/(2.0-c))-hi);
return 1.0 - ((x*c)/(c-2.0) - x);
y = 1.0-((lo-(x*c)/(2.0-c))-hi);
if (k < -1021)
return y*twopk*twom1000;
if (k == 1024)

9
src/math/expf.c
View File

@ -16,7 +16,6 @@
#include "libm.h"
static const float
one = 1.0,
halF[2] = {0.5,-0.5,},
huge = 1.0e+30,
o_threshold = 8.8721679688e+01, /* 0x42b17180 */
@ -72,8 +71,8 @@ float expf(float x)
STRICT_ASSIGN(float, x, hi - lo);
} else if(hx < 0x39000000) { /* |x|<2**-14 */
/* raise inexact */
if (huge+x > one)
return one + x;
if (huge+x > 1.0f)
return 1.0f + x;
} else
k = 0;
@ -85,8 +84,8 @@ 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 - 2.0f) - x);
y = one - ((lo - (x*c)/(2.0f - c)) - hi);
return 1.0f - ((x*c)/(c - 2.0f) - x);
y = 1.0f - ((lo - (x*c)/(2.0f - c)) - hi);
if (k < -125)
return y*twopk*twom100;
if (k == 128)

8
src/math/expl.c
View File

@ -102,13 +102,13 @@ long double expl(long double x)
if (x > MAXLOGL)
return INFINITY;
if (x < MINLOGL)
return 0.0L;
return 0.0;
/* Express e**x = e**g 2**n
* = e**g e**(n loge(2))
* = e**(g + n loge(2))
*/
px = floorl(LOG2EL * x + 0.5L); /* floor() truncates toward -infinity. */
px = floorl(LOG2EL * x + 0.5); /* floor() truncates toward -infinity. */
n = px;
x -= px * C1;
x -= px * C2;
@ -120,8 +120,8 @@ long double expl(long double x)
xx = x * x;
px = x * __polevll(xx, P, 2);
x = px/(__polevll(xx, Q, 3) - px);
x = 1.0L + ldexpl(x, 1);
x = ldexpl(x, n);
x = 1.0 + 2.0 * x;
x = scalbnl(x, n);
return x;
}
#endif

15
src/math/expm1.c
View File

@ -107,7 +107,6 @@
#include "libm.h"
static const double
one = 1.0,
huge = 1.0e+300,
tiny = 1.0e-300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
@ -148,7 +147,7 @@ double expm1(double x)
if (xsb != 0) { /* x < -56*ln2, return -1.0 with inexact */
/* raise inexact */
if(x+tiny<0.0)
return tiny-one; /* return -1 */
return tiny-1.0; /* return -1 */
}
}
@ -182,7 +181,7 @@ double expm1(double x)
/* x is now in primary range */
hfx = 0.5*x;
hxs = x*hfx;
r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
r1 = 1.0+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
t = 3.0-r1*hfx;
e = hxs*((r1-t)/(6.0 - x*t));
if (k == 0) /* c is 0 */
@ -195,17 +194,17 @@ double expm1(double x)
if (k == 1) {
if (x < -0.25)
return -2.0*(e-(x+0.5));
return one+2.0*(x-e);
return 1.0+2.0*(x-e);
}
if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
y = one - (e-x);
y = 1.0 - (e-x);
if (k == 1024)
y = y*2.0*0x1p1023;
else
y = y*twopk;
return y - one;
return y - 1.0;
}
t = one;
t = 1.0;
if (k < 20) {
SET_HIGH_WORD(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */
y = t-(e-x);
@ -213,7 +212,7 @@ double expm1(double x)
} else {
SET_HIGH_WORD(t, ((0x3ff-k)<<20)); /* 2^-k */
y = x-(e+t);
y += one;
y += 1.0;
y = y*twopk;
}
return y;

15
src/math/expm1f.c
View File

@ -16,7 +16,6 @@
#include "libm.h"
static const float
one = 1.0,
huge = 1.0e+30,
tiny = 1.0e-30,
o_threshold = 8.8721679688e+01, /* 0x42b17180 */
@ -54,7 +53,7 @@ float expm1f(float x)
if (xsb != 0) { /* x < -27*ln2 */
/* raise inexact */
if (x+tiny < 0.0f)
return tiny-one; /* return -1 */
return tiny-1.0f; /* return -1 */
}
}
@ -87,7 +86,7 @@ float expm1f(float x)
/* x is now in primary range */
hfx = 0.5f*x;
hxs = x*hfx;
r1 = one+hxs*(Q1+hxs*Q2);
r1 = 1.0f+hxs*(Q1+hxs*Q2);
t = 3.0f - r1*hfx;
e = hxs*((r1-t)/(6.0f - x*t));
if (k == 0) /* c is 0 */
@ -100,17 +99,17 @@ float expm1f(float x)
if (k == 1) {
if (x < -0.25f)
return -2.0f*(e-(x+0.5f));
return one + 2.0f*(x-e);
return 1.0f + 2.0f*(x-e);
}
if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
y = one - (e - x);
y = 1.0f - (e - x);
if (k == 128)
y = y*2.0f*0x1p127f;
else
y = y*twopk;
return y - one;
return y - 1.0f;
}
t = one;
t = 1.0f;
if (k < 23) {
SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
y = t - (e - x);
@ -118,7 +117,7 @@ float expm1f(float x)
} else {
SET_FLOAT_WORD(t, (0x7f-k)<<23); /* 2^-k */
y = x - (e + t);
y += one;
y += 1.0f;
y = y*twopk;
}
return y;

6
src/math/expm1l.c
View File

@ -97,11 +97,11 @@ long double expm1l(long double x)
return x;
/* Minimum value.*/
if (x < minarg)
return -1.0L;
return -1.0;
xx = C1 + C2;
/* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
px = floorl (0.5 + x / xx);
px = floorl(0.5 + x / xx);
k = px;
/* remainder times ln 2 */
x -= px * C1;
@ -116,7 +116,7 @@ long double expm1l(long double x)
/* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
We have qx = exp(remainder ln 2) - 1, so
exp(x) - 1 = 2^k (qx + 1) - 1 = 2^k qx + 2^k - 1. */
px = ldexpl(1.0L, k);
px = scalbnl(1.0, k);
x = px * qx + (px - 1.0);
return x;
}

12
src/math/fma.c
View File

@ -247,7 +247,7 @@ static inline double add_and_denormalize(double a, double b, int scale)
INSERT_WORD64(sum.hi, hibits);
}
}
return (ldexp(sum.hi, scale));
return scalbn(sum.hi, scale);
}
/*
@ -364,7 +364,7 @@ double fma(double x, double y, double z)
}
}
if (spread <= DBL_MANT_DIG * 2)
zs = ldexp(zs, -spread);
zs = scalbn(zs, -spread);
else
zs = copysign(DBL_MIN, zs);
@ -390,7 +390,7 @@ double fma(double x, double y, double z)
*/
fesetround(oround);
volatile double vzs = zs; /* XXX gcc CSE bug workaround */
return (xy.hi + vzs + ldexp(xy.lo, spread));
return xy.hi + vzs + scalbn(xy.lo, spread);
}
if (oround != FE_TONEAREST) {
@ -400,13 +400,13 @@ double fma(double x, double y, double z)
*/
fesetround(oround);
adj = r.lo + xy.lo;
return (ldexp(r.hi + adj, spread));
return scalbn(r.hi + adj, spread);
}
adj = add_adjusted(r.lo, xy.lo);
if (spread + ilogb(r.hi) > -1023)
return (ldexp(r.hi + adj, spread));
return scalbn(r.hi + adj, spread);
else
return (add_and_denormalize(r.hi, adj, spread));
return add_and_denormalize(r.hi, adj, spread);
}
#endif

12
src/math/fmal.c
View File

@ -115,7 +115,7 @@ static inline long double add_and_denormalize(long double a, long double b, int
if (bits_lost != 1 ^ (int)(u.bits.manl & 1))
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
return (ldexp(sum.hi, scale));
return scalbnl(sum.hi, scale);
}
/*
@ -228,7 +228,7 @@ long double fmal(long double x, long double y, long double z)
}
}
if (spread <= LDBL_MANT_DIG * 2)
zs = ldexpl(zs, -spread);
zs = scalbnl(zs, -spread);
else
zs = copysignl(LDBL_MIN, zs);
@ -254,7 +254,7 @@ long double fmal(long double x, long double y, long double z)
*/
fesetround(oround);
volatile long double vzs = zs; /* XXX gcc CSE bug workaround */
return (xy.hi + vzs + ldexpl(xy.lo, spread));
return xy.hi + vzs + scalbnl(xy.lo, spread);
}
if (oround != FE_TONEAREST) {
@ -264,13 +264,13 @@ long double fmal(long double x, long double y, long double z)
*/
fesetround(oround);
adj = r.lo + xy.lo;
return (ldexpl(r.hi + adj, spread));
return scalbnl(r.hi + adj, spread);
}
adj = add_adjusted(r.lo, xy.lo);
if (spread + ilogbl(r.hi) > -16383)
return (ldexpl(r.hi + adj, spread));
return scalbnl(r.hi + adj, spread);
else
return (add_and_denormalize(r.hi, adj, spread));
return add_and_denormalize(r.hi, adj, spread);
}
#endif

3
src/math/fmod.c
View File

@ -17,7 +17,7 @@
#include "libm.h"
static const double one = 1.0, Zero[] = {0.0, -0.0,};
static const double Zero[] = {0.0, -0.0,};
double fmod(double x, double y)
{
@ -140,7 +140,6 @@ double fmod(double x, double y)
lx = hx>>(n-32); hx = sx;
}
INSERT_WORDS(x, hx|sx, lx);
x *= one; /* create necessary signal */
}
return x; /* exact output */
}

3
src/math/fmodf.c
View File

@ -20,7 +20,7 @@
#include "libm.h"
static const float one = 1.0, Zero[] = {0.0, -0.0,};
static const float Zero[] = {0.0, -0.0,};
float fmodf(float x, float y)
{
@ -99,7 +99,6 @@ float fmodf(float x, float y)
n = -126 - iy;
hx >>= n;
SET_FLOAT_WORD(x, hx|sx);
x *= one; /* create necessary signal */
}
return x; /* exact output */
}

5
src/math/fmodl.c
View File

@ -48,7 +48,7 @@ typedef uint32_t manh_t;
#define MANL_SHIFT (LDBL_MANL_SIZE - 1)
static const long double one = 1.0, Zero[] = {0.0, -0.0,};
static const long double Zero[] = {0.0, -0.0,};
/*
* fmodl(x,y)
@ -153,7 +153,6 @@ long double fmodl(long double x, long double y)
} else {
ux.bits.exp = iy + BIAS;
}
x = ux.e * one; /* create necessary signal */
return x; /* exact output */
return ux.e; /* exact output */
}
#endif

39
src/math/j0.c
View File

@ -60,7 +60,6 @@ static double pzero(double), qzero(double);
static const double
huge = 1e300,
one = 1.0,
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0, 2.00] */
@ -73,8 +72,6 @@ S02 = 1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
S03 = 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
static const double zero = 0.0;
double j0(double x)
{
double z, s,c,ss,cc,r,u,v;
@ -83,7 +80,7 @@ double j0(double x)
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000)
return one/(x*x);
return 1.0/(x*x);
x = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
@ -92,7 +89,7 @@ double j0(double x)
cc = s+c;
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x);
if ((s*c) < zero)
if (s*c < 0.0)
cc = z/ss;
else
ss = z/cc;
@ -112,20 +109,20 @@ double j0(double x)
}
if (ix < 0x3f200000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */
if (huge+x > one) {
if (huge+x > 1.0) {
if (ix < 0x3e400000) /* |x| < 2**-27 */
return one;
return one - 0.25*x*x;
return 1.0;
return 1.0 - 0.25*x*x;
}
}
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
s = one+z*(S01+z*(S02+z*(S03+z*S04)));
s = 1.0+z*(S01+z*(S02+z*(S03+z*S04)));
if (ix < 0x3FF00000) { /* |x| < 1.00 */
return one + z*(-0.25+(r/s));
return 1.0 + z*(-0.25+(r/s));
} else {
u = 0.5*x;
return (one+u)*(one-u) + z*(r/s);
return (1.0+u)*(1.0-u) + z*(r/s);
}
}
@ -151,11 +148,11 @@ double y0(double x)
ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if (ix >= 0x7ff00000)
return one/(x+x*x);
return 1.0/(x+x*x);
if ((ix|lx) == 0)
return -one/zero;
return -1.0/0.0;
if (hx < 0)
return zero/zero;
return 0.0/0.0;
if (ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
@ -178,7 +175,7 @@ double y0(double x)
*/
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x);
if (s*c < zero)
if (s*c < 0.0)
cc = z/ss;
else
ss = z/cc;
@ -197,7 +194,7 @@ double y0(double x)
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = one+z*(v01+z*(v02+z*(v03+z*v04)));
v = 1.0+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0(x)*log(x));
}
@ -286,10 +283,10 @@ static double pzero(double x)
else if (ix >= 0x40122E8B){p = pR5; q = pS5;}
else if (ix >= 0x4006DB6D){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;}
z = one/(x*x);
z = 1.0/(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]))));
return one + r/s;
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0 + r/s;
}
@ -382,8 +379,8 @@ static double qzero(double x)
else if (ix >= 0x40122E8B){p = qR5; q = qS5;}
else if (ix >= 0x4006DB6D){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;}
z = one/(x*x);
z = 1.0/(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])))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-.125 + r/s)/x;
}

39
src/math/j0f.c
View File

@ -19,7 +19,6 @@ static float pzerof(float), qzerof(float);
static const float
huge = 1e30,
one = 1.0,
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */
/* R0/S0 on [0, 2.00] */
@ -32,8 +31,6 @@ S02 = 1.1692678527e-04, /* 0x38f53697 */
S03 = 5.1354652442e-07, /* 0x3509daa6 */
S04 = 1.1661400734e-09; /* 0x30a045e8 */
static const float zero = 0.0;
float j0f(float x)
{
float z, s,c,ss,cc,r,u,v;
@ -42,7 +39,7 @@ float j0f(float x)
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000)
return one/(x*x);
return 1.0f/(x*x);
x = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
@ -51,7 +48,7 @@ float j0f(float x)
cc = s+c;
if (ix < 0x7f000000) { /* make sure x+x does not overflow */
z = -cosf(x+x);
if (s*c < zero)
if (s*c < 0.0f)
cc = z/ss;
else
ss = z/cc;
@ -71,20 +68,20 @@ float j0f(float x)
}
if (ix < 0x39000000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */
if (huge+x > one) {
if (huge+x > 1.0f) {
if (ix < 0x32000000) /* |x| < 2**-27 */
return one;
return one - 0.25f*x*x;
return 1.0f;
return 1.0f - 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)));
s = 1.0f+z*(S01+z*(S02+z*(S03+z*S04)));
if(ix < 0x3F800000) { /* |x| < 1.00 */
return one + z*(-0.25f + (r/s));
return 1.0f + z*(-0.25f + (r/s));
} else {
u = 0.5f*x;
return (one+u)*(one-u) + z*(r/s);
return (1.0f+u)*(1.0f-u) + z*(r/s);
}
}
@ -110,11 +107,11 @@ float y0f(float x)
ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if (ix >= 0x7f800000)
return one/(x+x*x);
return 1.0f/(x+x*x);
if (ix == 0)
return -one/zero;
return -1.0f/0.0f;
if (hx < 0)
return zero/zero;
return 0.0f/0.0f;
if (ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
@ -137,7 +134,7 @@ float y0f(float x)
*/
if (ix < 0x7f000000) { /* make sure x+x not overflow */
z = -cosf(x+x);
if (s*c < zero)
if (s*c < 0.0f)
cc = z/ss;
else
ss = z/cc;
@ -156,7 +153,7 @@ float y0f(float x)
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = one+z*(v01+z*(v02+z*(v03+z*v04)));
v = 1.0f+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0f(x)*logf(x));
}
@ -244,10 +241,10 @@ static float pzerof(float x)
else if (ix >= 0x40f71c58){p = pR5; q = pS5;}
else if (ix >= 0x4036db68){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;}
z = one/(x*x);
z = 1.0f/(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]))));
return one + r/s;
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0f + r/s;
}
@ -340,8 +337,8 @@ static float qzerof(float x)
else if (ix >= 0x40f71c58){p = qR5; q = qS5;}
else if (ix >= 0x4036db68){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;}
z = one/(x*x);
z = 1.0f/(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])))));
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-.125f + r/s)/x;
}

33
src/math/j1.c
View File

@ -60,7 +60,6 @@ static double pone(double), qone(double);
static const double
huge = 1e300,
one = 1.0,
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0,2] */
@ -74,8 +73,6 @@ s03 = 1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
s04 = 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
static const double zero = 0.0;
double j1(double x)
{
double z,s,c,ss,cc,r,u,v,y;
@ -84,7 +81,7 @@ double j1(double x)
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000)
return one/x;
return 1.0/x;
y = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(y);
@ -93,7 +90,7 @@ double j1(double x)
cc = s-c;
if (ix < 0x7fe00000) { /* make sure y+y not overflow */
z = cos(y+y);
if (s*c > zero)
if (s*c > 0.0)
cc = z/ss;
else
ss = z/cc;
@ -116,12 +113,12 @@ double j1(double x)
}
if (ix < 0x3e400000) { /* |x| < 2**-27 */
/* raise inexact if x!=0 */
if (huge+x > one)
if (huge+x > 1.0)
return 0.5*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))));
s = 1.0+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
r *= x;
return x*0.5 + r/s;
}
@ -151,11 +148,11 @@ double y1(double x)
ix = 0x7fffffff & hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (ix >= 0x7ff00000)
return one/(x+x*x);
return 1.0/(x+x*x);
if ((ix|lx) == 0)
return -one/zero;
return -1.0/0.0;
if (hx < 0)
return zero/zero;
return 0.0/0.0;
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
c = cos(x);
@ -163,7 +160,7 @@ double y1(double x)
cc = s-c;
if (ix < 0x7fe00000) { /* make sure x+x not overflow */
z = cos(x+x);
if (s*c > zero)
if (s*c > 0.0)
cc = z/ss;
else
ss = z/cc;
@ -192,8 +189,8 @@ double y1(double x)
return -tpi/x;
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + tpi*(j1(x)*log(x)-one/x);
v = 1.0+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + tpi*(j1(x)*log(x)-1.0/x);
}
/* For x >= 8, the asymptotic expansions of pone is
@ -282,10 +279,10 @@ static double pone(double x)
else if (ix >= 0x40122E8B){p = pr5; q = ps5;}
else if (ix >= 0x4006DB6D){p = pr3; q = ps3;}
else if (ix >= 0x40000000){p = pr2; q = ps2;}
z = one/(x*x);
z = 1.0/(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]))));
return one+ r/s;
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0+ r/s;
}
/* For x >= 8, the asymptotic expansions of qone is
@ -378,8 +375,8 @@ static double qone(double x)
else if (ix >= 0x40122E8B){p = qr5; q = qs5;}
else if (ix >= 0x4006DB6D){p = qr3; q = qs3;}
else if (ix >= 0x40000000){p = qr2; q = qs2;}
z = one/(x*x);
z = 1.0/(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])))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (.375 + r/s)/x;
}

33
src/math/j1f.c
View File

@ -19,7 +19,6 @@ static float ponef(float), qonef(float);
static const float
huge = 1e30,
one = 1.0,
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */
/* R0/S0 on [0,2] */
@ -33,8 +32,6 @@ s03 = 1.1771846857e-06, /* 0x359dffc2 */
s04 = 5.0463624390e-09, /* 0x31ad6446 */
s05 = 1.2354227016e-11; /* 0x2d59567e */
static const float zero = 0.0;
float j1f(float x)
{
float z,s,c,ss,cc,r,u,v,y;
@ -43,7 +40,7 @@ float j1f(float x)
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000)
return one/x;
return 1.0f/x;
y = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(y);
@ -52,7 +49,7 @@ float j1f(float x)
cc = s-c;
if (ix < 0x7f000000) { /* make sure y+y not overflow */
z = cosf(y+y);
if (s*c > zero)
if (s*c > 0.0f)
cc = z/ss;
else
ss = z/cc;
@ -74,12 +71,12 @@ float j1f(float x)
}
if (ix < 0x32000000) { /* |x| < 2**-27 */
/* raise inexact if x!=0 */
if (huge+x > one)
if (huge+x > 1.0f)
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))));
s = 1.0f+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
r *= x;
return 0.5f*x + r/s;
}
@ -108,11 +105,11 @@ float y1f(float x)
ix = 0x7fffffff & hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (ix >= 0x7f800000)
return one/(x+x*x);
return 1.0f/(x+x*x);
if (ix == 0)
return -one/zero;
return -1.0f/0.0f;
if (hx < 0)
return zero/zero;
return 0.0f/0.0f;
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
c = cosf(x);
@ -120,7 +117,7 @@ float y1f(float x)
cc = s-c;
if (ix < 0x7f000000) { /* make sure x+x not overflow */
z = cosf(x+x);
if (s*c > zero)
if (s*c > 0.0f)
cc = z/ss;
else
ss = z/cc;
@ -149,8 +146,8 @@ float y1f(float x)
return -tpi/x;
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + tpi*(j1f(x)*logf(x)-one/x);
v = 1.0f+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
return x*(u/v) + tpi*(j1f(x)*logf(x)-1.0f/x);
}
/* For x >= 8, the asymptotic expansions of pone is
@ -239,10 +236,10 @@ static float ponef(float x)
else if (ix >= 0x40f71c58){p = pr5; q = ps5;}
else if (ix >= 0x4036db68){p = pr3; q = ps3;}
else if (ix >= 0x40000000){p = pr2; q = ps2;}
z = one/(x*x);
z = 1.0f/(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]))));
return one + r/s;
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
return 1.0f + r/s;
}
/* For x >= 8, the asymptotic expansions of qone is
@ -335,8 +332,8 @@ static float qonef(float x)
else if (ix >= 0x40f71c58){p = qr5; q = qs5;}
else if (ix >= 0x4036db68){p = qr3; q = qs3;}
else if (ix >= 0x40000000){p = qr2; q = qs2;}
z = one/(x*x);
z = 1.0f/(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])))));
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (.375f + r/s)/x;
}

33
src/math/jn.c
View File

@ -37,12 +37,7 @@
#include "libm.h"
static const double
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
static const double zero = 0.00000000000000000000e+00;
static const double invsqrtpi = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
double jn(int n, double x)
{
@ -69,7 +64,7 @@ double jn(int n, double x)
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabs(x);
if ((ix|lx) == 0 || ix >= 0x7ff00000) /* if x is 0 or inf */
b = zero;
b = 0.0;
else if ((double)n <= x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if (ix >= 0x52D00000) { /* x > 2**302 */
@ -108,11 +103,11 @@ double jn(int n, double x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n > 33) /* underflow */
b = zero;
b = 0.0;
else {
temp = x*0.5;
b = temp;
for (a=one,i=2; i<=n; i++) {
for (a=1.0,i=2; i<=n; i++) {
a *= (double)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
@ -165,10 +160,10 @@ double jn(int n, double x)
q1 = tmp;
}
m = n+n;
for (t=zero, i = 2*(n+k); i>=m; i -= 2)
t = one/(i/x-t);
for (t=0.0, i = 2*(n+k); i>=m; i -= 2)
t = 1.0/(i/x-t);
a = t;
b = one;
b = 1.0;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
@ -178,7 +173,7 @@ double jn(int n, double x)
* likely underflow to zero
*/
tmp = n;
v = two/x;
v = 2.0/x;
tmp = tmp*log(fabs(v*tmp));
if (tmp < 7.09782712893383973096e+02) {
for (i=n-1,di=(double)(i+i); i>0; i--) {
@ -186,7 +181,7 @@ double jn(int n, double x)
b *= di;
b = b/x - a;
a = temp;
di -= two;
di -= 2.0;
}
} else {
for (i=n-1,di=(double)(i+i); i>0; i--) {
@ -194,12 +189,12 @@ double jn(int n, double x)
b *= di;
b = b/x - a;
a = temp;
di -= two;
di -= 2.0;
/* scale b to avoid spurious overflow */
if (b > 1e100) {
a /= b;
t /= b;
b = one;
b = 1.0;
}
}
}
@ -229,9 +224,9 @@ double yn(int n, double x)
if ((ix|((uint32_t)(lx|-lx))>>31) > 0x7ff00000)
return x+x;
if ((ix|lx) == 0)
return -one/zero;
return -1.0/0.0;
if (hx < 0)
return zero/zero;
return 0.0/0.0;
sign = 1;
if (n < 0) {
n = -n;
@ -242,7 +237,7 @@ double yn(int n, double x)
if (n == 1)
return sign*y1(x);
if (ix == 0x7ff00000)
return zero;
return 0.0;
if (ix >= 0x52D00000) { /* x > 2**302 */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)

32
src/math/jnf.c
View File

@ -16,12 +16,6 @@
#define _GNU_SOURCE
#include "libm.h"
static const float
two = 2.0000000000e+00, /* 0x40000000 */
one = 1.0000000000e+00; /* 0x3F800000 */
static const float zero = 0.0000000000e+00;
float jnf(int n, float x)
{
int32_t i,hx,ix, sgn;
@ -47,7 +41,7 @@ float jnf(int n, float x)
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x);
if (ix == 0 || ix >= 0x7f800000) /* if x is 0 or inf */
b = zero;
b = 0.0f;
else if((float)n <= x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = j0f(x);
@ -63,11 +57,11 @@ float jnf(int n, float x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n > 33) /* underflow */
b = zero;
b = 0.0f;
else {
temp = 0.5f * x;
b = temp;
for (a=one,i=2; i<=n; i++) {
for (a=1.0f,i=2; i<=n; i++) {
a *= (float)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
@ -121,10 +115,10 @@ float jnf(int n, float x)
q1 = tmp;
}
m = n+n;
for (t=zero, i = 2*(n+k); i>=m; i -= 2)
t = one/(i/x-t);
for (t=0.0f, i = 2*(n+k); i>=m; i -= 2)
t = 1.0f/(i/x-t);
a = t;
b = one;
b = 1.0f;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
@ -134,7 +128,7 @@ float jnf(int n, float x)
* likely underflow to zero
*/
tmp = n;
v = two/x;
v = 2.0f/x;
tmp = tmp*logf(fabsf(v*tmp));
if (tmp < 88.721679688f) {
for (i=n-1,di=(float)(i+i); i>0; i--) {
@ -142,7 +136,7 @@ float jnf(int n, float x)
b *= di;
b = b/x - a;
a = temp;
di -= two;
di -= 2.0f;
}
} else {
for (i=n-1,di=(float)(i+i); i>0; i--){
@ -150,12 +144,12 @@ float jnf(int n, float x)
b *= di;
b = b/x - a;
a = temp;
di -= two;
di -= 2.0f;
/* scale b to avoid spurious overflow */
if (b > 1e10f) {
a /= b;
t /= b;
b = one;
b = 1.0f;
}
}
}
@ -183,9 +177,9 @@ float ynf(int n, float x)
if (ix > 0x7f800000)
return x+x;
if (ix == 0)
return -one/zero;
return -1.0f/0.0f;
if (hx < 0)
return zero/zero;
return 0.0f/0.0f;
sign = 1;
if (n < 0) {
n = -n;
@ -196,7 +190,7 @@ float ynf(int n, float x)
if (n == 1)
return sign*y1f(x);
if (ix == 0x7f800000)
return zero;
return 0.0f;
a = y0f(x);
b = y1f(x);

50
src/math/lgamma_r.c
View File

@ -82,8 +82,6 @@
static const double
two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
@ -148,8 +146,6 @@ w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
static const double zero = 0.00000000000000000000e+00;
static double sin_pi(double x)
{
double y,z;
@ -159,7 +155,7 @@ static double sin_pi(double x)
ix &= 0x7fffffff;
if (ix < 0x3fd00000)
return __sin(pi*x, zero, 0);
return __sin(pi*x, 0.0, 0);
y = -x; /* negative x is assumed */
@ -174,7 +170,7 @@ static double sin_pi(double x)
n = (int)(y*4.0);
} else {
if (ix >= 0x43400000) {
y = zero; /* y must be even */
y = 0.0; /* y must be even */
n = 0;
} else {
if (ix < 0x43300000)
@ -186,14 +182,14 @@ static double sin_pi(double x)
}
}
switch (n) {
case 0: y = __sin(pi*y, zero, 0); break;
case 0: y = __sin(pi*y, 0.0, 0); break;
case 1:
case 2: y = __cos(pi*(0.5-y), zero); break;
case 2: y = __cos(pi*(0.5-y), 0.0); break;
case 3:
case 4: y = __sin(pi*(one-y), zero, 0); break;
case 4: y = __sin(pi*(1.0-y), 0.0, 0); break;
case 5:
case 6: y = -__cos(pi*(y-1.5), zero); break;
default: y = __sin(pi*(y-2.0), zero, 0); break;
case 6: y = -__cos(pi*(y-1.5), 0.0); break;
default: y = __sin(pi*(y-2.0), 0.0, 0); break;
}
return -y;
}
@ -213,7 +209,7 @@ double __lgamma_r(double x, int *signgamp)
if (ix >= 0x7ff00000)
return x*x;
if ((ix|lx) == 0)
return one/zero;
return 1.0/0.0;
if (ix < 0x3b900000) { /* |x|<2**-70, return -log(|x|) */
if(hx < 0) {
*signgamp = -1;
@ -223,12 +219,12 @@ double __lgamma_r(double x, int *signgamp)
}
if (hx < 0) {
if (ix >= 0x43300000) /* |x|>=2**52, must be -integer */
return one/zero;
return 1.0/0.0;
t = sin_pi(x);
if (t == zero) /* -integer */
return one/zero;
if (t == 0.0) /* -integer */
return 1.0/0.0;
nadj = log(pi/fabs(t*x));
if (t < zero)
if (t < 0.0)
*signgamp = -1;
x = -x;
}
@ -241,17 +237,17 @@ double __lgamma_r(double x, int *signgamp)
if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -log(x);
if (ix >= 0x3FE76944) {
y = one - x;
y = 1.0 - x;
i = 0;
} else if (ix >= 0x3FCDA661) {
y = x - (tc-one);
y = x - (tc-1.0);
i = 1;
} else {
y = x;
i = 2;
}
} else {
r = zero;
r = 0.0;
if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */
y = 2.0 - x;
i = 0;
@ -259,7 +255,7 @@ double __lgamma_r(double x, int *signgamp)
y = x - tc;
i = 1;
} else {
y = x - one;
y = x - 1.0;
i = 2;
}
}
@ -282,16 +278,16 @@ double __lgamma_r(double x, int *signgamp)
break;
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))));
p2 = 1.0+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += -0.5*y + p1/p2;
}
} else if (ix < 0x40200000) { /* x < 8.0 */
i = (int)x;
y = x - (double)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) */
q = 1.0+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
r = 0.5*y+p/q;
z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) {
case 7: z *= y + 6.0; /* FALLTHRU */
case 6: z *= y + 5.0; /* FALLTHRU */
@ -303,12 +299,12 @@ double __lgamma_r(double x, int *signgamp)
}
} else if (ix < 0x43900000) { /* 8.0 <= x < 2**58 */
t = log(x);
z = one/x;
z = 1.0/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
r = (x-half)*(t-one)+w;
r = (x-0.5)*(t-1.0)+w;
} else /* 2**58 <= x <= inf */
r = x*(log(x)-one);
r = x*(log(x)-1.0);
if (hx < 0)
r = nadj - r;
return r;

40
src/math/lgammaf_r.c
View File

@ -17,8 +17,6 @@
static const float
two23= 8.3886080000e+06, /* 0x4b000000 */
half= 5.0000000000e-01, /* 0x3f000000 */
one = 1.0000000000e+00, /* 0x3f800000 */
pi = 3.1415927410e+00, /* 0x40490fdb */
a0 = 7.7215664089e-02, /* 0x3d9e233f */
a1 = 3.2246702909e-01, /* 0x3ea51a66 */
@ -83,8 +81,6 @@ w4 = -5.9518753551e-04, /* 0xba1c065c */
w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */
w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
static const float zero = 0.0000000000e+00;
static float sin_pif(float x)
{
float y,z;
@ -109,7 +105,7 @@ static float sin_pif(float x)
n = (int)(y*4.0f);
} else {
if (ix >= 0x4b800000) {
y = zero; /* y must be even */
y = 0.0f; /* y must be even */
n = 0;
} else {
if (ix < 0x4b000000)
@ -125,7 +121,7 @@ static float sin_pif(float x)
case 1:
case 2: y = __cosdf(pi*(0.5f - y)); break;
case 3:
case 4: y = __sindf(pi*(one - y)); break;
case 4: y = __sindf(pi*(1.0f - y)); break;
case 5:
case 6: y = -__cosdf(pi*(y - 1.5f)); break;
default: y = __sindf(pi*(y - 2.0f)); break;
@ -148,7 +144,7 @@ float __lgammaf_r(float x, int *signgamp)
if (ix >= 0x7f800000)
return x*x;
if (ix == 0)
return one/zero;
return 1.0f/0.0f;
if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */
if (hx < 0) {
*signgamp = -1;
@ -158,12 +154,12 @@ float __lgammaf_r(float x, int *signgamp)
}
if (hx < 0) {
if (ix >= 0x4b000000) /* |x| >= 2**23, must be -integer */
return one/zero;
return 1.0f/0.0f;
t = sin_pif(x);
if (t == zero) /* -integer */
return one/zero;
if (t == 0.0f) /* -integer */
return 1.0f/0.0f;
nadj = logf(pi/fabsf(t*x));
if (t < zero)
if (t < 0.0f)
*signgamp = -1;
x = -x;
}
@ -176,17 +172,17 @@ float __lgammaf_r(float x, int *signgamp)
if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -logf(x);
if (ix >= 0x3f3b4a20) {
y = one - x;
y = 1.0f - x;
i = 0;
} else if (ix >= 0x3e6d3308) {
y = x - (tc-one);
y = x - (tc-1.0f);
i = 1;
} else {
y = x;
i = 2;
}
} else {
r = zero;
r = 0.0f;
if (ix >= 0x3fdda618) { /* [1.7316,2] */
y = 2.0f - x;
i = 0;
@ -194,7 +190,7 @@ float __lgammaf_r(float x, int *signgamp)
y = x - tc;
i = 1;
} else {
y = x - one;
y = x - 1.0f;
i = 2;
}
}
@ -217,16 +213,16 @@ float __lgammaf_r(float x, int *signgamp)
break;
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))));
p2 = 1.0f+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += -0.5f*y + p1/p2;
}
} else if (ix < 0x41000000) { /* x < 8.0 */
i = (int)x;
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) */
q = 1.0f+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
r = 0.5f*y+p/q;
z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) {
case 7: z *= y + 6.0f; /* FALLTHRU */
case 6: z *= y + 5.0f; /* FALLTHRU */
@ -238,12 +234,12 @@ float __lgammaf_r(float x, int *signgamp)
}
} else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */
t = logf(x);
z = one/x;
z = 1.0f/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
r = (x-half)*(t-one)+w;
r = (x-0.5f)*(t-1.0f)+w;
} else /* 2**58 <= x <= inf */
r = x*(logf(x)-one);
r = x*(logf(x)-1.0f);
if (hx < 0)
r = nadj - r;
return r;

48
src/math/lgammal.c
View File

@ -95,8 +95,6 @@ long double __lgammal_r(long double x, int *sg)
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
half = 0.5L,
one = 1.0L,
pi = 3.14159265358979323846264L,
two63 = 9.223372036854775808e18L,
@ -200,8 +198,6 @@ w5 = 8.412723297322498080632E-4L,
w6 = -1.880801938119376907179E-3L,
w7 = 4.885026142432270781165E-3L;
static const long double zero = 0.0L;
static long double sin_pi(long double x)
{
long double y, z;
@ -226,7 +222,7 @@ static long double sin_pi(long double x)
n = (int) (y*4.0);
} else {
if (ix >= 0x403f8000) { /* 2^64 */
y = zero; /* y must be even */
y = 0.0; /* y must be even */
n = 0;
} else {
if (ix < 0x403e8000) /* 2^63 */
@ -244,11 +240,11 @@ static long double sin_pi(long double x)
break;
case 1:
case 2:
y = cosl(pi * (half - y));
y = cosl(pi * (0.5 - y));
break;
case 3:
case 4:
y = sinl(pi * (one - y));
y = sinl(pi * (1.0 - y));
break;
case 5:
case 6:
@ -273,7 +269,7 @@ long double __lgammal_r(long double x, int *sg) {
if ((ix | i0 | i1) == 0) {
if (se & 0x8000)
*sg = -1;
return one / fabsl(x);
return 1.0 / fabsl(x);
}
ix = (ix << 16) | (i0 >> 16);
@ -291,10 +287,10 @@ long double __lgammal_r(long double x, int *sg) {
}
if (se & 0x8000) {
t = sin_pi (x);
if (t == zero)
return one / fabsl(t); /* -integer */
if (t == 0.0)
return 1.0 / fabsl(t); /* -integer */
nadj = logl(pi / fabsl(t * x));
if (t < zero)
if (t < 0.0)
*sg = -1;
x = -x;
}
@ -306,19 +302,19 @@ long double __lgammal_r(long double x, int *sg) {
else if (ix < 0x40008000) { /* x < 2.0 */
if (ix <= 0x3ffee666) { /* 8.99993896484375e-1 */
/* lgamma(x) = lgamma(x+1) - log(x) */
r = -logl (x);
r = -logl(x);
if (ix >= 0x3ffebb4a) { /* 7.31597900390625e-1 */
y = x - one;
y = x - 1.0;
i = 0;
} else if (ix >= 0x3ffced33) { /* 2.31639862060546875e-1 */
y = x - (tc - one);
y = x - (tc - 1.0);
i = 1;
} else { /* x < 0.23 */
y = x;
i = 2;
}
} else {
r = zero;
r = 0.0;
if (ix >= 0x3fffdda6) { /* 1.73162841796875 */
/* [1.7316,2] */
y = x - 2.0;
@ -329,7 +325,7 @@ long double __lgammal_r(long double x, int *sg) {
i = 1;
} else {
/* [0.9, 1.23] */
y = x - one;
y = x - 1.0;
i = 2;
}
}
@ -337,7 +333,7 @@ long double __lgammal_r(long double x, int *sg) {
case 0:
p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
r += half * y + y * p1/p2;
r += 0.5 * y + y * p1/p2;
break;
case 1:
p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
@ -348,17 +344,17 @@ long double __lgammal_r(long double x, int *sg) {
case 2:
p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
r += (-half * y + p1 / p2);
r += (-0.5 * y + p1 / p2);
}
} else if (ix < 0x40028000) { /* 8.0 */
/* x < 8.0 */
i = (int)x;
t = zero;
t = 0.0;
y = x - (double)i;
p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
r = half * y + p / q;
z = one;/* lgamma(1+s) = log(s) + lgamma(s) */
r = 0.5 * y + p / q;
z = 1.0;/* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) {
case 7:
z *= (y + 6.0); /* FALLTHRU */
@ -370,18 +366,18 @@ long double __lgammal_r(long double x, int *sg) {
z *= (y + 3.0); /* FALLTHRU */
case 3:
z *= (y + 2.0); /* FALLTHRU */
r += logl (z);
r += logl(z);
break;
}
} else if (ix < 0x40418000) { /* 2^66 */
/* 8.0 <= x < 2**66 */
t = logl (x);
z = one / x;
t = logl(x);
z = 1.0 / x;
y = z * z;
w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
r = (x - half) * (t - one) + w;
r = (x - 0.5) * (t - 1.0) + w;
} else /* 2**66 <= x <= inf */
r = x * (logl (x) - one);
r = x * (logl(x) - 1.0);
if (se & 0x8000)
r = nadj - r;
return r;

10
src/math/log.c
View File

@ -74,8 +74,6 @@ Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
double log(double x)
{
double hfsq,f,s,z,R,w,t1,t2,dk;
@ -87,9 +85,9 @@ double log(double x)
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
return -two54/zero; /* log(+-0)=-inf */
return -two54/0.0; /* log(+-0)=-inf */
if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */
return (x-x)/0.0; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 54;
x *= two54;
@ -104,9 +102,9 @@ double log(double x)
k += i>>20;
f = x - 1.0;
if ((0x000fffff&(2+hx)) < 3) { /* -2**-20 <= f < 2**-20 */
if (f == zero) {
if (f == 0.0) {
if (k == 0) {
return zero;
return 0.0;
}
dk = (double)k;
return dk*ln2_hi + dk*ln2_lo;

8
src/math/log10.c
View File

@ -27,8 +27,6 @@ ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
static const double zero = 0.0;
double log10(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
@ -40,9 +38,9 @@ double log10(double x)
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
return -two54/zero; /* log(+-0)=-inf */
return -two54/0.0; /* log(+-0)=-inf */
if (hx<0)
return (x-x)/zero; /* log(-#) = NaN */
return (x-x)/0.0; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 54;
x *= two54;
@ -51,7 +49,7 @@ double log10(double x)
if (hx >= 0x7ff00000)
return x+x;
if (hx == 0x3ff00000 && lx == 0)
return zero; /* log(1) = +0 */
return 0.0; /* log(1) = +0 */
k += (hx>>20) - 1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;

8
src/math/log10f.c
View File

@ -23,8 +23,6 @@ ivln10lo = -3.1689971365e-05, /* 0xb804ead9 */
log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */
log10_2lo = 7.9034151668e-07; /* 0x355427db */
static const float zero = 0.0;
float log10f(float x)
{
float f,hfsq,hi,lo,r,y;
@ -35,9 +33,9 @@ float log10f(float x)
k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0)
return -two25/zero; /* log(+-0)=-inf */
return -two25/0.0f; /* log(+-0)=-inf */
if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */
return (x-x)/0.0f; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 25;
x *= two25;
@ -46,7 +44,7 @@ float log10f(float x)
if (hx >= 0x7f800000)
return x+x;
if (hx == 0x3f800000)
return zero; /* log(1) = +0 */
return 0.0f; /* log(1) = +0 */
k += (hx>>23) - 127;
hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000;

22
src/math/log10l.c
View File

@ -123,9 +123,9 @@ long double log10l(long double x)
if (isnan(x))
return x;
if(x <= 0.0L) {
if(x == 0.0L)
return -1.0L / (x - x);
if(x <= 0.0) {
if(x == 0.0)
return -1.0 / (x - x);
return (x - x) / (x - x);
}
if (x == INFINITY)
@ -142,12 +142,12 @@ long double log10l(long double x)
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
z = x - 0.5;
y = 0.5 * z + 0.5;
} else { /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
z = x*x;
@ -158,13 +158,13 @@ long double log10l(long double x)
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
x = ldexpl(x, 1) - 1.0L; /* 2x - 1 */
x = 2.0*x - 1.0;
} else {
x = x - 1.0L;
x = x - 1.0;
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
y = y - ldexpl(z, -1); /* -0.5x^2 + ... */
y = y - 0.5*z;
done:
/* Multiply log of fraction by log10(e)

10
src/math/log1p.c
View File

@ -88,8 +88,6 @@ Lp5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
double log1p(double x)
{
double hfsq,f,c,s,z,R,u;
@ -102,12 +100,12 @@ double log1p(double x)
if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */
if (ax >= 0x3ff00000) { /* x <= -1.0 */
if (x == -1.0)
return -two54/zero; /* log1p(-1)=+inf */
return -two54/0.0; /* log1p(-1)=+inf */
return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if (ax < 0x3e200000) { /* |x| < 2**-29 */
/* raise inexact */
if (two54 + x > zero && ax < 0x3c900000) /* |x| < 2**-54 */
if (two54 + x > 0.0 && ax < 0x3c900000) /* |x| < 2**-54 */
return x;
return x - x*x*0.5;
}
@ -151,9 +149,9 @@ double log1p(double x)
}
hfsq = 0.5*f*f;
if (hu == 0) { /* |f| < 2**-20 */
if (f == zero) {
if (f == 0.0) {
if(k == 0)
return zero;
return 0.0;
c += k*ln2_lo;
return k*ln2_hi + c;
}

10
src/math/log1pf.c
View File

@ -27,8 +27,6 @@ Lp5 = 1.8183572590e-01, /* 3E3A3325 */
Lp6 = 1.5313838422e-01, /* 3E1CD04F */
Lp7 = 1.4798198640e-01; /* 3E178897 */
static const float zero = 0.0;
float log1pf(float x)
{
float hfsq,f,c,s,z,R,u;
@ -41,12 +39,12 @@ float log1pf(float x)
if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */
if (ax >= 0x3f800000) { /* x <= -1.0 */
if (x == -1.0f)
return -two25/zero; /* log1p(-1)=+inf */
return -two25/0.0f; /* log1p(-1)=+inf */
return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if (ax < 0x38000000) { /* |x| < 2**-15 */
/* raise inexact */
if (two25 + x > zero && ax < 0x33800000) /* |x| < 2**-24 */
if (two25 + x > 0.0f && ax < 0x33800000) /* |x| < 2**-24 */
return x;
return x - x*x*0.5f;
}
@ -91,9 +89,9 @@ float log1pf(float x)
}
hfsq = 0.5f * f * f;
if (hu == 0) { /* |f| < 2**-20 */
if (f == zero) {
if (f == 0.0f) {
if (k == 0)
return zero;
return 0.0f;
c += k*ln2_lo;
return k*ln2_hi+c;
}

20
src/math/log1pl.c
View File

@ -118,11 +118,11 @@ long double log1pl(long double xm1)
if (xm1 == 0.0)
return xm1;
x = xm1 + 1.0L;
x = xm1 + 1.0;
/* Test for domain errors. */
if (x <= 0.0L) {
if (x == 0.0L)
if (x <= 0.0) {
if (x == 0.0)
return -INFINITY;
return NAN;
}
@ -136,12 +136,12 @@ long double log1pl(long double xm1)
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
z = x - 0.5;
y = 0.5 * z + 0.5;
} else { /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
z = x*x;
@ -156,12 +156,12 @@ long double log1pl(long double xm1)
if (x < SQRTH) {
e -= 1;
if (e != 0)
x = 2.0 * x - 1.0L;
x = 2.0 * x - 1.0;
else
x = xm1;
} else {
if (e != 0)
x = x - 1.0L;
x = x - 1.0;
else
x = xm1;
}

8
src/math/log2.c
View File

@ -27,8 +27,6 @@ two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
static const double zero = 0.0;
double log2(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
@ -40,9 +38,9 @@ double log2(double x)
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
return -two54/zero; /* log(+-0)=-inf */
return -two54/0.0; /* log(+-0)=-inf */
if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */
return (x-x)/0.0; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 54;
x *= two54;
@ -51,7 +49,7 @@ double log2(double x)
if (hx >= 0x7ff00000)
return x+x;
if (hx == 0x3ff00000 && lx == 0)
return zero; /* log(1) = +0 */
return 0.0; /* log(1) = +0 */
k += (hx>>20) - 1023;
hx &= 0x000fffff;
i = (hx+0x95f64) & 0x100000;

8
src/math/log2f.c
View File

@ -21,8 +21,6 @@ two25 = 3.3554432000e+07, /* 0x4c000000 */
ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */
ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
static const float zero = 0.0;
float log2f(float x)
{
float f,hfsq,hi,lo,r,y;
@ -33,9 +31,9 @@ float log2f(float x)
k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0)
return -two25/zero; /* log(+-0)=-inf */
return -two25/0.0f; /* log(+-0)=-inf */
if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */
return (x-x)/0.0f; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 25;
x *= two25;
@ -44,7 +42,7 @@ float log2f(float x)
if (hx >= 0x7f800000)
return x+x;
if (hx == 0x3f800000)
return zero; /* log(1) = +0 */
return 0.0f; /* log(1) = +0 */
k += (hx>>23) - 127;
hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000;

20
src/math/log2l.c
View File

@ -121,8 +121,8 @@ long double log2l(long double x)
return x;
if (x == INFINITY)
return x;
if (x <= 0.0L) {
if (x == 0.0L)
if (x <= 0.0) {
if (x == 0.0)
return -INFINITY;
return NAN;
}
@ -139,12 +139,12 @@ long double log2l(long double x)
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
z = x - 0.5;
y = 0.5 * z + 0.5;
} else { /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
z = x*x;
@ -155,13 +155,13 @@ long double log2l(long double x)
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
x = ldexpl(x, 1) - 1.0L; /* 2x - 1 */
x = 2.0*x - 1.0;
} else {
x = x - 1.0L;
x = x - 1.0;
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
y = y - ldexpl(z, -1); /* -0.5x^2 + ... */
y = y - 0.5*z;
done:
/* Multiply log of fraction by log2(e)

10
src/math/logf.c
View File

@ -25,8 +25,6 @@ Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
static const float zero = 0.0;
float logf(float x)
{
float hfsq,f,s,z,R,w,t1,t2,dk;
@ -37,9 +35,9 @@ float logf(float x)
k = 0;
if (ix < 0x00800000) { /* x < 2**-126 */
if ((ix & 0x7fffffff) == 0)
return -two25/zero; /* log(+-0)=-inf */
return -two25/0.0f; /* log(+-0)=-inf */
if (ix < 0)
return (x-x)/zero; /* log(-#) = NaN */
return (x-x)/0.0f; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 25;
x *= two25;
@ -54,9 +52,9 @@ float logf(float x)
k += i>>23;
f = x - 1.0f;
if ((0x007fffff & (0x8000 + ix)) < 0xc000) { /* -2**-9 <= f < 2**-9 */
if (f == zero) {
if (f == 0.0f) {
if (k == 0)
return zero;
return 0.0f;
dk = (float)k;
return dk*ln2_hi + dk*ln2_lo;
}

20
src/math/logl.c
View File

@ -119,8 +119,8 @@ long double logl(long double x)
return x;
if (x == INFINITY)
return x;
if (x <= 0.0L) {
if (x == 0.0L)
if (x <= 0.0) {
if (x == 0.0)
return -INFINITY;
return NAN;
}
@ -137,12 +137,12 @@ long double logl(long double x)
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
z = x - 0.5;
y = 0.5 * z + 0.5;
} else { /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
z = x*x;
@ -156,14 +156,14 @@ long double logl(long double x)
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
x = ldexpl(x, 1) - 1.0L; /* 2x - 1 */
x = 2.0*x - 1.0;
} else {
x = x - 1.0L;
x = x - 1.0;
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
y = y + e * C2;
z = y - ldexpl(z, -1); /* y - 0.5 * z */
z = y - 0.5*z;
/* Note, the sum of above terms does not exceed x/4,
* so it contributes at most about 1/4 lsb to the error.
*/

4
src/math/modf.c
View File

@ -21,8 +21,6 @@
#include "libm.h"
static const double one = 1.0;
double modf(double x, double *iptr)
{
int32_t i0,i1,j0;
@ -51,7 +49,7 @@ double modf(double x, double *iptr)
*iptr = x;
return 0.0 / x;
}
*iptr = x*one;
*iptr = x;
GET_HIGH_WORD(high, x);
INSERT_WORDS(x, high & 0x80000000, 0); /* return +-0 */
return x;

66
src/math/modff.c
View File

@ -1,51 +1,33 @@
/* origin: FreeBSD /usr/src/lib/msun/src/s_modff.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* 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"
static const float one = 1.0;
float modff(float x, float *iptr)
{
int32_t i0,j0;
uint32_t i;
uint32_t u, mask;
int e;
GET_FLOAT_WORD(i0, x);
j0 = ((i0>>23) & 0xff) - 0x7f; /* exponent of x */
if (j0 < 23) { /* integer part in x */
if (j0 < 0) { /* |x| < 1 */
SET_FLOAT_WORD(*iptr, i0 & 0x80000000); /* *iptr = +-0 */
GET_FLOAT_WORD(u, x);
e = (int)(u>>23 & 0xff) - 0x7f;
/* no fractional part */
if (e >= 23) {
*iptr = x;
if (e == 0x80 && u<<9 != 0) /* nan */
return x;
}
i = 0x007fffff >> j0;
if ((i0&i) == 0) { /* x is integral */
uint32_t ix;
*iptr = x;
GET_FLOAT_WORD(ix, x);
SET_FLOAT_WORD(x, ix & 0x80000000); /* return +-0 */
return x;
}
SET_FLOAT_WORD(*iptr, i0&~i);
return x - *iptr;
} else { /* no fraction part */
uint32_t ix;
*iptr = x*one;
if (x != x) /* NaN */
return x;
GET_FLOAT_WORD(ix, x);
SET_FLOAT_WORD(x, ix & 0x80000000); /* return +-0 */
SET_FLOAT_WORD(x, u & 0x80000000);
return x;
}
/* no integral part */
if (e < 0) {
SET_FLOAT_WORD(*iptr, u & 0x80000000);
return x;
}
mask = 0x007fffff>>e;
if ((u & mask) == 0) {
*iptr = x;
SET_FLOAT_WORD(x, u & 0x80000000);
return x;
}
SET_FLOAT_WORD(*iptr, u & ~mask);
return x - *iptr;
}

4
src/math/modfl.c
View File

@ -54,7 +54,7 @@ long double modfl(long double x, long double *iptr)
/* The number of fraction bits in manh, not counting the integer bit */
#define HIBITS (LDBL_MANT_DIG - LDBL_MANL_SIZE)
static const long double zero[] = { 0.0L, -0.0L };
static const long double zero[] = { 0.0, -0.0 };
long double modfl(long double x, long double *iptr)
{
@ -81,7 +81,7 @@ long double modfl(long double x, long double *iptr)
return x - ux.e;
} else if (e >= LDBL_MANT_DIG - 1) { /* x has no fraction part. */
*iptr = x;
if (x != x) /* Handle NaNs. */
if (e == LDBL_MAX_EXP && (ux.bits.manh|ux.bits.manl)) /* nan */
return x;
return zero[ux.bits.sign];
} else { /* Fraction part is in manl. */

37
src/math/pow.c
View File

@ -59,9 +59,6 @@ static const double
bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
zero = 0.0,
one = 1.0,
two = 2.0,
two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
huge = 1.0e300,
tiny = 1.0e-300,
@ -101,15 +98,15 @@ double pow(double x, double y)
ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
/* y == zero: x**0 = 1 */
/* y == 0.0: x**0 = 1 */
if ((iy|ly) == 0)
return one;
return 1.0;
/* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3ff00000 && lx == 0)
return one;
return 1.0;
/* y != zero: result is NaN if either arg is NaN */
/* y != 0.0: 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);
@ -141,15 +138,15 @@ double pow(double x, double y)
if (ly == 0) {
if (iy == 0x7ff00000) { /* y is +-inf */
if (((ix-0x3ff00000)|lx) == 0) /* (-1)**+-inf is 1 */
return one;
return 1.0;
else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
return hy >= 0 ? y : zero;
return hy >= 0 ? y : 0.0;
else /* (|x|<1)**+-inf = 0,inf */
return hy < 0 ? -y : zero;
return hy < 0 ? -y : 0.0;
}
if (iy == 0x3ff00000) { /* y is +-1 */
if (hy < 0)
return one/x;
return 1.0/x;
return x;
}
if (hy == 0x40000000) /* y is 2 */
@ -166,7 +163,7 @@ double pow(double x, double y)
if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) { /* x is +-0,+-inf,+-1 */
z = ax;
if (hy < 0) /* z = (1/|x|) */
z = one/z;
z = 1.0/z;
if (hx < 0) {
if (((ix-0x3ff00000)|yisint) == 0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
@ -187,9 +184,9 @@ double pow(double x, double y)
if ((n|yisint) == 0)
return (x-x)/(x-x);
s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
s = 1.0; /* s (sign of result -ve**odd) = -1 else = 1 */
if ((n|(yisint-1)) == 0)
s = -one;/* (-ve)**(odd int) */
s = -1.0;/* (-ve)**(odd int) */
/* |y| is huge */
if (iy > 0x41e00000) { /* if |y| > 2**31 */
@ -206,7 +203,7 @@ double pow(double x, double y)
return hy > 0 ? s*huge*huge : s*tiny*tiny;
/* 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 - one; /* t has 20 trailing zeros */
t = ax - 1.0; /* t has 20 trailing zeros */
w = (t*t)*(0.5 - t*(0.3333333333333333333333-t*0.25));
u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l - w*ivln2;
@ -239,12 +236,12 @@ double pow(double x, double y)
/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one/(ax+bp[k]);
v = 1.0/(ax+bp[k]);
ss = u*v;
s_h = ss;
SET_LOW_WORD(s_h, 0);
/* t_h=ax+bp[k] High */
t_h = zero;
t_h = 0.0;
SET_HIGH_WORD(t_h, ((ix>>1)|0x20000000) + 0x00080000 + (k<<18));
t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l);
@ -299,7 +296,7 @@ double pow(double x, double y)
if (i > 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
n = j + (0x00100000>>(k+1));
k = ((n&0x7fffffff)>>20) - 0x3ff; /* new k for n */
t = zero;
t = 0.0;
SET_HIGH_WORD(t, n & ~(0x000fffff>>k));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if (j < 0)
@ -314,8 +311,8 @@ double pow(double x, double y)
w = v - (z-u);
t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
r = (z*t1)/(t1-two) - (w + z*w);
z = one - (r-z);
r = (z*t1)/(t1-2.0) - (w + z*w);
z = 1.0 - (r-z);
GET_HIGH_WORD(j, z);
j += n<<20;
if ((j>>20) <= 0) /* subnormal output */

31
src/math/powf.c
View File

@ -19,9 +19,6 @@ static const float
bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
zero = 0.0,
one = 1.0,
two = 2.0,
two24 = 16777216.0, /* 0x4b800000 */
huge = 1.0e30,
tiny = 1.0e-30,
@ -60,15 +57,15 @@ float powf(float x, float y)
ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
/* y == zero: x**0 = 1 */
/* y == 0: x**0 = 1 */
if (iy == 0)
return one;
return 1.0f;
/* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3f800000)
return one;
return 1.0f;
/* y != zero: result is NaN if either arg is NaN */
/* y != 0: result is NaN if either arg is NaN */
if (ix > 0x7f800000 || iy > 0x7f800000)
return (x+0.0f) + (y+0.0f);
@ -92,15 +89,15 @@ float powf(float x, float y)
/* special value of y */
if (iy == 0x7f800000) { /* y is +-inf */
if (ix == 0x3f800000) /* (-1)**+-inf is 1 */
return one;
return 1.0f;
else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */
return hy >= 0 ? y : zero;
return hy >= 0 ? y : 0.0f;
else /* (|x|<1)**+-inf = 0,inf */
return hy < 0 ? -y : zero;
return hy < 0 ? -y : 0.0f;
}
if (iy == 0x3f800000) { /* y is +-1 */
if (hy < 0)
return one/x;
return 1.0f/x;
return x;
}
if (hy == 0x40000000) /* y is 2 */
@ -115,7 +112,7 @@ float powf(float x, float y)
if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) { /* x is +-0,+-inf,+-1 */
z = ax;
if (hy < 0) /* z = (1/|x|) */
z = one/z;
z = 1.0f/z;
if (hx < 0) {
if (((ix-0x3f800000)|yisint) == 0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
@ -131,9 +128,9 @@ float powf(float x, float y)
if ((n|yisint) == 0)
return (x-x)/(x-x);
sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */
sn = 1.0f; /* s (sign of result -ve**odd) = -1 else = 1 */
if ((n|(yisint-1)) == 0) /* (-ve)**(odd int) */
sn = -one;
sn = -1.0f;
/* |y| is huge */
if (iy > 0x4d000000) { /* if |y| > 2**27 */
@ -178,7 +175,7 @@ float powf(float x, float y)
/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one/(ax+bp[k]);
v = 1.0f/(ax+bp[k]);
s = u*v;
s_h = s;
GET_FLOAT_WORD(is, s_h);
@ -257,8 +254,8 @@ float powf(float x, float y)
w = v - (z - u);
t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
r = (z*t1)/(t1-two) - (w+z*w);
z = one - (r - z);
r = (z*t1)/(t1-2.0f) - (w+z*w);
z = 1.0f - (r - z);
GET_FLOAT_WORD(j, z);
j += n<<23;
if ((j>>23) <= 0) /* subnormal output */

103
src/math/powl.c
View File

@ -203,44 +203,44 @@ long double powl(long double x, long double y)
volatile long double z=0;
long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0;
if (y == 0.0L)
return 1.0L;
if (y == 0.0)
return 1.0;
if (isnan(x))
return x;
if (isnan(y))
return y;
if (y == 1.0L)
if (y == 1.0)
return x;
// FIXME: this is wrong, see pow special cases in c99 F.9.4.4
if (!isfinite(y) && (x == -1.0L || x == 1.0L) )
if (!isfinite(y) && (x == -1.0 || x == 1.0) )
return y - y; /* +-1**inf is NaN */
if (x == 1.0L)
return 1.0L;
if (x == 1.0)
return 1.0;
if (y >= LDBL_MAX) {
if (x > 1.0L)
if (x > 1.0)
return INFINITY;
if (x > 0.0L && x < 1.0L)
return 0.0L;
if (x < -1.0L)
if (x > 0.0 && x < 1.0)
return 0.0;
if (x < -1.0)
return INFINITY;
if (x > -1.0L && x < 0.0L)
return 0.0L;
if (x > -1.0 && x < 0.0)
return 0.0;
}
if (y <= -LDBL_MAX) {
if (x > 1.0L)
return 0.0L;
if (x > 0.0L && x < 1.0L)
if (x > 1.0)
return 0.0;
if (x > 0.0 && x < 1.0)
return INFINITY;
if (x < -1.0L)
return 0.0L;
if (x > -1.0L && x < 0.0L)
if (x < -1.0)
return 0.0;
if (x > -1.0 && x < 0.0)
return INFINITY;
}
if (x >= LDBL_MAX) {
if (y > 0.0L)
if (y > 0.0)
return INFINITY;
return 0.0L;
return 0.0;
}
w = floorl(y);
@ -253,29 +253,29 @@ long double powl(long double x, long double y)
yoddint = 0;
if (iyflg) {
ya = fabsl(y);
ya = floorl(0.5L * ya);
yb = 0.5L * fabsl(w);
ya = floorl(0.5 * ya);
yb = 0.5 * fabsl(w);
if( ya != yb )
yoddint = 1;
}
if (x <= -LDBL_MAX) {
if (y > 0.0L) {
if (y > 0.0) {
if (yoddint)
return -INFINITY;
return INFINITY;
}
if (y < 0.0L) {
if (y < 0.0) {
if (yoddint)
return -0.0L;
return -0.0;
return 0.0;
}
}
nflg = 0; /* flag = 1 if x<0 raised to integer power */
if (x <= 0.0L) {
if (x == 0.0L) {
if (x <= 0.0) {
if (x == 0.0) {
if (y < 0.0) {
if (signbit(x) && yoddint)
return -INFINITY;
@ -283,12 +283,12 @@ long double powl(long double x, long double y)
}
if (y > 0.0) {
if (signbit(x) && yoddint)
return -0.0L;
return -0.0;
return 0.0;
}
if (y == 0.0L)
return 1.0L; /* 0**0 */
return 0.0L; /* 0**y */
if (y == 0.0)
return 1.0; /* 0**0 */
return 0.0; /* 0**y */
}
if (iyflg == 0)
return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
@ -343,7 +343,7 @@ long double powl(long double x, long double y)
*/
z = x*x;
w = x * (z * __polevll(x, P, 3) / __p1evll(x, Q, 3));
w = w - ldexpl(z, -1); /* w - 0.5 * z */
w = w - 0.5*z;
/* Convert to base 2 logarithm:
* multiply by log2(e) = 1 + LOG2EA
@ -355,7 +355,8 @@ long double powl(long double x, long double y)
/* Compute exponent term of the base 2 logarithm. */
w = -i;
w = ldexpl(w, -LNXT); /* divide by NXT */
// TODO: use w * 0x1p-5;
w = scalbnl(w, -LNXT); /* divide by NXT */
w += e;
/* Now base 2 log of x is w + z. */
@ -380,7 +381,7 @@ long double powl(long double x, long double y)
H = Fb + Gb;
Ha = reducl(H);
w = ldexpl( Ga+Ha, LNXT );
w = scalbnl( Ga+Ha, LNXT );
/* Test the power of 2 for overflow */
if (w > MEXP)
@ -391,9 +392,9 @@ long double powl(long double x, long double y)
e = w;
Hb = H - Ha;
if (Hb > 0.0L) {
if (Hb > 0.0) {
e += 1;
Hb -= 1.0L/NXT; /*0.0625L;*/
Hb -= 1.0/NXT; /*0.0625L;*/
}
/* Now the product y * log2(x) = Hb + e/NXT.
@ -415,16 +416,16 @@ long double powl(long double x, long double y)
w = douba(e);
z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
z = z + w;
z = ldexpl(z, i); /* multiply by integer power of 2 */
z = scalbnl(z, i); /* multiply by integer power of 2 */
if (nflg) {
/* For negative x,
* find out if the integer exponent
* is odd or even.
*/
w = ldexpl(y, -1);
w = 0.5*y;
w = floorl(w);
w = ldexpl(w, 1);
w = 2.0*w;
if (w != y)
z = -z; /* odd exponent */
}
@ -438,9 +439,9 @@ static long double reducl(long double x)
{
long double t;
t = ldexpl(x, LNXT);
t = scalbnl(x, LNXT);
t = floorl(t);
t = ldexpl(t, -LNXT);
t = scalbnl(t, -LNXT);
return t;
}
@ -483,18 +484,18 @@ static long double powil(long double x, int nn)
long double s;
int n, e, sign, asign, lx;
if (x == 0.0L) {
if (x == 0.0) {
if (nn == 0)
return 1.0L;
return 1.0;
else if (nn < 0)
return LDBL_MAX;
return 0.0L;
return 0.0;
}
if (nn == 0)
return 1.0L;
return 1.0;
if (x < 0.0L) {
if (x < 0.0) {
asign = -1;
x = -x;
} else
@ -516,7 +517,7 @@ static long double powil(long double x, int nn)
e = (lx - 1)*n;
if ((e == 0) || (e > 64) || (e < -64)) {
s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L);
s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
s = (2.9142135623730950L * s - 0.5 + lx) * nn * LOGE2L;
} else {
s = LOGE2L * e;
}
@ -530,8 +531,8 @@ static long double powil(long double x, int nn)
* since roundoff error in 1.0/x will be amplified.
* The precise demarcation should be the gradual underflow threshold.
*/
if (s < -MAXLOGL+2.0L) {
x = 1.0L/x;
if (s < -MAXLOGL+2.0) {
x = 1.0/x;
sign = -sign;
}
@ -539,7 +540,7 @@ static long double powil(long double x, int nn)
if (n & 1)
y = x;
else {
y = 1.0L;
y = 1.0;
asign = 0;
}
@ -555,7 +556,7 @@ static long double powil(long double x, int nn)
if (asign)
y = -y; /* odd power of negative number */
if (sign < 0)
y = 1.0L/y;
y = 1.0/y;
return y;
}

12
src/math/remainder.c
View File

@ -20,8 +20,6 @@
#include "libm.h"
static const double zero = 0.0;
double remainder(double x, double p)
{
int32_t hx,hp;
@ -35,17 +33,15 @@ double remainder(double x, double p)
hx &= 0x7fffffff;
/* purge off exception values */
if ((hp|lp) == 0) /* p = 0 */
return (x*p)/(x*p);
if (hx >= 0x7ff00000 || /* x not finite */
if ((hp|lp) == 0 || /* p = 0 */
hx >= 0x7ff00000 || /* x not finite */
(hp >= 0x7ff00000 && (hp-0x7ff00000 | lp) != 0)) /* p is NaN */
// FIXME: why long double?
return ((long double)x*p)/((long double)x*p);
return (x*p)/(x*p);
if (hp <= 0x7fdfffff)
x = fmod(x, p+p); /* now x < 2p */
if (((hx-hp)|(lx-lp)) == 0)
return zero*x;
return 0.0*x;
x = fabs(x);
p = fabs(p);
if (hp < 0x00200000) {

9
src/math/remainderf.c
View File

@ -15,8 +15,6 @@
#include "libm.h"
static const float zero = 0.0;
float remainderf(float x, float p)
{
int32_t hx,hp;
@ -30,16 +28,13 @@ float remainderf(float x, float p)
hx &= 0x7fffffff;
/* purge off exception values */
if (hp == 0) /* p = 0 */
if (hp == 0 || hx >= 0x7f800000 || hp > 0x7f800000) /* p = 0, x not finite, p is NaN */
return (x*p)/(x*p);
if (hx >= 0x7f800000 || hp > 0x7f800000) /* x not finite, p is NaN */
// FIXME: why long double?
return ((long double)x*p)/((long double)x*p);
if (hp <= 0x7effffff)
x = fmodf(x, p + p); /* now x < 2p */
if (hx - hp == 0)
return zero*x;
return 0.0f*x;
x = fabsf(x);
p = fabsf(p);
if (hp < 0x01000000) {

2
src/math/remquol.c
View File

@ -48,7 +48,7 @@ typedef uint32_t manh_t;
#define MANL_SHIFT (LDBL_MANL_SIZE - 1)
static const long double Zero[] = {0.0L, -0.0L};
static const long double Zero[] = {0.0, -0.0};
/*
* Return the IEEE remainder and set *quo to the last n bits of the

2
src/math/rintl.c
View File

@ -79,7 +79,7 @@ long double rintl(long double x)
* If the result is +-0, then it must have the same sign as x, but
* the above calculation doesn't always give this. Fix up the sign.
*/
if (ex < BIAS && x == 0.0L)
if (ex < BIAS && x == 0.0)
return zero[sign];
return x;

2
src/math/scalbnl.c
View File

@ -29,7 +29,7 @@ long double scalbnl(long double x, int n)
return x * 0x1p-16382L;
}
}
scale.e = 1.0L;
scale.e = 1.0;
scale.bits.exp = 0x3fff + n;
return x * scale.e;
}

2
src/math/sincosl.c
View File

@ -9,8 +9,6 @@ void sincosl(long double x, long double *sin, long double *cos)
*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;

8
src/math/sinh.c
View File

@ -29,7 +29,7 @@
#include "libm.h"
static const double one = 1.0, huge = 1.0e307;
static const double huge = 1.0e307;
double sinh(double x)
{
@ -50,12 +50,12 @@ double sinh(double x)
if (ix < 0x40360000) { /* |x|<22 */
if (ix < 0x3e300000) /* |x|<2**-28 */
/* raise inexact, return x */
if (huge+x > one)
if (huge+x > 1.0)
return x;
t = expm1(fabs(x));
if (ix < 0x3ff00000)
return h*(2.0*t - t*t/(t+one));
return h*(t + t/(t+one));
return h*(2.0*t - t*t/(t+1.0));
return h*(t + t/(t+1.0));
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */

8
src/math/sinhf.c
View File

@ -15,7 +15,7 @@
#include "libm.h"
static const float one = 1.0, huge = 1.0e37;
static const float huge = 1.0e37;
float sinhf(float x)
{
@ -36,12 +36,12 @@ float sinhf(float x)
if (ix < 0x41100000) { /* |x|<9 */
if (ix < 0x39800000) /* |x|<2**-12 */
/* raise inexact, return x */
if (huge+x > one)
if (huge+x > 1.0f)
return x;
t = expm1f(fabsf(x));
if (ix < 0x3f800000)
return h*(2.0f*t - t*t/(t+one));
return h*(t + t/(t+one));
return h*(2.0f*t - t*t/(t+1.0f));
return h*(t + t/(t+1.0f));
}
/* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */

8
src/math/sinhl.c
View File

@ -35,7 +35,7 @@ long double sinhl(long double x)
return sinh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double one = 1.0, huge = 1.0e4931L;
static const long double huge = 1.0e4931L;
long double sinhl(long double x)
{
@ -55,12 +55,12 @@ long double sinhl(long double x)
/* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x| < 25 */
if (ix < 0x3fdf) /* |x|<2**-32 */
if (huge + x > one)
if (huge + x > 1.0)
return x;/* sinh(tiny) = tiny with inexact */
t = expm1l(fabsl(x));
if (ix < 0x3fff)
return h*(2.0*t - t*t/(t + one));
return h*(t + t/(t + one));
return h*(2.0*t - t*t/(t + 1.0));
return h*(t + t/(t + 1.0));
}
/* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */

2
src/math/sinl.c
View File

@ -34,8 +34,6 @@ long double sinl(long double x)
return sin(x);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__rem_pio2l.h"
long double sinl(long double x)
{
union IEEEl2bits z;

10
src/math/sqrt.c
View File

@ -78,7 +78,7 @@
#include "libm.h"
static const double one = 1.0, tiny = 1.0e-300;
static const double tiny = 1.0e-300;
double sqrt(double x)
{
@ -161,13 +161,13 @@ double sqrt(double x)
/* use floating add to find out rounding direction */
if ((ix0|ix1) != 0) {
z = one - tiny; /* raise inexact flag */
if (z >= one) {
z = one + tiny;
z = 1.0 - tiny; /* raise inexact flag */
if (z >= 1.0) {
z = 1.0 + tiny;
if (q1 == (uint32_t)0xffffffff) {
q1 = 0;
q++;
} else if (z > one) {
} else if (z > 1.0) {
if (q1 == (uint32_t)0xfffffffe)
q++;
q1 += 2;

10
src/math/sqrtf.c
View File

@ -15,7 +15,7 @@
#include "libm.h"
static const float one = 1.0, tiny = 1.0e-30;
static const float tiny = 1.0e-30;
float sqrtf(float x)
{
@ -68,10 +68,10 @@ float sqrtf(float x)
/* use floating add to find out rounding direction */
if (ix != 0) {
z = one - tiny; /* raise inexact flag */
if (z >= one) {
z = one + tiny;
if (z > one)
z = 1.0f - tiny; /* raise inexact flag */
if (z >= 1.0f) {
z = 1.0f + tiny;
if (z > 1.0f)
q += 2;
else
q += q & 1;

18
src/math/tanh.c
View File

@ -35,7 +35,7 @@
#include "libm.h"
static const double one = 1.0, two = 2.0, tiny = 1.0e-300, huge = 1.0e300;
static const double tiny = 1.0e-300, huge = 1.0e300;
double tanh(double x)
{
@ -48,26 +48,26 @@ double tanh(double x)
/* x is INF or NaN */
if (ix >= 0x7ff00000) {
if (jx >= 0)
return one/x + one; /* tanh(+-inf)=+-1 */
return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */
else
return one/x - one; /* tanh(NaN) = NaN */
return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */
}
if (ix < 0x40360000) { /* |x| < 22 */
if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* tanh(tiny) = tiny with inexact */
if (huge+x > one)
if (huge+x > 1.0f)
return x;
}
if (ix >= 0x3ff00000) { /* |x| >= 1 */
t = expm1(two*fabs(x));
z = one - two/(t+two);
t = expm1(2.0f*fabs(x));
z = 1.0f - 2.0f/(t+2.0f);
} else {
t = expm1(-two*fabs(x));
z= -t/(t+two);
t = expm1(-2.0f*fabs(x));
z= -t/(t+2.0f);
}
} else { /* |x| >= 22, return +-1 */
z = one - tiny; /* raise inexact */
z = 1.0f - tiny; /* raise inexact */
}
return jx >= 0 ? z : -z;
}

20
src/math/tanhf.c
View File

@ -15,7 +15,9 @@
#include "libm.h"
static const float one = 1.0, two = 2.0, tiny = 1.0e-30, huge = 1.0e30;
static const float
tiny = 1.0e-30,
huge = 1.0e30;
float tanhf(float x)
{
@ -28,26 +30,26 @@ float tanhf(float x)
/* x is INF or NaN */
if(ix >= 0x7f800000) {
if (jx >= 0)
return one/x + one; /* tanh(+-inf)=+-1 */
return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */
else
return one/x - one; /* tanh(NaN) = NaN */
return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */
}
if (ix < 0x41100000) { /* |x| < 9 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
/* tanh(tiny) = tiny with inexact */
if (huge+x > one)
if (huge+x > 1.0f)
return x;
}
if (ix >= 0x3f800000) { /* |x|>=1 */
t = expm1f(two*fabsf(x));
z = one - two/(t+two);
t = expm1f(2.0f*fabsf(x));
z = 1.0f - 2.0f/(t+2.0f);
} else {
t = expm1f(-two*fabsf(x));
z = -t/(t+two);
t = expm1f(-2.0f*fabsf(x));
z = -t/(t+2.0f);
}
} else { /* |x| >= 9, return +-1 */
z = one - tiny; /* raise inexact */
z = 1.0f - tiny; /* raise inexact */
}
return jx >= 0 ? z : -z;
}

18
src/math/tanhl.c
View File

@ -41,7 +41,7 @@ long double tanhl(long double x)
return tanh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double one=1.0, two=2.0, tiny = 1.0e-4900L;
static const long double tiny = 1.0e-4900L;
long double tanhl(long double x)
{
@ -57,8 +57,8 @@ long double tanhl(long double x)
if (ix == 0x7fff) {
/* for NaN it's not important which branch: tanhl(NaN) = NaN */
if (se & 0x8000)
return one/x-one; /* tanhl(-inf)= -1; */
return one/x+one; /* tanhl(+inf)= +1 */
return 1.0/x-1.0; /* tanhl(-inf)= -1; */
return 1.0/x+1.0; /* tanhl(+inf)= +1 */
}
/* |x| < 23 */
@ -66,17 +66,17 @@ long double tanhl(long double x)
if ((ix|jj0|jj1) == 0) /* x == +- 0 */
return x;
if (ix < 0x3fc8) /* |x| < 2**-55 */
return x*(one+tiny); /* tanh(small) = small */
return x*(1.0+tiny); /* tanh(small) = small */
if (ix >= 0x3fff) { /* |x| >= 1 */
t = expm1l(two*fabsl(x));
z = one - two/(t+two);
t = expm1l(2.0*fabsl(x));
z = 1.0 - 2.0/(t+2.0);
} else {
t = expm1l(-two*fabsl(x));
z = -t/(t+two);
t = expm1l(-2.0*fabsl(x));
z = -t/(t+2.0);
}
/* |x| > 23, return +-1 */
} else {
z = one - tiny; /* raise inexact flag */
z = 1.0 - tiny; /* raise inexact flag */
}
return se & 0x8000 ? -z : z;
}

2
src/math/tanl.c
View File

@ -38,8 +38,6 @@ long double tanl(long double x)
return tan(x);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__rem_pio2l.h"
long double tanl(long double x)
{
union IEEEl2bits z;

34
src/math/tgammal.c
View File

@ -183,18 +183,18 @@ static long double stirf(long double x)
{
long double y, w, v;
w = 1.0L/x;
w = 1.0/x;
/* For large x, use rational coefficients from the analytical expansion. */
if (x > 1024.0L)
if (x > 1024.0)
w = (((((6.97281375836585777429E-5L * w
+ 7.84039221720066627474E-4L) * w
- 2.29472093621399176955E-4L) * w
- 2.68132716049382716049E-3L) * w
+ 3.47222222222222222222E-3L) * w
+ 8.33333333333333333333E-2L) * w
+ 1.0L;
+ 1.0;
else
w = 1.0L + w * __polevll(w, STIR, 8);
w = 1.0 + w * __polevll(w, STIR, 8);
y = expl(x);
if (x > MAXSTIR) { /* Avoid overflow in pow() */
v = powl(x, 0.5L * x - 0.25L);
@ -219,10 +219,10 @@ long double tgammal(long double x)
if (x == -INFINITY)
return x - x;
q = fabsl(x);
if (q > 13.0L) {
if (q > 13.0) {
if (q > MAXGAML)
goto goverf;
if (x < 0.0L) {
if (x < 0.0) {
p = floorl(q);
if (p == q)
return (x - x) / (x - x);
@ -231,7 +231,7 @@ long double tgammal(long double x)
signgam = -1;
z = q - p;
if (z > 0.5L) {
p += 1.0L;
p += 1.0;
z = q - p;
}
z = q * sinl(PIL * z);
@ -247,25 +247,25 @@ goverf:
return signgam * z;
}
z = 1.0L;
while (x >= 3.0L) {
x -= 1.0L;
z = 1.0;
while (x >= 3.0) {
x -= 1.0;
z *= x;
}
while (x < -0.03125L) {
z /= x;
x += 1.0L;
x += 1.0;
}
if (x <= 0.03125L)
goto small;
while (x < 2.0L) {
while (x < 2.0) {
z /= x;
x += 1.0L;
x += 1.0;
}
if (x == 2.0L)
if (x == 2.0)
return z;
x -= 2.0L;
x -= 2.0;
p = __polevll(x, P, 7);
q = __polevll(x, Q, 8);
z = z * p / q;
@ -274,9 +274,9 @@ goverf:
return z;
small:
if (x == 0.0L)
if (x == 0.0)
return (x - x) / (x - x);
if (x < 0.0L) {
if (x < 0.0) {
x = -x;
q = z / (x * __polevll(x, SN, 8));
signgam = -1;