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