1 /* sf_cbrt.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 
17 #include "fdlibm.h"
18 
19 /* cbrtf(x)
20  * Return cube root of x
21  */
22 #ifdef __STDC__
23 static const __uint32_t
24 #else
25 static __uint32_t
26 #endif
27 	B1 = 709958130, /* B1 = (84+2/3-0.03306235651)*2**23 */
28 	B2 = 642849266; /* B2 = (76+2/3-0.03306235651)*2**23 */
29 
30 #ifdef __STDC__
31 static const float
32 #else
33 static float
34 #endif
35 C =  5.4285717010e-01, /* 19/35     = 0x3f0af8b0 */
36 D = -7.0530611277e-01, /* -864/1225 = 0xbf348ef1 */
37 E =  1.4142856598e+00, /* 99/70     = 0x3fb50750 */
38 F =  1.6071428061e+00, /* 45/28     = 0x3fcdb6db */
39 G =  3.5714286566e-01; /* 5/14      = 0x3eb6db6e */
40 
41 #ifdef __STDC__
cbrtf(float x)42 	float cbrtf(float x)
43 #else
44 	float cbrtf(x)
45 	float x;
46 #endif
47 {
48 	__int32_t	hx;
49 	float r,s,t;
50 	__uint32_t sign;
51 	__uint32_t high;
52 
53 	GET_FLOAT_WORD(hx,x);
54 	sign=hx&0x80000000; 		/* sign= sign(x) */
55 	hx  ^=sign;
56 	if(!FLT_UWORD_IS_FINITE(hx))
57 	    return(x+x);		/* cbrt(NaN,INF) is itself */
58 	if(FLT_UWORD_IS_ZERO(hx))
59 	    return(x);			/* cbrt(0) is itself */
60 
61 	SET_FLOAT_WORD(x,hx);	/* x <- |x| */
62     /* rough cbrt to 5 bits */
63 	if(FLT_UWORD_IS_SUBNORMAL(hx)) 		/* subnormal number */
64 	  {SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */
65 	   t*=x; GET_FLOAT_WORD(high,t); SET_FLOAT_WORD(t,high/3+B2);
66 	  }
67 	else
68 	  SET_FLOAT_WORD(t,hx/3+B1);
69 
70 
71     /* new cbrt to 23 bits */
72 	r=t*t/x;
73 	s=C+r*t;
74 	t*=G+F/(s+E+D/s);
75 
76     /* retore the sign bit */
77 	GET_FLOAT_WORD(high,t);
78 	SET_FLOAT_WORD(t,high|sign);
79 	return(t);
80 }
81 
82 _MATH_ALIAS_f_f(cbrt)
83