libquadmath/math/cbrtq.c

   1 #include "quadmath-imp.h"
   2 #include <math.h>
   3 #include <float.h>
   4
   5 __float128
   6 cbrtq (const __float128 x)
   7 {
   8   __float128 y;
   9   int exp, i;
  10
  11   if (x == 0)
  12     return x;
  13
  14   if (isnanq (x))
  15     return x;
  16
  17   if (x <= DBL_MAX && x >= DBL_MIN)
  18   {
  19     /* Use double result as starting point.  */
  20     y = cbrt ((double) x);
  21
  22     /* Two Newton iterations.  */
  23     y -= 0.333333333333333333333333333333333333333333333333333Q
  24           * (y - x / (y * y));
  25     y -= 0.333333333333333333333333333333333333333333333333333Q
  26           * (y - x / (y * y));
  27     return y;
  28   }
  29
  30 #ifdef HAVE_CBRTL
  31   if (x <= LDBL_MAX && x >= LDBL_MIN)
  32   {
  33     /* Use long double result as starting point.  */
  34     y = cbrtl ((long double) x);
  35
  36     /* One Newton iteration.  */
  37     y -= 0.333333333333333333333333333333333333333333333333333Q
  38           * (y - x / (y * y));
  39     return y;
  40   }
  41 #endif
  42
  43   /* If we're outside of the range of C types, we have to compute
  44      the initial guess the hard way.  */
  45   y = frexpq (x, &exp);
  46
  47   i = exp % 3;
  48   y = (i >= 0 ? i : -i);
  49   if (i == 1)
  50     y *= 2, exp--;
  51   else if (i == 2)
  52     y *= 4, exp -= 2;
  53
  54   y = cbrt (y);
  55   y = scalbnq (y, exp / 3);
  56
  57   /* Two Newton iterations.  */
  58   y -= 0.333333333333333333333333333333333333333333333333333Q
  59          * (y - x / (y * y));
  60   y -= 0.333333333333333333333333333333333333333333333333333Q
  61          * (y - x / (y * y));
  62   return y;
  63 }
  64