From 1cd0ec03d6696c84c743f05a105ab1537e35ac65 Mon Sep 17 00:00:00 2001 From: Steve Kargl Date: Sat, 12 Mar 2011 19:37:35 +0000 Subject: [PATCH] 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 --- lib/msun/Makefile | 2 +- lib/msun/Symbol.map | 1 + lib/msun/src/s_cbrtl.c | 157 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 159 insertions(+), 1 deletion(-) create mode 100644 lib/msun/src/s_cbrtl.c diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 674ebf06754f..bb99e45a25b2 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -93,7 +93,7 @@ COMMON_SRCS+= s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c 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 \ diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map index cb60ecf325a2..bb53ac39f9f4 100644 --- a/lib/msun/Symbol.map +++ b/lib/msun/Symbol.map @@ -222,6 +222,7 @@ FBSD_1.1 { /* First added in 9.0-CURRENT */ FBSD_1.2 { __isnanf; + cbrtl; cexp; cexpf; log2; diff --git a/lib/msun/src/s_cbrtl.c b/lib/msun/src/s_cbrtl.c new file mode 100644 index 000000000000..23c9184ae18b --- /dev/null +++ b/lib/msun/src/s_cbrtl.c @@ -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 +__FBSDID("$FreeBSD$"); + +#include +#include + +#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); +}