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  
      27  static const long double CBRT2 = 1.259921049894873164767210607278228350570251L;
      28  static const long double CBRT4 = 1.587401051968199474751705639272308260391493L;
      29  static const long double CBRT2I = 0.7937005259840997373758528196361541301957467L;
      30  static const long double CBRT4I = 0.6299605249474365823836053036391141752851257L;
      31  
      32  
      33  long 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, &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, 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  
     102  long_double_symbol (libm, __cbrtl, cbrtl);