| ... | ... | @@ -1,66 +0,0 @@ |
| 1 | | /* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtf.c */ |
| 2 | | /* |
| 3 | | * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
| 4 | | * Debugged and optimized by Bruce D. Evans. |
| 5 | | */ |
| 6 | | /* |
| 7 | | * ==================================================== |
| 8 | | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 9 | | * |
| 10 | | * Developed at SunPro, a Sun Microsystems, Inc. business. |
| 11 | | * Permission to use, copy, modify, and distribute this |
| 12 | | * software is freely granted, provided that this notice |
| 13 | | * is preserved. |
| 14 | | * ==================================================== |
| 15 | | */ |
| 16 | | /* cbrtf(x) |
| 17 | | * Return cube root of x |
| 18 | | */ |
| 19 | | |
| 20 | | #include <math.h> |
| 21 | | #include <stdint.h> |
| 22 | | |
| 23 | | static const unsigned |
| 24 | | B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */ |
| 25 | | B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */ |
| 26 | | |
| 27 | | float cbrtf(float x) |
| 28 | | { |
| 29 | | 	double_t r,T; |
| 30 | | 	union {float f; uint32_t i;} u = {x}; |
| 31 | | 	uint32_t hx = u.i & 0x7fffffff; |
| 32 | | |
| 33 | | 	if (hx >= 0x7f800000) /* cbrt(NaN,INF) is itself */ |
| 34 | | 		return x + x; |
| 35 | | |
| 36 | | 	/* rough cbrt to 5 bits */ |
| 37 | | 	if (hx < 0x00800000) { /* zero or subnormal? */ |
| 38 | | 		if (hx == 0) |
| 39 | | 			return x; /* cbrt(+-0) is itself */ |
| 40 | | 		u.f = x*0x1p24f; |
| 41 | | 		hx = u.i & 0x7fffffff; |
| 42 | | 		hx = hx/3 + B2; |
| 43 | | 	} else |
| 44 | | 		hx = hx/3 + B1; |
| 45 | | 	u.i &= 0x80000000; |
| 46 | | 	u.i |= hx; |
| 47 | | |
| 48 | | 	/* |
| 49 | | 	 * First step Newton iteration (solving t*t-x/t == 0) to 16 bits. In |
| 50 | | 	 * double precision so that its terms can be arranged for efficiency |
| 51 | | 	 * without causing overflow or underflow. |
| 52 | | 	 */ |
| 53 | | 	T = u.f; |
| 54 | | 	r = T*T*T; |
| 55 | | 	T = T*((double_t)x+x+r)/(x+r+r); |
| 56 | | |
| 57 | | 	/* |
| 58 | | 	 * Second step Newton iteration to 47 bits. In double precision for |
| 59 | | 	 * efficiency and accuracy. |
| 60 | | 	 */ |
| 61 | | 	r = T*T*T; |
| 62 | | 	T = T*((double_t)x+x+r)/(x+r+r); |
| 63 | | |
| 64 | | 	/* rounding to 24 bits is perfect in round-to-nearest mode */ |
| 65 | | 	return T; |
| 66 | | } |