| Index: fusl/src/math/cbrtl.c
|
| diff --git a/fusl/src/math/cbrtl.c b/fusl/src/math/cbrtl.c
|
| index ceff9136ebb507a604688ce07f3e6eae0fcd888c..525f4edcf6635267748e1d1894991baef7138774 100644
|
| --- a/fusl/src/math/cbrtl.c
|
| +++ b/fusl/src/math/cbrtl.c
|
| @@ -18,107 +18,109 @@
|
| #include "libm.h"
|
|
|
| #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
|
| -long double cbrtl(long double x)
|
| -{
|
| - return cbrt(x);
|
| +long double cbrtl(long double x) {
|
| + return cbrt(x);
|
| }
|
| #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
|
| -static const unsigned B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
|
| +static const unsigned B1 =
|
| + 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
|
|
|
| -long double cbrtl(long double x)
|
| -{
|
| - union ldshape u = {x}, v;
|
| - union {float f; uint32_t i;} uft;
|
| - long double r, s, t, w;
|
| - double_t dr, dt, dx;
|
| - float_t ft;
|
| - int e = u.i.se & 0x7fff;
|
| - int sign = u.i.se & 0x8000;
|
| +long double cbrtl(long double x) {
|
| + union ldshape u = {x}, v;
|
| + union {
|
| + float f;
|
| + uint32_t i;
|
| + } uft;
|
| + long double r, s, t, w;
|
| + double_t dr, dt, dx;
|
| + float_t ft;
|
| + int e = u.i.se & 0x7fff;
|
| + int sign = u.i.se & 0x8000;
|
|
|
| - /*
|
| - * If x = +-Inf, then cbrt(x) = +-Inf.
|
| - * If x = NaN, then cbrt(x) = NaN.
|
| - */
|
| - if (e == 0x7fff)
|
| - return x + x;
|
| - if (e == 0) {
|
| - /* Adjust subnormal numbers. */
|
| - u.f *= 0x1p120;
|
| - e = u.i.se & 0x7fff;
|
| - /* If x = +-0, then cbrt(x) = +-0. */
|
| - if (e == 0)
|
| - return x;
|
| - e -= 120;
|
| - }
|
| - e -= 0x3fff;
|
| - u.i.se = 0x3fff;
|
| - x = u.f;
|
| - switch (e % 3) {
|
| - case 1:
|
| - case -2:
|
| - x *= 2;
|
| - e--;
|
| - break;
|
| - case 2:
|
| - case -1:
|
| - x *= 4;
|
| - e -= 2;
|
| - break;
|
| - }
|
| - v.f = 1.0;
|
| - v.i.se = sign | (0x3fff + e/3);
|
| + /*
|
| + * If x = +-Inf, then cbrt(x) = +-Inf.
|
| + * If x = NaN, then cbrt(x) = NaN.
|
| + */
|
| + if (e == 0x7fff)
|
| + return x + x;
|
| + if (e == 0) {
|
| + /* Adjust subnormal numbers. */
|
| + u.f *= 0x1p120;
|
| + e = u.i.se & 0x7fff;
|
| + /* If x = +-0, then cbrt(x) = +-0. */
|
| + if (e == 0)
|
| + return x;
|
| + e -= 120;
|
| + }
|
| + e -= 0x3fff;
|
| + u.i.se = 0x3fff;
|
| + x = u.f;
|
| + switch (e % 3) {
|
| + case 1:
|
| + case -2:
|
| + x *= 2;
|
| + e--;
|
| + break;
|
| + case 2:
|
| + case -1:
|
| + x *= 4;
|
| + e -= 2;
|
| + break;
|
| + }
|
| + v.f = 1.0;
|
| + v.i.se = sign | (0x3fff + e / 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.
|
| - */
|
| + /*
|
| + * 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: */
|
| - uft.f = x;
|
| - uft.i = (uft.i & 0x7fffffff)/3 + B1;
|
| - ft = uft.f;
|
| + /* ~5-bit estimate: */
|
| + uft.f = x;
|
| + uft.i = (uft.i & 0x7fffffff) / 3 + B1;
|
| + ft = uft.f;
|
|
|
| - /* ~16-bit estimate: */
|
| - dx = x;
|
| - dt = ft;
|
| - dr = dt * dt * dt;
|
| - dt = dt * (dx + dx + dr) / (dx + dr + dr);
|
| + /* ~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);
|
| + /* ~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).
|
| - */
|
| - t = dt + (0x1.0p32L + 0x1.0p-31L) - 0x1.0p32;
|
| + /*
|
| + * 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).
|
| + */
|
| + t = dt + (0x1.0p32L + 0x1.0p-31L) - 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;
|
| + /*
|
| + * 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;
|
| #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 */
|
| + /*
|
| + * 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.f;
|
| - return t;
|
| + t *= v.f;
|
| + return t;
|
| }
|
| #endif
|
|
|