Skip site navigation (1)Skip section navigation (2)
Date:      Mon, 17 Oct 2011 05:38:22 +0000 (UTC)
From:      David Schultz <das@FreeBSD.org>
To:        src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-stable@freebsd.org, svn-src-stable-7@freebsd.org
Subject:   svn commit: r226457 - in stable/7/lib/msun: . man src
Message-ID:  <201110170538.p9H5cM5U096770@svn.freebsd.org>

next in thread | raw e-mail | index | archive | help
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 <das@FreeBSD.org>
+.\" Copyright (c) 2008-2010 David Schultz <das@FreeBSD.org>
 .\" 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 <sys/cdefs.h>
 __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 <sys/cdefs.h>
 __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 <sys/cdefs.h>
 __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 <sys/cdefs.h>
 __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



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201110170538.p9H5cM5U096770>