geemili
/
SDL
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Updated math code from the uClibc 0.9.33 release

Sam Lantinga 2017-11-04 15:53:19 -07:00
parent 34502143d9
commit 6cf065753c
18 changed files with 945 additions and 1198 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/libm/e_atan2.c
View File

@ -57,8 +57,8 @@ double attribute_hidden __ieee754_atan2(double y, double x)
ix = hx&0x7fffffff; ix = hx&0x7fffffff;
EXTRACT_WORDS(hy,ly,y); EXTRACT_WORDS(hy,ly,y);
iy = hy&0x7fffffff; iy = hy&0x7fffffff;
if(((ix|((lx|-(int32_t)lx)>>31))>0x7ff00000)|| if(((ix|((lx|-lx)>>31))>0x7ff00000)||
((iy|((ly|-(int32_t)ly)>>31))>0x7ff00000)) /* x or y is NaN */ ((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */
return x+y; return x+y;
if(((hx-0x3ff00000)|lx)==0) return atan(y); /* x=1.0 */ if(((hx-0x3ff00000)|lx)==0) return atan(y); /* x=1.0 */
m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
@ -81,8 +81,8 @@ double attribute_hidden __ieee754_atan2(double y, double x)
switch(m) { switch(m) {
case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
case 2: return 3.0*pi_o_4+tiny;/* atan(+INF,-INF) */ case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
case 3: return -3.0*pi_o_4-tiny;/* atan(-INF,-INF) */ case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
} }
} else { } else {
switch(m) { switch(m) {
@ -114,3 +114,21 @@ double attribute_hidden __ieee754_atan2(double y, double x)
return (z-pi_lo)-pi;/* atan(-,-) */ return (z-pi_lo)-pi;/* atan(-,-) */
} }
} }
/*
* wrapper atan2(y,x)
*/
#ifndef _IEEE_LIBM
double atan2(double y, double x)
{
double z = __ieee754_atan2(y, x);
if (_LIB_VERSION == _IEEE_ || isnan(x) || isnan(y))
return z;
if (x == 0.0 && y == 0.0)
return __kernel_standard(y,x,3); /* atan2(+-0,+-0) */
return z;
}
#else
strong_alias(__ieee754_atan2, atan2)
#endif
libm_hidden_def(atan2)

170
src/libm/e_log.c
View File

@ -1,4 +1,3 @@
/* @(#)e_log.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: e_log.c,v 1.8 1995/05/10 20:45:49 jtc Exp $";
#endif
/* __ieee754_log(x) /* __ieee754_log(x)
* Return the logrithm of x * Return the logrithm of x
* *
@ -69,99 +63,85 @@ static const char rcsid[] =
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
#ifdef __STDC__
static const double static const double
#else ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
static double ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
#endif two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
#ifdef __STDC__ static const double zero = 0.0;
static const double zero = 0.0;
#else
static double zero = 0.0;
#endif
#ifdef __STDC__ double attribute_hidden __ieee754_log(double x)
double attribute_hidden
__ieee754_log(double x)
#else
double attribute_hidden
__ieee754_log(x)
double x;
#endif
{ {
double hfsq, f, s, z, R, w, t1, t2, dk; double hfsq,f,s,z,R,w,t1,t2,dk;
int32_t k, hx, i, j; int32_t k,hx,i,j;
u_int32_t lx; u_int32_t lx;
EXTRACT_WORDS(hx, lx, x); EXTRACT_WORDS(hx,lx,x);
k = 0; k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */ if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx & 0x7fffffff) | lx) == 0) if (((hx&0x7fffffff)|lx)==0)
return -two54 / zero; /* log(+-0)=-inf */ return -two54/zero; /* log(+-0)=-inf */
if (hx < 0) if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
return (x - x) / zero; /* log(-#) = NaN */ k -= 54; x *= two54; /* subnormal number, scale up x */
k -= 54; GET_HIGH_WORD(hx,x);
x *= two54; /* subnormal number, scale up x */ }
GET_HIGH_WORD(hx, x); if (hx >= 0x7ff00000) return x+x;
} k += (hx>>20)-1023;
if (hx >= 0x7ff00000) hx &= 0x000fffff;
return x + x; i = (hx+0x95f64)&0x100000;
k += (hx >> 20) - 1023; SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */
hx &= 0x000fffff; k += (i>>20);
i = (hx + 0x95f64) & 0x100000; f = x-1.0;
SET_HIGH_WORD(x, hx | (i ^ 0x3ff00000)); /* normalize x or x/2 */ if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
k += (i >> 20); if(f==zero) {if(k==0) return zero; else {dk=(double)k;
f = x - 1.0; return dk*ln2_hi+dk*ln2_lo;}
if ((0x000fffff & (2 + hx)) < 3) { /* |f| < 2**-20 */ }
if (f == zero) { R = f*f*(0.5-0.33333333333333333*f);
if (k == 0) if(k==0) return f-R; else {dk=(double)k;
return zero; return dk*ln2_hi-((R-dk*ln2_lo)-f);}
else { }
dk = (double) k; s = f/(2.0+f);
return dk * ln2_hi + dk * ln2_lo; dk = (double)k;
} z = s*s;
} i = hx-0x6147a;
R = f * f * (0.5 - 0.33333333333333333 * f); w = z*z;
if (k == 0) j = 0x6b851-hx;
return f - R; t1= w*(Lg2+w*(Lg4+w*Lg6));
else { t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
dk = (double) k; i |= j;
return dk * ln2_hi - ((R - dk * ln2_lo) - f); R = t2+t1;
} if(i>0) {
} hfsq=0.5*f*f;
s = f / (2.0 + f); if(k==0) return f-(hfsq-s*(hfsq+R)); else
dk = (double) k; return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
z = s * s; } else {
i = hx - 0x6147a; if(k==0) return f-s*(f-R); else
w = z * z; return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
j = 0x6b851 - hx; }
t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
i |= j;
R = t2 + t1;
if (i > 0) {
hfsq = 0.5 * f * f;
if (k == 0)
return f - (hfsq - s * (hfsq + R));
else
return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) -
f);
} else {
if (k == 0)
return f - s * (f - R);
else
return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f);
}
} }
/*
* wrapper log(x)
*/
#ifndef _IEEE_LIBM
double log(double x)
{
double z = __ieee754_log(x);
if (_LIB_VERSION == _IEEE_ || isnan(x) || x > 0.0)
return z;
if (x == 0.0)
return __kernel_standard(x, x, 16); /* log(0) */
return __kernel_standard(x, x, 17); /* log(x<0) */
}
#else
strong_alias(__ieee754_log, log)
#endif
libm_hidden_def(log)

596
src/libm/e_pow.c
View File

@ -1,4 +1,3 @@
/* @(#)e_pow.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,10 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
#endif
/* __ieee754_pow(x,y) return x**y /* __ieee754_pow(x,y) return x**y
* *
* n * n
@ -26,25 +21,26 @@ static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
* 3. Return x**y = 2**n*exp(y'*log2) * 3. Return x**y = 2**n*exp(y'*log2)
* *
* Special cases: * Special cases:
* 1. (anything) ** 0 is 1 * 1. +-1 ** anything is 1.0
* 2. (anything) ** 1 is itself * 2. +-1 ** +-INF is 1.0
* 3. (anything) ** NAN is NAN * 3. (anything) ** 0 is 1
* 4. NAN ** (anything except 0) is NAN * 4. (anything) ** 1 is itself
* 5. +-(|x| > 1) ** +INF is +INF * 5. (anything) ** NAN is NAN
* 6. +-(|x| > 1) ** -INF is +0 * 6. NAN ** (anything except 0) is NAN
* 7. +-(|x| < 1) ** +INF is +0 * 7. +-(|x| > 1) ** +INF is +INF
* 8. +-(|x| < 1) ** -INF is +INF * 8. +-(|x| > 1) ** -INF is +0
* 9. +-1 ** +-INF is NAN * 9. +-(|x| < 1) ** +INF is +0
* 10. +0 ** (+anything except 0, NAN) is +0 * 10 +-(|x| < 1) ** -INF is +INF
* 11. -0 ** (+anything except 0, NAN, odd integer) is +0 * 11. +0 ** (+anything except 0, NAN) is +0
* 12. +0 ** (-anything except 0, NAN) is +INF * 12. -0 ** (+anything except 0, NAN, odd integer) is +0
* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF * 13. +0 ** (-anything except 0, NAN) is +INF
* 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) * 14. -0 ** (-anything except 0, NAN, odd integer) is +INF
* 15. +INF ** (+anything except 0,NAN) is +INF * 15. -0 ** (odd integer) = -( +0 ** (odd integer) )
* 16. +INF ** (-anything except 0,NAN) is +0 * 16. +INF ** (+anything except 0,NAN) is +INF
* 17. -INF ** (anything) = -0 ** (-anything) * 17. +INF ** (-anything except 0,NAN) is +0
* 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) * 18. -INF ** (anything) = -0 ** (-anything)
* 19. (-anything except 0 and inf) ** (non-integer) is NAN * 19. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
* 20. (-anything except 0 and inf) ** (non-integer) is NAN
* *
* Accuracy: * Accuracy:
* pow(x,y) returns x**y nearly rounded. In particular * pow(x,y) returns x**y nearly rounded. In particular
@ -62,281 +58,281 @@ static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
libm_hidden_proto(scalbn) static const double
libm_hidden_proto(fabs) bp[] = {1.0, 1.5,},
#ifdef __STDC__ dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
static const double 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,
/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
double attribute_hidden __ieee754_pow(double x, double y)
{
double z,ax,z_h,z_l,p_h,p_l;
double y1,t1,t2,r,s,t,u,v,w;
int32_t i,j,k,yisint,n;
int32_t hx,hy,ix,iy;
u_int32_t lx,ly;
EXTRACT_WORDS(hx,lx,x);
/* x==1: 1**y = 1 (even if y is NaN) */
if (hx==0x3ff00000 && lx==0) {
return x;
}
ix = hx&0x7fffffff;
EXTRACT_WORDS(hy,ly,y);
iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */
if((iy|ly)==0) return one;
/* +-NaN return x+y */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
return x+y;
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 1 ... y is an odd int
* yisint = 2 ... y is an even int
*/
yisint = 0;
if(hx<0) {
if(iy>=0x43400000) yisint = 2; /* even integer y */
else if(iy>=0x3ff00000) {
k = (iy>>20)-0x3ff; /* exponent */
if(k>20) {
j = ly>>(52-k);
if((j<<(52-k))==ly) yisint = 2-(j&1);
} else if(ly==0) {
j = iy>>(20-k);
if((j<<(20-k))==iy) yisint = 2-(j&1);
}
}
}
/* special value of y */
if(ly==0) {
if (iy==0x7ff00000) { /* y is +-inf */
if (((ix-0x3ff00000)|lx)==0)
return one; /* +-1**+-inf is 1 (yes, weird rule) */
if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
return (hy>=0) ? y : zero;
/* (|x|<1)**-,+inf = inf,0 */
return (hy<0) ? -y : zero;
}
if(iy==0x3ff00000) { /* y is +-1 */
if(hy<0) return one/x; else return x;
}
if(hy==0x40000000) return x*x; /* y is 2 */
if(hy==0x3fe00000) { /* y is 0.5 */
if(hx>=0) /* x >= +0 */
return __ieee754_sqrt(x);
}
}
ax = fabs(x);
/* special value of x */
if(lx==0) {
if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
z = ax; /*x is +-0,+-inf,+-1*/
if(hy<0) z = one/z; /* z = (1/|x|) */
if(hx<0) {
if(((ix-0x3ff00000)|yisint)==0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
} else if(yisint==1)
z = -z; /* (x<0)**odd = -(|x|**odd) */
}
return z;
}
}
/* (x<0)**(non-int) is NaN */
if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
/* |y| is huge */
if(iy>0x41e00000) { /* if |y| > 2**31 */
if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
}
/* over/underflow if x is not close to one */
if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
if(ix>0x3ff00000) return (hy>0)? huge*huge: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 = x-1; /* 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;
t1 = u+v;
SET_LOW_WORD(t1,0);
t2 = v-(t1-u);
} else {
double s2,s_h,s_l,t_h,t_l;
n = 0;
/* take care subnormal number */
if(ix<0x00100000)
{ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
n += ((ix)>>20)-0x3ff;
j = ix&0x000fffff;
/* determine interval */
ix = j|0x3ff00000; /* normalize ix */
if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
else {k=0;n+=1;ix -= 0x00100000;}
SET_HIGH_WORD(ax,ix);
/* 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]);
s = u*v;
s_h = s;
SET_LOW_WORD(s_h,0);
/* t_h=ax+bp[k] High */
t_h = zero;
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);
/* compute log(ax) */
s2 = s*s;
r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
r += s_l*(s_h+s);
s2 = s_h*s_h;
t_h = 3.0+s2+r;
SET_LOW_WORD(t_h,0);
t_l = r-((t_h-3.0)-s2);
/* u+v = s*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*s;
/* 2/(3log2)*(s+...) */
p_h = u+v;
SET_LOW_WORD(p_h,0);
p_l = v-(p_h-u);
z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l*p_h+p_l*cp+dp_l[k];
/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
t = (double)n;
t1 = (((z_h+z_l)+dp_h[k])+t);
SET_LOW_WORD(t1,0);
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
}
s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
s = -one;/* (-ve)**(odd int) */
/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
y1 = y;
SET_LOW_WORD(y1,0);
p_l = (y-y1)*t1+y*t2;
p_h = y1*t1;
z = p_l+p_h;
EXTRACT_WORDS(j,i,z);
if (j>=0x40900000) { /* z >= 1024 */
if(((j-0x40900000)|i)!=0) /* if z > 1024 */
return s*huge*huge; /* overflow */
else {
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
}
} else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
return s*tiny*tiny; /* underflow */
else {
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
}
}
/*
* compute 2**(p_h+p_l)
*/
i = j&0x7fffffff;
k = (i>>20)-0x3ff;
n = 0;
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;
SET_HIGH_WORD(t,n&~(0x000fffff>>k));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if(j<0) n = -n;
p_h -= t;
}
t = p_l+p_h;
SET_LOW_WORD(t,0);
u = t*lg2_h;
v = (p_l-(t-p_h))*lg2+t*lg2_l;
z = u+v;
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);
GET_HIGH_WORD(j,z);
j += (n<<20);
if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
else SET_HIGH_WORD(z,j);
return s*z;
}
/*
* wrapper pow(x,y) return x**y
*/
#ifndef _IEEE_LIBM
double pow(double x, double y)
{
double z = __ieee754_pow(x, y);
if (_LIB_VERSION == _IEEE_|| isnan(y))
return z;
if (isnan(x)) {
if (y == 0.0)
return __kernel_standard(x, y, 42); /* pow(NaN,0.0) */
return z;
}
if (x == 0.0) {
if (y == 0.0)
return __kernel_standard(x, y, 20); /* pow(0.0,0.0) */
if (isfinite(y) && y < 0.0)
return __kernel_standard(x,y,23); /* pow(0.0,negative) */
return z;
}
if (!isfinite(z)) {
if (isfinite(x) && isfinite(y)) {
if (isnan(z))
return __kernel_standard(x, y, 24); /* pow neg**non-int */
return __kernel_standard(x, y, 21); /* pow overflow */
}
}
if (z == 0.0 && isfinite(x) && isfinite(y))
return __kernel_standard(x, y, 22); /* pow underflow */
return z;
}
#else #else
static double strong_alias(__ieee754_pow, pow)
#endif #endif
bp[] = { 1.0, 1.5, }, dp_h[] = { libm_hidden_def(pow)
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_val = 1.0e300, tiny = 1.0e-300,
/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h */
ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2 */
ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail */
#ifdef __STDC__
double attribute_hidden __ieee754_pow(double x, double y)
#else
double attribute_hidden __ieee754_pow(x, y)
double x, y;
#endif
{
double z, ax, z_h, z_l, p_h, p_l;
double y1, t1, t2, r, s, t, u, v, w;
int32_t i, j, k, yisint, n;
int32_t hx, hy, ix, iy;
u_int32_t lx, ly;
EXTRACT_WORDS(hx, lx, x);
EXTRACT_WORDS(hy, ly, y);
ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
/* y==zero: x**0 = 1 */
if ((iy | ly) == 0)
return one;
/* +-NaN return x+y */
if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) ||
iy > 0x7ff00000 || ((iy == 0x7ff00000) && (ly != 0)))
return x + y;
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 1 ... y is an odd int
* yisint = 2 ... y is an even int
*/
yisint = 0;
if (hx < 0) {
if (iy >= 0x43400000)
yisint = 2; /* even integer y */
else if (iy >= 0x3ff00000) {
k = (iy >> 20) - 0x3ff; /* exponent */
if (k > 20) {
j = ly >> (52 - k);
if ((j << (52 - k)) == ly)
yisint = 2 - (j & 1);
} else if (ly == 0) {
j = iy >> (20 - k);
if ((j << (20 - k)) == iy)
yisint = 2 - (j & 1);
}
}
}
/* special value of y */
if (ly == 0) {
if (iy == 0x7ff00000) { /* y is +-inf */
if (((ix - 0x3ff00000) | lx) == 0)
return y - y; /* inf**+-1 is NaN */
else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
return (hy >= 0) ? y : zero;
else /* (|x|<1)**-,+inf = inf,0 */
return (hy < 0) ? -y : zero;
}
if (iy == 0x3ff00000) { /* y is +-1 */
if (hy < 0)
return one / x;
else
return x;
}
if (hy == 0x40000000)
return x * x; /* y is 2 */
if (hy == 0x3fe00000) { /* y is 0.5 */
if (hx >= 0) /* x >= +0 */
return __ieee754_sqrt(x);
}
}
ax = fabs(x);
/* special value of x */
if (lx == 0) {
if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) {
z = ax; /* x is +-0,+-inf,+-1 */
if (hy < 0)
z = one / z; /* z = (1/|x|) */
if (hx < 0) {
if (((ix - 0x3ff00000) | yisint) == 0) {
z = (z - z) / (z - z); /* (-1)**non-int is NaN */
} else if (yisint == 1)
z = -z; /* (x<0)**odd = -(|x|**odd) */
}
return z;
}
}
/* (x<0)**(non-int) is NaN */
if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
return (x - x) / (x - x);
/* |y| is huge */
if (iy > 0x41e00000) { /* if |y| > 2**31 */
if (iy > 0x43f00000) { /* if |y| > 2**64, must o/uflow */
if (ix <= 0x3fefffff)
return (hy < 0) ? huge_val * huge_val : tiny * tiny;
if (ix >= 0x3ff00000)
return (hy > 0) ? huge_val * huge_val : tiny * tiny;
}
/* over/underflow if x is not close to one */
if (ix < 0x3fefffff)
return (hy < 0) ? huge_val * huge_val : tiny * tiny;
if (ix > 0x3ff00000)
return (hy > 0) ? huge_val * huge_val : 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 = x - 1; /* 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;
t1 = u + v;
SET_LOW_WORD(t1, 0);
t2 = v - (t1 - u);
} else {
double s2, s_h, s_l, t_h, t_l;
n = 0;
/* take care subnormal number */
if (ix < 0x00100000) {
ax *= two53;
n -= 53;
GET_HIGH_WORD(ix, ax);
}
n += ((ix) >> 20) - 0x3ff;
j = ix & 0x000fffff;
/* determine interval */
ix = j | 0x3ff00000; /* normalize ix */
if (j <= 0x3988E)
k = 0; /* |x|<sqrt(3/2) */
else if (j < 0xBB67A)
k = 1; /* |x|<sqrt(3) */
else {
k = 0;
n += 1;
ix -= 0x00100000;
}
SET_HIGH_WORD(ax, ix);
/* 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]);
s = u * v;
s_h = s;
SET_LOW_WORD(s_h, 0);
/* t_h=ax+bp[k] High */
t_h = zero;
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);
/* compute log(ax) */
s2 = s * s;
r = s2 * s2 * (L1 +
s2 * (L2 +
s2 * (L3 +
s2 * (L4 + s2 * (L5 + s2 * L6)))));
r += s_l * (s_h + s);
s2 = s_h * s_h;
t_h = 3.0 + s2 + r;
SET_LOW_WORD(t_h, 0);
t_l = r - ((t_h - 3.0) - s2);
/* u+v = s*(1+...) */
u = s_h * t_h;
v = s_l * t_h + t_l * s;
/* 2/(3log2)*(s+...) */
p_h = u + v;
SET_LOW_WORD(p_h, 0);
p_l = v - (p_h - u);
z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l * p_h + p_l * cp + dp_l[k];
/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
t = (double) n;
t1 = (((z_h + z_l) + dp_h[k]) + t);
SET_LOW_WORD(t1, 0);
t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
}
s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
s = -one; /* (-ve)**(odd int) */
/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
y1 = y;
SET_LOW_WORD(y1, 0);
p_l = (y - y1) * t1 + y * t2;
p_h = y1 * t1;
z = p_l + p_h;
EXTRACT_WORDS(j, i, z);
if (j >= 0x40900000) { /* z >= 1024 */
if (((j - 0x40900000) | i) != 0) /* if z > 1024 */
return s * huge_val * huge_val; /* overflow */
else {
if (p_l + ovt > z - p_h)
return s * huge_val * huge_val; /* overflow */
}
} else if ((j & 0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */
if (((j - 0xc090cc00) | i) != 0) /* z < -1075 */
return s * tiny * tiny; /* underflow */
else {
if (p_l <= z - p_h)
return s * tiny * tiny; /* underflow */
}
}
/*
* compute 2**(p_h+p_l)
*/
i = j & 0x7fffffff;
k = (i >> 20) - 0x3ff;
n = 0;
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;
SET_HIGH_WORD(t, n & ~(0x000fffff >> k));
n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
if (j < 0)
n = -n;
p_h -= t;
}
t = p_l + p_h;
SET_LOW_WORD(t, 0);
u = t * lg2_h;
v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
z = u + v;
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);
GET_HIGH_WORD(j, z);
j += (n << 20);
if ((j >> 20) <= 0)
z = scalbn(z, n); /* subnormal output */
else
SET_HIGH_WORD(z, j);
return s * z;
}

277
src/libm/e_rem_pio2.c
View File

@ -1,4 +1,3 @@
/* @(#)e_rem_pio2.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $";
#endif
/* __ieee754_rem_pio2(x,y) /* __ieee754_rem_pio2(x,y)
* *
* return the remainder of x rem pi/2 in y[0]+y[1] * return the remainder of x rem pi/2 in y[0]+y[1]
@ -24,42 +18,30 @@ static const char rcsid[] =
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
#include "SDL_assert.h"
libm_hidden_proto(fabs)
/* /*
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
*/ */
#ifdef __STDC__ static const int32_t two_over_pi[] = {
static const int32_t two_over_pi[] = { 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
#else 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
static int32_t two_over_pi[] = { 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
#endif 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, };
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
};
#ifdef __STDC__
static const int32_t npio2_hw[] = { static const int32_t npio2_hw[] = {
#else 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
static int32_t npio2_hw[] = { 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
#endif 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 0x404858EB, 0x404921FB,
0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
0x404858EB, 0x404921FB,
}; };
/* /*
@ -72,137 +54,108 @@ static int32_t npio2_hw[] = {
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/ */
#ifdef __STDC__
static const double static const double
#else zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
static double half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
#endif two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
#ifdef __STDC__ int32_t attribute_hidden __ieee754_rem_pio2(double x, double *y)
int32_t attribute_hidden
__ieee754_rem_pio2(double x, double *y)
#else
int32_t attribute_hidden
__ieee754_rem_pio2(x, y)
double x, y[];
#endif
{ {
double z = 0.0, w, t, r, fn; double z=0.0,w,t,r,fn;
double tx[3]; double tx[3];
int32_t e0, i, j, nx, n, ix, hx; int32_t e0,i,j,nx,n,ix,hx;
u_int32_t low; u_int32_t low;
GET_HIGH_WORD(hx, x); /* high word of x */ GET_HIGH_WORD(hx,x); /* high word of x */
ix = hx & 0x7fffffff; ix = hx&0x7fffffff;
if (ix <= 0x3fe921fb) { /* |x| ~<= pi/4 , no need for reduction */ if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
y[0] = x; {y[0] = x; y[1] = 0; return 0;}
y[1] = 0; if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */
return 0; if(hx>0) {
} z = x - pio2_1;
if (ix < 0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
if (hx > 0) { y[0] = z - pio2_1t;
z = x - pio2_1; y[1] = (z-y[0])-pio2_1t;
if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */ } else { /* near pi/2, use 33+33+53 bit pi */
y[0] = z - pio2_1t; z -= pio2_2;
y[1] = (z - y[0]) - pio2_1t; y[0] = z - pio2_2t;
} else { /* near pi/2, use 33+33+53 bit pi */ y[1] = (z-y[0])-pio2_2t;
z -= pio2_2; }
y[0] = z - pio2_2t; return 1;
y[1] = (z - y[0]) - pio2_2t; } else { /* negative x */
} z = x + pio2_1;
return 1; if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
} else { /* negative x */ y[0] = z + pio2_1t;
z = x + pio2_1; y[1] = (z-y[0])+pio2_1t;
if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */ } else { /* near pi/2, use 33+33+53 bit pi */
y[0] = z + pio2_1t; z += pio2_2;
y[1] = (z - y[0]) + pio2_1t; y[0] = z + pio2_2t;
} else { /* near pi/2, use 33+33+53 bit pi */ y[1] = (z-y[0])+pio2_2t;
z += pio2_2; }
y[0] = z + pio2_2t; return -1;
y[1] = (z - y[0]) + pio2_2t; }
} }
return -1; if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
} t = fabs(x);
} n = (int32_t) (t*invpio2+half);
if (ix <= 0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ fn = (double)n;
t = fabs(x); r = t-fn*pio2_1;
n = (int32_t) (t * invpio2 + half); w = fn*pio2_1t; /* 1st round good to 85 bit */
fn = (double) n; if(n<32&&ix!=npio2_hw[n-1]) {
r = t - fn * pio2_1; y[0] = r-w; /* quick check no cancellation */
w = fn * pio2_1t; /* 1st round good to 85 bit */ } else {
if (n < 32 && ix != npio2_hw[n - 1]) { u_int32_t high;
y[0] = r - w; /* quick check no cancellation */ j = ix>>20;
} else { y[0] = r-w;
u_int32_t high; GET_HIGH_WORD(high,y[0]);
j = ix >> 20; i = j-((high>>20)&0x7ff);
y[0] = r - w; if(i>16) { /* 2nd iteration needed, good to 118 */
GET_HIGH_WORD(high, y[0]); t = r;
i = j - ((high >> 20) & 0x7ff); w = fn*pio2_2;
if (i > 16) { /* 2nd iteration needed, good to 118 */ r = t-w;
t = r; w = fn*pio2_2t-((t-r)-w);
w = fn * pio2_2; y[0] = r-w;
r = t - w; GET_HIGH_WORD(high,y[0]);
w = fn * pio2_2t - ((t - r) - w); i = j-((high>>20)&0x7ff);
y[0] = r - w; if(i>49) { /* 3rd iteration need, 151 bits acc */
GET_HIGH_WORD(high, y[0]); t = r; /* will cover all possible cases */
i = j - ((high >> 20) & 0x7ff); w = fn*pio2_3;
if (i > 49) { /* 3rd iteration need, 151 bits acc */ r = t-w;
t = r; /* will cover all possible cases */ w = fn*pio2_3t-((t-r)-w);
w = fn * pio2_3; y[0] = r-w;
r = t - w; }
w = fn * pio2_3t - ((t - r) - w); }
y[0] = r - w; }
} y[1] = (r-y[0])-w;
} if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
} else return n;
y[1] = (r - y[0]) - w; }
if (hx < 0) {
y[0] = -y[0];
y[1] = -y[1];
return -n;
} else
return n;
}
/* /*
* all other (large) arguments * all other (large) arguments
*/ */
if (ix >= 0x7ff00000) { /* x is inf or NaN */ if(ix>=0x7ff00000) { /* x is inf or NaN */
y[0] = y[1] = x - x; y[0]=y[1]=x-x; return 0;
return 0; }
}
/* set z = scalbn(|x|,ilogb(x)-23) */ /* set z = scalbn(|x|,ilogb(x)-23) */
GET_LOW_WORD(low, x); GET_LOW_WORD(low,x);
SET_LOW_WORD(z, low); SET_LOW_WORD(z,low);
e0 = (ix >> 20) - 1046; /* e0 = ilogb(z)-23; */ e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */
SET_HIGH_WORD(z, ix - ((int32_t) (e0 << 20))); SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20)));
for (i = 0; i < 2; i++) { for(i=0;i<2;i++) {
tx[i] = (double) ((int32_t) (z)); tx[i] = (double)((int32_t)(z));
z = (z - tx[i]) * two24; z = (z-tx[i])*two24;
} }
tx[2] = z; tx[2] = z;
nx = 3; nx = 3;
while(tx[nx-1]==zero) nx--; /* skip zero term */
/* If this assertion ever fires, here's the static analysis that warned about it: n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
http://buildbot.libsdl.org/sdl-static-analysis/sdl-macosx-static-analysis/sdl-macosx-static-analysis-1101/report-8c6ccb.html#EndPath */ if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
SDL_assert((tx[0] != zero) || (tx[1] != zero) || (tx[2] != zero)); return n;
while (nx && tx[nx - 1] == zero)
nx--; /* skip zero term */
n = __kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi);
if (hx < 0) {
y[0] = -y[0];
y[1] = -y[1];
return -n;
}
return n;
} }

209
src/libm/e_sqrt.c
View File

@ -1,4 +1,3 @@
/* @(#)e_sqrt.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
#endif
/* __ieee754_sqrt(x) /* __ieee754_sqrt(x)
* Return correctly rounded sqrt. * Return correctly rounded sqrt.
* ------------------------------------------ * ------------------------------------------
@ -88,124 +82,123 @@ static const char rcsid[] =
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
#ifdef __STDC__
static const double one = 1.0, tiny = 1.0e-300; static const double one = 1.0, tiny = 1.0e-300;
#else
static double one = 1.0, tiny = 1.0e-300;
#endif
#ifdef __STDC__ double attribute_hidden __ieee754_sqrt(double x)
double attribute_hidden
__ieee754_sqrt(double x)
#else
double attribute_hidden
__ieee754_sqrt(x)
double x;
#endif
{ {
double z; double z;
int32_t sign = (int) 0x80000000; int32_t sign = (int)0x80000000;
int32_t ix0, s0, q, m, t, i; int32_t ix0,s0,q,m,t,i;
u_int32_t r, t1, s1, ix1, q1; u_int32_t r,t1,s1,ix1,q1;
EXTRACT_WORDS(ix0, ix1, x); EXTRACT_WORDS(ix0,ix1,x);
/* take care of Inf and NaN */ /* take care of Inf and NaN */
if ((ix0 & 0x7ff00000) == 0x7ff00000) { if((ix0&0x7ff00000)==0x7ff00000) {
return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
sqrt(-inf)=sNaN */ sqrt(-inf)=sNaN */
} }
/* take care of zero */ /* take care of zero */
if (ix0 <= 0) { if(ix0<=0) {
if (((ix0 & (~sign)) | ix1) == 0) if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
return x; /* sqrt(+-0) = +-0 */ else if(ix0<0)
else if (ix0 < 0) return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
return (x - x) / (x - x); /* sqrt(-ve) = sNaN */ }
}
/* normalize x */ /* normalize x */
m = (ix0 >> 20); m = (ix0>>20);
if (m == 0) { /* subnormal x */ if(m==0) { /* subnormal x */
while (ix0 == 0) { while(ix0==0) {
m -= 21; m -= 21;
ix0 |= (ix1 >> 11); ix0 |= (ix1>>11); ix1 <<= 21;
ix1 <<= 21; }
} for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
for (i = 0; (ix0 & 0x00100000) == 0; i++) m -= i-1;
ix0 <<= 1; ix0 |= (ix1>>(32-i));
m -= i - 1; ix1 <<= i;
ix0 |= (ix1 >> (32 - i)); }
ix1 <<= i; m -= 1023; /* unbias exponent */
} ix0 = (ix0&0x000fffff)|0x00100000;
m -= 1023; /* unbias exponent */ if(m&1){ /* odd m, double x to make it even */
ix0 = (ix0 & 0x000fffff) | 0x00100000; ix0 += ix0 + ((ix1&sign)>>31);
if (m & 1) { /* odd m, double x to make it even */ ix1 += ix1;
ix0 += ix0 + ((ix1 & sign) >> 31); }
ix1 += ix1; m >>= 1; /* m = [m/2] */
}
m >>= 1; /* m = [m/2] */
/* generate sqrt(x) bit by bit */ /* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1 & sign) >> 31); ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1; ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */ q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x00200000; /* r = moving bit from right to left */ r = 0x00200000; /* r = moving bit from right to left */
while (r != 0) { while(r!=0) {
t = s0 + r; t = s0+r;
if (t <= ix0) { if(t<=ix0) {
s0 = t + r; s0 = t+r;
ix0 -= t; ix0 -= t;
q += r; q += r;
} }
ix0 += ix0 + ((ix1 & sign) >> 31); ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1; ix1 += ix1;
r >>= 1; r>>=1;
} }
r = sign; r = sign;
while (r != 0) { while(r!=0) {
t1 = s1 + r; t1 = s1+r;
t = s0; t = s0;
if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) { if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
s1 = t1 + r; s1 = t1+r;
if (((t1 & sign) == sign) && (s1 & sign) == 0) if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
s0 += 1; ix0 -= t;
ix0 -= t; if (ix1 < t1) ix0 -= 1;
if (ix1 < t1) ix1 -= t1;
ix0 -= 1; q1 += r;
ix1 -= t1; }
q1 += r; ix0 += ix0 + ((ix1&sign)>>31);
} ix1 += ix1;
ix0 += ix0 + ((ix1 & sign) >> 31); r>>=1;
ix1 += ix1; }
r >>= 1;
}
/* use floating add to find out rounding direction */ /* use floating add to find out rounding direction */
if ((ix0 | ix1) != 0) { if((ix0|ix1)!=0) {
z = one - tiny; /* trigger inexact flag */ z = one-tiny; /* trigger inexact flag */
if (z >= one) { if (z>=one) {
z = one + tiny; z = one+tiny;
if (q1 == (u_int32_t) 0xffffffff) { if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
q1 = 0; else if (z>one) {
q += 1; if (q1==(u_int32_t)0xfffffffe) q+=1;
} else if (z > one) { q1+=2;
if (q1 == (u_int32_t) 0xfffffffe) } else
q += 1; q1 += (q1&1);
q1 += 2; }
} else }
q1 += (q1 & 1); ix0 = (q>>1)+0x3fe00000;
} ix1 = q1>>1;
} if ((q&1)==1) ix1 |= sign;
ix0 = (q >> 1) + 0x3fe00000; ix0 += (m <<20);
ix1 = q1 >> 1; INSERT_WORDS(z,ix0,ix1);
if ((q & 1) == 1) return z;
ix1 |= sign;
ix0 += (m << 20);
INSERT_WORDS(z, ix0, ix1);
return z;
} }
/*
* wrapper sqrt(x)
*/
#ifndef _IEEE_LIBM
double sqrt(double x)
{
double z = __ieee754_sqrt(x);
if (_LIB_VERSION == _IEEE_ || isnan(x))
return z;
if (x < 0.0)
return __kernel_standard(x, x, 26); /* sqrt(negative) */
return z;
}
#else
strong_alias(__ieee754_sqrt, sqrt)
#endif
libm_hidden_def(sqrt)
/* /*
Other methods (use floating-point arithmetic) Other methods (use floating-point arithmetic)
------------- -------------

76
src/libm/k_cos.c
View File

@ -1,4 +1,3 @@
/* @(#)k_cos.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: k_cos.c,v 1.8 1995/05/10 20:46:22 jtc Exp $";
#endif
/* /*
* __kernel_cos( x, y ) * __kernel_cos( x, y )
* kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
@ -53,48 +47,36 @@ static const char rcsid[] =
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
#ifdef __STDC__
static const double static const double
#else one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
static double C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
#endif C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
#ifdef __STDC__ double attribute_hidden __kernel_cos(double x, double y)
double attribute_hidden
__kernel_cos(double x, double y)
#else
double attribute_hidden
__kernel_cos(x, y)
double x, y;
#endif
{ {
double a, hz, z, r, qx; double a,hz,z,r,qx;
int32_t ix; int32_t ix;
GET_HIGH_WORD(ix, x); GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff; /* ix = |x|'s high word */ ix &= 0x7fffffff; /* ix = |x|'s high word*/
if (ix < 0x3e400000) { /* if x < 2**27 */ if(ix<0x3e400000) { /* if x < 2**27 */
if (((int) x) == 0) if(((int)x)==0) return one; /* generate inexact */
return one; /* generate inexact */ }
} z = x*x;
z = x * x; r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6))))); if(ix < 0x3FD33333) /* if |x| < 0.3 */
if (ix < 0x3FD33333) /* if |x| < 0.3 */ return one - (0.5*z - (z*r - x*y));
return one - (0.5 * z - (z * r - x * y)); else {
else { if(ix > 0x3fe90000) { /* x > 0.78125 */
if (ix > 0x3fe90000) { /* x > 0.78125 */ qx = 0.28125;
qx = 0.28125; } else {
} else { INSERT_WORDS(qx,ix-0x00200000,0); /* x/4 */
INSERT_WORDS(qx, ix - 0x00200000, 0); /* x/4 */ }
} hz = 0.5*z-qx;
hz = 0.5 * z - qx; a = one-qx;
a = one - qx; return a - (hz - (z*r-x*y));
return a - (hz - (z * r - x * y)); }
}
} }

344
src/libm/k_rem_pio2.c
View File

@ -1,4 +1,3 @@
/* @(#)k_rem_pio2.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $";
#endif
/* /*
* __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
* double x[],y[]; int e0,nx,prec; int ipio2[]; * double x[],y[]; int e0,nx,prec; int ipio2[];
@ -134,235 +128,171 @@ static const char rcsid[] =
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
#include "SDL_assert.h" static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
libm_hidden_proto(scalbn)
libm_hidden_proto(floor)
#ifdef __STDC__
static const int init_jk[] = { 2, 3, 4, 6 }; /* initial value for jk */
#else
static int init_jk[] = { 2, 3, 4, 6 };
#endif
#ifdef __STDC__
static const double PIo2[] = { static const double PIo2[] = {
#else 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
static double PIo2[] = { 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
#endif 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
}; };
#ifdef __STDC__
static const double static const double
#else zero = 0.0,
static double one = 1.0,
#endif two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
#ifdef __STDC__ int attribute_hidden __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
int attribute_hidden
__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec,
const int32_t * ipio2)
#else
int attribute_hidden
__kernel_rem_pio2(x, y, e0, nx, prec, ipio2)
double x[], y[];
int e0, nx, prec;
int32_t ipio2[];
#endif
{ {
int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih; int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
double z, fw, f[20], fq[20], q[20]; double z,fw,f[20],fq[20],q[20];
/* initialize jk */ /* initialize jk*/
SDL_assert((prec >= 0) && (prec < SDL_arraysize(init_jk))); jk = init_jk[prec];
jk = init_jk[prec]; jp = jk;
SDL_assert((jk >= 2) && (jk <= 6));
jp = jk;
/* determine jx,jv,q0, note that 3>q0 */ /* determine jx,jv,q0, note that 3>q0 */
SDL_assert(nx > 0); jx = nx-1;
jx = nx - 1; jv = (e0-3)/24; if(jv<0) jv=0;
jv = (e0 - 3) / 24; q0 = e0-24*(jv+1);
if (jv < 0)
jv = 0;
q0 = e0 - 24 * (jv + 1);
/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
j = jv - jx; j = jv-jx; m = jx+jk;
m = jx + jk; for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
for (i = 0; i <= m; i++, j++)
f[i] = (j < 0) ? zero : (double) ipio2[j];
/* compute q[0],q[1],...q[jk] */ /* compute q[0],q[1],...q[jk] */
for (i = 0; i <= jk; i++) { for (i=0;i<=jk;i++) {
for (j = 0, fw = 0.0; j <= jx; j++) { for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
const int32_t idx = jx + i - j; }
SDL_assert(idx >= 0);
SDL_assert(idx < 20);
SDL_assert(idx <= m);
fw += x[j] * f[idx];
}
q[i] = fw;
}
jz = jk; jz = jk;
recompute: recompute:
/* distill q[] into iq[] reversingly */ /* distill q[] into iq[] reversingly */
for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) { for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
fw = (double) ((int32_t) (twon24 * z)); fw = (double)((int32_t)(twon24* z));
iq[i] = (int32_t) (z - two24 * fw); iq[i] = (int32_t)(z-two24*fw);
z = q[j - 1] + fw; z = q[j-1]+fw;
} }
/* compute n */ /* compute n */
z = scalbn(z, q0); /* actual value of z */ z = scalbn(z,q0); /* actual value of z */
z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */ z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
n = (int32_t) z; n = (int32_t) z;
z -= (double) n; z -= (double)n;
ih = 0; ih = 0;
if (q0 > 0) { /* need iq[jz-1] to determine n */ if(q0>0) { /* need iq[jz-1] to determine n */
i = (iq[jz - 1] >> (24 - q0)); i = (iq[jz-1]>>(24-q0)); n += i;
n += i; iq[jz-1] -= i<<(24-q0);
iq[jz - 1] -= i << (24 - q0); ih = iq[jz-1]>>(23-q0);
ih = iq[jz - 1] >> (23 - q0); }
} else if (q0 == 0) else if(q0==0) ih = iq[jz-1]>>23;
ih = iq[jz - 1] >> 23; else if(z>=0.5) ih=2;
else if (z >= 0.5)
ih = 2;
if (ih > 0) { /* q > 0.5 */ if(ih>0) { /* q > 0.5 */
n += 1; n += 1; carry = 0;
carry = 0; for(i=0;i<jz ;i++) { /* compute 1-q */
for (i = 0; i < jz; i++) { /* compute 1-q */ j = iq[i];
j = iq[i]; if(carry==0) {
if (carry == 0) { if(j!=0) {
if (j != 0) { carry = 1; iq[i] = 0x1000000- j;
carry = 1; }
iq[i] = 0x1000000 - j; } else iq[i] = 0xffffff - j;
} }
} else if(q0>0) { /* rare case: chance is 1 in 12 */
iq[i] = 0xffffff - j; switch(q0) {
} case 1:
if (q0 > 0) { /* rare case: chance is 1 in 12 */ iq[jz-1] &= 0x7fffff; break;
switch (q0) { case 2:
case 1: iq[jz-1] &= 0x3fffff; break;
iq[jz - 1] &= 0x7fffff; }
break; }
case 2: if(ih==2) {
iq[jz - 1] &= 0x3fffff; z = one - z;
break; if(carry!=0) z -= scalbn(one,q0);
} }
} }
if (ih == 2) {
z = one - z;
if (carry != 0)
z -= scalbn(one, q0);
}
}
/* check if recomputation is needed */ /* check if recomputation is needed */
if (z == zero) { if(z==zero) {
j = 0; j = 0;
for (i = jz - 1; i >= jk; i--) for (i=jz-1;i>=jk;i--) j |= iq[i];
j |= iq[i]; if(j==0) { /* need recomputation */
if (j == 0) { /* need recomputation */ for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */
for (k = 1; iq[jk - k] == 0; k++); /* k = no. of terms needed */
for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */ for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */
f[jx + i] = (double) ipio2[jv + i]; f[jx+i] = (double) ipio2[jv+i];
for (j = 0, fw = 0.0; j <= jx; j++) for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
fw += x[j] * f[jx + i - j]; q[i] = fw;
q[i] = fw; }
} jz += k;
jz += k; goto recompute;
goto recompute; }
} }
}
/* chop off zero terms */ /* chop off zero terms */
if (z == 0.0) { if(z==0.0) {
jz -= 1; jz -= 1; q0 -= 24;
q0 -= 24; while(iq[jz]==0) { jz--; q0-=24;}
while (iq[jz] == 0) { } else { /* break z into 24-bit if necessary */
jz--; z = scalbn(z,-q0);
q0 -= 24; if(z>=two24) {
} fw = (double)((int32_t)(twon24*z));
} else { /* break z into 24-bit if necessary */ iq[jz] = (int32_t)(z-two24*fw);
z = scalbn(z, -q0); jz += 1; q0 += 24;
if (z >= two24) { iq[jz] = (int32_t) fw;
fw = (double) ((int32_t) (twon24 * z)); } else iq[jz] = (int32_t) z ;
iq[jz] = (int32_t) (z - two24 * fw); }
jz += 1;
q0 += 24;
iq[jz] = (int32_t) fw;
} else
iq[jz] = (int32_t) z;
}
/* convert integer "bit" chunk to floating-point value */ /* convert integer "bit" chunk to floating-point value */
fw = scalbn(one, q0); fw = scalbn(one,q0);
for (i = jz; i >= 0; i--) { for(i=jz;i>=0;i--) {
q[i] = fw * (double) iq[i]; q[i] = fw*(double)iq[i]; fw*=twon24;
fw *= twon24; }
}
/* compute PIo2[0,...,jp]*q[jz,...,0] */ /* compute PIo2[0,...,jp]*q[jz,...,0] */
for (i = jz; i >= 0; i--) { for(i=jz;i>=0;i--) {
for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++) for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
fw += PIo2[k] * q[i + k]; fq[jz-i] = fw;
fq[jz - i] = fw; }
}
/* compress fq[] into y[] */ /* compress fq[] into y[] */
switch (prec) { switch(prec) {
case 0: case 0:
fw = 0.0; fw = 0.0;
for (i = jz; i >= 0; i--) for (i=jz;i>=0;i--) fw += fq[i];
fw += fq[i]; y[0] = (ih==0)? fw: -fw;
y[0] = (ih == 0) ? fw : -fw; break;
break; case 1:
case 1: case 2:
case 2: fw = 0.0;
fw = 0.0; for (i=jz;i>=0;i--) fw += fq[i];
for (i = jz; i >= 0; i--) y[0] = (ih==0)? fw: -fw;
fw += fq[i]; fw = fq[0]-fw;
y[0] = (ih == 0) ? fw : -fw; for (i=1;i<=jz;i++) fw += fq[i];
fw = fq[0] - fw; y[1] = (ih==0)? fw: -fw;
for (i = 1; i <= jz; i++) break;
fw += fq[i]; case 3: /* painful */
y[1] = (ih == 0) ? fw : -fw; for (i=jz;i>0;i--) {
break; fw = fq[i-1]+fq[i];
case 3: /* painful */ fq[i] += fq[i-1]-fw;
for (i = jz; i > 0; i--) { fq[i-1] = fw;
fw = fq[i - 1] + fq[i]; }
fq[i] += fq[i - 1] - fw; for (i=jz;i>1;i--) {
fq[i - 1] = fw; fw = fq[i-1]+fq[i];
} fq[i] += fq[i-1]-fw;
for (i = jz; i > 1; i--) { fq[i-1] = fw;
fw = fq[i - 1] + fq[i]; }
fq[i] += fq[i - 1] - fw; for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
fq[i - 1] = fw; if(ih==0) {
} y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
for (fw = 0.0, i = jz; i >= 2; i--) } else {
fw += fq[i]; y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
if (ih == 0) { }
y[0] = fq[0]; }
y[1] = fq[1]; return n&7;
y[2] = fw;
} else {
y[0] = -fq[0];
y[1] = -fq[1];
y[2] = -fw;
}
}
return n & 7;
} }

60
src/libm/k_sin.c
View File

@ -1,4 +1,3 @@
/* @(#)k_sin.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: k_sin.c,v 1.8 1995/05/10 20:46:31 jtc Exp $";
#endif
/* __kernel_sin( x, y, iy) /* __kernel_sin( x, y, iy)
* kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
* Input x is assumed to be bounded by ~pi/4 in magnitude. * Input x is assumed to be bounded by ~pi/4 in magnitude.
@ -46,42 +40,26 @@ static const char rcsid[] =
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
#ifdef __STDC__
static const double static const double
#else half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
static double S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
#endif S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
#ifdef __STDC__ double attribute_hidden __kernel_sin(double x, double y, int iy)
double attribute_hidden
__kernel_sin(double x, double y, int iy)
#else
double attribute_hidden
__kernel_sin(x, y, iy)
double x, y;
int iy; /* iy=0 if y is zero */
#endif
{ {
double z, r, v; double z,r,v;
int32_t ix; int32_t ix;
GET_HIGH_WORD(ix, x); GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff; /* high word of x */ ix &= 0x7fffffff; /* high word of x */
if (ix < 0x3e400000) { /* |x| < 2**-27 */ if(ix<0x3e400000) /* |x| < 2**-27 */
if ((int) x == 0) {if((int)x==0) return x;} /* generate inexact */
return x; z = x*x;
} /* generate inexact */ v = z*x;
z = x * x; r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
v = z * x; if(iy==0) return x+v*(S1+z*r);
r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6))); else return x-((z*(half*y-v*r)-y)-v*S1);
if (iy == 0)
return x + v * (S1 + z * r);
else
return x - ((z * (half * y - v * r) - y) - v * S1);
} }

2
src/libm/k_tan.c
View File

@ -66,7 +66,7 @@ T[] = {
2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */ 2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */
}; };
double __kernel_tan(double x, double y, int iy) double attribute_hidden __kernel_tan(double x, double y, int iy)
{ {
double z,r,v,w,s; double z,r,v,w,s;
int32_t ix,hx; int32_t ix,hx;

1
src/libm/math_libm.h
View File

@ -28,6 +28,7 @@ double SDL_uclibc_copysign(double x, double y);
double SDL_uclibc_cos(double x); double SDL_uclibc_cos(double x);
double SDL_uclibc_fabs(double x); double SDL_uclibc_fabs(double x);
double SDL_uclibc_floor(double x); double SDL_uclibc_floor(double x);
double SDL_uclibc_fmod(double x, double y);
double SDL_uclibc_log(double x); double SDL_uclibc_log(double x);
double SDL_uclibc_pow(double x, double y); double SDL_uclibc_pow(double x, double y);
double SDL_uclibc_scalbn(double x, int n); double SDL_uclibc_scalbn(double x, int n);

3
src/libm/math_private.h
View File

@ -21,9 +21,11 @@
#include "SDL_endian.h" #include "SDL_endian.h"
/* #include <sys/types.h> */ /* #include <sys/types.h> */
#define _IEEE_LIBM
#define attribute_hidden #define attribute_hidden
#define libm_hidden_proto(x) #define libm_hidden_proto(x)
#define libm_hidden_def(x) #define libm_hidden_def(x)
#define strong_alias(x, y)
#ifndef __HAIKU__ /* already defined in a system header. */ #ifndef __HAIKU__ /* already defined in a system header. */
typedef unsigned int u_int32_t; typedef unsigned int u_int32_t;
@ -35,6 +37,7 @@ typedef unsigned int u_int32_t;
#define cos SDL_uclibc_cos #define cos SDL_uclibc_cos
#define fabs SDL_uclibc_fabs #define fabs SDL_uclibc_fabs
#define floor SDL_uclibc_floor #define floor SDL_uclibc_floor
#define __ieee754_fmod SDL_uclibc_fmod
#define __ieee754_log SDL_uclibc_log #define __ieee754_log SDL_uclibc_log
#define __ieee754_pow SDL_uclibc_pow #define __ieee754_pow SDL_uclibc_pow
#define scalbn SDL_uclibc_scalbn #define scalbn SDL_uclibc_scalbn

1
src/libm/s_atan.c
View File

@ -112,4 +112,3 @@ double atan(double x)
} }
} }
libm_hidden_def(atan) libm_hidden_def(atan)

25
src/libm/s_copysign.c
View File

@ -1,4 +1,3 @@
/* @(#)s_copysign.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: s_copysign.c,v 1.8 1995/05/10 20:46:57 jtc Exp $";
#endif
/* /*
* copysign(double x, double y) * copysign(double x, double y)
* copysign(x,y) returns a value with the magnitude of x and * copysign(x,y) returns a value with the magnitude of x and
@ -24,19 +18,12 @@ static const char rcsid[] =
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
libm_hidden_proto(copysign) double copysign(double x, double y)
#ifdef __STDC__
double copysign(double x, double y)
#else
double copysign(x, y)
double x, y;
#endif
{ {
u_int32_t hx, hy; u_int32_t hx,hy;
GET_HIGH_WORD(hx, x); GET_HIGH_WORD(hx,x);
GET_HIGH_WORD(hy, y); GET_HIGH_WORD(hy,y);
SET_HIGH_WORD(x, (hx & 0x7fffffff) | (hy & 0x80000000)); SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
return x; return x;
} }
libm_hidden_def(copysign) libm_hidden_def(copysign)

52
src/libm/s_cos.c
View File

@ -1,4 +1,3 @@
/* @(#)s_cos.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: s_cos.c,v 1.7 1995/05/10 20:47:02 jtc Exp $";
#endif
/* cos(x) /* cos(x)
* Return cosine function of x. * Return cosine function of x.
* *
@ -49,43 +43,31 @@ static const char rcsid[] =
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
libm_hidden_proto(cos) double cos(double x)
#ifdef __STDC__
double cos(double x)
#else
double cos(x)
double x;
#endif
{ {
double y[2], z = 0.0; double y[2],z=0.0;
int32_t n, ix; int32_t n, ix;
/* High word of x. */ /* High word of x. */
GET_HIGH_WORD(ix, x); GET_HIGH_WORD(ix,x);
/* |x| ~< pi/4 */ /* |x| ~< pi/4 */
ix &= 0x7fffffff; ix &= 0x7fffffff;
if (ix <= 0x3fe921fb) if(ix <= 0x3fe921fb) return __kernel_cos(x,z);
return __kernel_cos(x, z);
/* cos(Inf or NaN) is NaN */ /* cos(Inf or NaN) is NaN */
else if (ix >= 0x7ff00000) else if (ix>=0x7ff00000) return x-x;
return x - x;
/* argument reduction needed */ /* argument reduction needed */
else { else {
n = __ieee754_rem_pio2(x, y); n = __ieee754_rem_pio2(x,y);
switch (n & 3) { switch(n&3) {
case 0: case 0: return __kernel_cos(y[0],y[1]);
return __kernel_cos(y[0], y[1]); case 1: return -__kernel_sin(y[0],y[1],1);
case 1: case 2: return -__kernel_cos(y[0],y[1]);
return -__kernel_sin(y[0], y[1], 1); default:
case 2: return __kernel_sin(y[0],y[1],1);
return -__kernel_cos(y[0], y[1]); }
default: }
return __kernel_sin(y[0], y[1], 1);
}
}
} }
libm_hidden_def(cos) libm_hidden_def(cos)

26
src/libm/s_fabs.c
View File

@ -1,4 +1,3 @@
/* @(#)s_fabs.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,30 +9,21 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: s_fabs.c,v 1.7 1995/05/10 20:47:13 jtc Exp $";
#endif
/* /*
* fabs(x) returns the absolute value of x. * fabs(x) returns the absolute value of x.
*/ */
/*#include <features.h>*/
/* Prevent math.h from defining a colliding inline */
#undef __USE_EXTERN_INLINES
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
libm_hidden_proto(fabs) double fabs(double x)
#ifdef __STDC__
double fabs(double x)
#else
double fabs(x)
double x;
#endif
{ {
u_int32_t high; u_int32_t high;
GET_HIGH_WORD(high, x); GET_HIGH_WORD(high,x);
SET_HIGH_WORD(x, high & 0x7fffffff); SET_HIGH_WORD(x,high&0x7fffffff);
return x; return x;
} }
libm_hidden_def(fabs) libm_hidden_def(fabs)

113
src/libm/s_floor.c
View File

@ -1,4 +1,3 @@
/* @(#)s_floor.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: s_floor.c,v 1.8 1995/05/10 20:47:20 jtc Exp $";
#endif
/* /*
* floor(x) * floor(x)
* Return x rounded toward -inf to integral value * Return x rounded toward -inf to integral value
@ -24,73 +18,54 @@ static const char rcsid[] =
* Inexact flag raised if x not equal to floor(x). * Inexact flag raised if x not equal to floor(x).
*/ */
/*#include <features.h>*/
/* Prevent math.h from defining a colliding inline */
#undef __USE_EXTERN_INLINES
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
#ifdef __STDC__ static const double huge = 1.0e300;
static const double huge_val = 1.0e300;
#else
static double huge_val = 1.0e300;
#endif
libm_hidden_proto(floor) double floor(double x)
#ifdef __STDC__
double floor(double x)
#else
double floor(x)
double x;
#endif
{ {
int32_t i0, i1, j0; int32_t i0,i1,j0;
u_int32_t i, j; u_int32_t i,j;
EXTRACT_WORDS(i0, i1, x); EXTRACT_WORDS(i0,i1,x);
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; j0 = ((i0>>20)&0x7ff)-0x3ff;
if (j0 < 20) { if(j0<20) {
if (j0 < 0) { /* raise inexact if x != 0 */ if(j0<0) { /* raise inexact if x != 0 */
if (huge_val + x > 0.0) { /* return 0*sign(x) if |x|<1 */ if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
if (i0 >= 0) { if(i0>=0) {i0=i1=0;}
i0 = i1 = 0; else if(((i0&0x7fffffff)|i1)!=0)
} else if (((i0 & 0x7fffffff) | i1) != 0) { { i0=0xbff00000;i1=0;}
i0 = 0xbff00000; }
i1 = 0; } else {
} i = (0x000fffff)>>j0;
} if(((i0&i)|i1)==0) return x; /* x is integral */
} else { if(huge+x>0.0) { /* raise inexact flag */
i = (0x000fffff) >> j0; if(i0<0) i0 += (0x00100000)>>j0;
if (((i0 & i) | i1) == 0) i0 &= (~i); i1=0;
return x; /* x is integral */ }
if (huge_val + x > 0.0) { /* raise inexact flag */ }
if (i0 < 0) } else if (j0>51) {
i0 += (0x00100000) >> j0; if(j0==0x400) return x+x; /* inf or NaN */
i0 &= (~i); else return x; /* x is integral */
i1 = 0; } else {
} i = ((u_int32_t)(0xffffffff))>>(j0-20);
} if((i1&i)==0) return x; /* x is integral */
} else if (j0 > 51) { if(huge+x>0.0) { /* raise inexact flag */
if (j0 == 0x400) if(i0<0) {
return x + x; /* inf or NaN */ if(j0==20) i0+=1;
else else {
return x; /* x is integral */ j = i1+(1<<(52-j0));
} else { if(j<i1) i0 +=1 ; /* got a carry */
i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); i1=j;
if ((i1 & i) == 0) }
return x; /* x is integral */ }
if (huge_val + x > 0.0) { /* raise inexact flag */ i1 &= (~i);
if (i0 < 0) { }
if (j0 == 20) }
i0 += 1; INSERT_WORDS(x,i0,i1);
else { return x;
j = i1 + (1 << (52 - j0));
if (j < (u_int32_t) i1)
i0 += 1; /* got a carry */
i1 = j;
}
}
i1 &= (~i);
}
}
INSERT_WORDS(x, i0, i1);
return x;
} }
libm_hidden_def(floor) libm_hidden_def(floor)

110
src/libm/s_scalbn.c
View File

@ -1,4 +1,3 @@
/* @(#)s_scalbn.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,70 +9,69 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $";
#endif
/* /*
* scalbn (double x, int n) * scalbln(double x, long n)
* scalbn(x,n) returns x* 2**n computed by exponent * scalbln(x,n) returns x * 2**n computed by exponent
* manipulation rather than by actually performing an * manipulation rather than by actually performing an
* exponentiation or a multiplication. * exponentiation or a multiplication.
*/ */
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
#include <limits.h>
libm_hidden_proto(copysign) static const double
#ifdef __STDC__ two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
static const double twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
#else huge = 1.0e+300,
static double tiny = 1.0e-300;
#endif
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
huge_val = 1.0e+300, tiny = 1.0e-300;
libm_hidden_proto(scalbn) double scalbln(double x, long n)
#ifdef __STDC__
double scalbn(double x, int n)
#else
double scalbn(x, n)
double x;
int n;
#endif
{ {
int32_t k, hx, lx; int32_t k, hx, lx;
EXTRACT_WORDS(hx, lx, x);
k = (hx & 0x7ff00000) >> 20; /* extract exponent */
if (k == 0) { /* 0 or subnormal x */
if ((lx | (hx & 0x7fffffff)) == 0)
return x; /* +-0 */
x *= two54;
GET_HIGH_WORD(hx, x);
k = ((hx & 0x7ff00000) >> 20) - 54;
if (n < -50000)
return tiny * x; /* underflow */
}
if (k == 0x7ff)
return x + x; /* NaN or Inf */
k = k + n;
if (k > 0x7fe)
return huge_val * copysign(huge_val, x); /* overflow */
if (k > 0) { /* normal result */
SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
return x;
}
if (k <= -54) {
if (n > 50000) /* in case integer overflow in n+k */
return huge_val * copysign(huge_val, x); /* overflow */
else
return tiny * copysign(tiny, x); /* underflow */
}
k += 54; /* subnormal result */
SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
return x * twom54;
}
EXTRACT_WORDS(hx, lx, x);
k = (hx & 0x7ff00000) >> 20; /* extract exponent */
if (k == 0) { /* 0 or subnormal x */
if ((lx | (hx & 0x7fffffff)) == 0)
return x; /* +-0 */
x *= two54;
GET_HIGH_WORD(hx, x);
k = ((hx & 0x7ff00000) >> 20) - 54;
}
if (k == 0x7ff)
return x + x; /* NaN or Inf */
k = k + n;
if (k > 0x7fe)
return huge * copysign(huge, x); /* overflow */
if (n < -50000)
return tiny * copysign(tiny, x); /* underflow */
if (k > 0) { /* normal result */
SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
return x;
}
if (k <= -54) {
if (n > 50000) /* in case integer overflow in n+k */
return huge * copysign(huge, x); /* overflow */
return tiny * copysign(tiny, x); /* underflow */
}
k += 54; /* subnormal result */
SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
return x * twom54;
}
libm_hidden_def(scalbln)
#if LONG_MAX == INT_MAX
/* strong_alias(scalbln, scalbn) - "error: conflicting types for 'scalbn'"
* because it tries to declare "typeof(scalbln) scalbn;"
* which tries to give "long" parameter to scalbn.
* Doing it by hand:
*/
__typeof(scalbn) scalbn __attribute__((alias("scalbln")));
#else
double scalbn(double x, int n)
{
return scalbln(x, n);
}
#endif
libm_hidden_def(scalbn) libm_hidden_def(scalbn)

52
src/libm/s_sin.c
View File

@ -1,4 +1,3 @@
/* @(#)s_sin.c 5.1 93/09/24 */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@ -10,11 +9,6 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static const char rcsid[] =
"$NetBSD: s_sin.c,v 1.7 1995/05/10 20:48:15 jtc Exp $";
#endif
/* sin(x) /* sin(x)
* Return sine function of x. * Return sine function of x.
* *
@ -49,43 +43,31 @@ static const char rcsid[] =
#include "math_libm.h" #include "math_libm.h"
#include "math_private.h" #include "math_private.h"
libm_hidden_proto(sin) double sin(double x)
#ifdef __STDC__
double sin(double x)
#else
double sin(x)
double x;
#endif
{ {
double y[2], z = 0.0; double y[2],z=0.0;
int32_t n, ix; int32_t n, ix;
/* High word of x. */ /* High word of x. */
GET_HIGH_WORD(ix, x); GET_HIGH_WORD(ix,x);
/* |x| ~< pi/4 */ /* |x| ~< pi/4 */
ix &= 0x7fffffff; ix &= 0x7fffffff;
if (ix <= 0x3fe921fb) if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0);
return __kernel_sin(x, z, 0);
/* sin(Inf or NaN) is NaN */ /* sin(Inf or NaN) is NaN */
else if (ix >= 0x7ff00000) else if (ix>=0x7ff00000) return x-x;
return x - x;
/* argument reduction needed */ /* argument reduction needed */
else { else {
n = __ieee754_rem_pio2(x, y); n = __ieee754_rem_pio2(x,y);
switch (n & 3) { switch(n&3) {
case 0: case 0: return __kernel_sin(y[0],y[1],1);
return __kernel_sin(y[0], y[1], 1); case 1: return __kernel_cos(y[0],y[1]);
case 1: case 2: return -__kernel_sin(y[0],y[1],1);
return __kernel_cos(y[0], y[1]); default:
case 2: return -__kernel_cos(y[0],y[1]);
return -__kernel_sin(y[0], y[1], 1); }
default: }
return -__kernel_cos(y[0], y[1]);
}
}
} }
libm_hidden_def(sin) libm_hidden_def(sin)