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   {
  36     long double xl = (long double) x;
  37     if (xl <= LDBL_MAX && xl >= LDBL_MIN)
  38       {
  39         /* Use long double result as starting point.  */
  40         y = (__float128) sqrtl (xl);
  41
  42         /* One Newton iteration.  */
  43         y -= 0.5q * (y - x / y);
  44         return y;
  45       }
  46   }
  47 #endif
  48
  49   /* If we're outside of the range of C types, we have to compute
  50      the initial guess the hard way.  */
  51   y = frexpq (x, &exp);
  52   if (exp % 2)
  53     y *= 2, exp--;
  54
  55   y = sqrt (y);
  56   y = scalbnq (y, exp / 2);
  57
  58   /* Two Newton iterations.  */
  59   y -= 0.5q * (y - x / y);
  60   y -= 0.5q * (y - x / y);
  61   return y;
  62 }
  63