1/* Implementation of cbrtl. IBM Extended Precision version.
2 Cephes Math Library Release 2.2: January, 1991
3 Copyright 1984, 1991 by Stephen L. Moshier
4 Adapted for glibc October, 2001.
5
6 This library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
10
11 This library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with this library; if not, see
18 <https://www.gnu.org/licenses/>. */
19
20/* This file was copied from sysdeps/ieee754/ldbl-128/e_j0l.c. */
21
22
23#include <math_ldbl_opt.h>
24#include <math.h>
25#include <math_private.h>
26
27static const long double CBRT2 = 1.259921049894873164767210607278228350570251L;
28static const long double CBRT4 = 1.587401051968199474751705639272308260391493L;
29static const long double CBRT2I = 0.7937005259840997373758528196361541301957467L;
30static const long double CBRT4I = 0.6299605249474365823836053036391141752851257L;
31
32
33long double
34__cbrtl (long double x)
35{
36 int e, rem, sign;
37 long double z;
38
39 if (!isfinite (x))
40 return x + x;
41
42 if (x == 0)
43 return (x);
44
45 if (x > 0)
46 sign = 1;
47 else
48 {
49 sign = -1;
50 x = -x;
51 }
52
53 z = x;
54 /* extract power of 2, leaving mantissa between 0.5 and 1 */
55 x = __frexpl (x: x, exponent: &e);
56
57 /* Approximate cube root of number between .5 and 1,
58 peak relative error = 1.2e-6 */
59 x = ((((1.3584464340920900529734e-1L * x
60 - 6.3986917220457538402318e-1L) * x
61 + 1.2875551670318751538055e0L) * x
62 - 1.4897083391357284957891e0L) * x
63 + 1.3304961236013647092521e0L) * x + 3.7568280825958912391243e-1L;
64
65 /* exponent divided by 3 */
66 if (e >= 0)
67 {
68 rem = e;
69 e /= 3;
70 rem -= 3 * e;
71 if (rem == 1)
72 x *= CBRT2;
73 else if (rem == 2)
74 x *= CBRT4;
75 }
76 else
77 { /* argument less than 1 */
78 e = -e;
79 rem = e;
80 e /= 3;
81 rem -= 3 * e;
82 if (rem == 1)
83 x *= CBRT2I;
84 else if (rem == 2)
85 x *= CBRT4I;
86 e = -e;
87 }
88
89 /* multiply by power of 2 */
90 x = __ldexpl (x: x, exponent: e);
91
92 /* Newton iteration */
93 x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
94 x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
95 x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
96
97 if (sign < 0)
98 x = -x;
99 return (x);
100}
101
102long_double_symbol (libm, __cbrtl, cbrtl);
103

source code of glibc/sysdeps/ieee754/ldbl-128ibm/s_cbrtl.c