diff options
Diffstat (limited to 'src/math/cbrtf.c')
| -rw-r--r-- | src/math/cbrtf.c | 31 |
1 files changed, 14 insertions, 17 deletions
diff --git a/src/math/cbrtf.c b/src/math/cbrtf.c index 4a984b10..89c2c865 100644 --- a/src/math/cbrtf.c +++ b/src/math/cbrtf.c @@ -17,7 +17,8 @@ * Return cube root of x */ -#include "libm.h" +#include <math.h> +#include <stdint.h> static const unsigned B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */ @@ -25,15 +26,10 @@ B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */ float cbrtf(float x) { - double r,T; - float t; - int32_t hx; - uint32_t sign; - uint32_t high; + double_t r,T; + union {float f; uint32_t i;} u = {x}; + uint32_t hx = u.i & 0x7fffffff; - GET_FLOAT_WORD(hx, x); - sign = hx & 0x80000000; - hx ^= sign; if (hx >= 0x7f800000) /* cbrt(NaN,INF) is itself */ return x + x; @@ -41,28 +37,29 @@ float cbrtf(float x) if (hx < 0x00800000) { /* zero or subnormal? */ if (hx == 0) return x; /* cbrt(+-0) is itself */ - SET_FLOAT_WORD(t, 0x4b800000); /* set t = 2**24 */ - t *= x; - GET_FLOAT_WORD(high, t); - SET_FLOAT_WORD(t, sign|((high&0x7fffffff)/3+B2)); + u.f = x*0x1p24f; + hx = u.i & 0x7fffffff; + hx = hx/3 + B2; } else - SET_FLOAT_WORD(t, sign|(hx/3+B1)); + hx = hx/3 + B1; + u.i &= 0x80000000; + u.i |= hx; /* * First step Newton iteration (solving t*t-x/t == 0) to 16 bits. In * double precision so that its terms can be arranged for efficiency * without causing overflow or underflow. */ - T = t; + T = u.f; r = T*T*T; - T = T*((double)x+x+r)/(x+r+r); + T = T*((double_t)x+x+r)/(x+r+r); /* * Second step Newton iteration to 47 bits. In double precision for * efficiency and accuracy. */ r = T*T*T; - T = T*((double)x+x+r)/(x+r+r); + T = T*((double_t)x+x+r)/(x+r+r); /* rounding to 24 bits is perfect in round-to-nearest mode */ return T; |