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).

2 files changed, 22 insertions, 8 deletions

diff --git a/lib/msun/src/s_log1p.c b/lib/msun/src/s_log1p.c
index 4ce4c3a..0ee69a6 100644
--- a/lib/msun/src/s_log1p.c
+++ b/lib/msun/src/s_log1p.c
@@ -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;
diff --git a/lib/msun/src/s_log1pf.c b/lib/msun/src/s_log1pf.c
index 76fb33c..eb0f4af 100644
--- a/lib/msun/src/s_log1pf.c
+++ b/lib/msun/src/s_log1pf.c
@@ -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;