133965Sjdp/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtl.c */ 277298Sobrien/*- 389857Sobrien * ==================================================== 4218822Sdim * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5218822Sdim * Copyright (c) 2009-2011, Bruce D. Evans, Steven G. Kargl, David Schultz. 668765Sobrien * 7218822Sdim * Developed at SunPro, a Sun Microsystems, Inc. business. 868765Sobrien * Permission to use, copy, modify, and distribute this 933965Sjdp * software is freely granted, provided that this notice 1033965Sjdp * is preserved. 1133965Sjdp * ==================================================== 1233965Sjdp * 1333965Sjdp * The argument reduction and testing for exceptional cases was 1433965Sjdp * written by Steven G. Kargl with input from Bruce D. Evans 1533965Sjdp * and David A. Schultz. 1633965Sjdp */ 1733965Sjdp 1833965Sjdp#include "libm.h" 1933965Sjdp 2033965Sjdp#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 2133965Sjdplong double cbrtl(long double x) 2233965Sjdp{ 2333965Sjdp return cbrt(x); 2433965Sjdp} 25218822Sdim#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 26218822Sdimstatic const unsigned B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */ 27218822Sdim 2833965Sjdplong double cbrtl(long double x) 2933965Sjdp{ 3033965Sjdp union ldshape u = {x}, v; 3133965Sjdp union {float f; uint32_t i;} uft; 3233965Sjdp long double r, s, t, w; 33218822Sdim double_t dr, dt, dx; 3489857Sobrien float_t ft; 3589857Sobrien int e = u.i.se & 0x7fff; 3660484Sobrien int sign = u.i.se & 0x8000; 3733965Sjdp 3833965Sjdp /* 3933965Sjdp * If x = +-Inf, then cbrt(x) = +-Inf. 4033965Sjdp * If x = NaN, then cbrt(x) = NaN. 4133965Sjdp */ 4233965Sjdp if (e == 0x7fff) 4333965Sjdp return x + x; 4433965Sjdp if (e == 0) { 4533965Sjdp /* Adjust subnormal numbers. */ 4633965Sjdp u.f *= 0x1p120; 4733965Sjdp e = u.i.se & 0x7fff; 4833965Sjdp /* If x = +-0, then cbrt(x) = +-0. */ 4933965Sjdp if (e == 0) 5033965Sjdp return x; 5133965Sjdp e -= 120; 5233965Sjdp } 5333965Sjdp e -= 0x3fff; 5433965Sjdp u.i.se = 0x3fff; 5533965Sjdp x = u.f; 5633965Sjdp switch (e % 3) { 5768765Sobrien case 1: 5833965Sjdp case -2: 5968765Sobrien x *= 2; 6068765Sobrien e--; 6168765Sobrien break; 6268765Sobrien case 2: 6368765Sobrien case -1: 6468765Sobrien x *= 4; 6568765Sobrien e -= 2; 6677298Sobrien break; 6777298Sobrien } 6877298Sobrien v.f = 1.0; 6968765Sobrien v.i.se = sign | (0x3fff + e/3); 7077298Sobrien 7177298Sobrien /* 7277298Sobrien * The following is the guts of s_cbrtf, with the handling of 7377298Sobrien * special values removed and extra care for accuracy not taken, 7477298Sobrien * but with most of the extra accuracy not discarded. 75218822Sdim */ 7677298Sobrien 7777298Sobrien /* ~5-bit estimate: */ 7877298Sobrien uft.f = x; 7977298Sobrien uft.i = (uft.i & 0x7fffffff)/3 + B1; 8077298Sobrien ft = uft.f; 8168765Sobrien 8268765Sobrien /* ~16-bit estimate: */ 8368765Sobrien dx = x; 8468765Sobrien dt = ft; 8568765Sobrien dr = dt * dt * dt; 8677298Sobrien dt = dt * (dx + dx + dr) / (dx + dr + dr); 8777298Sobrien 88218822Sdim /* ~47-bit estimate: */ 8977298Sobrien dr = dt * dt * dt; 90218822Sdim dt = dt * (dx + dx + dr) / (dx + dr + dr); 9168765Sobrien 92218822Sdim#if LDBL_MANT_DIG == 64 9368765Sobrien /* 9468765Sobrien * dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8). 9568765Sobrien * Round it away from zero to 32 bits (32 so that t*t is exact, and 9668765Sobrien * away from zero for technical reasons). 9768765Sobrien */ 9877298Sobrien t = dt + (0x1.0p32L + 0x1.0p-31L) - 0x1.0p32; 9968765Sobrien#elif LDBL_MANT_DIG == 113 10068765Sobrien /* 10168765Sobrien * Round dt away from zero to 47 bits. Since we don't trust the 47, 10268765Sobrien * add 2 47-bit ulps instead of 1 to round up. Rounding is slow and 10368765Sobrien * might be avoidable in this case, since on most machines dt will 104218822Sdim * have been evaluated in 53-bit precision and the technical reasons 10568765Sobrien * for rounding up might not apply to either case in cbrtl() since 10668765Sobrien * dt is much more accurate than needed. 10768765Sobrien */ 10868765Sobrien t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60; 10968765Sobrien#endif 11068765Sobrien 11168765Sobrien /* 11268765Sobrien * Final step Newton iteration to 64 or 113 bits with 11368765Sobrien * error < 0.667 ulps 11468765Sobrien */ 11568765Sobrien s = t*t; /* t*t is exact */ 11668765Sobrien r = x/s; /* error <= 0.5 ulps; |r| < |t| */ 11733965Sjdp w = t+t; /* t+t is exact */ 11833965Sjdp r = (r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */ 11933965Sjdp t = t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */ 12033965Sjdp 12133965Sjdp t *= v.f; 12233965Sjdp return t; 123218822Sdim} 124218822Sdim#endif 125218822Sdim