1 /* cbrt.c
2 *
3 * Cube root
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, cbrt();
10 *
11 * y = cbrt( x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns the cube root of the argument, which may be negative.
18 *
19 * Range reduction involves determining the power of 2 of
20 * the argument. A polynomial of degree 2 applied to the
21 * mantissa, and multiplication by the cube root of 1, 2, or 4
22 * approximates the root to within about 0.1%. Then Newton's
23 * iteration is used three times to converge to an accurate
24 * result.
25 *
26 *
27 *
28 * ACCURACY:
29 *
30 * Relative error:
31 * arithmetic domain # trials peak rms
32 * DEC -10,10 200000 1.8e-17 6.2e-18
33 * IEEE 0,1e308 30000 1.5e-16 5.0e-17
34 *
35 */
36 /* cbrt.c */
37
38 /*
39 Cephes Math Library Release 2.2: January, 1991
40 Copyright 1984, 1991 by Stephen L. Moshier
41 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
42 */
43
44
45 #include "mconf.h"
46 #include "cephes.h"
47
48 #ifndef HAVE_CBRT
49
50 static double CBRT2 = 1.2599210498948731647672;
51 static double CBRT4 = 1.5874010519681994747517;
52 static double CBRT2I = 0.79370052598409973737585;
53 static double CBRT4I = 0.62996052494743658238361;
54
cbrt(x)55 double cbrt(x)
56 double x;
57 {
58 int e, rem, sign;
59 double z;
60
61 #ifdef NANS
62 if( isnan(x) )
63 return x;
64 #endif
65 #ifdef INFINITIES
66 if( !finite(x) )
67 return x;
68 #endif
69 if( x == 0 )
70 return( x );
71 if( x > 0 )
72 sign = 1;
73 else
74 {
75 sign = -1;
76 x = -x;
77 }
78
79 z = x;
80 /* extract power of 2, leaving
81 * mantissa between 0.5 and 1
82 */
83 x = frexp( x, &e );
84
85 /* Approximate cube root of number between .5 and 1,
86 * peak relative error = 9.2e-6
87 */
88 x = (((-1.3466110473359520655053e-1 * x
89 + 5.4664601366395524503440e-1) * x
90 - 9.5438224771509446525043e-1) * x
91 + 1.1399983354717293273738e0 ) * x
92 + 4.0238979564544752126924e-1;
93
94 /* exponent divided by 3 */
95 if( e >= 0 )
96 {
97 rem = e;
98 e /= 3;
99 rem -= 3*e;
100 if( rem == 1 )
101 x *= CBRT2;
102 else if( rem == 2 )
103 x *= CBRT4;
104 }
105
106
107 /* argument less than 1 */
108
109 else
110 {
111 e = -e;
112 rem = e;
113 e /= 3;
114 rem -= 3*e;
115 if( rem == 1 )
116 x *= CBRT2I;
117 else if( rem == 2 )
118 x *= CBRT4I;
119 e = -e;
120 }
121
122 /* multiply by power of 2 */
123 x = ldexp( x, e );
124
125 /* Newton iteration */
126 x -= ( x - (z/(x*x)) )*0.33333333333333333333;
127 #ifdef DEC
128 x -= ( x - (z/(x*x)) )/3.0;
129 #else
130 x -= ( x - (z/(x*x)) )*0.33333333333333333333;
131 #endif
132
133 if( sign < 0 )
134 x = -x;
135 return(x);
136 }
137
138 #endif
139
140