libquadmath/math/sqrtq.c

   1 #include "quadmath-imp.h"
   2 #include <math.h>
   3 #include <float.h>
   4
   5 __float128
   6 sqrtq (const __float128 x)
   7 {
   8   __float128 y;
   9   int exp;
  10
  11   if (x == 0)
  12     return x;
  13
  14   if (isnanq (x))
  15     return x;
  16
  17   if (x < 0)
  18     return nanq ("");
  19
  20   if (x <= DBL_MAX && x >= DBL_MIN)
  21   {
  22     /* Use double result as starting point.  */
  23     y = sqrt ((double) x);
  24
  25     /* Two Newton iterations.  */
  26     y -= 0.5q * (y - x / y);
  27     y -= 0.5q * (y - x / y);
  28     return y;
  29   }
  30
  31 #ifdef HAVE_SQRTL
  32   if (x <= LDBL_MAX && x >= LDBL_MIN)
  33   {
  34     /* Use long double result as starting point.  */
  35     y = sqrtl ((long double) x);
  36
  37     /* One Newton iteration.  */
  38     y -= 0.5q * (y - x / 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   if (exp % 2)
  47     y *= 2, exp--;
  48
  49   y = sqrt (y);
  50   y = scalbnq (y, exp / 2);
  51
  52   /* Two Newton iterations.  */
  53   y -= 0.5q * (y - x / y);
  54   y -= 0.5q * (y - x / y);
  55   return y;
  56 }
  57