d/freebsd-skq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

When j0() and j1() were converted to j0f() and j1f(), the threshold

Browse Source 
values for the different invervals were not converted correctly.
Adjust the threshold values to values, which should agree with the
comments.

Reported by:	cognet (j1f only)
Discussed with: pfg, bde
Reviewed by:	bde
This commit is contained in:
Steve Kargl 2015-03-01 20:26:03 +00:00
parent 7190a55c3e
commit 8192c85cb3
2 changed files with 16 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
lib/msun/src/e_j0f.c
View File

@ -62,7 +62,7 @@ __ieee754_j0f(float x)
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if(ix>0x54000000) z = (invsqrtpi*cc)/sqrtf(x);
if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */
else {
u = pzerof(x); v = qzerof(x);
z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
@ -136,14 +136,14 @@ __ieee754_y0f(float x)
if ((s*c)<zero) cc = z/ss;
else ss = z/cc;
}
if(ix>0x54800000) z = (invsqrtpi*ss)/sqrtf(x);
if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
else {
u = pzerof(x); v = qzerof(x);
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
}
return z;
}
if(ix<=0x32000000) { /* x < 2**-27 */
if(ix<=0x39000000) { /* x < 2**-13 */
return(u00 + tpi*__ieee754_logf(x));
}
z = x*x;
@ -232,8 +232,8 @@ static const float pS2[5] = {
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x41000000) {p = pR8; q= pS8;}
else if(ix>=0x40f71c58){p = pR5; q= pS5;}
else if(ix>=0x4036db68){p = pR3; q= pS3;}
else if(ix>=0x409173eb){p = pR5; q= pS5;}
else if(ix>=0x4036d917){p = pR3; q= pS3;}
else {p = pR2; q= pS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
@ -327,8 +327,8 @@ static const float qS2[6] = {
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x41000000) {p = qR8; q= qS8;}
else if(ix>=0x40f71c58){p = qR5; q= qS5;}
else if(ix>=0x4036db68){p = qR3; q= qS3;}
else if(ix>=0x409173eb){p = qR5; q= qS5;}
else if(ix>=0x4036d917){p = qR3; q= qS3;}
else {p = qR2; q= qS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));

18
lib/msun/src/e_j1f.c
View File

@ -63,7 +63,7 @@ __ieee754_j1f(float x)
* j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
* y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
*/
if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */
else {
u = ponef(y); v = qonef(y);
z = invsqrtpi*(u*cc-v*ss)/sqrtf(y);
@ -71,7 +71,7 @@ __ieee754_j1f(float x)
if(hx<0) return -z;
else return z;
}
if(ix<0x32000000) { /* |x|<2**-27 */
if(ix<0x39000000) { /* |x|<2**-13 */
if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */
}
z = x*x;
@ -129,14 +129,14 @@ __ieee754_y1f(float x)
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x);
if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
else {
u = ponef(x); v = qonef(x);
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
}
return z;
}
if(ix<=0x24800000) { /* x < 2**-54 */
if(ix<=0x33000000) { /* x < 2**-25 */
return(-tpi/x);
}
z = x*x;
@ -227,8 +227,8 @@ static const float ps2[5] = {
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x41000000) {p = pr8; q= ps8;}
else if(ix>=0x40f71c58){p = pr5; q= ps5;}
else if(ix>=0x4036db68){p = pr3; q= ps3;}
else if(ix>=0x409173eb){p = pr5; q= ps5;}
else if(ix>=0x4036d917){p = pr3; q= ps3;}
else {p = pr2; q= ps2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
@ -322,9 +322,9 @@ static const float qs2[6] = {
int32_t ix;
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x40200000) {p = qr8; q= qs8;}
else if(ix>=0x40f71c58){p = qr5; q= qs5;}
else if(ix>=0x4036db68){p = qr3; q= qs3;}
if(ix>=0x41000000) {p = qr8; q= qs8;}
else if(ix>=0x409173eb){p = qr5; q= qs5;}
else if(ix>=0x4036d917){p = qr3; q= qs3;}
else {p = qr2; q= qs2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));