1/* Compute cubic root of double value. 2 Copyright (C) 1997 Free Software Foundation, Inc. 3 This file is part of the GNU C Library. 4 Contributed by Dirk Alboth <dirka@uni-paderborn.de> and 5 Ulrich Drepper <drepper@cygnus.com>, 1997. 6 7 The GNU C Library is free software; you can redistribute it and/or 8 modify it under the terms of the GNU Lesser General Public 9 License as published by the Free Software Foundation; either 10 version 2.1 of the License, or (at your option) any later version. 11 12 The GNU C Library 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 GNU 15 Lesser General Public License for more details. 16 17 You should have received a copy of the GNU Lesser General Public 18 License along with the GNU C Library; if not, write to the Free 19 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 20 02111-1307 USA. */ 21 22#include "math.h" 23#include "math_private.h" 24 25 26#define CBRT2 1.2599210498948731648 /* 2^(1/3) */ 27#define SQR_CBRT2 1.5874010519681994748 /* 2^(2/3) */ 28 29static const double factor[5] = 30{ 31 1.0 / SQR_CBRT2, 32 1.0 / CBRT2, 33 1.0, 34 CBRT2, 35 SQR_CBRT2 36}; 37 38 39double 40__cbrt (double x) 41{ 42 double xm, ym, u, t2; 43 int xe; 44 45 /* Reduce X. XM now is an range 1.0 to 0.5. */ 46 xm = __frexp (fabs (x), &xe); 47 48 /* If X is not finite or is null return it (with raising exceptions 49 if necessary. 50 Note: *Our* version of `frexp' sets XE to zero if the argument is 51 Inf or NaN. This is not portable but faster. */ 52 if (xe == 0 && fpclassify (x) <= FP_ZERO) 53 return x + x; 54 55 u = (0.354895765043919860 56 + ((1.50819193781584896 57 + ((-2.11499494167371287 58 + ((2.44693122563534430 59 + ((-1.83469277483613086 60 + (0.784932344976639262 - 0.145263899385486377 * xm) * xm) 61 * xm)) 62 * xm)) 63 * xm)) 64 * xm)); 65 66 t2 = u * u * u; 67 68 ym = u * (t2 + 2.0 * xm) / (2.0 * t2 + xm) * factor[2 + xe % 3]; 69 70 return __ldexp (x > 0.0 ? ym : -ym, xe / 3); 71} 72weak_alias (__cbrt, cbrt) 73#ifdef NO_LONG_DOUBLE 74strong_alias (__cbrt, __cbrtl) 75weak_alias (__cbrt, cbrtl) 76#endif 77