1*181254a7Smrg /* cbrtq.c
2*181254a7Smrg *
3*181254a7Smrg * Cube root, long double precision
4*181254a7Smrg *
5*181254a7Smrg *
6*181254a7Smrg *
7*181254a7Smrg * SYNOPSIS:
8*181254a7Smrg *
9*181254a7Smrg * long double x, y, cbrtq();
10*181254a7Smrg *
11*181254a7Smrg * y = cbrtq( x );
12*181254a7Smrg *
13*181254a7Smrg *
14*181254a7Smrg *
15*181254a7Smrg * DESCRIPTION:
16*181254a7Smrg *
17*181254a7Smrg * Returns the cube root of the argument, which may be negative.
18*181254a7Smrg *
19*181254a7Smrg * Range reduction involves determining the power of 2 of
20*181254a7Smrg * the argument. A polynomial of degree 2 applied to the
21*181254a7Smrg * mantissa, and multiplication by the cube root of 1, 2, or 4
22*181254a7Smrg * approximates the root to within about 0.1%. Then Newton's
23*181254a7Smrg * iteration is used three times to converge to an accurate
24*181254a7Smrg * result.
25*181254a7Smrg *
26*181254a7Smrg *
27*181254a7Smrg *
28*181254a7Smrg * ACCURACY:
29*181254a7Smrg *
30*181254a7Smrg * Relative error:
31*181254a7Smrg * arithmetic domain # trials peak rms
32*181254a7Smrg * IEEE -8,8 100000 1.3e-34 3.9e-35
33*181254a7Smrg * IEEE exp(+-707) 100000 1.3e-34 4.3e-35
34*181254a7Smrg *
35*181254a7Smrg */
36*181254a7Smrg
37*181254a7Smrg /*
38*181254a7Smrg Cephes Math Library Release 2.2: January, 1991
39*181254a7Smrg Copyright 1984, 1991 by Stephen L. Moshier
40*181254a7Smrg Adapted for glibc October, 2001.
41*181254a7Smrg
42*181254a7Smrg This library is free software; you can redistribute it and/or
43*181254a7Smrg modify it under the terms of the GNU Lesser General Public
44*181254a7Smrg License as published by the Free Software Foundation; either
45*181254a7Smrg version 2.1 of the License, or (at your option) any later version.
46*181254a7Smrg
47*181254a7Smrg This library is distributed in the hope that it will be useful,
48*181254a7Smrg but WITHOUT ANY WARRANTY; without even the implied warranty of
49*181254a7Smrg MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
50*181254a7Smrg Lesser General Public License for more details.
51*181254a7Smrg
52*181254a7Smrg You should have received a copy of the GNU Lesser General Public
53*181254a7Smrg License along with this library; if not, see
54*181254a7Smrg <http://www.gnu.org/licenses/>. */
55*181254a7Smrg
56*181254a7Smrg #include "quadmath-imp.h"
57*181254a7Smrg
58*181254a7Smrg static const __float128 CBRT2 = 1.259921049894873164767210607278228350570251Q;
59*181254a7Smrg static const __float128 CBRT4 = 1.587401051968199474751705639272308260391493Q;
60*181254a7Smrg static const __float128 CBRT2I = 0.7937005259840997373758528196361541301957467Q;
61*181254a7Smrg static const __float128 CBRT4I = 0.6299605249474365823836053036391141752851257Q;
62*181254a7Smrg
63*181254a7Smrg
64*181254a7Smrg __float128
cbrtq(__float128 x)65*181254a7Smrg cbrtq (__float128 x)
66*181254a7Smrg {
67*181254a7Smrg int e, rem, sign;
68*181254a7Smrg __float128 z;
69*181254a7Smrg
70*181254a7Smrg if (!finiteq (x))
71*181254a7Smrg return x + x;
72*181254a7Smrg
73*181254a7Smrg if (x == 0)
74*181254a7Smrg return (x);
75*181254a7Smrg
76*181254a7Smrg if (x > 0)
77*181254a7Smrg sign = 1;
78*181254a7Smrg else
79*181254a7Smrg {
80*181254a7Smrg sign = -1;
81*181254a7Smrg x = -x;
82*181254a7Smrg }
83*181254a7Smrg
84*181254a7Smrg z = x;
85*181254a7Smrg /* extract power of 2, leaving mantissa between 0.5 and 1 */
86*181254a7Smrg x = frexpq (x, &e);
87*181254a7Smrg
88*181254a7Smrg /* Approximate cube root of number between .5 and 1,
89*181254a7Smrg peak relative error = 1.2e-6 */
90*181254a7Smrg x = ((((1.3584464340920900529734e-1Q * x
91*181254a7Smrg - 6.3986917220457538402318e-1Q) * x
92*181254a7Smrg + 1.2875551670318751538055e0Q) * x
93*181254a7Smrg - 1.4897083391357284957891e0Q) * x
94*181254a7Smrg + 1.3304961236013647092521e0Q) * x + 3.7568280825958912391243e-1Q;
95*181254a7Smrg
96*181254a7Smrg /* exponent divided by 3 */
97*181254a7Smrg if (e >= 0)
98*181254a7Smrg {
99*181254a7Smrg rem = e;
100*181254a7Smrg e /= 3;
101*181254a7Smrg rem -= 3 * e;
102*181254a7Smrg if (rem == 1)
103*181254a7Smrg x *= CBRT2;
104*181254a7Smrg else if (rem == 2)
105*181254a7Smrg x *= CBRT4;
106*181254a7Smrg }
107*181254a7Smrg else
108*181254a7Smrg { /* argument less than 1 */
109*181254a7Smrg e = -e;
110*181254a7Smrg rem = e;
111*181254a7Smrg e /= 3;
112*181254a7Smrg rem -= 3 * e;
113*181254a7Smrg if (rem == 1)
114*181254a7Smrg x *= CBRT2I;
115*181254a7Smrg else if (rem == 2)
116*181254a7Smrg x *= CBRT4I;
117*181254a7Smrg e = -e;
118*181254a7Smrg }
119*181254a7Smrg
120*181254a7Smrg /* multiply by power of 2 */
121*181254a7Smrg x = ldexpq (x, e);
122*181254a7Smrg
123*181254a7Smrg /* Newton iteration */
124*181254a7Smrg x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333Q;
125*181254a7Smrg x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333Q;
126*181254a7Smrg x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333Q;
127*181254a7Smrg
128*181254a7Smrg if (sign < 0)
129*181254a7Smrg x = -x;
130*181254a7Smrg return (x);
131*181254a7Smrg }
132