mirror of
https://git.FreeBSD.org/src.git
synced 2024-12-16 10:20:30 +00:00
For log1pf(), fixed the approximations to sqrt(2), sqrt(2)-1 and
sqrt(2)/2-1. For log1p(), fixed the approximation to sqrt(2)/2-1. The end result is to fix an error of 1.293 ulps in log1pf(0.41421395540 (hex 0x3ed413da)) and an error of 1.783 ulps in log1p(-0.292893409729003961761) (hex 0x12bec4 00000001)). The former was the only error of > 1 ulp for log1pf() and the latter is the only such error that I know of for log1p(). The approximations don't need to be very accurate, but the last 2 need to be related to the first and be rounded up a little (even more than 1 ulp for sqrt(2)/2-1) for the following implementation-detail reason: when the arg (x) is not between (the approximations to) sqrt(2)/2-1 and sqrt(2)-1, we commit to using a correction term, but we only actually use it if 1+x is between sqrt(2)/2 and sqrt(2) according to the first approximation. Thus we must ensure that !(sqrt(2)/2-1 < x < sqrt(2)-1) implies !(sqrt(2)/2 < x+1 < sqrt(2)), where all the sqrt(2)'s are really slightly different approximations to sqrt(2) and some of the "<"'s are really "<="'s. This was not done. In log1pf(), the last 2 approximations were rounded up by about 6 ulps more than needed relative to a good approximation to sqrt(2), but the actual approximation to sqrt(2) was off by 3 ulps. The approximation to sqrt(2)-1 ended up being 4 ulps too small, so the algoritm was broken in 4 cases. The result happened to be broken in 1 case. This is fixed by using a natural approximation to sqrt(2) and derived approximations for the others. In logf(), all the approximations made sense, but the approximation to sqrt(2)/2-1 was 2 ulps too small (a tiny amount, since we compare with a granularity of 2**32 ulps), so the algorithm was broken in 2 cases. The result was broken in 1 case. This is fixed by rounding up the approximation to sqrt(2)/2-1 by 2**32 ulps, so 2**32 cases are now handled a little differently (still correctly according to my assertion that the approximations don't need to be very accurate, but this has not been checked).
This commit is contained in:
parent
69a7db097e
commit
d48ea9753c
Notes:
svn2git
2020-12-20 02:59:44 +00:00
svn path=/head/; revision=153086
@ -106,7 +106,7 @@ log1p(double x)
|
||||
ax = hx&0x7fffffff;
|
||||
|
||||
k = 1;
|
||||
if (hx < 0x3FDA827A) { /* x < 0.41422 */
|
||||
if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */
|
||||
if(ax>=0x3ff00000) { /* x <= -1.0 */
|
||||
if(x==-1.0) return -two54/zero; /* log1p(-1)=+inf */
|
||||
else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
|
||||
@ -118,8 +118,8 @@ log1p(double x)
|
||||
else
|
||||
return x - x*x*0.5;
|
||||
}
|
||||
if(hx>0||hx<=((int32_t)0xbfd2bec3)) {
|
||||
k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */
|
||||
if(hx>0||hx<=((int32_t)0xbfd2bec4)) {
|
||||
k=0;f=x;hu=1;} /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
|
||||
}
|
||||
if (hx >= 0x7ff00000) return x+x;
|
||||
if(k!=0) {
|
||||
@ -136,7 +136,14 @@ log1p(double x)
|
||||
c = 0;
|
||||
}
|
||||
hu &= 0x000fffff;
|
||||
if(hu<0x6a09e) {
|
||||
/*
|
||||
* The approximation to sqrt(2) used in thresholds is not
|
||||
* critical. However, the ones used above must give less
|
||||
* strict bounds than the one here so that the k==0 case is
|
||||
* never reached from here, since here we have committed to
|
||||
* using the correction term but don't use it if k==0.
|
||||
*/
|
||||
if(hu<0x6a09e) { /* u ~< sqrt(2) */
|
||||
SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */
|
||||
} else {
|
||||
k += 1;
|
||||
|
@ -44,7 +44,7 @@ log1pf(float x)
|
||||
ax = hx&0x7fffffff;
|
||||
|
||||
k = 1;
|
||||
if (hx < 0x3ed413d7) { /* x < 0.41422 */
|
||||
if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */
|
||||
if(ax>=0x3f800000) { /* x <= -1.0 */
|
||||
if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */
|
||||
else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
|
||||
@ -56,8 +56,8 @@ log1pf(float x)
|
||||
else
|
||||
return x - x*x*(float)0.5;
|
||||
}
|
||||
if(hx>0||hx<=((int32_t)0xbe95f61f)) {
|
||||
k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */
|
||||
if(hx>0||hx<=((int32_t)0xbe95f619)) {
|
||||
k=0;f=x;hu=1;} /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
|
||||
}
|
||||
if (hx >= 0x7f800000) return x+x;
|
||||
if(k!=0) {
|
||||
@ -75,7 +75,14 @@ log1pf(float x)
|
||||
c = 0;
|
||||
}
|
||||
hu &= 0x007fffff;
|
||||
if(hu<0x3504f7) {
|
||||
/*
|
||||
* The approximation to sqrt(2) used in thresholds is not
|
||||
* critical. However, the ones used above must give less
|
||||
* strict bounds than the one here so that the k==0 case is
|
||||
* never reached from here, since here we have committed to
|
||||
* using the correction term but don't use it if k==0.
|
||||
*/
|
||||
if(hu<0x3504f4) { /* u < sqrt(2) */
|
||||
SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */
|
||||
} else {
|
||||
k += 1;
|
||||
|
Loading…
Reference in New Issue
Block a user