1 #include "quadmath-imp.h"
2 #include <math.h>
3 #include <float.h>
4 
5 __float128
sqrtq(const __float128 x)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 
64