diff options
Diffstat (limited to 'lib/msun')
| -rw-r--r-- | lib/msun/Makefile | 11 | ||||
| -rw-r--r-- | lib/msun/Symbol.map | 2 | ||||
| -rw-r--r-- | lib/msun/man/log.3 | 18 | ||||
| -rw-r--r-- | lib/msun/man/math.3 | 21 | ||||
| -rw-r--r-- | lib/msun/src/e_log2.c | 60 | ||||
| -rw-r--r-- | lib/msun/src/e_log2f.c | 58 | ||||
| -rw-r--r-- | lib/msun/src/k_log.h | 116 | ||||
| -rw-r--r-- | lib/msun/src/k_logf.h | 55 | ||||
| -rw-r--r-- | lib/msun/src/math.h | 2 | ||||
| -rw-r--r-- | lib/msun/src/math_private.h | 2 |
10 files changed, 331 insertions, 14 deletions
diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 87bd954..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 \ @@ -113,6 +114,12 @@ COMMON_SRCS:= ${COMMON_SRCS:N${i:R}.c} .endfor .endif +# Some files need certain gcc built-in functions to be disabled, since gcc's +# model of the functions bogusly assumes -fno-trapping-math. +XRINT_CFLAGS= -fno-builtin-rint -fno-builtin-rintf -fno-builtin-rintl +CFLAGS+= ${XRINT_CFLAGS} +XRINT_CFLAGS:= ${.IMPSRC:M*/s_nearbyint.c:C/^.+$/${XRINT_CFLAGS}/:C/^$//} + SRCS= ${COMMON_SRCS} ${ARCH_SRCS} INCS= fenv.h math.h @@ -169,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..39aa8ab 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,8 +194,8 @@ 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 -in this section may not produce a result that is correctly rounded, +The routines +in this section might 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 rarely, and then only in very-close-to-halfway cases. @@ -221,18 +221,19 @@ 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 -.Em ulp -for large arguments to +On the i386 platform, trigonometric argument reduction is not +performed accurately for huge arguments, resulting in +large errors +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..6cf3dbc --- /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$"); + +/* + * Return the base 2 logarithm of x. See k_log.c for details on the algorithm. + */ + +#include "math.h" +#include "math_private.h" +#include "k_log.h" + +static const double +two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ +ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */ +ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */ + +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..bb308d3 --- /dev/null +++ b/lib/msun/src/e_log2f.c @@ -0,0 +1,58 @@ +/* + * ==================================================== + * 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$"); + +/* + * Return the base 2 logarithm of x. See k_log.c for details on the algorithm. + */ + +#include "math.h" +#include "math_private.h" +#include "k_logf.h" + +static const float +two25 = 3.3554432000e+07, /* 0x4c000000 */ +ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */ +ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */ + +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 - (float)1.0; + 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/k_log.h b/lib/msun/src/k_log.h new file mode 100644 index 0000000..206355c --- /dev/null +++ b/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); + } +} diff --git a/lib/msun/src/k_logf.h b/lib/msun/src/k_logf.h new file mode 100644 index 0000000..d9f0f3d --- /dev/null +++ b/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$"); + +/* + * float version of __kernel_log(x). See k_log.c for details. + */ + +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.0f) return 0.0f; + 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); + } +} 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 |
