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  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
61