68
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
1 #include "quadmath-imp.h"
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
2 #include <math.h>
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
3 #include <float.h>
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
4
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
5 __float128
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
6 sqrtq (const __float128 x)
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
7 {
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
8 __float128 y;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
9 int exp;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
10
|
111
|
11 if (isnanq (x) || (isinfq (x) && x > 0))
|
|
12 return x;
|
|
13
|
68
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
14 if (x == 0)
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
15 return x;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
16
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
17 if (x < 0)
|
111
|
18 {
|
|
19 /* Return NaN with invalid signal. */
|
|
20 return (x - x) / (x - x);
|
|
21 }
|
68
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
22
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
23 if (x <= DBL_MAX && x >= DBL_MIN)
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
24 {
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
25 /* Use double result as starting point. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
26 y = sqrt ((double) x);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
27
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
28 /* Two Newton iterations. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
29 y -= 0.5q * (y - x / y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
30 y -= 0.5q * (y - x / y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
31 return y;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
32 }
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
33
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
34 #ifdef HAVE_SQRTL
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
35 {
|
111
|
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);
|
68
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
41
|
111
|
42 /* One Newton iteration. */
|
|
43 y -= 0.5q * (y - x / y);
|
|
44 return y;
|
|
45 }
|
68
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
46 }
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
47 #endif
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
48
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
49 /* If we're outside of the range of C types, we have to compute
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
50 the initial guess the hard way. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
51 y = frexpq (x, &exp);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
52 if (exp % 2)
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
53 y *= 2, exp--;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
54
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
55 y = sqrt (y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
56 y = scalbnq (y, exp / 2);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
57
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
58 /* Two Newton iterations. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
59 y -= 0.5q * (y - x / y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
60 y -= 0.5q * (y - x / y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
61 return y;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
62 }
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
63
|