summaryrefslogtreecommitdiffstats
path: root/lib/msun/src
diff options
context:
space:
mode:
Diffstat (limited to 'lib/msun/src')
-rw-r--r--lib/msun/src/e_log2.c60
-rw-r--r--lib/msun/src/e_log2f.c58
-rw-r--r--lib/msun/src/k_log.h116
-rw-r--r--lib/msun/src/k_logf.h55
-rw-r--r--lib/msun/src/math.h2
-rw-r--r--lib/msun/src/math_private.h2
6 files changed, 293 insertions, 0 deletions
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
OpenPOWER on IntegriCloud