2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2024 Free Software Foundation, Inc.
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.
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.
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/>.
19 /*********************************************************************/
20 /* MODULE_NAME: uroot.c */
24 /* FILES NEEDED: dla.h endian.h mydefs.h uroot.h */
27 /* An ultimate sqrt routine. Given an IEEE double machine number x */
28 /* it computes the correctly rounded (to nearest) value of square */
30 /* Assumption: Machine arithmetic operations are performed in */
31 /* round to nearest mode of IEEE 754 standard. */
33 /*********************************************************************/
35 #include <math_private.h>
36 #include <libm-alias-finite.h>
38 typedef union {int64_t i
[2]; long double x
; double d
[2]; } mynumber
;
43 two54
= 0x1p
54, /* 0x4350000000000000 */
44 twom54
= 0x1p
-54; /* 0x3C90000000000000 */
46 /*********************************************************************/
47 /* An ultimate sqrt routine. Given an IEEE double machine number x */
48 /* it computes the correctly rounded (to nearest) value of square */
50 /*********************************************************************/
51 long double __ieee754_sqrtl(long double x
)
53 static const long double big
= 134217728.0, big1
= 134217729.0;
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;
71 m
= ((a
.i
[1] >> 52) & 0x7ff) - 54;
75 a
.i
[1] = (a
.i
[1] & INT64_C(0x800fffffffffffff)) | (m
<< 52);
77 a
.i
[1] &= INT64_C(0x8000000000000000);
80 a
.i
[1] = (a
.i
[1] & INT64_C(0x800fffffffffffff)) | (m
<< 52);
86 d
= __ieee754_sqrt (a
.d
[0]);
87 c
.i
[0] = INT64_C(0x2000000000000000)+((k
&INT64_C(0x7fe0000000000000))>>1);
90 t
= 0.5L * (i
+ s
/ i
);
91 i
= 0.5L * (t
+ s
/ t
);
95 if (k
>=INT64_C(0x7ff0000000000000))
96 /* sqrt (-Inf) = NaN, sqrt (NaN) = NaN, sqrt (+Inf) = +Inf. */
99 if (x
< 0) return (big1
-big1
)/(big
-big
);
100 return tm256
*__ieee754_sqrtl(x
*t512
);
103 libm_alias_finite (__ieee754_sqrtl
, __sqrtl
)