Mercurial > mxe-octave

comparison src/build-msvctools/math/cbrt.c @ 3061:f8299bb6c872

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial support for native MSVC compilation. * add MSVC support files: compiler wrappers and support libraries * adapt libiconv to work with MSVC * adapt gettext to work with MSVC
author Michael Goffioul <michael.goffioul@gmail.com>
date Mon, 17 Jun 2013 22:43:11 -0400
parents
children
comparison
equal deleted inserted replaced
3060:cbdf4575016d 3061:f8299bb6c872
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 Modified for mingwex.a
46 2002-07-01 Danny Smith <dannysmith@users.sourceforge.net>
47 */
48 #ifdef __MINGW32__
49 #include <math.h>
50 #include "cephes_mconf.h"
51 #else
52 #include "mconf.h"
53 #endif
54
55
56 static const double CBRT2 = 1.2599210498948731647672;
57 static const double CBRT4 = 1.5874010519681994747517;
58 static const double CBRT2I = 0.79370052598409973737585;
59 static const double CBRT4I = 0.62996052494743658238361;
60
61 #ifndef __MINGW32__
62 #ifdef ANSIPROT
63 extern double frexp ( double, int * );
64 extern double ldexp ( double, int );
65 extern int isnan ( double );
66 extern int isfinite ( double );
67 #else
68 double frexp(), ldexp();
69 int isnan(), isfinite();
70 #endif
71 #endif
72
73 double cbrt(x)
74 double x;
75 {
76 int e, rem, sign;
77 double z;
78
79 #ifdef __MINGW32__
80 if (!isfinite (x) || x == 0 )
81 return x;
82 #else
83
84 #ifdef NANS
85 if( isnan(x) )
86 return x;
87 #endif
88 #ifdef INFINITIES
89 if( !isfinite(x) )
90 return x;
91 #endif
92 if( x == 0 )
93 return( x );
94
95 #endif /* __MINGW32__ */
96
97 if( x > 0 )
98 sign = 1;
99 else
100 {
101 sign = -1;
102 x = -x;
103 }
104
105 z = x;
106 /* extract power of 2, leaving
107 * mantissa between 0.5 and 1
108 */
109 x = frexp( x, &e );
110
111 /* Approximate cube root of number between .5 and 1,
112 * peak relative error = 9.2e-6
113 */
114 x = (((-1.3466110473359520655053e-1 * x
115 + 5.4664601366395524503440e-1) * x
116 - 9.5438224771509446525043e-1) * x
117 + 1.1399983354717293273738e0 ) * x
118 + 4.0238979564544752126924e-1;
119
120 /* exponent divided by 3 */
121 if( e >= 0 )
122 {
123 rem = e;
124 e /= 3;
125 rem -= 3*e;
126 if( rem == 1 )
127 x *= CBRT2;
128 else if( rem == 2 )
129 x *= CBRT4;
130 }
131
132
133 /* argument less than 1 */
134
135 else
136 {
137 e = -e;
138 rem = e;
139 e /= 3;
140 rem -= 3*e;
141 if( rem == 1 )
142 x *= CBRT2I;
143 else if( rem == 2 )
144 x *= CBRT4I;
145 e = -e;
146 }
147
148 /* multiply by power of 2 */
149 x = ldexp( x, e );
150
151 /* Newton iteration */
152 x -= ( x - (z/(x*x)) )*0.33333333333333333333;
153 #ifdef DEC
154 x -= ( x - (z/(x*x)) )/3.0;
155 #else
156 x -= ( x - (z/(x*x)) )*0.33333333333333333333;
157 #endif
158
159 if( sign < 0 )
160 x = -x;
161 return(x);
162 }