1 /* s_cbrtf.c -- float version of s_cbrt.c. 2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3 */ 4 5 /* 6 * ==================================================== 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8 * 9 * Developed at SunPro, a Sun Microsystems, Inc. business. 10 * Permission to use, copy, modify, and distribute this 11 * software is freely granted, provided that this notice 12 * is preserved. 13 * ==================================================== 14 */ 15 16 #include <sys/cdefs.h> 17 #if defined(LIBM_SCCS) && !defined(lint) 18 __RCSID("$NetBSD: s_cbrtf.c,v 1.6 1999/07/02 15:37:42 simonb 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 #ifdef __STDC__ 28 static const unsigned 29 #else 30 static unsigned 31 #endif 32 B1 = 709958130, /* B1 = (84+2/3-0.03306235651)*2**23 */ 33 B2 = 642849266; /* B2 = (76+2/3-0.03306235651)*2**23 */ 34 35 #ifdef __STDC__ 36 static const float 37 #else 38 static float 39 #endif 40 C = 5.4285717010e-01, /* 19/35 = 0x3f0af8b0 */ 41 D = -7.0530611277e-01, /* -864/1225 = 0xbf348ef1 */ 42 E = 1.4142856598e+00, /* 99/70 = 0x3fb50750 */ 43 F = 1.6071428061e+00, /* 45/28 = 0x3fcdb6db */ 44 G = 3.5714286566e-01; /* 5/14 = 0x3eb6db6e */ 45 46 #ifdef __STDC__ 47 float cbrtf(float x) 48 #else 49 float cbrtf(x) 50 float x; 51 #endif 52 { 53 float r,s,t; 54 int32_t hx; 55 u_int32_t sign; 56 u_int32_t high; 57 58 GET_FLOAT_WORD(hx,x); 59 sign=hx&0x80000000; /* sign= sign(x) */ 60 hx ^=sign; 61 if(hx>=0x7f800000) return(x+x); /* cbrt(NaN,INF) is itself */ 62 if(hx==0) 63 return(x); /* cbrt(0) is itself */ 64 65 SET_FLOAT_WORD(x,hx); /* x <- |x| */ 66 /* rough cbrt to 5 bits */ 67 if(hx<0x00800000) /* subnormal number */ 68 {SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */ 69 t*=x; GET_FLOAT_WORD(high,t); SET_FLOAT_WORD(t,high/3+B2); 70 } 71 else 72 SET_FLOAT_WORD(t,hx/3+B1); 73 74 75 /* new cbrt to 23 bits */ 76 r=t*t/x; 77 s=C+r*t; 78 t*=G+F/(s+E+D/s); 79 80 /* retore the sign bit */ 81 GET_FLOAT_WORD(high,t); 82 SET_FLOAT_WORD(t,high|sign); 83 return(t); 84 } 85