winsup/cygwin/math/cbrtl.c

   1 /**
   2  * This file has no copyright assigned and is placed in the Public Domain.
   3  * This file is part of the mingw-w64 runtime package.
   4  * No warranty is given; refer to the file DISCLAIMER.PD within this package.
   5  */
   6 #include <math.h>
   7
   8 static const long double CBRT2  = 1.2599210498948731647672L;
   9 static const long double CBRT4  = 1.5874010519681994747517L;
  10 static const long double CBRT2I = 0.79370052598409973737585L;
  11 static const long double CBRT4I = 0.62996052494743658238361L;
  12
  13 long double cbrtl(long double x)
  14 {
  15         int e, rem, sign;
  16         long double z;
  17
  18         if (!isfinite (x) || x == 0.0L)
  19                 return (x);
  20
  21         if (x > 0)
  22                 sign = 1;
  23         else
  24         {
  25                 sign = -1;
  26                 x = -x;
  27         }
  28
  29         z = x;
  30         /* extract power of 2, leaving
  31          * mantissa between 0.5 and 1
  32          */
  33         x = frexpl(x, &e);
  34
  35         /* Approximate cube root of number between .5 and 1,
  36          * peak relative error = 1.2e-6
  37          */
  38         x = (((( 1.3584464340920900529734e-1L * x
  39                - 6.3986917220457538402318e-1L) * x
  40                + 1.2875551670318751538055e0L) * x
  41                - 1.4897083391357284957891e0L) * x
  42                + 1.3304961236013647092521e0L) * x
  43                + 3.7568280825958912391243e-1L;
  44
  45         /* exponent divided by 3 */
  46         if (e >= 0)
  47         {
  48                 rem = e;
  49                 e /= 3;
  50                 rem -= 3*e;
  51                 if (rem == 1)
  52                         x *= CBRT2;
  53                 else if (rem == 2)
  54                         x *= CBRT4;
  55         }
  56         else
  57         { /* argument less than 1 */
  58                 e = -e;
  59                 rem = e;
  60                 e /= 3;
  61                 rem -= 3*e;
  62                 if (rem == 1)
  63                         x *= CBRT2I;
  64                 else if (rem == 2)
  65                         x *= CBRT4I;
  66                 e = -e;
  67         }
  68
  69         /* multiply by power of 2 */
  70         x = ldexpl(x, e);
  71
  72         /* Newton iteration */
  73
  74         x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
  75         x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
  76
  77         if (sign < 0)
  78                 x = -x;
  79         return (x);
  80 }