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
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
11 if (x == 0)
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
12 return x;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
13
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
14 if (isnanq (x))
|
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)
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
18 return nanq ("");
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
19
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
20 if (x <= DBL_MAX && x >= DBL_MIN)
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
21 {
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
22 /* Use double result as starting point. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
23 y = sqrt ((double) x);
|
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 /* Two Newton iterations. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
26 y -= 0.5q * (y - x / y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
27 y -= 0.5q * (y - x / y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
28 return y;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
29 }
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
30
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
31 #ifdef HAVE_SQRTL
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
32 if (x <= LDBL_MAX && x >= LDBL_MIN)
|
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 /* Use long double result as starting point. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
35 y = sqrtl ((long double) x);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
36
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
37 /* One Newton iteration. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
38 y -= 0.5q * (y - x / y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
39 return y;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
40 }
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
41 #endif
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
42
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
43 /* 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
|
44 the initial guess the hard way. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
45 y = frexpq (x, &exp);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
46 if (exp % 2)
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
47 y *= 2, exp--;
|
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 y = sqrt (y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
50 y = scalbnq (y, exp / 2);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
51
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
52 /* Two Newton iterations. */
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
53 y -= 0.5q * (y - x / y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
54 y -= 0.5q * (y - x / y);
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
55 return y;
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
56 }
|
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff
changeset
|
57
|