1 /* Compute cubic root of float 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/flt-32/s_cbrtf.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 float 47 cbrtf (float x) /* */ 48 { 49 if (isfinite (x) && x != 0.0f) 50 { 51 float xm, ym, u, t2; 52 int xe; 53 54 /* Reduce X. XM now is an range 1.0 to 0.5. */ 55 xm = frexpf (fabsf (x), &xe); 56 57 u = (0.492659620528969547 58 + (0.697570460207922770 - 0.191502161678719066 * xm) * xm); 59 60 t2 = u * u * u; 61 62 ym = u * (t2 + 2.0 * xm) / (2.0 * t2 + xm) * factor[2 + xe % 3]; 63 64 return ldexpf (x > 0.0 ? ym : -ym, xe / 3); 65 } 66 else 67 return x + x; 68 }