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