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 (isnanq (x) || (isinfq (x) && x > 0))
  12     return x;
  13
  14   if (x == 0)
  15     return x;
  16
  17   if (x < 0)
  18     {
  19       /* Return NaN with invalid signal.  */
  20       return (x - x) / (x - x);
  21     }
  22
  23   if (x <= DBL_MAX && x >= DBL_MIN)
  24   {
  25     /* Use double result as starting point.  */
  26     y = sqrt ((double) x);
  27
  28     /* Two Newton iterations.  */
  29     y -= 0.5q * (y - x / y);
  30     y -= 0.5q * (y - x / y);
  31     return y;
  32   }
  33
  34 #ifdef HAVE_SQRTL
  35   if (x <= LDBL_MAX && x >= LDBL_MIN)
  36   {
  37     /* Use long double result as starting point.  */
  38     y = sqrtl ((long double) x);
  39
  40     /* One Newton iteration.  */
  41     y -= 0.5q * (y - x / y);
  42     return y;
  43   }
  44 #endif
  45
  46   /* If we're outside of the range of C types, we have to compute
  47      the initial guess the hard way.  */
  48   y = frexpq (x, &exp);
  49   if (exp % 2)
  50     y *= 2, exp--;
  51
  52   y = sqrt (y);
  53   y = scalbnq (y, exp / 2);
  54
  55   /* Two Newton iterations.  */
  56   y -= 0.5q * (y - x / y);
  57   y -= 0.5q * (y - x / y);
  58   return y;
  59 }
  60