Add log2 and log2f.
authorPeter Avalos <pavalos@dragonflybsd.org>
Sat, 16 Jun 2007 22:26:53 +0000 (22:26 +0000)
committerPeter Avalos <pavalos@dragonflybsd.org>
Sat, 16 Jun 2007 22:26:53 +0000 (22:26 +0000)
Dragonfly-bug:  <http://bugs.dragonflybsd.org/issue632>
Obtained-from:  NetBSD

include/math.h
lib/libm/arch/i386/Makefile.inc
lib/libm/arch/i386/e_log2.S [new file with mode: 0644]
lib/libm/arch/i386/e_log2f.S [new file with mode: 0644]
lib/libm/man/exp.3
lib/libm/src/Makefile.inc
lib/libm/src/e_log2.c [new file with mode: 0644]
lib/libm/src/e_log2f.c [new file with mode: 0644]

index 0d0913b..04d399f 100644 (file)
@@ -1,5 +1,3 @@
-/*     $NetBSD: math.h,v 1.40 2005/02/03 04:39:32 perry Exp $  */
-
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -10,8 +8,8 @@
  * is preserved.
  * ====================================================
  *
- * $NetBSD: math.h,v 1.40 2005/02/03 04:39:32 perry Exp $
- * $DragonFly: src/include/math.h,v 1.8 2006/05/17 14:06:36 swildner Exp $
+ * $NetBSD: math.h,v 1.41 2005/07/21 12:56:29 christos Exp $
+ * $DragonFly: src/include/math.h,v 1.9 2007/06/16 22:26:53 pavalos Exp $
  */
 
 /*
@@ -180,6 +178,7 @@ double      exp(double);
 double frexp(double, int *);
 double ldexp(double, int);
 double log(double);
+double log2(double);
 double log10(double);
 double modf(double, double *);
 
@@ -264,6 +263,7 @@ float       frexpf(float, int *);
 int    ilogbf(float);
 float  ldexpf(float, int);
 float  logf(float);
+float  log2f(float);
 float  log10f(float);
 float  log1pf(float);
 float  logbf(float);
index e832ea8..efe236f 100644 (file)
@@ -1,9 +1,9 @@
-# $DragonFly: src/lib/libm/arch/i386/Makefile.inc,v 1.1 2005/07/26 21:15:19 joerg Exp $
+# $DragonFly: src/lib/libm/arch/i386/Makefile.inc,v 1.2 2007/06/16 22:26:53 pavalos Exp $
 
 .PATH: ${.CURDIR}/arch/i386
 
 SRCS+= e_acos.S e_asin.S e_atan2.S e_atan2f.S e_exp.S e_expf.S e_fmod.S \
-       e_log.S e_log10.S e_log10f.S e_logf.S e_remainder.S e_remainderf.S \
+       e_log.S e_log10.S e_log10f.S e_log2.S e_log2f.S e_logf.S e_remainder.S e_remainderf.S \
        e_scalb.S e_scalbf.S e_sqrt.S e_sqrtf.S lrint.S s_atan.S s_atanf.S \
        s_ceil.S s_ceilf.S s_copysign.S s_copysignf.S s_cos.S s_cosf.S \
        s_finite.S s_finitef.S s_floor.S s_floorf.S s_ilogb.S s_ilogbf.S \
diff --git a/lib/libm/arch/i386/e_log2.S b/lib/libm/arch/i386/e_log2.S
new file mode 100644 (file)
index 0000000..2214fa2
--- /dev/null
@@ -0,0 +1,18 @@
+/*
+ * Written by Rui Paulo <rpaulo@NetBSD.org>, based on e_log.S.
+ * Public domain.
+ * $NetBSD: e_log2.S,v 1.1 2005/07/21 20:58:21 rpaulo Exp $
+ * $DragonFly: src/lib/libm/arch/i386/e_log2.S,v 1.1 2007/06/16 22:26:53 pavalos Exp $
+ */
+
+#include <machine/asm.h>
+
+#include "abi.h"
+
+ENTRY(log2)
+       XMM_ONE_ARG_DOUBLE_PROLOGUE
+       fld1
+       fldl    ARG_DOUBLE_ONE
+       fyl2x
+       XMM_DOUBLE_EPILOGUE
+       ret
diff --git a/lib/libm/arch/i386/e_log2f.S b/lib/libm/arch/i386/e_log2f.S
new file mode 100644 (file)
index 0000000..9db9041
--- /dev/null
@@ -0,0 +1,18 @@
+/*
+ * Written by Rui Paulo <rpaulo@NetBSD.org>, based on e_logf.S.
+ * Public domain.
+ * $NetBSD: e_log2f.S,v 1.1 2005/07/21 20:58:21 rpaulo Exp $
+ * $DragonFly: src/lib/libm/arch/i386/e_log2f.S,v 1.1 2007/06/16 22:26:53 pavalos Exp $
+ */
+
+#include <machine/asm.h>
+
+#include "abi.h"
+
+ENTRY(log2f)
+       XMM_ONE_ARG_FLOAT_PROLOGUE
+       fld1
+       flds    ARG_FLOAT_ONE
+       fyl2x
+       XMM_FLOAT_EPILOGUE
+       ret
index ad5a297..0faf395 100644 (file)
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)exp.3        6.12 (Berkeley) 7/31/91
-.\"    $NetBSD: exp.3,v 1.21 2003/08/07 16:44:47 agc Exp $
-.\"    $DragonFly: src/lib/libm/man/exp.3,v 1.1 2005/07/26 21:15:20 joerg Exp $
+.\"    $NetBSD: exp.3,v 1.23 2005/07/21 12:58:22 wiz Exp $
+.\"    $DragonFly: src/lib/libm/man/exp.3,v 1.2 2007/06/16 22:26:53 pavalos Exp $
 .\"
-.Dd July 31, 1991
+.Dd July 21, 2005
 .Dt EXP 3
 .Os
 .Sh NAME
@@ -39,6 +39,8 @@
 .Nm expm1f ,
 .Nm log ,
 .Nm logf ,
+.Nm log2 ,
+.Nm log2f ,
 .Nm log10 ,
 .Nm log10f ,
 .Nm log1p ,
 .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"
@@ -103,6 +109,14 @@ the value of log(1+x) accurately even for tiny 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 pow
 computes the value
 of
index 775e7d0..505c648 100644 (file)
@@ -1,4 +1,4 @@
-# $DragonFly: src/lib/libm/src/Makefile.inc,v 1.1 2005/07/26 21:15:20 joerg Exp $
+# $DragonFly: src/lib/libm/src/Makefile.inc,v 1.2 2007/06/16 22:26:53 pavalos Exp $
 
 .PATH: ${.CURDIR}/src
 
@@ -17,7 +17,7 @@ MI_FUNCS+=    \
        e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c \
        e_exp.c e_expf.c e_fmod.c e_fmodf.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_r.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 \
        s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_cbrt.c s_cbrtf.c s_ceil.c \
diff --git a/lib/libm/src/e_log2.c b/lib/libm/src/e_log2.c
new file mode 100644 (file)
index 0000000..45d3460
--- /dev/null
@@ -0,0 +1,78 @@
+
+/* @(#)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.
+ * ====================================================
+ *
+ * $NetBSD: e_log2.c,v 1.1 2005/07/21 12:55:58 christos Exp $
+ * $DragonFly: src/lib/libm/src/e_log2.c,v 1.1 2007/06/16 22:26:53 pavalos Exp $
+ */
+
+#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
+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 --git a/lib/libm/src/e_log2f.c b/lib/libm/src/e_log2f.c
new file mode 100644 (file)
index 0000000..549db93
--- /dev/null
@@ -0,0 +1,79 @@
+/* e_logf.c -- float version of e_log.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@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.
+ * ====================================================
+ *
+ * $NetBSD: e_log2f.c,v 1.1 2005/07/21 12:55:58 christos Exp $
+ * $DragonFly: src/lib/libm/src/e_log2f.c,v 1.1 2007/06/16 22:26:53 pavalos Exp $
+ */
+
+#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
+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);
+}