1 /* s_cbrtf.c -- float version of s_cbrt.c. 2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian (at) cygnus.com. 3 * Debugged and optimized by Bruce D. Evans. 4 */ 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 17 #ifndef lint 18 static char rcsid[] = "$FreeBSD: src/lib/msun/src/s_cbrtf.c,v 1.12 2005/12/13 20:17:23 bde Exp $"; 19 #endif 20 21 #include "math.h" 22 #include "math_private.h" 23 24 /* cbrtf(x) 25 * Return cube root of x 26 */ 27 static const unsigned 28 B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */ 29 B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */ 30 31 static const float 32 C = 5.4285717010e-01, /* 19/35 = 0x3f0af8b0 */ 33 D = -7.0530611277e-01, /* -864/1225 = 0xbf348ef1 */ 34 E = 1.4142856598e+00, /* 99/70 = 0x3fb50750 */ 35 F = 1.6071428061e+00, /* 45/28 = 0x3fcdb6db */ 36 G = 3.5714286566e-01; /* 5/14 = 0x3eb6db6e */ 37 38 float 39 cbrtf(float x) 40 { 41 float r,s,t,w; 42 int32_t hx; 43 u_int32_t sign; 44 u_int32_t high; 45 46 GET_FLOAT_WORD(hx,x); 47 sign=hx&0x80000000; /* sign= sign(x) */ 48 hx ^=sign; 49 if(hx>=0x7f800000) return(x+x); /* cbrt(NaN,INF) is itself */ 50 if(hx==0) 51 return(x); /* cbrt(0) is itself */ 52 53 /* rough cbrt to 5 bits */ 54 if(hx<0x00800000) { /* subnormal number */ 55 SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */ 56 t*=x; 57 GET_FLOAT_WORD(high,t); 58 SET_FLOAT_WORD(t,sign|((high&0x7fffffff)/3+B2)); 59 } else 60 SET_FLOAT_WORD(t,sign|(hx/3+B1)); 61 62 /* new cbrt to 23 bits */ 63 r=t*t/x; 64 s=C+r*t; 65 t*=G+F/(s+E+D/s); 66 67 /* chop t to 12 bits and make it larger in magnitude than cbrt(x) */ 68 GET_FLOAT_WORD(high,t); 69 SET_FLOAT_WORD(t,(high&0xfffff000)+0x00001000); 70 71 /* one step Newton iteration to 24 bits with error less than 0.667 ulps */ 72 s=t*t; /* t*t is exact */ 73 r=x/s; 74 w=t+t; 75 r=(r-t)/(w+r); /* r-t is exact */ 76 t=t+t*r; 77 78 return(t); 79 } 80
