1 | /* cbrtl.c |
2 | * |
3 | * Cube root, long double precision |
4 | * |
5 | * |
6 | * |
7 | * SYNOPSIS: |
8 | * |
9 | * long double x, y, cbrtl(); |
10 | * |
11 | * y = cbrtl( 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 -8,8 100000 1.3e-34 3.9e-35 |
33 | * IEEE exp(+-707) 100000 1.3e-34 4.3e-35 |
34 | * |
35 | */ |
36 | |
37 | /* |
38 | Cephes Math Library Release 2.2: January, 1991 |
39 | Copyright 1984, 1991 by Stephen L. Moshier |
40 | Adapted for glibc October, 2001. |
41 | |
42 | This library is free software; you can redistribute it and/or |
43 | modify it under the terms of the GNU Lesser General Public |
44 | License as published by the Free Software Foundation; either |
45 | version 2.1 of the License, or (at your option) any later version. |
46 | |
47 | This library is distributed in the hope that it will be useful, |
48 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
49 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
50 | Lesser General Public License for more details. |
51 | |
52 | You should have received a copy of the GNU Lesser General Public |
53 | License along with this library; if not, see |
54 | <https://www.gnu.org/licenses/>. */ |
55 | |
56 | |
57 | #include <math.h> |
58 | #include <math_private.h> |
59 | #include <libm-alias-ldouble.h> |
60 | |
61 | static const _Float128 CBRT2 = L(1.259921049894873164767210607278228350570251); |
62 | static const _Float128 CBRT4 = L(1.587401051968199474751705639272308260391493); |
63 | static const _Float128 CBRT2I = L(0.7937005259840997373758528196361541301957467); |
64 | static const _Float128 CBRT4I = L(0.6299605249474365823836053036391141752851257); |
65 | |
66 | |
67 | _Float128 |
68 | __cbrtl (_Float128 x) |
69 | { |
70 | int e, rem, sign; |
71 | _Float128 z; |
72 | |
73 | if (!isfinite (x)) |
74 | return x + x; |
75 | |
76 | if (x == 0) |
77 | return (x); |
78 | |
79 | if (x > 0) |
80 | sign = 1; |
81 | else |
82 | { |
83 | sign = -1; |
84 | x = -x; |
85 | } |
86 | |
87 | z = x; |
88 | /* extract power of 2, leaving mantissa between 0.5 and 1 */ |
89 | x = __frexpl (x, &e); |
90 | |
91 | /* Approximate cube root of number between .5 and 1, |
92 | peak relative error = 1.2e-6 */ |
93 | x = ((((L(1.3584464340920900529734e-1) * x |
94 | - L(6.3986917220457538402318e-1)) * x |
95 | + L(1.2875551670318751538055e0)) * x |
96 | - L(1.4897083391357284957891e0)) * x |
97 | + L(1.3304961236013647092521e0)) * x + L(3.7568280825958912391243e-1); |
98 | |
99 | /* exponent divided by 3 */ |
100 | if (e >= 0) |
101 | { |
102 | rem = e; |
103 | e /= 3; |
104 | rem -= 3 * e; |
105 | if (rem == 1) |
106 | x *= CBRT2; |
107 | else if (rem == 2) |
108 | x *= CBRT4; |
109 | } |
110 | else |
111 | { /* argument less than 1 */ |
112 | e = -e; |
113 | rem = e; |
114 | e /= 3; |
115 | rem -= 3 * e; |
116 | if (rem == 1) |
117 | x *= CBRT2I; |
118 | else if (rem == 2) |
119 | x *= CBRT4I; |
120 | e = -e; |
121 | } |
122 | |
123 | /* multiply by power of 2 */ |
124 | x = __ldexpl (x, e); |
125 | |
126 | /* Newton iteration */ |
127 | x -= (x - (z / (x * x))) * L(0.3333333333333333333333333333333333333333); |
128 | x -= (x - (z / (x * x))) * L(0.3333333333333333333333333333333333333333); |
129 | x -= (x - (z / (x * x))) * L(0.3333333333333333333333333333333333333333); |
130 | |
131 | if (sign < 0) |
132 | x = -x; |
133 | return (x); |
134 | } |
135 | |
136 | libm_alias_ldouble (__cbrt, cbrt) |
137 | |