1 /*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2023 Free Software Foundation, Inc.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <https://www.gnu.org/licenses/>.
18 */
19 /*********************************************************************/
20 /* MODULE_NAME: uroot.c */
21 /* */
22 /* FUNCTION: usqrt */
23 /* */
24 /* FILES NEEDED: dla.h endian.h mydefs.h uroot.h */
25 /* uroot.tbl */
26 /* */
27 /* An ultimate sqrt routine. Given an IEEE double machine number x */
28 /* it computes the correctly rounded (to nearest) value of square */
29 /* root of x. */
30 /* Assumption: Machine arithmetic operations are performed in */
31 /* round to nearest mode of IEEE 754 standard. */
32 /* */
33 /*********************************************************************/
34
35 #include <math_private.h>
36 #include <libm-alias-finite.h>
37
38 typedef union {int64_t i[2]; long double x; double d[2]; } mynumber;
39
40 static const double
41 t512 = 0x1p512,
42 tm256 = 0x1p-256,
43 two54 = 0x1p54, /* 0x4350000000000000 */
44 twom54 = 0x1p-54; /* 0x3C90000000000000 */
45
46 /*********************************************************************/
47 /* An ultimate sqrt routine. Given an IEEE double machine number x */
48 /* it computes the correctly rounded (to nearest) value of square */
49 /* root of x. */
50 /*********************************************************************/
51 long double __ieee754_sqrtl(long double x)
52 {
53 static const long double big = 134217728.0, big1 = 134217729.0;
54 long double t,s,i;
55 mynumber a,c;
56 uint64_t k, l;
57 int64_t m, n;
58 double d;
59
60 a.x=x;
61 k=a.i[0] & INT64_C(0x7fffffffffffffff);
62 /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/
63 if (k>INT64_C(0x000fffff00000000) && k<INT64_C(0x7ff0000000000000)) {
64 if (x < 0) return (big1-big1)/(big-big);
65 l = (k&INT64_C(0x001fffffffffffff))|INT64_C(0x3fe0000000000000);
66 if ((a.i[1] & INT64_C(0x7fffffffffffffff)) != 0) {
67 n = (int64_t) ((l - k) * 2) >> 53;
68 m = (a.i[1] >> 52) & 0x7ff;
69 if (m == 0) {
70 a.d[1] *= two54;
71 m = ((a.i[1] >> 52) & 0x7ff) - 54;
72 }
73 m += n;
74 if (m > 0)
75 a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52);
76 else if (m <= -54) {
77 a.i[1] &= INT64_C(0x8000000000000000);
78 } else {
79 m += 54;
80 a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52);
81 a.d[1] *= twom54;
82 }
83 }
84 a.i[0] = l;
85 s = a.x;
86 d = __ieee754_sqrt (a.d[0]);
87 c.i[0] = INT64_C(0x2000000000000000)+((k&INT64_C(0x7fe0000000000000))>>1);
88 c.i[1] = 0;
89 i = d;
90 t = 0.5L * (i + s / i);
91 i = 0.5L * (t + s / t);
92 return c.x * i;
93 }
94 else {
95 if (k>=INT64_C(0x7ff0000000000000))
96 /* sqrt (-Inf) = NaN, sqrt (NaN) = NaN, sqrt (+Inf) = +Inf. */
97 return x * x + x;
98 if (x == 0) return x;
99 if (x < 0) return (big1-big1)/(big-big);
100 return tm256*__ieee754_sqrtl(x*t512);
101 }
102 }
103 libm_alias_finite (__ieee754_sqrtl, __sqrtl)