cbrt.c
Go to the documentation of this file.
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  * IEEE 0,1e308 30000 1.5e-16 5.0e-17
33  *
34  */
35  /* cbrt.c */
36 
37 /*
38  * Cephes Math Library Release 2.2: January, 1991
39  * Copyright 1984, 1991 by Stephen L. Moshier
40  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
41  */
42 
43 
44 #include "mconf.h"
45 
46 static double CBRT2 = 1.2599210498948731647672;
47 static double CBRT4 = 1.5874010519681994747517;
48 static double CBRT2I = 0.79370052598409973737585;
49 static double CBRT4I = 0.62996052494743658238361;
50 
51 double cbrt(double x)
52 {
53  int e, rem, sign;
54  double z;
55 
56  if (!cephes_isfinite(x))
57  return x;
58  if (x == 0)
59  return (x);
60  if (x > 0)
61  sign = 1;
62  else {
63  sign = -1;
64  x = -x;
65  }
66 
67  z = x;
68  /* extract power of 2, leaving
69  * mantissa between 0.5 and 1
70  */
71  x = frexp(x, &e);
72 
73  /* Approximate cube root of number between .5 and 1,
74  * peak relative error = 9.2e-6
75  */
76  x = (((-1.3466110473359520655053e-1 * x
77  + 5.4664601366395524503440e-1) * x
78  - 9.5438224771509446525043e-1) * x
79  + 1.1399983354717293273738e0) * x + 4.0238979564544752126924e-1;
80 
81  /* exponent divided by 3 */
82  if (e >= 0) {
83  rem = e;
84  e /= 3;
85  rem -= 3 * e;
86  if (rem == 1)
87  x *= CBRT2;
88  else if (rem == 2)
89  x *= CBRT4;
90  }
91 
92 
93  /* argument less than 1 */
94 
95  else {
96  e = -e;
97  rem = e;
98  e /= 3;
99  rem -= 3 * e;
100  if (rem == 1)
101  x *= CBRT2I;
102  else if (rem == 2)
103  x *= CBRT4I;
104  e = -e;
105  }
106 
107  /* multiply by power of 2 */
108  x = ldexp(x, e);
109 
110  /* Newton iteration */
111  x -= (x - (z / (x * x))) * 0.33333333333333333333;
112  x -= (x - (z / (x * x))) * 0.33333333333333333333;
113 
114  if (sign < 0)
115  x = -x;
116  return (x);
117 }
CBRT4I
static double CBRT4I
Definition: cbrt.c:49
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
sign
const EIGEN_DEVICE_FUNC SignReturnType sign() const
Definition: ArrayCwiseUnaryOps.h:219
CBRT2I
static double CBRT2I
Definition: cbrt.c:48
CBRT4
static double CBRT4
Definition: cbrt.c:47
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
cbrt
double cbrt(double x)
Definition: cbrt.c:51
mconf.h
cephes_isfinite
#define cephes_isfinite(x)
Definition: mconf.h:101
CBRT2
static double CBRT2
Definition: cbrt.c:46


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:01:57