From 54f91c64155d6a83781ff9b7b42b456c84869b96 Mon Sep 17 00:00:00 2001 From: Peter Avalos Date: Sat, 16 Jun 2007 22:26:53 +0000 Subject: [PATCH] Add log2 and log2f. Dragonfly-bug: Obtained-from: NetBSD --- include/math.h | 8 ++-- lib/libm/arch/i386/Makefile.inc | 4 +- lib/libm/arch/i386/e_log2.S | 18 ++++++++ lib/libm/arch/i386/e_log2f.S | 18 ++++++++ lib/libm/man/exp.3 | 20 +++++++-- lib/libm/src/Makefile.inc | 4 +- lib/libm/src/e_log2.c | 78 ++++++++++++++++++++++++++++++++ lib/libm/src/e_log2f.c | 79 +++++++++++++++++++++++++++++++++ 8 files changed, 218 insertions(+), 11 deletions(-) create mode 100644 lib/libm/arch/i386/e_log2.S create mode 100644 lib/libm/arch/i386/e_log2f.S create mode 100644 lib/libm/src/e_log2.c create mode 100644 lib/libm/src/e_log2f.c diff --git a/include/math.h b/include/math.h index 0d0913b855..04d399fc9a 100644 --- a/include/math.h +++ b/include/math.h @@ -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); diff --git a/lib/libm/arch/i386/Makefile.inc b/lib/libm/arch/i386/Makefile.inc index e832ea87ab..efe236f037 100644 --- a/lib/libm/arch/i386/Makefile.inc +++ b/lib/libm/arch/i386/Makefile.inc @@ -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 index 0000000000..2214fa2560 --- /dev/null +++ b/lib/libm/arch/i386/e_log2.S @@ -0,0 +1,18 @@ +/* + * Written by Rui Paulo , 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 + +#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 index 0000000000..9db904168b --- /dev/null +++ b/lib/libm/arch/i386/e_log2f.S @@ -0,0 +1,18 @@ +/* + * Written by Rui Paulo , 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 + +#include "abi.h" + +ENTRY(log2f) + XMM_ONE_ARG_FLOAT_PROLOGUE + fld1 + flds ARG_FLOAT_ONE + fyl2x + XMM_FLOAT_EPILOGUE + ret diff --git a/lib/libm/man/exp.3 b/lib/libm/man/exp.3 index ad5a2971a0..0faf395d68 100644 --- a/lib/libm/man/exp.3 +++ b/lib/libm/man/exp.3 @@ -26,10 +26,10 @@ .\" 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 , @@ -63,6 +65,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" @@ -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 diff --git a/lib/libm/src/Makefile.inc b/lib/libm/src/Makefile.inc index 775e7d0c9e..505c6489b3 100644 --- a/lib/libm/src/Makefile.inc +++ b/lib/libm/src/Makefile.inc @@ -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 index 0000000000..45d346057b --- /dev/null +++ b/lib/libm/src/e_log2.c @@ -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 +#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 index 0000000000..549db93bd5 --- /dev/null +++ b/lib/libm/src/e_log2f.c @@ -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 +#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); +} -- 2.41.0