diff --git a/lib/msun/src/e_log10.c b/lib/msun/src/e_log10.c index ad0aed5f1de2..104d2574c659 100644 --- a/lib/msun/src/e_log10.c +++ b/lib/msun/src/e_log10.c @@ -17,6 +17,9 @@ __FBSDID("$FreeBSD$"); /* * Return the base 10 logarithm of x. See e_log.c and k_log.h for most * comments. + * + * log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2) + * in not-quite-routine extra precision. */ #include "math.h" @@ -35,7 +38,7 @@ static const double zero = 0.0; double __ieee754_log10(double x) { - double f,hi,lo,y,z; + double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2; int32_t i,k,hx; u_int32_t lx; @@ -50,16 +53,35 @@ __ieee754_log10(double x) GET_HIGH_WORD(hx,x); } if (hx >= 0x7ff00000) return x+x; + if (hx == 0x3ff00000 && lx == 0) + return zero; /* log(1) = +0 */ k += (hx>>20)-1023; hx &= 0x000fffff; i = (hx+0x95f64)&0x100000; SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ k += (i>>20); y = (double)k; - f = __kernel_log(x); - hi = x = x - 1; + f = x - 1.0; + hfsq = 0.5*f*f; + r = k_log1p(f); + + /* See e_log2.c for most details. */ + hi = f - hfsq; SET_LOW_WORD(hi,0); - lo = x - hi; - z = y*log10_2lo + (x+f)*ivln10lo + (lo+f)*ivln10hi + hi*ivln10hi; - return z+y*log10_2hi; + lo = (f - hi) - hfsq + r; + val_hi = hi*ivln10hi; + y2 = y*log10_2hi; + val_lo = y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi; + + /* + * Extra precision in for adding y*log10_2hi is not strictly needed + * since there is no very large cancellation near x = sqrt(2) or + * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs + * with some parallelism and it reduces the error for many args. + */ + w = y2 + val_hi; + val_lo += (y2 - w) + val_hi; + val_hi = w; + + return val_lo + val_hi; } diff --git a/lib/msun/src/e_log10f.c b/lib/msun/src/e_log10f.c index 960b608bf2b4..46658986eeea 100644 --- a/lib/msun/src/e_log10f.c +++ b/lib/msun/src/e_log10f.c @@ -32,7 +32,7 @@ static const float zero = 0.0; float __ieee754_log10f(float x) { - float f,hi,lo,y,z; + float f,hfsq,hi,lo,r,y,y2; int32_t i,k,hx; GET_FLOAT_WORD(hx,x); @@ -46,17 +46,26 @@ __ieee754_log10f(float x) GET_FLOAT_WORD(hx,x); } if (hx >= 0x7f800000) return x+x; + if (hx == 0x3f800000) + return zero; /* log(1) = +0 */ k += (hx>>23)-127; hx &= 0x007fffff; i = (hx+(0x4afb0d))&0x800000; SET_FLOAT_WORD(x,hx|(i^0x3f800000)); /* normalize x or x/2 */ k += (i>>23); y = (float)k; - f = __kernel_logf(x); - x = x - (float)1.0; - GET_FLOAT_WORD(hx,x); + f = x - (float)1.0; + hfsq = (float)0.5*f*f; + r = k_log1pf(f); + + /* See e_log2f.c and e_log2.c for details. */ + if (sizeof(float_t) > sizeof(float)) + return (r - hfsq + f) * ((float_t)ivln10lo + ivln10hi) + + y * ((float_t)log10_2lo + log10_2hi); + hi = f - hfsq; + GET_FLOAT_WORD(hx,hi); SET_FLOAT_WORD(hi,hx&0xfffff000); - lo = x - hi; - z = y*log10_2lo + (x+f)*ivln10lo + (lo+f)*ivln10hi + hi*ivln10hi; - return z+y*log10_2hi; + lo = (f - hi) - hfsq + r; + return y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi + hi*ivln10hi + + y*log10_2hi; } diff --git a/lib/msun/src/e_log2.c b/lib/msun/src/e_log2.c index d237759d8a41..1fc44a5fedec 100644 --- a/lib/msun/src/e_log2.c +++ b/lib/msun/src/e_log2.c @@ -17,6 +17,11 @@ __FBSDID("$FreeBSD$"); /* * Return the base 2 logarithm of x. See e_log.c and k_log.h for most * comments. + * + * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel, + * then does the combining and scaling steps + * log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k + * in not-quite-routine extra precision. */ #include "math.h" @@ -33,7 +38,7 @@ static const double zero = 0.0; double __ieee754_log2(double x) { - double f,hi,lo; + double f,hfsq,hi,lo,r,val_hi,val_lo,w,y; int32_t i,k,hx; u_int32_t lx; @@ -48,14 +53,58 @@ __ieee754_log2(double x) GET_HIGH_WORD(hx,x); } if (hx >= 0x7ff00000) return x+x; + if (hx == 0x3ff00000 && lx == 0) + return zero; /* log(1) = +0 */ k += (hx>>20)-1023; hx &= 0x000fffff; i = (hx+0x95f64)&0x100000; SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ k += (i>>20); - f = __kernel_log(x); - hi = x = x - 1; + y = (double)k; + f = x - 1.0; + hfsq = 0.5*f*f; + r = k_log1p(f); + + /* + * f-hfsq must (for args near 1) be evaluated in extra precision + * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2). + * This is fairly efficient since f-hfsq only depends on f, so can + * be evaluated in parallel with R. Not combining hfsq with R also + * keeps R small (though not as small as a true `lo' term would be), + * so that extra precision is not needed for terms involving R. + * + * Compiler bugs involving extra precision used to break Dekker's + * theorem for spitting f-hfsq as hi+lo, unless double_t was used + * or the multi-precision calculations were avoided when double_t + * has extra precision. These problems are now automatically + * avoided as a side effect of the optimization of combining the + * Dekker splitting step with the clear-low-bits step. + * + * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra + * precision to avoid a very large cancellation when x is very near + * these values. Unlike the above cancellations, this problem is + * specific to base 2. It is strange that adding +-1 is so much + * harder than adding +-ln2 or +-log10_2. + * + * This uses Dekker's theorem to normalize y+val_hi, so the + * compiler bugs are back in some configurations, sigh. And I + * don't want to used double_t to avoid them, since that gives a + * pessimization and the support for avoiding the pessimization + * is not yet available. + * + * The multi-precision calculations for the multiplications are + * routine. + */ + hi = f - hfsq; SET_LOW_WORD(hi,0); - lo = x - hi; - return (x+f)*ivln2lo + (lo+f)*ivln2hi + hi*ivln2hi + k; + lo = (f - hi) - hfsq + r; + val_hi = hi*ivln2hi; + val_lo = (lo+hi)*ivln2lo + lo*ivln2hi; + + /* spadd(val_hi, val_lo, y), except for not using double_t: */ + w = y + val_hi; + val_lo += (y - w) + val_hi; + val_hi = w; + + return val_lo + val_hi; } diff --git a/lib/msun/src/e_log2f.c b/lib/msun/src/e_log2f.c index ba63bc3d255c..716634664a42 100644 --- a/lib/msun/src/e_log2f.c +++ b/lib/msun/src/e_log2f.c @@ -30,7 +30,7 @@ static const float zero = 0.0; float __ieee754_log2f(float x) { - float f,hi,lo; + float f,hfsq,hi,lo,r,y; int32_t i,k,hx; GET_FLOAT_WORD(hx,x); @@ -44,15 +44,38 @@ __ieee754_log2f(float x) GET_FLOAT_WORD(hx,x); } if (hx >= 0x7f800000) return x+x; + if (hx == 0x3f800000) + return zero; /* log(1) = +0 */ k += (hx>>23)-127; hx &= 0x007fffff; i = (hx+(0x4afb0d))&0x800000; SET_FLOAT_WORD(x,hx|(i^0x3f800000)); /* normalize x or x/2 */ k += (i>>23); - f = __kernel_logf(x); - x = x - (float)1.0; - GET_FLOAT_WORD(hx,x); + y = (float)k; + f = x - (float)1.0; + hfsq = (float)0.5*f*f; + r = k_log1pf(f); + + /* + * We no longer need to avoid falling into the multi-precision + * calculations due to compiler bugs breaking Dekker's theorem. + * Keep avoiding this as an optimization. See e_log2.c for more + * details (some details are here only because the optimization + * is not yet available in double precision). + * + * Another compiler bug turned up. With gcc on i386, + * (ivln2lo + ivln2hi) would be evaluated in float precision + * despite runtime evaluations using double precision. So we + * must cast one of its terms to float_t. This makes the whole + * expression have type float_t, so return is forced to waste + * time clobbering its extra precision. + */ + if (sizeof(float_t) > sizeof(float)) + return (r - hfsq + f) * ((float_t)ivln2lo + ivln2hi) + y; + + hi = f - hfsq; + GET_FLOAT_WORD(hx,hi); SET_FLOAT_WORD(hi,hx&0xfffff000); - lo = x - hi; - return (x+f)*ivln2lo + (lo+f)*ivln2hi + hi*ivln2hi + k; + lo = (f - hi) - hfsq + r; + return (lo+hi)*ivln2lo + lo*ivln2hi + hi*ivln2hi + y; } diff --git a/lib/msun/src/k_log.h b/lib/msun/src/k_log.h index 206355c272ac..aaff8bd9000e 100644 --- a/lib/msun/src/k_log.h +++ b/lib/msun/src/k_log.h @@ -14,8 +14,9 @@ #include __FBSDID("$FreeBSD$"); -/* __kernel_log(x) - * Return log(x) - (x-1) for x in ~[sqrt(2)/2, sqrt(2)]. +/* + * k_log1p(f): + * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)]. * * The following describes the overall strategy for computing * logarithms in base e. The argument reduction and adding the final @@ -80,37 +81,20 @@ Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ /* - * We always inline __kernel_log(), since doing so produces a + * We always inline k_log1p(), since doing so produces a * substantial performance improvement (~40% on amd64). */ static inline double -__kernel_log(double x) +k_log1p(double f) { - double hfsq,f,s,z,R,w,t1,t2; - int32_t hx,i,j; - u_int32_t lx; + double hfsq,s,z,R,w,t1,t2; - EXTRACT_WORDS(hx,lx,x); - - f = x-1.0; - if((0x000fffff&(2+hx))<3) { /* -2**-20 <= f < 2**-20 */ - if(f==0.0) return 0.0; - return f*f*(0.33333333333333333*f-0.5); - } - s = f/(2.0+f); + s = f/(2.0+f); z = s*s; - hx &= 0x000fffff; - i = hx-0x6147a; w = z*z; - j = 0x6b851-hx; - t1= w*(Lg2+w*(Lg4+w*Lg6)); - t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); - i |= j; + t1= w*(Lg2+w*(Lg4+w*Lg6)); + t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); R = t2+t1; - if (i>0) { - hfsq=0.5*f*f; - return s*(hfsq+R) - hfsq; - } else { - return s*(R-f); - } + hfsq=0.5*f*f; + return s*(hfsq+R); } diff --git a/lib/msun/src/k_logf.h b/lib/msun/src/k_logf.h index d9f0f3ddfde1..71c547e888ae 100644 --- a/lib/msun/src/k_logf.h +++ b/lib/msun/src/k_logf.h @@ -13,7 +13,7 @@ __FBSDID("$FreeBSD$"); /* - * float version of __kernel_log(x). See k_log.c for details. + * Float version of k_log.h. See the latter for most comments. */ static const float @@ -24,32 +24,16 @@ Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ static inline float -__kernel_logf(float x) +k_log1pf(float f) { - float hfsq,f,s,z,R,w,t1,t2; - int32_t ix,i,j; + float hfsq,s,z,R,w,t1,t2; - GET_FLOAT_WORD(ix,x); - - f = x-(float)1.0; - if((0x007fffff&(0x8000+ix))<0xc000) { /* -2**-9 <= f < 2**-9 */ - if(f==0.0f) return 0.0f; - return f*f*((float)0.33333333333333333*f-(float)0.5); - } s = f/((float)2.0+f); z = s*s; - ix &= 0x007fffff; - i = ix-(0x6147a<<3); w = z*z; - j = (0x6b851<<3)-ix; t1= w*(Lg2+w*Lg4); t2= z*(Lg1+w*Lg3); - i |= j; R = t2+t1; - if(i>0) { - hfsq=(float)0.5*f*f; - return s*(hfsq+R) - hfsq; - } else { - return s*(R-f); - } + hfsq=(float)0.5*f*f; + return s*(hfsq+R); }