s_cbrtf.c (4bb9780353ed9c7311a608508fea74bc29262d92) s_cbrtf.c (fd2891004df86c868f77e042a840510c496c81aa)
1/* s_cbrtf.c -- float version of s_cbrt.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@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.

--- 39 unchanged lines hidden (view full) ---

48 if(hx<0x00800000) { /* subnormal number */
49 SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */
50 t*=x;
51 GET_FLOAT_WORD(high,t);
52 SET_FLOAT_WORD(t,sign|((high&0x7fffffff)/3+B2));
53 } else
54 SET_FLOAT_WORD(t,sign|(hx/3+B1));
55
1/* s_cbrtf.c -- float version of s_cbrt.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@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.

--- 39 unchanged lines hidden (view full) ---

48 if(hx<0x00800000) { /* subnormal number */
49 SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */
50 t*=x;
51 GET_FLOAT_WORD(high,t);
52 SET_FLOAT_WORD(t,sign|((high&0x7fffffff)/3+B2));
53 } else
54 SET_FLOAT_WORD(t,sign|(hx/3+B1));
55
56 /* first step Newton iteration (solving t*t-x/t == 0) to 16 bits */
57 /* in double precision to avoid problems with denormals */
56 /*
57 * First step Newton iteration (solving t*t-x/t == 0) to 16 bits. In
58 * double precision so that its terms can be arranged for efficiency
59 * without causing overflow or underflow.
60 */
58 T=t;
59 r=T*T*T;
61 T=t;
62 r=T*T*T;
60 T=T*(x+x+r)/(x+r+r);
63 T=T*((double)x+x+r)/(x+r+r);
61
64
62 /* second step Newton iteration to 47 bits */
63 /* in double precision for accuracy */
65 /*
66 * Second step Newton iteration to 47 bits. In double precision for
67 * efficiency and accuracy.
68 */
64 r=T*T*T;
69 r=T*T*T;
65 T=T*(x+x+r)/(x+r+r);
70 T=T*((double)x+x+r)/(x+r+r);
66
67 /* rounding to 24 bits is perfect in round-to-nearest mode */
68 return(T);
69}
71
72 /* rounding to 24 bits is perfect in round-to-nearest mode */
73 return(T);
74}