libm/ldouble/cbrtl.c

   1 /*                                                      cbrtl.c
   2  *
   3  *      Cube root, long double precision
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * long double x, y, cbrtl();
  10  *
  11  * y = cbrtl( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns the cube root of the argument, which may be negative.
  18  *
  19  * Range reduction involves determining the power of 2 of
  20  * the argument.  A polynomial of degree 2 applied to the
  21  * mantissa, and multiplication by the cube root of 1, 2, or 4
  22  * approximates the root to within about 0.1%.  Then Newton's
  23  * iteration is used three times to converge to an accurate
  24  * result.
  25  *
  26  *
  27  *
  28  * ACCURACY:
  29  *
  30  *                      Relative error:
  31  * arithmetic   domain     # trials      peak         rms
  32  *    IEEE     .125,8        80000      7.0e-20     2.2e-20
  33  *    IEEE    exp(+-707)    100000      7.0e-20     2.4e-20
  34  *
  35  */
  36 \f
  37
  38 /*
  39 Cephes Math Library Release 2.2: January, 1991
  40 Copyright 1984, 1991 by Stephen L. Moshier
  41 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  42 */
  43
  44
  45 #include <math.h>
  46
  47 static long double CBRT2  = 1.2599210498948731647672L;
  48 static long double CBRT4  = 1.5874010519681994747517L;
  49 static long double CBRT2I = 0.79370052598409973737585L;
  50 static long double CBRT4I = 0.62996052494743658238361L;
  51
  52 #ifdef ANSIPROT
  53 extern long double frexpl ( long double, int * );
  54 extern long double ldexpl ( long double, int );
  55 extern int isnanl ( long double );
  56 #else
  57 long double frexpl(), ldexpl();
  58 extern int isnanl();
  59 #endif
  60
  61 #ifdef INFINITIES
  62 extern long double INFINITYL;
  63 #endif
  64
  65 long double cbrtl(x)
  66 long double x;
  67 {
  68 int e, rem, sign;
  69 long double z;
  70
  71
  72 #ifdef NANS
  73 if(isnanl(x))
  74         return(x);
  75 #endif
  76 #ifdef INFINITIES
  77 if( x == INFINITYL)
  78         return(x);
  79 if( x == -INFINITYL)
  80         return(x);
  81 #endif
  82 if( x == 0 )
  83         return( x );
  84 if( x > 0 )
  85         sign = 1;
  86 else
  87         {
  88         sign = -1;
  89         x = -x;
  90         }
  91
  92 z = x;
  93 /* extract power of 2, leaving
  94  * mantissa between 0.5 and 1
  95  */
  96 x = frexpl( x, &e );
  97
  98 /* Approximate cube root of number between .5 and 1,
  99  * peak relative error = 1.2e-6
 100  */
 101 x = (((( 1.3584464340920900529734e-1L * x
 102        - 6.3986917220457538402318e-1L) * x
 103        + 1.2875551670318751538055e0L) * x
 104        - 1.4897083391357284957891e0L) * x
 105        + 1.3304961236013647092521e0L) * x
 106        + 3.7568280825958912391243e-1L;
 107
 108 /* exponent divided by 3 */
 109 if( e >= 0 )
 110         {
 111         rem = e;
 112         e /= 3;
 113         rem -= 3*e;
 114         if( rem == 1 )
 115                 x *= CBRT2;
 116         else if( rem == 2 )
 117                 x *= CBRT4;
 118         }
 119 else
 120         { /* argument less than 1 */
 121         e = -e;
 122         rem = e;
 123         e /= 3;
 124         rem -= 3*e;
 125         if( rem == 1 )
 126                 x *= CBRT2I;
 127         else if( rem == 2 )
 128                 x *= CBRT4I;
 129         e = -e;
 130         }
 131
 132 /* multiply by power of 2 */
 133 x = ldexpl( x, e );
 134
 135 /* Newton iteration */
 136
 137 x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
 138 x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
 139
 140 if( sign < 0 )
 141         x = -x;
 142 return(x);
 143 }