From a2f75f88a983eec32b81ea942988bb3703a86f8b Mon Sep 17 00:00:00 2001 From: das Date: Mon, 17 Oct 2011 05:38:07 +0000 Subject: [PATCH] MFC various log* improvements. r216247 - log2f style r216248 - log2f insignificant bug r219360 - log10 converted to use k_log r219361 - log10f converted to use k_log r226375 - log2/log10 style r226376 - log2/log10 bde's improvements; fix log(1) with FE_DOWNWARD rounding git-svn-id: svn://svn.freebsd.org/base/stable/8@226456 ccf9f872-aa2e-dd11-9fc8-001c23d0bc1f --- lib/msun/man/math.3 | 7 +++--- lib/msun/src/e_log10.c | 54 ++++++++++++----------------------------- lib/msun/src/e_log10f.c | 30 ++++++++++++++--------- lib/msun/src/e_log2.c | 8 +++--- lib/msun/src/e_log2f.c | 10 +++++--- lib/msun/src/k_logf.h | 6 ++--- 6 files changed, 52 insertions(+), 63 deletions(-) diff --git a/lib/msun/man/math.3 b/lib/msun/man/math.3 index 366d18553..f8521a17b 100644 --- a/lib/msun/man/math.3 +++ b/lib/msun/man/math.3 @@ -198,7 +198,7 @@ yn Bessel function of the second kind of the order n .El .Pp The routines -in this section may not produce a result that is correctly rounded, +in this section might not produce a result that is correctly rounded, so reproducible results cannot be guaranteed across platforms. For most of these functions, however, incorrect rounding occurs rarely, and then only in very-close-to-halfway cases. @@ -234,9 +234,8 @@ Many of the routines to compute transcendental functions produce inaccurate results in other than the default rounding mode. .Pp On the i386 platform, trigonometric argument reduction is not -performed accurately for very large arguments, resulting in errors -greater than 1 -.Em ulp +performed accurately for huge arguments, resulting in +large errors for such arguments to .Fn cos , .Fn sin , diff --git a/lib/msun/src/e_log10.c b/lib/msun/src/e_log10.c index 7871b5e63..135f0dcf6 100644 --- a/lib/msun/src/e_log10.c +++ b/lib/msun/src/e_log10.c @@ -14,45 +14,18 @@ #include __FBSDID("$FreeBSD$"); -/* __ieee754_log10(x) - * Return the base 10 logarithm of x - * - * Method : - * Let log10_2hi = leading 40 bits of log10(2) and - * log10_2lo = log10(2) - log10_2hi, - * ivln10 = 1/log(10) rounded. - * Then - * n = ilogb(x), - * if(n<0) n = n+1; - * x = scalbn(x,-n); - * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x)) - * - * Note 1: - * To guarantee log10(10**n)=n, where 10**n is normal, the rounding - * mode must set to Round-to-Nearest. - * Note 2: - * [1/log(10)] rounded to 53 bits has error .198 ulps; - * log10 is monotonic at all binary break points. - * - * Special cases: - * log10(x) is NaN with signal if x < 0; - * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; - * log10(NaN) is that NaN with no signal; - * log10(10**N) = N for N=0,1,...,22. - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. +/* + * Return the base 10 logarithm of x. See k_log.c for details on the algorithm. */ #include "math.h" #include "math_private.h" +#include "k_log.h" static const double two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ -ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */ +ivln10hi = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */ +ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */ log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ @@ -61,7 +34,7 @@ static const double zero = 0.0; double __ieee754_log10(double x) { - double y,z; + double f,hi,lo,y,z; int32_t i,k,hx; u_int32_t lx; @@ -77,10 +50,15 @@ __ieee754_log10(double x) } if (hx >= 0x7ff00000) return x+x; k += (hx>>20)-1023; - i = ((u_int32_t)k&0x80000000)>>31; - hx = (hx&0x000fffff)|((0x3ff-i)<<20); - y = (double)(k+i); - SET_HIGH_WORD(x,hx); - z = y*log10_2lo + ivln10*__ieee754_log(x); + 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; + 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; } diff --git a/lib/msun/src/e_log10f.c b/lib/msun/src/e_log10f.c index a0b1618a4..940b8319f 100644 --- a/lib/msun/src/e_log10f.c +++ b/lib/msun/src/e_log10f.c @@ -1,7 +1,3 @@ -/* e_log10f.c -- float version of e_log10.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -16,12 +12,18 @@ #include __FBSDID("$FreeBSD$"); +/* + * Return the base 10 logarithm of x. See k_log.c for details on the algorithm. + */ + #include "math.h" #include "math_private.h" +#include "k_logf.h" static const float two25 = 3.3554432000e+07, /* 0x4c000000 */ -ivln10 = 4.3429449201e-01, /* 0x3ede5bd9 */ +ivln10hi = 4.3432617188e-01, /* 0x3ede6000 */ +ivln10lo = -3.1689971365e-05, /* 0xb804ead9 */ log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ log10_2lo = 7.9034151668e-07; /* 0x355427db */ @@ -30,7 +32,7 @@ static const float zero = 0.0; float __ieee754_log10f(float x) { - float y,z; + float f,hi,lo,y,z; int32_t i,k,hx; GET_FLOAT_WORD(hx,x); @@ -45,10 +47,16 @@ __ieee754_log10f(float x) } if (hx >= 0x7f800000) return x+x; k += (hx>>23)-127; - i = ((u_int32_t)k&0x80000000)>>31; - hx = (hx&0x007fffff)|((0x7f-i)<<23); - y = (float)(k+i); - SET_FLOAT_WORD(x,hx); - z = y*log10_2lo + ivln10*__ieee754_logf(x); + 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); + 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; } diff --git a/lib/msun/src/e_log2.c b/lib/msun/src/e_log2.c index 8a729d5c3..6cf3dbc19 100644 --- a/lib/msun/src/e_log2.c +++ b/lib/msun/src/e_log2.c @@ -14,8 +14,8 @@ #include __FBSDID("$FreeBSD$"); -/* log2(x) - * Return the base 2 logarithm of x. +/* + * Return the base 2 logarithm of x. See k_log.c for details on the algorithm. */ #include "math.h" @@ -24,8 +24,8 @@ __FBSDID("$FreeBSD$"); static const double two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ -ivln2hi = 0x1.71547652000p+0, -ivln2lo = 0x1.705fc2eefa2p-33; +ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */ +ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */ static const double zero = 0.0; diff --git a/lib/msun/src/e_log2f.c b/lib/msun/src/e_log2f.c index 6b2b966fe..bb308d3e8 100644 --- a/lib/msun/src/e_log2f.c +++ b/lib/msun/src/e_log2f.c @@ -12,14 +12,18 @@ #include __FBSDID("$FreeBSD$"); +/* + * Return the base 2 logarithm of x. See k_log.c for details on the algorithm. + */ + #include "math.h" #include "math_private.h" #include "k_logf.h" static const float two25 = 3.3554432000e+07, /* 0x4c000000 */ -ivln2hi = 0x1.716p+0f, -ivln2lo = -0x1.7135a8fa03d11p-13; +ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */ +ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */ static const float zero = 0.0; @@ -46,7 +50,7 @@ __ieee754_log2f(float x) SET_FLOAT_WORD(x,hx|(i^0x3f800000)); /* normalize x or x/2 */ k += (i>>23); f = __kernel_logf(x); - x = x - 1; + x = x - (float)1.0; GET_FLOAT_WORD(hx,x); SET_FLOAT_WORD(hi,hx&0xfffff000); lo = x - hi; diff --git a/lib/msun/src/k_logf.h b/lib/msun/src/k_logf.h index d11a56306..d9f0f3ddf 100644 --- a/lib/msun/src/k_logf.h +++ b/lib/msun/src/k_logf.h @@ -12,8 +12,8 @@ #include __FBSDID("$FreeBSD$"); -/* __kernel_logf(x) - * Return log(x) - (x-1) for x in ~[sqrt(2)/2, sqrt(2)]. +/* + * float version of __kernel_log(x). See k_log.c for details. */ static const float @@ -33,7 +33,7 @@ __kernel_logf(float x) f = x-(float)1.0; if((0x007fffff&(0x8000+ix))<0xc000) { /* -2**-9 <= f < 2**-9 */ - if(f==0.0) return 0.0; + if(f==0.0f) return 0.0f; return f*f*((float)0.33333333333333333*f-(float)0.5); } s = f/((float)2.0+f); -- 2.45.0