Add log2() and log2f().

8 files changed, 149 insertions, 12 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 929af97..0048e25 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -45,7 +45,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 \
@@ -175,7 +176,7 @@ 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 lgamma.3 tgammaf.3
-MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log1p.3 log.3 log1pf.3 log.3 logf.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 llrintl.3 \
 	lrint.3 lrintf.3 lrint.3 lrintl.3
 MLINKS+=lround.3 llround.3 lround.3 llroundf.3 lround.3 llroundl.3 \
diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map
index 429a76f..69618e6 100644
--- a/lib/msun/Symbol.map
+++ b/lib/msun/Symbol.map
@@ -222,4 +222,6 @@ FBSD_1.1 {
 /* First added in 9.0-CURRENT */
 FBSD_1.2 {
 	__isnanf;
+	log2;
+	log2f;
 };
diff --git a/lib/msun/man/log.3 b/lib/msun/man/log.3
index 1934473..b9fd83c 100644
--- a/lib/msun/man/log.3
+++ b/lib/msun/man/log.3
@@ -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,7 +24,7 @@
 .\"
 .\" $FreeBSD$
 .\"
-.Dd January 17, 2008
+.Dd December 5, 2010
 .Dt LOG 3
 .Os
 .Sh NAME
@@ -33,6 +33,8 @@
 .Nm logl ,
 .Nm log10 ,
 .Nm log10f ,
+.Nm log2 ,
+.Nm log2f ,
 .Nm log1p ,
 .Nm log1pf
 .Nd logarithm functions
@@ -49,6 +51,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 +71,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 +109,8 @@ The
 .Fn logf ,
 .Fn log10 ,
 .Fn log10f ,
+.Fn log2 ,
+.Fn log2f ,
 .Fn log1p ,
 and
 .Fn log1pf
diff --git a/lib/msun/man/math.3 b/lib/msun/man/math.3
index 4690a9a..8e79ce5 100644
--- a/lib/msun/man/math.3
+++ b/lib/msun/man/math.3
@@ -28,7 +28,7 @@
 .\"	from: @(#)math.3	6.10 (Berkeley) 5/6/91
 .\" $FreeBSD$
 .\"
-.Dd December 16, 2007
+.Dd December 5, 2010
 .Dt MATH 3
 .Os
 .Sh NAME
@@ -182,7 +182,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
@@ -194,7 +194,7 @@ y1	Bessel function of the second kind of the order 1
 yn	Bessel function of the second kind of the order n
 .El
 .Pp
-Unlike the algebraic functions listed earlier, the routines
+The routines
 in this section may 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
@@ -221,18 +221,20 @@ and
 values, were written for or imported into subsequent versions of FreeBSD.
 .Sh BUGS
 The
-.Fn log2
-function is missing, and many functions are not available in their
+.Fn cbrt
+function and many of the transcendental 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
+On the i386 platform, trigonometric argument reduction is not
+performed accurately for very large arguments, resulting in errors
+greater than 1
 .Em ulp
-for large arguments to
+for such arguments to
 .Fn cos ,
 .Fn sin ,
 and
diff --git a/lib/msun/src/e_log2.c b/lib/msun/src/e_log2.c
new file mode 100644
index 0000000..8a729d5
--- /dev/null
+++ b/lib/msun/src/e_log2.c
@@ -0,0 +1,60 @@
+
+/* @(#)e_log10.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* log2(x)
+ * Return the base 2 logarithm of x.
+ */
+
+#include "math.h"
+#include "math_private.h"
+#include "k_log.h"
+
+static const double
+two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+ivln2hi    =  0x1.71547652000p+0,
+ivln2lo    =  0x1.705fc2eefa2p-33;
+
+static const double zero   =  0.0;
+
+double
+__ieee754_log2(double x)
+{
+	double f,hi,lo;
+	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 */
+	    GET_HIGH_WORD(hx,x);
+        }
+	if (hx >= 0x7ff00000) return x+x;
+	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;
+	SET_LOW_WORD(hi,0);
+	lo = x - hi;
+	return (x+f)*ivln2lo + (lo+f)*ivln2hi + hi*ivln2hi + k;
+}
diff --git a/lib/msun/src/e_log2f.c b/lib/msun/src/e_log2f.c
new file mode 100644
index 0000000..6b2b966
--- /dev/null
+++ b/lib/msun/src/e_log2f.c
@@ -0,0 +1,54 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#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;
+
+static const float zero   =  0.0;
+
+float
+__ieee754_log2f(float x)
+{
+	float f,hi,lo;
+	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 */
+	    GET_FLOAT_WORD(hx,x);
+        }
+	if (hx >= 0x7f800000) return x+x;
+	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);
+	SET_FLOAT_WORD(hi,hx&0xfffff000);
+	lo = x - hi;
+	return (x+f)*ivln2lo + (lo+f)*ivln2hi + hi*ivln2hi + k;
+}
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index 9109727..9aa2484 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -236,6 +236,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);
@@ -319,6 +320,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 */
 
diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
index 10c92f4..d79f808 100644
--- a/lib/msun/src/math_private.h
+++ b/lib/msun/src/math_private.h
@@ -292,6 +292,7 @@ irint(double x)
 #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
@@ -330,6 +331,7 @@ irint(double x)
 #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