From owner-svn-src-all@FreeBSD.ORG Mon Oct 17 05:38:22 2011 Return-Path: Delivered-To: svn-src-all@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id BEB9D1065673; Mon, 17 Oct 2011 05:38:22 +0000 (UTC) (envelope-from das@FreeBSD.org) Received: from svn.freebsd.org (svn.freebsd.org [IPv6:2001:4f8:fff6::2c]) by mx1.freebsd.org (Postfix) with ESMTP id AC07D8FC0A; Mon, 17 Oct 2011 05:38:22 +0000 (UTC) Received: from svn.freebsd.org (localhost [127.0.0.1]) by svn.freebsd.org (8.14.4/8.14.4) with ESMTP id p9H5cM7p096784; Mon, 17 Oct 2011 05:38:22 GMT (envelope-from das@svn.freebsd.org) Received: (from das@localhost) by svn.freebsd.org (8.14.4/8.14.4/Submit) id p9H5cM5U096770; Mon, 17 Oct 2011 05:38:22 GMT (envelope-from das@svn.freebsd.org) Message-Id: <201110170538.p9H5cM5U096770@svn.freebsd.org> From: David Schultz Date: Mon, 17 Oct 2011 05:38:22 +0000 (UTC) To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-stable@freebsd.org, svn-src-stable-7@freebsd.org X-SVN-Group: stable-7 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Cc: Subject: svn commit: r226457 - in stable/7/lib/msun: . man src X-BeenThere: svn-src-all@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "SVN commit messages for the entire src tree \(except for " user" and " projects" \)" List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 17 Oct 2011 05:38:22 -0000 Author: das Date: Mon Oct 17 05:38:22 2011 New Revision: 226457 URL: http://svn.freebsd.org/changeset/base/226457 Log: Belated MFC of log2 and log2f to 7-STABLE by popular demand, along with other log* improvements. r175461 - log(3) manpage (partial MFC) r216210 - add k_log (used by log*) r216211 - add log2 and log2f 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 Added: stable/7/lib/msun/man/log.3 - copied, changed from r175461, head/lib/msun/man/log.3 stable/7/lib/msun/src/e_log2.c - copied, changed from r216211, head/lib/msun/src/e_log2.c stable/7/lib/msun/src/e_log2f.c - copied, changed from r216211, head/lib/msun/src/e_log2f.c stable/7/lib/msun/src/k_log.h - copied, changed from r216210, head/lib/msun/src/k_log.h stable/7/lib/msun/src/k_logf.h - copied, changed from r216210, head/lib/msun/src/k_logf.h Modified: stable/7/lib/msun/Makefile stable/7/lib/msun/Symbol.map stable/7/lib/msun/man/exp.3 stable/7/lib/msun/man/math.3 stable/7/lib/msun/src/e_log10.c stable/7/lib/msun/src/e_log10f.c stable/7/lib/msun/src/math.h stable/7/lib/msun/src/math_private.h Directory Properties: stable/7/lib/msun/ (props changed) Modified: stable/7/lib/msun/Makefile ============================================================================== --- stable/7/lib/msun/Makefile Mon Oct 17 05:38:07 2011 (r226456) +++ stable/7/lib/msun/Makefile Mon Oct 17 05:38:22 2011 (r226457) @@ -33,7 +33,8 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c e_expf.c e_fmod.c e_fmodf.c e_gamma.c e_gamma_r.c e_gammaf.c \ e_gammaf_r.c e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \ e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \ - e_log.c e_log10.c e_log10f.c e_logf.c e_pow.c e_powf.c e_rem_pio2.c \ + e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \ + e_pow.c e_powf.c e_rem_pio2.c \ e_rem_pio2f.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \ e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \ k_cos.c k_cosf.c k_rem_pio2.c k_sin.c k_sinf.c \ @@ -97,7 +98,7 @@ MAN= acos.3 acosh.3 asin.3 asinh.3 atan. feclearexcept.3 feenableexcept.3 fegetenv.3 \ fegetround.3 fenv.3 floor.3 \ fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \ - lgamma.3 lrint.3 lround.3 math.3 nextafter.3 remainder.3 rint.3 \ + lgamma.3 log.3 lrint.3 lround.3 math.3 nextafter.3 remainder.3 rint.3 \ round.3 scalbn.3 signbit.3 sin.3 sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 MLINKS+=acos.3 acosf.3 @@ -115,7 +116,7 @@ MLINKS+=copysign.3 copysignf.3 copysign. MLINKS+=cos.3 cosf.3 MLINKS+=cosh.3 coshf.3 MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 -MLINKS+=exp.3 expm1.3 exp.3 log.3 exp.3 log10.3 exp.3 log1p.3 exp.3 pow.3 \ +MLINKS+=exp.3 expm1.3 exp.3 log10.3 exp.3 log1p.3 exp.3 pow.3 \ exp.3 exp2.3 exp.3 exp2f.3 exp.3 expf.3 \ exp.3 expm1f.3 exp.3 logf.3 exp.3 powf.3 \ exp.3 log10f.3 exp.3 log1pf.3 @@ -140,6 +141,7 @@ MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl. MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3 MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 lgammaf.3 lgamma.3 tgamma.3 +MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log1p.3 log.3 log1pf.3 log.3 logf.3 log.3 log2.3 log.3 log2f.3 MLINKS+=lrint.3 llrint.3 lrint.3 llrintf.3 lrint.3 lrintf.3 MLINKS+=lround.3 llround.3 lround.3 llroundf.3 lround.3 llroundl.3 \ lround.3 lroundf.3 lround.3 lroundl.3 Modified: stable/7/lib/msun/Symbol.map ============================================================================== --- stable/7/lib/msun/Symbol.map Mon Oct 17 05:38:07 2011 (r226456) +++ stable/7/lib/msun/Symbol.map Mon Oct 17 05:38:22 2011 (r226457) @@ -181,3 +181,9 @@ FBSD_1.0 { drem; dremf; }; + +/* First added in 9.0-CURRENT */ +FBSD_1.2 { + log2; + log2f; +}; Modified: stable/7/lib/msun/man/exp.3 ============================================================================== --- stable/7/lib/msun/man/exp.3 Mon Oct 17 05:38:07 2011 (r226456) +++ stable/7/lib/msun/man/exp.3 Mon Oct 17 05:38:22 2011 (r226457) @@ -28,7 +28,7 @@ .\" from: @(#)exp.3 6.12 (Berkeley) 7/31/91 .\" $FreeBSD$ .\" -.Dd April 5, 2005 +.Dd January 17, 2008 .Dt EXP 3 .Os .Sh NAME @@ -39,15 +39,9 @@ .Nm exp2f , .Nm expm1 , .Nm expm1f , -.Nm log , -.Nm logf , -.Nm log10 , -.Nm log10f , -.Nm log1p , -.Nm log1pf , .Nm pow , .Nm powf -.Nd exponential, logarithm, power functions +.Nd exponential and power functions .Sh LIBRARY .Lb libm .Sh SYNOPSIS @@ -65,18 +59,6 @@ .Ft float .Fn expm1f "float x" .Ft double -.Fn log "double x" -.Ft float -.Fn logf "float x" -.Ft double -.Fn log10 "double x" -.Ft float -.Fn log10f "float x" -.Ft double -.Fn log1p "double x" -.Ft float -.Fn log1pf "float x" -.Ft double .Fn pow "double x" "double y" .Ft float .Fn powf "float x" "float y" @@ -92,7 +74,7 @@ exponential value of the given argument .Pp The .Fn exp2 -and the +and .Fn exp2f functions compute the base 2 exponential of the given argument .Fa x . @@ -105,29 +87,6 @@ functions compute the value exp(x)\-1 ac .Fa x . .Pp The -.Fn log -and the -.Fn logf -functions compute the value of the natural logarithm of argument -.Fa x . -.Pp -The -.Fn log10 -and the -.Fn log10f -functions compute the value of the logarithm of argument -.Fa x -to base 10. -.Pp -The -.Fn log1p -and the -.Fn log1pf -functions compute -the value of log(1+x) accurately even for tiny argument -.Fa x . -.Pp -The .Fn pow and the .Fn powf @@ -159,30 +118,7 @@ raise an invalid exception and return an < 0 and .Fa y is not an integer. -An attempt to take the logarithm of \*(Pm0 will result in -a divide-by-zero exception, and an infinity will be returned. -An attempt to take the logarithm of a negative number will -result in an invalid exception, and an \*(Na will be generated. .Sh NOTES -The functions exp(x)\-1 and log(1+x) are called -expm1 and logp1 in -.Tn BASIC -on the Hewlett\-Packard -.Tn HP Ns \-71B -and -.Tn APPLE -Macintosh, -.Tn EXP1 -and -.Tn LN1 -in Pascal, exp1 and log1 in C -on -.Tn APPLE -Macintoshes, where they have been provided to make -sure financial calculations of ((1+x)**n\-1)/x, namely -expm1(n\(**log1p(x))/x, will be accurate when x is tiny. -They also provide accurate inverse hyperbolic functions. -.Pp The function .Fn pow x 0 returns x**0 = 1 for all x including x = 0, \*(If, and \*(Na . @@ -229,4 +165,9 @@ and infinite x, i.e., independently of x .El .Sh SEE ALSO .Xr fenv 3 , +.Xr ldexp 3 , +.Xr log 3 , .Xr math 3 +.Sh STANDARDS +These functions conform to +.St -isoC-99 . Copied and modified: stable/7/lib/msun/man/log.3 (from r175461, head/lib/msun/man/log.3) ============================================================================== --- head/lib/msun/man/log.3 Fri Jan 18 21:43:00 2008 (r175461, copy source) +++ stable/7/lib/msun/man/log.3 Mon Oct 17 05:38:22 2011 (r226457) @@ -1,4 +1,4 @@ -.\" Copyright (c) 2008 David Schultz +.\" Copyright (c) 2008-2010 David Schultz .\" All rights reserved. .\" .\" Redistribution and use in source and binary forms, with or without @@ -24,15 +24,16 @@ .\" .\" $FreeBSD$ .\" -.Dd January 17, 2008 +.Dd December 5, 2010 .Dt LOG 3 .Os .Sh NAME .Nm log , .Nm logf , -.Nm logl , .Nm log10 , .Nm log10f , +.Nm log2 , +.Nm log2f , .Nm log1p , .Nm log1pf .Nd logarithm functions @@ -49,6 +50,10 @@ .Ft float .Fn log10f "float x" .Ft double +.Fn log2 "double x" +.Ft float +.Fn log2f "float x" +.Ft double .Fn log1p "double x" .Ft float .Fn log1pf "float x" @@ -65,6 +70,12 @@ The and .Fn log10f functions compute the logarithm base 10 of +.Fa x , +while +.Fn log2 +and +.Fn log2f +compute the logarithm base 2 of .Fa x . .Pp The @@ -97,6 +108,8 @@ The .Fn logf , .Fn log10 , .Fn log10f , +.Fn log2 , +.Fn log2f , .Fn log1p , and .Fn log1pf Modified: stable/7/lib/msun/man/math.3 ============================================================================== --- stable/7/lib/msun/man/math.3 Mon Oct 17 05:38:07 2011 (r226456) +++ stable/7/lib/msun/man/math.3 Mon Oct 17 05:38:22 2011 (r226457) @@ -28,7 +28,7 @@ .\" from: @(#)math.3 6.10 (Berkeley) 5/6/91 .\" $FreeBSD$ .\" -.Dd November 6, 2005 +.Dd December 5, 2010 .Dt MATH 3 .Os .if n \{\ @@ -180,7 +180,7 @@ lgamma log gamma function log natural logarithm log10 logarithm to base 10 log1p log(1+x) -.\" log2 base 2 logarithm +log2 base 2 logarithm pow exponential x**y sin trigonometric function sinh hyperbolic function @@ -192,8 +192,8 @@ y1 Bessel function of the second kind of yn Bessel function of the second kind of the order n .El .Pp -Unlike the algebraic functions listed earlier, the routines -in this section may not produce a result that is correctly rounded, +The routines +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. @@ -218,20 +218,18 @@ and values, were written for or imported into subsequent versions of FreeBSD. .Sh BUGS The -.Fn log2 -and .Fn nan -functions are missing, and many functions are not available in their +function is missing, and many functions are not available in their .Vt "long double" variants. .Pp Many of the routines to compute transcendental functions produce inaccurate results in other than the default rounding mode. .Pp -On some architectures, trigonometric argument reduction is not -performed accurately, resulting in errors greater than 1 -.Em ulp -for large arguments to +On the i386 platform, trigonometric argument reduction is not +performed accurately for huge arguments, resulting in +large errors +for such arguments to .Fn cos , .Fn sin , and Modified: stable/7/lib/msun/src/e_log10.c ============================================================================== --- stable/7/lib/msun/src/e_log10.c Mon Oct 17 05:38:07 2011 (r226456) +++ stable/7/lib/msun/src/e_log10.c Mon Oct 17 05:38:22 2011 (r226457) @@ -15,45 +15,22 @@ static char rcsid[] = "$FreeBSD$"; #endif -/* __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)) +/* + * Return the base 10 logarithm of x. See e_log.c and k_log.h for most + * comments. * - * 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. + * log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2) + * in not-quite-routine extra precision. */ #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 */ @@ -62,26 +39,50 @@ static const double zero = 0.0; double __ieee754_log10(double x) { - double y,z; + double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2; int32_t i,k,hx; u_int32_t lx; EXTRACT_WORDS(hx,lx,x); - k=0; - if (hx < 0x00100000) { /* x < 2**-1022 */ - if (((hx&0x7fffffff)|lx)==0) - return -two54/zero; /* log(+-0)=-inf */ - if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ - k -= 54; x *= two54; /* subnormal number, scale up x */ + k=0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx&0x7fffffff)|lx)==0) + return -two54/zero; /* log(+-0)=-inf */ + if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ + k -= 54; x *= two54; /* subnormal number, scale up 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; - 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); - return z+y*log10_2hi; + 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 = 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 = (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; } Modified: stable/7/lib/msun/src/e_log10f.c ============================================================================== --- stable/7/lib/msun/src/e_log10f.c Mon Oct 17 05:38:07 2011 (r226456) +++ stable/7/lib/msun/src/e_log10f.c Mon Oct 17 05:38:22 2011 (r226457) @@ -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. @@ -17,12 +13,18 @@ static char rcsid[] = "$FreeBSD$"; #endif +/* + * Float version of e_log10.c. See the latter for most comments. + */ + #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 */ @@ -31,25 +33,40 @@ static const float zero = 0.0; float __ieee754_log10f(float x) { - float y,z; + float f,hfsq,hi,lo,r,y,y2; int32_t i,k,hx; GET_FLOAT_WORD(hx,x); - k=0; - if (hx < 0x00800000) { /* x < 2**-126 */ - if ((hx&0x7fffffff)==0) - return -two25/zero; /* log(+-0)=-inf */ - if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ - k -= 25; x *= two25; /* subnormal number, scale up x */ + k=0; + if (hx < 0x00800000) { /* x < 2**-126 */ + if ((hx&0x7fffffff)==0) + return -two25/zero; /* log(+-0)=-inf */ + if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ + k -= 25; x *= two25; /* subnormal number, scale up x */ GET_FLOAT_WORD(hx,x); - } + } if (hx >= 0x7f800000) return x+x; + if (hx == 0x3f800000) + return zero; /* log(1) = +0 */ 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); - return z+y*log10_2hi; + 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 = 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 = (f - hi) - hfsq + r; + return y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi + hi*ivln10hi + + y*log10_2hi; } Copied and modified: stable/7/lib/msun/src/e_log2.c (from r216211, head/lib/msun/src/e_log2.c) ============================================================================== --- head/lib/msun/src/e_log2.c Sun Dec 5 22:11:22 2010 (r216211, copy source) +++ stable/7/lib/msun/src/e_log2.c Mon Oct 17 05:38:22 2011 (r226457) @@ -14,8 +14,14 @@ #include __FBSDID("$FreeBSD$"); -/* log2(x) - * Return the base 2 logarithm of x. +/* + * 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" @@ -24,37 +30,81 @@ __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; 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; EXTRACT_WORDS(hx,lx,x); - k=0; - if (hx < 0x00100000) { /* x < 2**-1022 */ - if (((hx&0x7fffffff)|lx)==0) - return -two54/zero; /* log(+-0)=-inf */ - if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ - k -= 54; x *= two54; /* subnormal number, scale up x */ + k=0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx&0x7fffffff)|lx)==0) + return -two54/zero; /* log(+-0)=-inf */ + if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ + k -= 54; x *= two54; /* subnormal number, scale up 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; } Copied and modified: stable/7/lib/msun/src/e_log2f.c (from r216211, head/lib/msun/src/e_log2f.c) ============================================================================== --- head/lib/msun/src/e_log2f.c Sun Dec 5 22:11:22 2010 (r216211, copy source) +++ stable/7/lib/msun/src/e_log2f.c Mon Oct 17 05:38:22 2011 (r226457) @@ -12,43 +12,70 @@ #include __FBSDID("$FreeBSD$"); +/* + * Float version of e_log2.c. See the latter for most comments. + */ + #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; 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); - k=0; - if (hx < 0x00800000) { /* x < 2**-126 */ - if ((hx&0x7fffffff)==0) - return -two25/zero; /* log(+-0)=-inf */ - if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ - k -= 25; x *= two25; /* subnormal number, scale up x */ + k=0; + if (hx < 0x00800000) { /* x < 2**-126 */ + if ((hx&0x7fffffff)==0) + return -two25/zero; /* log(+-0)=-inf */ + if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ + k -= 25; x *= two25; /* subnormal number, scale up 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 - 1; - 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; } Copied and modified: stable/7/lib/msun/src/k_log.h (from r216210, head/lib/msun/src/k_log.h) ============================================================================== --- head/lib/msun/src/k_log.h Sun Dec 5 22:11:03 2010 (r216210, copy source) +++ stable/7/lib/msun/src/k_log.h Mon Oct 17 05:38:22 2011 (r226457) @@ -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, /* 3FC3 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; - - EXTRACT_WORDS(hx,lx,x); + double hfsq,s,z,R,w,t1,t2; - 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); } Copied and modified: stable/7/lib/msun/src/k_logf.h (from r216210, head/lib/msun/src/k_logf.h) ============================================================================== --- head/lib/msun/src/k_logf.h Sun Dec 5 22:11:03 2010 (r216210, copy source) +++ stable/7/lib/msun/src/k_logf.h Mon Oct 17 05:38:22 2011 (r226457) @@ -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 k_log.h. See the latter for most comments. */ static const float @@ -24,32 +24,16 @@ Lg3 = 0x91e9ee.0p-25, /* 0.28498786 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.0) return 0.0; - 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); } Modified: stable/7/lib/msun/src/math.h ============================================================================== --- stable/7/lib/msun/src/math.h Mon Oct 17 05:38:07 2011 (r226456) +++ stable/7/lib/msun/src/math.h Mon Oct 17 05:38:22 2011 (r226457) @@ -235,6 +235,7 @@ double lgamma(double); long long llrint(double); long long llround(double); double log1p(double); +double log2(double); double logb(double); long lrint(double); long lround(double); @@ -314,6 +315,7 @@ int ilogbf(float) __pure2; float ldexpf(float, int); float log10f(float); float log1pf(float); +float log2f(float); float logf(float); float modff(float, float *); /* fundamentally !__pure2 */ Modified: stable/7/lib/msun/src/math_private.h ============================================================================== --- stable/7/lib/msun/src/math_private.h Mon Oct 17 05:38:07 2011 (r226456) +++ stable/7/lib/msun/src/math_private.h Mon Oct 17 05:38:22 2011 (r226457) @@ -206,6 +206,7 @@ cpackl(long double x, long double y) #define __ieee754_acos acos #define __ieee754_acosh acosh #define __ieee754_log log +#define __ieee754_log2 log2 #define __ieee754_atanh atanh #define __ieee754_asin asin #define __ieee754_atan2 atan2 @@ -244,6 +245,7 @@ cpackl(long double x, long double y) #define __ieee754_lgammaf_r lgammaf_r #define __ieee754_gammaf_r gammaf_r #define __ieee754_log10f log10f +#define __ieee754_log2f log2f #define __ieee754_sinhf sinhf #define __ieee754_hypotf hypotf #define __ieee754_j0f j0f