Date: Sun, 6 Mar 2011 08:49:44 +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-8@freebsd.org Subject: svn commit: r219323 - in stable/8/lib/msun: . man src Message-ID: <201103060849.p268nitH091290@svn.freebsd.org>
next in thread | raw e-mail | index | archive | help
Author: das Date: Sun Mar 6 08:49:44 2011 New Revision: 219323 URL: http://svn.freebsd.org/changeset/base/219323 Log: MFC r216210: refactor log(3) r216211: add log2(3) and log2f(3) Added: stable/8/lib/msun/src/e_log2.c - copied unchanged from r216211, head/lib/msun/src/e_log2.c stable/8/lib/msun/src/e_log2f.c - copied unchanged from r216211, head/lib/msun/src/e_log2f.c stable/8/lib/msun/src/k_log.h - copied unchanged from r216210, head/lib/msun/src/k_log.h stable/8/lib/msun/src/k_logf.h - copied unchanged from r216210, head/lib/msun/src/k_logf.h Modified: stable/8/lib/msun/Makefile stable/8/lib/msun/Symbol.map stable/8/lib/msun/man/log.3 stable/8/lib/msun/man/math.3 stable/8/lib/msun/src/math.h stable/8/lib/msun/src/math_private.h Directory Properties: stable/8/lib/msun/ (props changed) Modified: stable/8/lib/msun/Makefile ============================================================================== --- stable/8/lib/msun/Makefile Sun Mar 6 08:35:50 2011 (r219322) +++ stable/8/lib/msun/Makefile Sun Mar 6 08:49:44 2011 (r219323) @@ -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 \ @@ -164,7 +165,7 @@ MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0 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 \ Modified: stable/8/lib/msun/Symbol.map ============================================================================== --- stable/8/lib/msun/Symbol.map Sun Mar 6 08:35:50 2011 (r219322) +++ stable/8/lib/msun/Symbol.map Sun Mar 6 08:49:44 2011 (r219323) @@ -218,3 +218,10 @@ FBSD_1.1 { cprojf; cprojl; }; + +/* First added in 9.0-CURRENT */ +FBSD_1.2 { + __isnanf; + log2; + log2f; +}; Modified: stable/8/lib/msun/man/log.3 ============================================================================== --- stable/8/lib/msun/man/log.3 Sun Mar 6 08:35:50 2011 (r219322) +++ stable/8/lib/msun/man/log.3 Sun Mar 6 08:49:44 2011 (r219323) @@ -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 Modified: stable/8/lib/msun/man/math.3 ============================================================================== --- stable/8/lib/msun/man/math.3 Sun Mar 6 08:35:50 2011 (r219322) +++ stable/8/lib/msun/man/math.3 Sun Mar 6 08:49:44 2011 (r219323) @@ -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 .if n \{\ @@ -185,7 +185,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 @@ -197,7 +197,7 @@ 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 +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 @@ -224,18 +224,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 Copied: stable/8/lib/msun/src/e_log2.c (from r216211, head/lib/msun/src/e_log2.c) ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ stable/8/lib/msun/src/e_log2.c Sun Mar 6 08:49:44 2011 (r219323, copy of r216211, head/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; +} Copied: stable/8/lib/msun/src/e_log2f.c (from r216211, head/lib/msun/src/e_log2f.c) ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ stable/8/lib/msun/src/e_log2f.c Sun Mar 6 08:49:44 2011 (r219323, copy of r216211, head/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; +} Copied: stable/8/lib/msun/src/k_log.h (from r216210, head/lib/msun/src/k_log.h) ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ stable/8/lib/msun/src/k_log.h Sun Mar 6 08:49:44 2011 (r219323, copy of r216210, head/lib/msun/src/k_log.h) @@ -0,0 +1,116 @@ + +/* @(#)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. + * ==================================================== + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +/* __kernel_log(x) + * Return log(x) - (x-1) for x 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 + * term of the polynomial are done by the caller for increased accuracy + * when different bases are used. + * + * Method : + * 1. Argument Reduction: find k and f such that + * x = 2^k * (1+f), + * where sqrt(2)/2 < 1+f < sqrt(2) . + * + * 2. Approximation of log(1+f). + * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) + * = 2s + 2/3 s**3 + 2/5 s**5 + ....., + * = 2s + s*R + * We use a special Reme algorithm on [0,0.1716] to generate + * a polynomial of degree 14 to approximate R The maximum error + * of this polynomial approximation is bounded by 2**-58.45. In + * other words, + * 2 4 6 8 10 12 14 + * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s + * (the values of Lg1 to Lg7 are listed in the program) + * and + * | 2 14 | -58.45 + * | Lg1*s +...+Lg7*s - R(z) | <= 2 + * | | + * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. + * In order to guarantee error in log below 1ulp, we compute log + * by + * log(1+f) = f - s*(f - R) (if f is not too large) + * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) + * + * 3. Finally, log(x) = k*ln2 + log(1+f). + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * Here ln2 is split into two floating point number: + * ln2_hi + ln2_lo, + * where n*ln2_hi is always exact for |n| < 2000. + * + * Special cases: + * log(x) is NaN with signal if x < 0 (including -INF) ; + * log(+INF) is +INF; log(0) is -INF with signal; + * log(NaN) is that NaN with no signal. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * 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. + */ + +static const double +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 */ + +/* + * We always inline __kernel_log(), since doing so produces a + * substantial performance improvement (~40% on amd64). + */ +static inline double +__kernel_log(double x) +{ + double hfsq,f,s,z,R,w,t1,t2; + int32_t hx,i,j; + u_int32_t lx; + + EXTRACT_WORDS(hx,lx,x); + + 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); + 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; + R = t2+t1; + if (i>0) { + hfsq=0.5*f*f; + return s*(hfsq+R) - hfsq; + } else { + return s*(R-f); + } +} Copied: stable/8/lib/msun/src/k_logf.h (from r216210, head/lib/msun/src/k_logf.h) ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ stable/8/lib/msun/src/k_logf.h Sun Mar 6 08:49:44 2011 (r219323, copy of r216210, head/lib/msun/src/k_logf.h) @@ -0,0 +1,55 @@ +/* + * ==================================================== + * 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$"); + +/* __kernel_logf(x) + * Return log(x) - (x-1) for x in ~[sqrt(2)/2, sqrt(2)]. + */ + +static const float +/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ +Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ +Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ +Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ +Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ + +static inline float +__kernel_logf(float x) +{ + float hfsq,f,s,z,R,w,t1,t2; + int32_t ix,i,j; + + 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); + } +} Modified: stable/8/lib/msun/src/math.h ============================================================================== --- stable/8/lib/msun/src/math.h Sun Mar 6 08:35:50 2011 (r219322) +++ stable/8/lib/msun/src/math.h Sun Mar 6 08:49:44 2011 (r219323) @@ -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); @@ -318,6 +319,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/8/lib/msun/src/math_private.h ============================================================================== --- stable/8/lib/msun/src/math_private.h Sun Mar 6 08:35:50 2011 (r219322) +++ stable/8/lib/msun/src/math_private.h Sun Mar 6 08:49:44 2011 (r219323) @@ -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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201103060849.p268nitH091290>