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