Fixed the threshold for using the simple Taylor approximation.

In e_log.c, there was just a off-by-1 (1 ulp) error in the comment
about the threshold.  The precision of the threshold is unimportant,
but the magic numbers in the code are easier to understand when the
threshold is described precisely.

In e_logf.c, mistranslation of the magic numbers gave an off-by-1
(1 * 16 ulps) error in the intended negative bound for the threshold
and an off-by-7 (7 * 16 ulps) error in the intended positive bound for
the threshold, and the intended bounds were not translated from the
double precision bounds so they were unnecessarily small by a factor
of about 2048.

The optimization of using the simple Taylor approximation for args
near a power of 2 is dubious since it only applies to a relatively
small proportion of args, but if it is done then doing it 2048 times
as often _may_ be more efficient.  (My benchmarks show unexplained
dependencies on the data that increase with further optimizations
in this area.)

2 files changed, 2 insertions, 2 deletions

diff --git a/lib/msun/src/e_log.c b/lib/msun/src/e_log.c
index eee7cf9..ed3c401 100644
--- a/lib/msun/src/e_log.c
+++ b/lib/msun/src/e_log.c
@@ -107,7 +107,7 @@ __ieee754_log(double x)
 	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
 	k += (i>>20);
 	f = x-1.0;
-	if((0x000fffff&(2+hx))<3) {	/* |f| < 2**-20 */
+	if((0x000fffff&(2+hx))<3) {	/* -2**-20 <= f < 2**-20 */
 	    if(f==zero) if(k==0) return zero;  else {dk=(double)k;
 				 return dk*ln2_hi+dk*ln2_lo;}
 	    R = f*f*(0.5-0.33333333333333333*f);
diff --git a/lib/msun/src/e_logf.c b/lib/msun/src/e_logf.c
index b47b77b..175c955 100644
--- a/lib/msun/src/e_logf.c
+++ b/lib/msun/src/e_logf.c
@@ -55,7 +55,7 @@ __ieee754_logf(float x)
 	SET_FLOAT_WORD(x,ix|(i^0x3f800000));	/* normalize x or x/2 */
 	k += (i>>23);
 	f = x-(float)1.0;
-	if((0x007fffff&(15+ix))<16) {	/* |f| < 2**-20 */
+	if((0x007fffff&(0x8000+ix))<0xc000) {	/* -2**-9 <= f < 2**-9 */
 	    if(f==zero) if(k==0) return zero;  else {dk=(float)k;
 				 return dk*ln2_hi+dk*ln2_lo;}
 	    R = f*f*((float)0.5-(float)0.33333333333333333*f);