standards/83845: [ patch ] add log2() and log2f() support for libm
Roman Bogorodskiy
bogorodskiy at gmail.com
Thu Jul 21 12:01:26 GMT 2005
>Number: 83845
>Category: standards
>Synopsis: [ patch ] add log2() and log2f() support for libm
>Confidential: no
>Severity: non-critical
>Priority: medium
>Responsible: freebsd-standards
>State: open
>Quarter:
>Keywords:
>Date-Required:
>Class: change-request
>Submitter-Id: current-users
>Arrival-Date: Thu Jul 21 12:00:34 GMT 2005
>Closed-Date:
>Last-Modified:
>Originator: Roman Bogorodskiy
>Release: FreeBSD 7.0-CURRENT i386
>Organization:
>Environment:
System: FreeBSD lame.novel.ru 7.0-CURRENT FreeBSD 7.0-CURRENT #17: Wed Jul 20 11:32:19 MSD 2005 root at lame.novel.ru:/usr/obj/usr/src/sys/NOVEL i386
>Description:
Add e_log2.c and e_log2f.c files based on e_log.c and e_logf.c
respectively which adds support of C99 funtions log2() and log2f().
The idea is simple:
log(N) to base A is equal to logN/logA, that way, log2(x) = log(x)/log(2).
So after some simplifications and deviding result by precalculated value of
log(2) we can compute the value of the logarithm of argument x to base 2.
Patch is attached. You can also pick it up here:
http://people.freebsd.org/~novel/patches/freebsd/msun_log2_and_log2f.diff
>How-To-Repeat:
>Fix:
--- msun_log2_and_log2f.diff begins here ---
diff -ruN /usr/src.bak/lib/msun/Makefile lib/msun/Makefile
--- /usr/src.bak/lib/msun/Makefile Sun Jul 17 17:47:52 2005
+++ lib/msun/Makefile Thu Jul 21 15:00:25 2005
@@ -33,7 +33,7 @@
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_log2.c e_log10.c e_log10f.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_rem_pio2f.c k_sin.c k_sinf.c \
@@ -114,7 +114,8 @@
MLINKS+=exp.3 expm1.3 exp.3 log.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
+ exp.3 log10f.3 exp.3 log1pf.3 \
+ exp.3 log2.3 exp.3 log2f.3
MLINKS+=fabs.3 fabsf.3 fabs.3 fabsl.3
MLINKS+=fdim.3 fdimf.3 fdim.3 fdiml.3
MLINKS+=feclearexcept.3 fegetexceptflag.3 feclearexcept.3 feraiseexcept.3 \
diff -ruN /usr/src.bak/lib/msun/man/exp.3 lib/msun/man/exp.3
--- /usr/src.bak/lib/msun/man/exp.3 Sun Jul 17 17:47:50 2005
+++ lib/msun/man/exp.3 Thu Jul 21 15:00:54 2005
@@ -45,6 +45,8 @@
.Nm expm1f ,
.Nm log ,
.Nm logf ,
+.Nm log2,
+.Nm log2f,
.Nm log10 ,
.Nm log10f ,
.Nm log1p ,
@@ -73,6 +75,10 @@
.Ft float
.Fn logf "float x"
.Ft double
+.Fn log2 "double x"
+.Ft float
+.Fn log2f "float x"
+.Ft double
.Fn log10 "double x"
.Ft float
.Fn log10f "float x"
@@ -114,6 +120,14 @@
.Fn logf
functions compute the value of the natural logarithm of argument
.Fa x .
+.Pp
+The
+.Fn log2
+and the
+.Fn log2f
+functions compute the value of the logarithm of argument
+.Fa x
+to base 2.
.Pp
The
.Fn log10
diff -ruN /usr/src.bak/lib/msun/src/e_log2.c lib/msun/src/e_log2.c
--- /usr/src.bak/lib/msun/src/e_log2.c Thu Jan 1 03:00:00 1970
+++ lib/msun/src/e_log2.c Thu Jul 21 15:12:49 2005
@@ -0,0 +1,79 @@
+
+/* @(#)e_log.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.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD$";
+#endif
+
+#include "math.h"
+#include "math_private.h"
+
+static const double
+ln2 = 0.6931471805599452862268,
+two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
+Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
+Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
+Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
+Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
+Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
+Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
+Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
+
+static const double zero = 0.0;
+
+double
+__ieee754_log2(double x)
+{
+ double hfsq,f,s,z,R,w,t1,t2,dk;
+ int32_t k,hx,i,j;
+ 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 = x-1.0;
+ dk = (double)k;
+ if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
+ if (f==zero)
+ return (dk);
+ R = f*f*(0.5-0.33333333333333333*f);
+ return (dk-(R-f)/ln2);
+ }
+ s = f/(2.0+f);
+ z = s*s;
+ 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;
+ R = t2+t1;
+ if(i>0) {
+ hfsq=0.5*f*f;
+ return (dk-(hfsq-s*(hfsq+R)-f)/ln2);
+ } else
+ return (dk-((s*(f-R))-f)/ln2);
+}
diff -ruN /usr/src.bak/lib/msun/src/e_log2f.c lib/msun/src/e_log2f.c
--- /usr/src.bak/lib/msun/src/e_log2f.c Thu Jan 1 03:00:00 1970
+++ lib/msun/src/e_log2f.c Thu Jul 21 14:43:57 2005
@@ -0,0 +1,80 @@
+/* e_logf.c -- float version of e_log.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian at cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD$";
+#endif
+
+#include "math.h"
+#include "math_private.h"
+
+static const float
+ln2 = 0.6931471805599452862268,
+two25 = 3.355443200e+07, /* 0x4c000000 */
+Lg1 = 6.6666668653e-01, /* 3F2AAAAB */
+Lg2 = 4.0000000596e-01, /* 3ECCCCCD */
+Lg3 = 2.8571429849e-01, /* 3E924925 */
+Lg4 = 2.2222198546e-01, /* 3E638E29 */
+Lg5 = 1.8183572590e-01, /* 3E3A3325 */
+Lg6 = 1.5313838422e-01, /* 3E1CD04F */
+Lg7 = 1.4798198640e-01; /* 3E178897 */
+
+static const float zero = 0.0;
+
+float
+__ieee754_log2f(float x)
+{
+ float hfsq,f,s,z,R,w,t1,t2,dk;
+ int32_t k,ix,i,j;
+
+ GET_FLOAT_WORD(ix,x);
+
+ k=0;
+ if (ix < 0x00800000) { /* x < 2**-126 */
+ if ((ix&0x7fffffff)==0)
+ return -two25/zero; /* log(+-0)=-inf */
+ if (ix<0) return (x-x)/zero; /* log(-#) = NaN */
+ k -= 25; x *= two25; /* subnormal number, scale up x */
+ GET_FLOAT_WORD(ix,x);
+ }
+ if (ix >= 0x7f800000) return x+x;
+ k += (ix>>23)-127;
+ ix &= 0x007fffff;
+ i = (ix+(0x95f64<<3))&0x800000;
+ SET_FLOAT_WORD(x,ix|(i^0x3f800000)); /* normalize x or x/2 */
+ k += (i>>23);
+ dk = (float)k;
+ f = x-(float)1.0;
+ if((0x007fffff&(15+ix))<16) { /* |f| < 2**-20 */
+ if (f==zero)
+ return (dk);
+ R = f*f*((float)0.5-(float)0.33333333333333333*f);
+ return (dk-(R-f)/ln2);
+ }
+ s = f/((float)2.0+f);
+ z = s*s;
+ i = ix-(0x6147a<<3);
+ w = z*z;
+ j = (0x6b851<<3)-ix;
+ t1= w*(Lg2+w*(Lg4+w*Lg6));
+ t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
+ i |= j;
+ R = t2+t1;
+ if(i>0) {
+ hfsq=(float)0.5*f*f;
+ return (dk-(hfsq-s*(hfsq+R)-f)/ln2);
+ } else
+ return (dk-((s*(f-R))-f)/ln2);
+}
diff -ruN /usr/src.bak/lib/msun/src/math.h lib/msun/src/math.h
--- /usr/src.bak/lib/msun/src/math.h Sun Jul 17 17:47:52 2005
+++ lib/msun/src/math.h Thu Jul 21 14:25:21 2005
@@ -203,6 +203,7 @@
double frexp(double, int *); /* fundamentally !__pure2 */
double ldexp(double, int);
double log(double);
+double log2(double);
double log10(double);
double modf(double, double *); /* fundamentally !__pure2 */
@@ -314,6 +315,7 @@
float ldexpf(float, int);
float log10f(float);
float log1pf(float);
+float log2f(float);
float logf(float);
float modff(float, float *); /* fundamentally !__pure2 */
diff -ruN /usr/src.bak/lib/msun/src/math_private.h lib/msun/src/math_private.h
--- /usr/src.bak/lib/msun/src/math_private.h Sun Jul 17 17:47:52 2005
+++ lib/msun/src/math_private.h Thu Jul 21 14:45:31 2005
@@ -176,6 +176,7 @@
#define __ieee754_lgamma_r lgamma_r
#define __ieee754_gamma_r gamma_r
#define __ieee754_log10 log10
+#define __ieee754_log2 log2
#define __ieee754_sinh sinh
#define __ieee754_hypot hypot
#define __ieee754_j0 j0
@@ -202,6 +203,7 @@
#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
--- msun_log2_and_log2f.diff ends here ---
>Release-Note:
>Audit-Trail:
>Unformatted:
More information about the freebsd-standards
mailing list