Author: kargl
Date: Sat Mar 12 19:37:35 2011
New Revision: 219576
URL: http://svn.freebsd.org/changeset/base/219576
Log:
Take two. Add the missing file that should have been committed
with r219571 and re-enable building of cbrtl.
Implement the long double version for the cube root function, cbrtl.
The algorithm uses Newton's iterations with a crude estimate of the
cube root to converge to a result.
Reviewed by: bde
Approved by: das
Added:
head/lib/msun/src/s_cbrtl.c (contents, props changed)
Modified:
head/lib/msun/Makefile
head/lib/msun/Symbol.map
Modified: head/lib/msun/Makefile
==============================================================================
--- head/lib/msun/Makefile Sat Mar 12 19:07:19 2011 (r219575)
+++ head/lib/msun/Makefile Sat Mar 12 19:37:35 2011 (r219576)
@@ -93,7 +93,7 @@ COMMON_SRCS+= s_copysignl.c s_fabsl.c s_
COMMON_SRCS+= e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \
e_hypotl.c e_remainderl.c e_sqrtl.c \
invtrig.c k_cosl.c k_sinl.c k_tanl.c \
- s_atanl.c s_ceill.c s_cosl.c s_cprojl.c \
+ s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
s_csqrtl.c s_exp2l.c s_floorl.c s_fmal.c \
s_frexpl.c s_logbl.c s_nanl.c s_nextafterl.c s_nexttoward.c \
s_remquol.c s_rintl.c s_scalbnl.c \
Modified: head/lib/msun/Symbol.map
==============================================================================
--- head/lib/msun/Symbol.map Sat Mar 12 19:07:19 2011 (r219575)
+++ head/lib/msun/Symbol.map Sat Mar 12 19:37:35 2011 (r219576)
@@ -222,6 +222,7 @@ FBSD_1.1 {
/* First added in 9.0-CURRENT */
FBSD_1.2 {
__isnanf;
+ cbrtl;
cexp;
cexpf;
log2;
Added: head/lib/msun/src/s_cbrtl.c
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ head/lib/msun/src/s_cbrtl.c Sat Mar 12 19:37:35 2011 (r219576)
@@ -0,0 +1,157 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2009-2011, Bruce D. Evans, Steven G. Kargl, David Schultz.
+ *
+ * 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.
+ * ====================================================
+ *
+ * The argument reduction and testing for exceptional cases was
+ * written by Steven G. Kargl with input from Bruce D. Evans
+ * and David A. Schultz.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+#include <ieeefp.h>
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+static const unsigned
+ B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
+
+long double
+cbrtl(long double x)
+{
+ union IEEEl2bits u, v;
+ long double r, s, t, w;
+ double dr, dt, dx;
+ float ft, fx;
+ uint32_t hx;
+ uint16_t expsign;
+ int k;
+
+ u.e = x;
+ expsign = u.xbits.expsign;
+ k = expsign & 0x7fff;
+
+ /*
+ * If x = +-Inf, then cbrt(x) = +-Inf.
+ * If x = NaN, then cbrt(x) = NaN.
+ */
+ if (k == BIAS + LDBL_MAX_EXP)
+ return (x + x);
+
+#ifdef __i386__
+ fp_prec_t oprec;
+
+ oprec = fpgetprec();
+ if (oprec != FP_PE)
+ fpsetprec(FP_PE);
+#endif
+
+ if (k == 0) {
+ /* If x = +-0, then cbrt(x) = +-0. */
+ if ((u.bits.manh | u.bits.manl) == 0) {
+#ifdef __i386__
+ if (oprec != FP_PE)
+ fpsetprec(oprec);
+#endif
+ return (x);
+ }
+ /* Adjust subnormal numbers. */
+ u.e *= 0x1.0p514;
+ k = u.bits.exp;
+ k -= BIAS + 514;
+ } else
+ k -= BIAS;
+ u.xbits.expsign = BIAS;
+ v.e = 1;
+
+ x = u.e;
+ switch (k % 3) {
+ case 1:
+ case -2:
+ x = 2*x;
+ k--;
+ break;
+ case 2:
+ case -1:
+ x = 4*x;
+ k -= 2;
+ break;
+ }
+ v.xbits.expsign = (expsign & 0x8000) | (BIAS + k / 3);
+
+ /*
+ * The following is the guts of s_cbrtf, with the handling of
+ * special values removed and extra care for accuracy not taken,
+ * but with most of the extra accuracy not discarded.
+ */
+
+ /* ~5-bit estimate: */
+ fx = x;
+ GET_FLOAT_WORD(hx, fx);
+ SET_FLOAT_WORD(ft, ((hx & 0x7fffffff) / 3 + B1));
+
+ /* ~16-bit estimate: */
+ dx = x;
+ dt = ft;
+ dr = dt * dt * dt;
+ dt = dt * (dx + dx + dr) / (dx + dr + dr);
+
+ /* ~47-bit estimate: */
+ dr = dt * dt * dt;
+ dt = dt * (dx + dx + dr) / (dx + dr + dr);
+
+#if LDBL_MANT_DIG == 64
+ /*
+ * dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8).
+ * Round it away from zero to 32 bits (32 so that t*t is exact, and
+ * away from zero for technical reasons).
+ */
+ volatile double vd2 = 0x1.0p32;
+ volatile double vd1 = 0x1.0p-31;
+ #define vd ((long double)vd2 + vd1)
+
+ t = dt + vd - 0x1.0p32;
+#elif LDBL_MANT_DIG == 113
+ /*
+ * Round dt away from zero to 47 bits. Since we don't trust the 47,
+ * add 2 47-bit ulps instead of 1 to round up. Rounding is slow and
+ * might be avoidable in this case, since on most machines dt will
+ * have been evaluated in 53-bit precision and the technical reasons
+ * for rounding up might not apply to either case in cbrtl() since
+ * dt is much more accurate than needed.
+ */
+ t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
+#else
+#error "Unsupported long double format"
+#endif
+
+ /*
+ * Final step Newton iteration to 64 or 113 bits with
+ * error < 0.667 ulps
+ */
+ s=t*t; /* t*t is exact */
+ r=x/s; /* error <= 0.5 ulps; |r| < |t| */
+ w=t+t; /* t+t is exact */
+ r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
+ t=t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */
+
+ t *= v.e;
+#ifdef __i386__
+ if (oprec != FP_PE)
+ fpsetprec(oprec);
+#endif
+ return (t);
+}
_______________________________________________
svn-src-all@freebsd.org mailing list
http://lists.freebsd.org/mailman/listinfo/svn-src-all
To unsubscribe, send any mail to "svn-src-all-unsubscr...@freebsd.org"
