1 /* Compute cubic root of double value. 2 Copyright (C) 1997, 2012-2021 Free Software Foundation, Inc. 3 4 Contributed by Dirk Alboth <dirka@uni-paderborn.de> and 5 Ulrich Drepper <drepper@cygnus.com>, 1997. 6 7 This file is free software: you can redistribute it and/or modify 8 it under the terms of the GNU Lesser General Public License as 9 published by the Free Software Foundation; either version 3 of the 10 License, or (at your option) any later version. 11 12 This file is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 GNU Lesser General Public License for more details. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along with this program. If not, see <https://www.gnu.org/licenses/>. */ 19 20 #include <config.h> 21 22 /* Specification. */ 23 #include <math.h> 24 25 /* MSVC with option -fp:strict refuses to compile constant initializers that 26 contain floating-point operations. Pacify this compiler. */ 27 #if defined _MSC_VER && !defined __clang__ 28 # pragma fenv_access (off) 29 #endif 30 31 /* Code based on glibc/sysdeps/ieee754/dbl-64/s_cbrt.c. */ 32 33 #define CBRT2 1.2599210498948731648 /* 2^(1/3) */ 34 #define SQR_CBRT2 1.5874010519681994748 /* 2^(2/3) */ 35 36 static const double factor[5] = 37 { 38 1.0 / SQR_CBRT2, 39 1.0 / CBRT2, 40 1.0, 41 CBRT2, 42 SQR_CBRT2 43 }; 44 45 46 double 47 cbrt (double x) /* */ 48 { 49 if (isfinite (x) && x != 0.0) 50 { 51 double xm, ym, u, t2; 52 int xe; 53 54 /* Reduce X. XM now is an range 1.0 to 0.5. */ 55 xm = frexp (fabs (x), &xe); 56 57 u = (0.354895765043919860 58 + ((1.50819193781584896 59 + ((-2.11499494167371287 60 + ((2.44693122563534430 61 + ((-1.83469277483613086 62 + (0.784932344976639262 - 0.145263899385486377 * xm) 63 * xm) 64 * xm)) 65 * xm)) 66 * xm)) 67 * xm)); 68 69 t2 = u * u * u; 70 71 ym = u * (t2 + 2.0 * xm) / (2.0 * t2 + xm) * factor[2 + xe % 3]; 72 73 return ldexp (x > 0.0 ? ym : -ym, xe / 3); 74 } 75 else 76 return x + x; 77 }