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